12 Interpreting the integration results
In the previous chapters, we have applied different integration tools to our example multi-omics dataset. In each of these chapters, we have shown some method-specific visualisations that could be used to interpret the integration results. With the moiraine package, it is possible to use method-agnostic functions to inspect and interpret these results. These functionalities make use of the metadata in the MultiDataSet object to enrich plots and other summary outputs with information about the samples and the omics features, and allow users to obtain consistent visualisations across the integration methods used.
_targets.R script
library(targets)
library(tarchetypes)
library(moiraine)
tar_option_set(
packages = c(
"moiraine",
"MOFA2",
"mixOmics",
"readr",
"tibble",
"tidyr",
"dplyr",
"ggplot2",
"patchwork"
)
)
list(
# Importing data ---------------------------------------------------------------
## Data import using a target factory
import_dataset_csv_factory(
files = c(
system.file("extdata/genomics_dataset.csv", package = "moiraine"),
system.file("extdata/transcriptomics_dataset.csv", package = "moiraine"),
system.file("extdata/metabolomics_dataset.csv", package = "moiraine")
),
col_ids = c("marker", "gene_id", "sample_id"),
features_as_rowss = c(TRUE, TRUE, FALSE),
target_name_suffixes = c("geno", "transcripto", "metabo")
),
## Genomics features metadata file
tar_target(
fmetadata_file_geno,
system.file("extdata/genomics_features_info.csv", package = "moiraine"),
format = "file"
),
## Genomics features metadata import
tar_target(
fmetadata_geno,
import_fmetadata_csv(
fmetadata_file_geno,
col_id = "marker",
col_types = c("chromosome" = "c")
)
),
## Metabolomics features metadata import
import_fmetadata_csv_factory(
files = c(
system.file("extdata/metabolomics_features_info.csv", package = "moiraine")
),
col_ids = c("feature_id"),
target_name_suffixes = c("metabo")
),
## Transcriptomics features metadata import
import_fmetadata_gff_factory(
files = system.file("extdata/bos_taurus_gene_model.gff3", package = "moiraine"),
feature_types = "genes",
add_fieldss = c("Name", "description"),
target_name_suffixes = "transcripto"
),
## Samples metadata import
import_smetadata_csv_factory(
files = system.file("extdata/samples_info.csv", package = "moiraine"),
col_ids = "animal_id",
target_name_suffixes = "all"
),
## Creating omics sets for each dataset
create_omics_set_factory(
datasets = c(data_geno, data_transcripto, data_metabo),
omics_types = c("genomics", "transcriptomics", "metabolomics"),
features_metadatas = c(fmetadata_geno, fmetadata_transcripto, fmetadata_metabo),
samples_metadatas = c(smetadata_all, smetadata_all, smetadata_all)
),
## Creating the MultiDataSet object
tar_target(
mo_set,
create_multiomics_set(
list(set_geno,
set_transcripto,
set_metabo)
)
),
# Inspecting the `MultiDataSet` object -----------------------------------------
## Creating a density plot for each dataset
tar_target(
density_plots,
plot_density_data(
mo_set,
combined = FALSE,
scales = "free"
)
),
## Plotting the relationship between features mean and standard deviation
## for each dataset
tar_target(
mean_sd_plots,
plot_meansd_data(mo_set)
),
## Assessing missing values
tar_target(
n_missing_values,
check_missing_values(mo_set)
),
# Modifying the `MultiDataSet` object ------------------------------------------
## RNAseq differential expression results file
tar_target(
rnaseq_de_res_file,
system.file(
"extdata/transcriptomics_de_results.csv",
package = "moiraine"
),
format = "file"
),
## Reading the RNAseq differential expression results
tar_target(
rnaseq_de_res_df,
read_csv(rnaseq_de_res_file) |>
rename(feature_id = gene_id) |>
mutate(dataset = "rnaseq")
),
## Adding the differential expression results to the MultiDataSet object
tar_target(
mo_set_de,
add_features_metadata(mo_set, rnaseq_de_res_df)
),
# Data pre-processing ----------------------------------------------------------
## Applying transformations to the datasets
transformation_datasets_factory(
mo_set_de,
c("rnaseq" = "vst-deseq2",
"metabolome" = "logx"),
log_bases = 2,
pre_log_functions = zero_to_half_min,
transformed_data_name = "mo_set_transformed"
),
## Density plot for each transformed dataset
tar_target(
density_plots_transformed,
plot_density_data(
mo_set_transformed,
combined = FALSE,
scales = "free"
)
),
## Plotting the mean-SD trend for transformed each dataset
tar_target(
mean_sd_plots_transformed,
plot_meansd_data(mo_set_transformed)
),
## Summary table of the transformations applied
tar_target(
transformation_summary,
get_table_transformations(transformations_runs_list)
),
## Running a PCA on each dataset
pca_complete_data_factory(
mo_set_transformed,
complete_data_name = "mo_set_complete"
),
## PCA screeplots
tar_target(
pca_screeplots,
plot_screeplot_pca(pca_runs_list)
),
## PCA sample plots
tar_target(
pca_sample_plots,
plot_samples_coordinates_pca(
pca_runs_list,
datasets = "snps",
pcs = 1:3,
mo_data = mo_set_de,
colour_upper = "geno_comp_cluster",
shape_upper = "status",
colour_lower = "feedlot"
)
),
# Dataset pre-filtering --------------------------------------------------------
## Unsupervised feature selection based on MAD score
feature_preselection_mad_factory(
mo_set_complete,
to_keep_ns = c("snps" = 1000, "rnaseq" = 1000),
with_ties = TRUE,
filtered_set_target_name = "mo_presel_unsupervised"
),
## Diagnostic plot for MAD-based feature selection
tar_target(
preselection_mad_plot,
plot_feature_preselection_mad(individual_mad_values)
),
## Supervised feature selection based on bruising groups
feature_preselection_splsda_factory(
mo_set_complete,
group = "status",
to_keep_ns = c("snps" = 1000, "rnaseq" = 1000),
filtered_set_target_name = "mo_presel_supervised"
),
## Diagnostic plot for sPLS-DA based feature selection
tar_target(
preselection_splsda_plot,
plot_feature_preselection_splsda(individual_splsda_perf)
),
# Integration with sPLS --------------------------------------------------------
## Creating sPLS input
tar_target(
spls_input,
get_input_spls(
mo_presel_supervised,
mode = "canonical",
datasets = c("rnaseq", "metabolome")
)
),
## Initial PLS run with no feature selection and large number of components
tar_target(
spls_novarsel,
spls_run(
spls_input,
ncomp = 4
)
),
## Cross-validation for number of components
tar_target(
spls_perf_res,
mixOmics::perf(
spls_novarsel,
validation = "Mfold",
folds = 10,
nrepeat = 10,
cpus = 3
)
),
## Plotting cross-validation results (for number of components)
## Can try criterion = 'Q2.total', 'cor.tpred', 'cor.upred', 'RSS.tpred',
## 'RSS.upred' (but avoid 'RSS' and 'PRESS')
tar_target(
spls_perf_plot,
plot(spls_perf_res, criterion = "Q2.total")
),
## Selected value for ncomp
tar_target(
spls_optim_ncomp,
spls_get_optim_ncomp(spls_perf_res, min_ncomp = 2)
),
## Cross-validation for number of features to retain
tar_target(
spls_tune_res,
spls_tune(
spls_input,
ncomp = spls_optim_ncomp,
keepX = seq(10, 100, 10),
keepY = seq(10, 100, 10),
validation = "Mfold",
folds = 10,
nrepeat = 5,
measure = "cor",
cpus = 3
)
),
## Plotting cross-validation results (for number of features)
tar_target(
spls_tune_plot,
spls_plot_tune(spls_tune_res)
),
## Final sPLS run
tar_target(
spls_final_run,
spls_run(
spls_input,
ncomp = spls_optim_ncomp,
keepX = spls_tune_res$choice.keepX,
keepY = spls_tune_res$choice.keepY
)
),
# Integration with sO2PLS ------------------------------------------------------
## Creating sO2PLS input
tar_target(
omicspls_input,
get_input_omicspls(
mo_presel_supervised,
datasets = c("rnaseq", "metabolome")
)
),
## Adjusted cross-validation for number of components
tar_target(
so2pls_cv_adj,
so2pls_crossval_o2m_adjR2(
omicspls_input,
a = 1:5,
ax = seq(0, 10, by = 2),
ay = seq(0, 10, by = 2),
nr_folds = 10,
nr_cores = 6,
seed = 127
)
),
tar_target(
so2pls_cv_adj_res,
so2pls_get_optim_ncomp_adj(so2pls_cv_adj)
),
## Plotting adjusted cross-validation results
tar_target(
so2pls_cv_adj_plot,
so2pls_plot_cv_adj(so2pls_cv_adj)
),
## Standard cross-validation for number of components
tar_target(
so2pls_cv,
so2pls_crossval_o2m(
omicspls_input,
so2pls_cv_adj,
nr_folds = 10,
nr_cores = 6,
seed = 356
)
),
tar_target(
so2pls_cv_res,
so2pls_get_optim_ncomp(so2pls_cv)
),
## Plotting standard cross-validation results
tar_target(
so2pls_cv_plot,
so2pls_plot_cv(so2pls_cv)
),
## Cross-validation for sparsity parameters
tar_target(
so2pls_cv_sparsity,
so2pls_crossval_sparsity(
omicspls_input,
n = so2pls_cv_res["n"],
nx = so2pls_cv_res["nx"],
ny = so2pls_cv_res["ny"],
nr_folds = 10,
keepx_seq = c(seq(5, 30, 5), seq(40, 100, 10)),
keepy_seq = c(seq(5, 40, 5))
)
),
tar_target(
so2pls_cv_sparsity_res,
so2pls_get_optim_keep(so2pls_cv_sparsity)
),
## Plotting the results of the cross-validation for the number of features
## to retain from each dataset for the different joint components
tar_target(
so2pls_cv_sparsity_plot,
so2pls_plot_cv_sparsity(so2pls_cv_sparsity)
),
## Extracting sparsity results in table format
tar_target(
so2pls_cv_sparsity_table,
so2pls_print_cv_sparsity(so2pls_cv_sparsity_res)
),
## Final sO2PLS run
tar_target(
so2pls_final_run,
so2pls_o2m(
omicspls_input,
so2pls_cv_res,
so2pls_cv_sparsity_res
)
),
## Summary plot of percentage of variance explained
tar_target(
so2pls_summary_plot,
so2pls_plot_summary(so2pls_final_run)
),
## Screeplot
tar_target(
so2pls_screeplot,
so2pls_screeplot(so2pls_final_run)
),
## Comparison of samples score for joint components
tar_target(
so2pls_joint_components_comparison_plot,
so2pls_compare_samples_joint_components(
so2pls_final_run,
mo_data = mo_set_de,
colour_by = "status",
shape_by = "feedlot"
)
),
## Coefficient plot for joint components
tar_target(
so2pls_joint_components_coefficients_plot,
so2pls_plot_joint_components_coefficients(so2pls_final_run)
),
## Joint component samples score plot
tar_target(
so2pls_joint_components_samples_score_plot,
so2pls_plot_samples_joint_components(
so2pls_final_run,
mo_data = mo_set_de,
colour_upper = "status",
scale_colour_upper = scale_colour_brewer(palette = "Paired"),
shape_upper = "feedlot"
) +
theme(legend.box = "vertical")
),
## Specific components samples score plot
tar_target(
so2pls_specific_components_samples_score_plot,
so2pls_plot_samples_specific_components(
so2pls_final_run,
mo_data = mo_set_de,
colour_upper = "feedlot",
scale_colour_upper = scale_colour_brewer(palette = "Paired"),
colour_lower = "rnaseq_batch",
shape_upper = "gender"
) |>
map(\(x) x + theme(legend.box = "vertical"))
),
# Integration with MOFA --------------------------------------------------------
## Creating MOFA input
tar_target(
mofa_input,
get_input_mofa(
mo_presel_supervised,
options_list = list(
data_options = list(scale_views = TRUE),
model_options = list(likelihoods = c(
"snps" = "poisson",
"rnaseq" = "gaussian",
"metabolome" = "gaussian")
),
training_options = list(seed = 43)
),
only_common_samples = FALSE
)
),
## Overview plot of the samples in each dataset
tar_target(
mofa_input_plot,
plot_data_overview(mofa_input)
),
## Training MOFA model
tar_target(
mofa_trained,
run_mofa(
mofa_input,
save_data = TRUE,
use_basilisk = TRUE
)
),
## Formatting MOFA output
tar_target(
mofa_output,
get_output(mofa_trained)
),
## Plots of variance explained
tar_target(
mofa_var_explained_plot,
plot_variance_explained(
mofa_trained,
x = "view", ## datasets on the x-axis
y = "factor" ## factors on the y-axis
)
),
tar_target(
mofa_total_var_explained_plot,
plot_variance_explained(
mofa_trained,
x = "view",
y = "factor",
plot_total = TRUE
)[[2]]
),
## Plot of factors correlation with covariates
tar_target(
mofa_factors_covariates_cor_plot,
mofa_plot_cor_covariates(mofa_trained)
),
# Integration with DIABLO ------------------------------------------------------
## Creating the DIABLO input
tar_target(
diablo_input,
get_input_mixomics_supervised(
mo_presel_supervised,
group = "status"
)
),
## Running sPLS on each dataset to construct the design matrix
diablo_pairwise_pls_factory(diablo_input),
## Initial DIABLO run with no feature selection and large number of components
tar_target(
diablo_novarsel,
diablo_run(
diablo_input,
diablo_design_matrix,
ncomp = 7
)
),
## Cross-validation for number of components
tar_target(
diablo_perf_res,
mixOmics::perf(
diablo_novarsel,
validation = "Mfold",
folds = 10,
nrepeat = 10,
cpus = 3
)
),
## Plotting cross-validation results (for number of components)
tar_target(
diablo_perf_plot,
diablo_plot_perf(diablo_perf_res)
),
## Selected value for ncomp
tar_target(
diablo_optim_ncomp,
diablo_get_optim_ncomp(diablo_perf_res)
),
## Cross-validation for number of features to retain
tar_target(
diablo_tune_res,
diablo_tune(
diablo_input,
diablo_design_matrix,
ncomp = diablo_optim_ncomp,
validation = "Mfold",
folds = 10,
nrepeat = 5,
dist = "centroids.dist",
cpus = 3
)
),
## Plotting cross-validation results (for number of features)
tar_target(
diablo_tune_plot,
diablo_plot_tune(diablo_tune_res)
),
## Final DIABLO run
tar_target(
diablo_final_run,
diablo_run(
diablo_input,
diablo_design_matrix,
ncomp = diablo_optim_ncomp,
keepX = diablo_tune_res$choice.keepX
)
)
)12.1 Dimension reduction methods
Despite relying on very different statistical approaches, the different integration methods included in the pipeline all perform dimension reduction of the omics datasets through feature extraction. That is, they construct a small number of latent components/variables/dimensions (that we refer to as latent dimensions in the moiraine package) that capture as much information from the original datasets as possible. A dimension reduction approach typically returns, for each latent dimension constructed, two sets of values:
Features weight: the contribution of the features from the different omics datasets to the latent dimension. All methods included in the pipeline construct latent dimensions as linear combinations of the original features, and therefore the features contribution is quantified by their weight in the linear combination for the corresponding latent dimension.
Samples score: the projection of the samples onto the latent dimension.
In addition, the fraction or percentage of variance that each latent dimension explains in the different omics datasets is usually calculated.
Interpreting the results of a dimension reduction method involves:
Understanding the source of the variation captured by each latent dimension: is a given latent dimension representing an important source of biological variation, such as effect of a treatment, or age of the samples? Or do they show a source of technical variation, for example highlighting a group of outlier samples with different omics profiles from the rest of the observations? Answering these questions allows us to identify which latent dimensions capture the biological phenomenon investigated, or whether there are some sources of noise that should be accounted for in follow-up experiments.
Investigating which omics features are driving the latent dimensions: once we have identified some latent dimensions of interest, we can look at the features that contribute the most to these dimensions, in order to understand the molecular mechanisms or pathways involved. This is typically done after looking into the phenomenon captured by the latent dimensions, but can also help to identify it.
12.2 Generating a standardised output
12.2.1 get_output function
In the moiraine package, the output of the different integration methods can be converted to a standardised output containing the pieces of information (features weight, samples score and percentage of variance explained) characteristic of dimension reduction methods, stored in a consistent format. This enables us to use functions for visualisation or analysis that can be applied to the results of any integration method, rather than having to implement one for each object type returned by the different integration packages.
The get_output() function transforms the output from any integration package included in moiraine into an output_dimension_reduction object, which is a list with three tibbles: features_weight, samples_score and variance_explained:
tar_target(
spls_output,
get_output(spls_final_run)
)tar_load(spls_output)
spls_output
#> $features_weight
#> # A tibble: 2,098 × 5
#> feature_id dataset latent_dimension weight importance
#> <chr> <fct> <fct> <dbl> <dbl>
#> 1 ENSBTAG00000000020 rnaseq Component 1 0 0
#> 2 ENSBTAG00000000046 rnaseq Component 1 0 0
#> 3 ENSBTAG00000000056 rnaseq Component 1 0 0
#> 4 ENSBTAG00000000061 rnaseq Component 1 0 0
#> 5 ENSBTAG00000000113 rnaseq Component 1 0 0
#> 6 ENSBTAG00000000149 rnaseq Component 1 0 0
#> 7 ENSBTAG00000000164 rnaseq Component 1 0 0
#> 8 ENSBTAG00000000205 rnaseq Component 1 0 0
#> 9 ENSBTAG00000000212 rnaseq Component 1 0 0
#> 10 ENSBTAG00000000289 rnaseq Component 1 0 0
#> # ℹ 2,088 more rows
#>
#> $samples_score
#> # A tibble: 278 × 3
#> sample_id latent_dimension score
#> <chr> <fct> <dbl>
#> 1 G1979 Component 1 -3.63
#> 2 G2500 Component 1 -2.46
#> 3 G3030 Component 1 -6.42
#> 4 G3068 Component 1 -6.16
#> 5 G3121 Component 1 -3.39
#> 6 G3315 Component 1 4.98
#> 7 G3473 Component 1 5.29
#> 8 G3474 Component 1 4.97
#> 9 G3550 Component 1 -4.17
#> 10 G3594 Component 1 4.71
#> # ℹ 268 more rows
#>
#> $variance_explained
#> # A tibble: 4 × 3
#> latent_dimension dataset prop_var_expl
#> <fct> <fct> <dbl>
#> 1 Component 1 rnaseq 0.205
#> 2 Component 1 metabolome 0.179
#> 3 Component 2 rnaseq 0.147
#> 4 Component 2 metabolome 0.108
#>
#> attr(,"class")
#> [1] "output_dimension_reduction"
#> attr(,"method")
#> [1] "sPLS"tar_target(
so2pls_output,
get_output(so2pls_final_run)
)tar_load(so2pls_output)
so2pls_output
#> $features_weight
#> # A tibble: 2,318 × 5
#> feature_id dataset latent_dimension weight importance
#> <chr> <fct> <fct> <dbl> <dbl>
#> 1 ENSBTAG00000000020 rnaseq joint component 1 0 0
#> 2 ENSBTAG00000000046 rnaseq joint component 1 0 0
#> 3 ENSBTAG00000000056 rnaseq joint component 1 0 0
#> 4 ENSBTAG00000000061 rnaseq joint component 1 0 0
#> 5 ENSBTAG00000000113 rnaseq joint component 1 0 0
#> 6 ENSBTAG00000000149 rnaseq joint component 1 0 0
#> 7 ENSBTAG00000000164 rnaseq joint component 1 0 0
#> 8 ENSBTAG00000000205 rnaseq joint component 1 0 0
#> 9 ENSBTAG00000000212 rnaseq joint component 1 0 0
#> 10 ENSBTAG00000000289 rnaseq joint component 1 0 0
#> # ℹ 2,308 more rows
#>
#> $samples_score
#> # A tibble: 973 × 3
#> sample_id latent_dimension score
#> <chr> <fct> <dbl>
#> 1 G1979 joint component 1 -8.12
#> 2 G2500 joint component 1 -9.81
#> 3 G3030 joint component 1 -9.93
#> 4 G3068 joint component 1 -14.7
#> 5 G3121 joint component 1 -5.01
#> 6 G3315 joint component 1 9.66
#> 7 G3473 joint component 1 9.79
#> 8 G3474 joint component 1 9.99
#> 9 G3550 joint component 1 -9.70
#> 10 G3594 joint component 1 8.37
#> # ℹ 963 more rows
#>
#> $variance_explained
#> # A tibble: 8 × 3
#> latent_dimension dataset prop_var_expl
#> <fct> <fct> <dbl>
#> 1 joint component 1 rnaseq 0.327
#> 2 joint component 1 metabolome 0.150
#> 3 rnaseq specific component 1 rnaseq 0.123
#> 4 metabolome specific component 1 metabolome 0.107
#> 5 metabolome specific component 2 metabolome 0.0767
#> 6 metabolome specific component 3 metabolome 0.0613
#> 7 metabolome specific component 4 metabolome 0.0605
#> 8 metabolome specific component 5 metabolome 0.0376
#>
#> attr(,"class")
#> [1] "output_dimension_reduction"
#> attr(,"method")
#> [1] "sO2PLS"tar_target(
mofa_output,
get_output(mofa_trained)
)tar_load(mofa_output)
mofa_output
#> $features_weight
#> # A tibble: 30,735 × 5
#> feature_id dataset latent_dimension weight importance
#> <chr> <fct> <fct> <dbl> <dbl>
#> 1 21-25977541-C-T-rs41974686 snps Factor 1 0.000374 0.139
#> 2 22-51403583-A-C-rs210306176 snps Factor 1 0.0000535 0.0199
#> 3 24-12959068-G-T-rs381471286 snps Factor 1 0.000268 0.0999
#> 4 8-85224224-T-C-rs43565287 snps Factor 1 -0.000492 0.183
#> 5 ARS-BFGL-BAC-16973 snps Factor 1 -0.000883 0.329
#> 6 ARS-BFGL-BAC-19403 snps Factor 1 -0.000500 0.186
#> 7 ARS-BFGL-BAC-2450 snps Factor 1 0.000555 0.207
#> 8 ARS-BFGL-BAC-2600 snps Factor 1 -0.0000325 0.0121
#> 9 ARS-BFGL-BAC-27911 snps Factor 1 -0.000568 0.212
#> 10 ARS-BFGL-BAC-35925 snps Factor 1 0.000803 0.299
#> # ℹ 30,725 more rows
#>
#> $samples_score
#> # A tibble: 2,160 × 3
#> sample_id latent_dimension score
#> <chr> <fct> <dbl>
#> 1 G1979 Factor 1 1.95
#> 2 G2500 Factor 1 1.98
#> 3 G3030 Factor 1 3.39
#> 4 G3068 Factor 1 3.65
#> 5 G3121 Factor 1 1.57
#> 6 G3315 Factor 1 -1.81
#> 7 G3473 Factor 1 -2.63
#> 8 G3474 Factor 1 -2.30
#> 9 G3550 Factor 1 2.98
#> 10 G3594 Factor 1 -1.88
#> # ℹ 2,150 more rows
#>
#> $variance_explained
#> # A tibble: 45 × 3
#> latent_dimension dataset prop_var_expl
#> <fct> <fct> <dbl>
#> 1 Factor 1 snps 0.000313
#> 2 Factor 1 rnaseq 0.496
#> 3 Factor 1 metabolome 0.137
#> 4 Factor 2 snps 0.239
#> 5 Factor 2 rnaseq 0.000197
#> 6 Factor 2 metabolome 0.0108
#> 7 Factor 3 snps 0.000108
#> 8 Factor 3 rnaseq 0.199
#> 9 Factor 3 metabolome 0.0000872
#> 10 Factor 4 snps 0.0000927
#> # ℹ 35 more rows
#>
#> attr(,"class")
#> [1] "output_dimension_reduction"
#> attr(,"method")
#> [1] "MOFA"tar_target(
diablo_output,
get_output(diablo_final_run)
)tar_load(diablo_output)
diablo_output
#> $features_weight
#> # A tibble: 8,196 × 5
#> feature_id dataset latent_dimension weight importance
#> <chr> <fct> <fct> <dbl> <dbl>
#> 1 21-25977541-C-T-rs41974686 snps Component 1 0 0
#> 2 22-51403583-A-C-rs210306176 snps Component 1 0 0
#> 3 24-12959068-G-T-rs381471286 snps Component 1 0 0
#> 4 8-85224224-T-C-rs43565287 snps Component 1 0 0
#> 5 ARS-BFGL-BAC-16973 snps Component 1 0 0
#> 6 ARS-BFGL-BAC-19403 snps Component 1 0 0
#> 7 ARS-BFGL-BAC-2450 snps Component 1 0 0
#> 8 ARS-BFGL-BAC-2600 snps Component 1 0 0
#> 9 ARS-BFGL-BAC-27911 snps Component 1 0 0
#> 10 ARS-BFGL-BAC-35925 snps Component 1 0 0
#> # ℹ 8,186 more rows
#>
#> $samples_score
#> # A tibble: 540 × 3
#> sample_id latent_dimension score
#> <chr> <fct> <dbl>
#> 1 G1979 Component 1 1.97
#> 2 G2500 Component 1 1.50
#> 3 G3030 Component 1 3.57
#> 4 G3068 Component 1 3.50
#> 5 G3121 Component 1 2.13
#> 6 G3315 Component 1 -3.26
#> 7 G3473 Component 1 -3.97
#> 8 G3474 Component 1 -3.20
#> 9 G3550 Component 1 2.46
#> 10 G3594 Component 1 -2.38
#> # ℹ 530 more rows
#>
#> $variance_explained
#> # A tibble: 12 × 3
#> latent_dimension dataset prop_var_expl
#> <fct> <fct> <dbl>
#> 1 Component 1 snps 0.0698
#> 2 Component 1 rnaseq 0.201
#> 3 Component 1 metabolome 0.171
#> 4 Component 2 snps 0.0213
#> 5 Component 2 rnaseq 0.0801
#> 6 Component 2 metabolome 0.0521
#> 7 Component 3 snps 0.00898
#> 8 Component 3 rnaseq 0.149
#> 9 Component 3 metabolome 0.0634
#> 10 Component 4 snps 0.0136
#> 11 Component 4 rnaseq 0.0601
#> 12 Component 4 metabolome 0.0399
#>
#> attr(,"class")
#> [1] "output_dimension_reduction"
#> attr(,"method")
#> [1] "DIABLO"The features_weight tibble contains one row per combination of feature and latent dimension. The ID of the features and the name of the dataset from which they originate are stored in the feature_id and dataset columns, respectively. The latent_dimension column gives the name of the latent dimension; this is a factor column. For each feature and latent dimension, the weight column shows the weight that was attributed to the feature for the corresponding latent dimension. In addition, the importance column contains the features importance score, which is computed as the absolute value of the features weight, divided by the maximum absolute weight across all features from the same omics dataset for the corresponding latent dimension. This importance score allows us to compare the contribution of the features across latent dimensions or integration methods, as the weight can be on different scales and thus cannot be directly compared. The importance scores range from 0 to 1. For any method performing feature selection (e.g. sPLS or DIABLO), features that were not selected for a given latent dimension are assigned a weight and importance score of 0.
The samples score tibble contains for each sample (sample_id) and latent dimension (latent_dimension) the sample’s coordinate for the corresponding latent dimension.
The variance_explained tibble gives for each latent dimension (latent_dimension) the proportion of variance explained (prop_var_expl) for each dataset (dataset). The values in prop_var_expl are between 0 and 1.
In addition, the name of the integration method used to obtain these results is stored in the method attribute of the output_dimension_reduction object:
For convenience, the get_latent_dimensions() function can be used on an output_dimension_reduction object to see the names of the latent dimensions (the levels used for the latent_dimension column in each tibble):
get_latent_dimensions(spls_output)
#> [1] "Component 1" "Component 2"get_latent_dimensions(so2pls_output)
#> [1] "joint component 1" "rnaseq specific component 1"
#> [3] "metabolome specific component 1" "metabolome specific component 2"
#> [5] "metabolome specific component 3" "metabolome specific component 4"
#> [7] "metabolome specific component 5"get_latent_dimensions(mofa_output)
#> [1] "Factor 1" "Factor 2" "Factor 3" "Factor 4" "Factor 5" "Factor 6"
#> [7] "Factor 7" "Factor 8" "Factor 9" "Factor 10" "Factor 11" "Factor 12"
#> [13] "Factor 13" "Factor 14" "Factor 15"get_latent_dimensions(diablo_output)
#> [1] "Component 1" "Component 2" "Component 3" "Component 4"get_output
Note that both PCA and sPLS-DA (the method used for supervised features preselection in Section 7.3.2) are also dimension reduction methods. Therefore, the get_output() function also converts pcaRes objects (from run_pca() or pcaMethods::pca()) and mixo_splsda objects (from run_splsda() or mixOmics::splsda()) into output_dimension_reduction objects.
12.2.2 Averaging latent dimensions over datasets
While MOFA computes one score per sample for each latent dimension created, sPLS, sO2PLS and DIABLO all compute one score per dataset for each sample and latent dimension. For each latent dimension, the samples score obtained for the different datasets are then compared, to assess the agreement or covariation between datasets. Ideally, these scores should be highly correlated across datasets, since the methods aim at maximising the variation between datasets, but it is not always the case. However, when they are highly correlated, it becomes redundant to interpret each dataset-specific version of the latent dimensions.
Instead, the mixOmics authors proposed a solution for DIABLO, which is to construct a weighted average space in which to represent the samples: for each latent component, the samples score are averaged over the different datasets. The weight attributed to each dataset is determined by how well the corresponding dataset discriminates between the samples group of interest. This way, rather than looking at samples score for each dataset for any given latent component, we can look at an average of them.
The get_output() function uses this idea to construct, for the output of sPLS, sO2PLS and DIABLO a set of average samples score for each latent dimension, rather than returning a set of samples score per dataset. For DIABLO, the average is weighted as explained above, while for sPLS and sO2PLS each dataset is given equal weight in the average. This calculation can be disabled in the get_output() function to extract the dataset-specific samples score, by setting the use_average_dimensions parameter to FALSE. Note that this only affects the samples_score tibble in terms of dimensions, but the name of the latent dimensions will change to reflect the dataset to which they refer.
tar_target(
spls_output_no_average,
get_output(spls_final_run, use_average_dimensions = FALSE)
)tar_load(spls_output_no_average)
get_latent_dimensions(spls_output_no_average)
#> [1] "rnaseq Component 1" "metabolome Component 1" "rnaseq Component 2"
#> [4] "metabolome Component 2"
nrow(spls_output$samples_score)
#> [1] 278
nrow(spls_output_no_average$samples_score)
#> [1] 556tar_target(
so2pls_output_no_average,
get_output(so2pls_final_run, use_average_dimensions = FALSE)
)tar_load(so2pls_output_no_average)
get_latent_dimensions(so2pls_output_no_average)
#> [1] "rnaseq joint component 1" "metabolome joint component 1"
#> [3] "rnaseq specific component 1" "metabolome specific component 1"
#> [5] "metabolome specific component 2" "metabolome specific component 3"
#> [7] "metabolome specific component 4" "metabolome specific component 5"
nrow(so2pls_output$samples_score)
#> [1] 973
nrow(so2pls_output_no_average$samples_score)
#> [1] 1112tar_target(
diablo_output_no_average,
get_output(diablo_final_run, use_average_dimensions = FALSE)
)tar_load(diablo_output_no_average)
get_latent_dimensions(diablo_output_no_average)
#> [1] "snps Component 1" "rnaseq Component 1" "metabolome Component 1"
#> [4] "snps Component 2" "rnaseq Component 2" "metabolome Component 2"
#> [7] "snps Component 3" "rnaseq Component 3" "metabolome Component 3"
#> [10] "snps Component 4" "rnaseq Component 4" "metabolome Component 4"
nrow(diablo_output$samples_score)
#> [1] 540
nrow(diablo_output_no_average$samples_score)
#> [1] 162012.3 Percentage of variance explained
The first step in interpreting the results of a dimension reduction method is to assess how much variance is explained by the different latent dimensions for each dataset. The function plot_variance_explained() takes as input a output_dimension_reduction object and constructs a dataset-specific screeplot:
plot_variance_explained(spls_output)
plot_variance_explained(so2pls_output, ncol = 1) +
theme(axis.text.x = element_text(size = 9, angle = 30, hjust = 1))
plot_variance_explained(mofa_output, ncol = 1) +
theme(axis.text.x = element_text(size = 9, angle = 30, hjust = 1))
plot_variance_explained(diablo_output, ncol = 2)
12.4 Samples plot
To understand the phenomenon driving the variation represented by the different latent components, we can investigate the coordinates of the samples in the space spanned by these latent components, in combination with information that we have about the samples.
We will start by loading the MultiDataSet object and generating some custom colour palettes for different samples covariates of interest, that we will use in the rest of the section.
tar_load(mo_set_complete)
## Choosing colour palettes for the samples covariates
palette_status <- scale_colour_brewer(palette = "Set1")
palette_feedlot <- scale_colour_brewer(palette = "Set2")
palette_geno_comp <- scale_colour_brewer(palette = "Dark2")
palette_rnaseq_batch <- scale_colour_brewer(palette = "Paired")
## For some plots will need both colour and fill
palette_status_fill <- scale_colour_brewer(
palette = "Set1",
aesthetics = c("colour", "fill")
)12.4.1 Matrix of scatter plots
The plot_samples_score() function shows the samples score in a matrix of scatter plots that represent all possible two-by-two combinations of the latent dimensions. The plots in the resulting matrix are redundant: the lower plots (below the diagonal) are just a rotated version of the upper plots (above the diagonal). The function takes as input an output_dimension_reduction object. In addition, it accepts a MultiDataSet object to extract information about the samples that can then be used to customise the plot. Since the plots are redundant, different properties of the samples can be shown in the upper and lower plots. For example, we will colour the samples according to disease status in the upper plots, while representing some other covariate in the lower plots. We will use the shape of the points to illustrate gender in both upper and lower plots. By default, all latent dimensions are included in the plots, but it is possible to select some of them by passing their names to the latent_dimensions argument. There are more options to further customise the plot, which are shown in the function’s help.
plot_samples_score(
spls_output,
mo_data = mo_set_complete,
colour_upper = "status",
scale_colour_upper = palette_status,
shape_upper = "gender",
colour_lower = "feedlot",
scale_colour_lower = palette_feedlot
) +
theme(legend.box = "vertical")
Component 1 clearly splits the control and BRD animals, while component 2 seems to separate two samples from the rest. There is no clear separation of the samples by feedlot or gender.
Here we will focus on the joint components. Since there is only one, instead of a matrix of scatterplot, the function will display the samples score with a violin plot, and the properties selected for the upper plots will be used to customise the points:
plot_samples_score(
so2pls_output,
latent_dimensions = "joint component 1",
mo_data = mo_set_complete,
colour_upper = "status",
scale_colour_upper = palette_status,
shape_upper = "gender",
colour_lower = "feedlot",
scale_colour_lower = palette_feedlot
)
#> Warning in plot_samples_score(so2pls_output, latent_dimensions = "joint
#> component 1", : Only one latent dimension to plot; 'colour_diag',
#> 'colour_lower' and 'shape_lower' argument will be ignored.
Joint component 1 separates the BRD and control animals, although some of them overlap.
For clarity we will focus here on the first four factors:
plot_samples_score(
mofa_output,
latent_dimensions = paste("Factor", 1:4),
mo_data = mo_set_complete,
colour_upper = "status",
scale_colour_upper = palette_status,
shape_upper = "gender",
colour_lower = "geno_comp_cluster",
scale_colour_lower = palette_geno_comp
) +
theme(legend.box = "vertical")
#> Warning: Removed 2 rows containing missing values (`geom_point()`).
#> Removed 2 rows containing missing values (`geom_point()`).
#> Removed 2 rows containing missing values (`geom_point()`).
#> Removed 2 rows containing missing values (`geom_point()`).
#> Removed 2 rows containing missing values (`geom_point()`).
#> Removed 2 rows containing missing values (`geom_point()`).
Factor 1 clearly separates the control and BRD animals. Factor 2 separates samples from the genomics composition cluster K1 from the other two clusters. None of the factors are separating the animals by gender, indicating that this covariate does not have a strong effect in the omics measurements.
plot_samples_score(
diablo_output,
mo_data = mo_set_complete,
colour_upper = "status",
scale_colour_upper = palette_status,
shape_upper = "gender",
colour_lower = "rnaseq_batch",
scale_colour_lower = palette_rnaseq_batch
) +
theme(legend.box = "vertical")
As DIABLO is a supervised method, it is no surprise that the first component clearly separates the two disease status groups. The second latent component separates samples according to the RNAseq batch. From this plot it is not clear what components 3 and 4 represent.
12.4.2 Scatter plot for pair of latent dimensions
We can also focus on a specific pair of latent dimensions and represent the samples in the space spanned by these two dimensions, with the plot_samples_score_pair() function. The two latent dimensions to plot are selected by passing their names to the function. This function also accepts a MultiDataSet object, whose samples metadata will be used to customise the plot.
plot_samples_score_pair(
spls_output,
paste("Component", 1:2),
mo_data = mo_set_complete,
colour_by = "status",
shape_by = "gender"
) +
palette_status
The two disease status groups are well separated in this space.
To showcase this function, we will look at the first two metabolomics-specific latent components:
plot_samples_score_pair(
so2pls_output,
paste("metabolome specific component", 1:2),
mo_data = mo_set_complete,
colour_by = "status",
shape_by = "gender"
) +
palette_status
These two specific components are not related to disease status or gender.
plot_samples_score_pair(
mofa_output,
paste("Factor", 1:2),
mo_data = mo_set_complete,
colour_by = "status",
shape_by = "geno_comp_cluster"
) +
palette_status
#> Warning: Removed 2 rows containing missing values (`geom_point()`).
We see that these first two factors discriminate the disease status groups and the genomics composition clusters.
plot_samples_score_pair(
diablo_output,
paste("Component", 1:2),
mo_data = mo_set_complete,
colour_by = "status",
shape_by = "rnaseq_batch"
) +
palette_status
We see that these first two components discriminate the disease status groups and the RNAseq batches.
12.4.3 Samples score vs covariate
Lastly, we can plot the samples score for each latent dimension against a samples covariate of interest, with the plot_samples_score_covariate() function. This is useful to focus on one dimension at a time rather than pairs of them. By default, the function displays all latent dimensions, but we can focus on a subset of them through the latent_dimensions arguments. The function needs a MultiDataSet object to extract information about the samples. The plot generated will depend on whether the covariate of interest is categorical or continuous. We will show both options below. In addition to the main covariate of interest that will be displayed on the x-axis, the colour and shape of the points can be further customised according to other samples properties.
First, with a categorical covariate:
plot_samples_score_covariate(
spls_output,
mo_set_complete,
"status",
colour_by = "status"
) +
palette_status_fill
For clarity we will focus on the first joint component only.
plot_samples_score_covariate(
so2pls_output,
mo_set_complete,
"status",
colour_by = "status",
latent_dimensions = "joint component 1"
) +
palette_status_fill
For clarity we will focus on the first two factors only.
plot_samples_score_covariate(
mofa_output,
mo_set_complete,
"status",
colour_by = "status",
shape_by = "geno_comp_cluster",
latent_dimensions = paste("Factor", 1:2)
) +
palette_status_fill
#> Warning: Removed 4 rows containing missing values (`geom_point()`).
For clarity we will focus on the first two components only.
plot_samples_score_covariate(
diablo_output,
mo_set_complete,
"status",
colour_by = "status",
shape_by = "rnaseq_batch",
latent_dimensions = paste("Component", 1:2)
) +
palette_status_fill
Now, with a continuous covariate:
plot_samples_score_covariate(
spls_output,
mo_set_complete,
"day_on_feed",
colour_by = "status"
) +
palette_status_fill
For clarity we will focus on the first joint component only.
plot_samples_score_covariate(
so2pls_output,
mo_set_complete,
"day_on_feed",
colour_by = "status",
latent_dimensions = "joint component 1"
) +
palette_status_fill
For clarity we will focus on the first two factors only.
plot_samples_score_covariate(
mofa_output,
mo_set_complete,
"day_on_feed",
colour_by = "status",
shape_by = "geno_comp_cluster",
latent_dimensions = paste("Factor", 1:2)
) +
palette_status_fill
#> Warning: Removed 4 rows containing non-finite values (`stat_smooth()`).
#> Warning: Removed 4 rows containing missing values (`geom_point()`).
For clarity we will focus on the first two components only.
plot_samples_score_covariate(
diablo_output,
mo_set_complete,
"day_on_feed",
colour_by = "status",
shape_by = "rnaseq_batch",
latent_dimensions = paste("Component", 1:2)
) +
palette_status_fill
12.5 Feature plots
Next, we can have a look at the features contributing to each latent dimension.
12.5.1 Features weight distribution
We can start by looking at the distribution of features weight for each latent dimension across the datasets, with the plot_features_weight_distr() function. Again, the function takes an output_dimension_reduction object as first input; by default it will show all latent dimensions and datasets, but they can be specified through the latent_dimensions and datasets arguments. The function shows by default the distribution of the “signed” importance scores, i.e. the features weights normalised by the highest absolute weight value for the corresponding latent dimension and dataset, but we can also show the distribution of importance scores (which are absolute values) and non-transformed features weight (through the features_metric argument). Note that the function returns a patchwork of ggplots.
plot_features_weight_distr(spls_output)
As sPLS performs feature selection, most of the weights are equal to 0.
plot_features_weight_distr(so2pls_output) +
## showing two plots per column (from patchwork package)
plot_layout(ncol = 2)
As sO2PLS performs feature selection for the joint components, most of the weights are equal to 0 for joint component 1.
We will only look at factors 1 to 4 for clarity:
plot_features_weight_distr(
mofa_output,
latent_dimensions = paste("Factor", 1:4)
) +
## showing one plot per column (from patchwork package)
plot_layout(ncol = 1)
We can see for example for factor 1 that only a couple of metabolites have an importance score above 0.5, this small number of features are driving the variation captured by factor 1.
plot_features_weight_distr(
diablo_output
) +
## showing one plot per column (from patchwork package)
plot_layout(ncol = 1)
As DIABLO performs feature selection, most of the weights are equal to zero.
12.5.2 Comparing features importance between latent components
We can also compare the importance score given to the features between any two latent dimensions with the plot_features_weight_pair() function.
plot_features_weight_pair(
spls_output,
paste("Component", 1:2),
mo_data = mo_set_complete,
label_cols = list(
"rnaseq" = "Name",
"metabolome" = "name"
)
)
All genes and most metabolites selected for the two components are different.
plot_features_weight_pair(
so2pls_output,
paste("metabolome specific component", 1:2),
mo_data = mo_set_complete,
label_cols = list(
"rnaseq" = "Name",
"metabolome" = "name"
)
)
The plot shows an interesting pattern: the weights of the compound for these two specific components are negatively correlated, except for a couple of outliers.
plot_features_weight_pair(
mofa_output,
paste("Factor", 1:2),
mo_data = mo_set_complete,
label_cols = list(
"rnaseq" = "Name",
"metabolome" = "name"
)
)
We can see that the genomic markers are mostly given weights of opposite signs for the first two MOFA factors (e.g. positive weight for factor 1 and negative weight for factor 2). For the transcriptomics and metabolomics datasets, the importance score are mostly uncorrelated between the first two factors.
plot_features_weight_pair(
diablo_output,
paste("Component", 1:2),
mo_data = mo_set_complete,
label_cols = list(
"rnaseq" = "Name",
"metabolome" = "name"
)
)
We can see that different features were selected for the first two latent components from each omics dataset.
These plots allow us to quickly check whether different latent dimensions capture different or related aspects of the data. The top 5 features according to their consensus importance score (see Section 14.9) are highlighted in red.
By default, the function shows the signed importance score of the features, i.e. their importance score to which the sign of their weight was added. This can be changed through the features_metric argument of the function.
12.5.3 Plotting top features weight
We can then investigate which features are assigned the highest weights for each latent dimension. These are the features driving the variation captured by the latent components. This is done with the plot_top_features() function.
By default, the feature IDs will be used as labels in the plot, but by passing a MultiDataSet object to the function we can use a column from the features metadata table instead. To do so, we need to pass a named list to the label_cols argument: the names of the elements in the list should correspond to dataset names in the MultiDataSet object, and the value should be the name of the column in the features metadata table of the corresponding dataset to use as feature label. The number of features displayed is controlled through the top_n argument (the default value is 20). Note that if a lower number of features is selected for a certain dataset and latent dimension, only the selected features will be displayed. Again, the function returns a patchwork of ggplots.
plot_top_features(
spls_output,
mo_data = mo_set_complete,
label_cols = list(
"rnaseq" = "Name",
"metabolome" = "name"
)
)
We can see that citric acid and D-mannose have an importance score close to 1 for component 1. Citric acid has a negative weight, meaning that its abundance is negatively correlated with component 1, while D-mannose has a positive weight, so its abundance is positively correlated with component 1. For the second component, only two genes have an importance score above 0.5.
We will focus on joint component 1:
plot_top_features(
so2pls_output,
mo_data = mo_set_complete,
label_cols = list(
"rnaseq" = "Name",
"metabolome" = "name"
),
latent_dimensions = "joint component 1"
)
While all top 20 genes have an importance score above 0.5 for joint component 1, only citric acid is given a high importance score for joint component 1 in the metabolomics dataset.
We will focus on the first two factors here:
plot_top_features(
mofa_output,
mo_data = mo_set_complete,
label_cols = list(
"rnaseq" = "Name",
"metabolome" = "name"
),
latent_dimensions = paste("Factor", 1:2)
)
For factor 1, citric acid is the only metabolite with an importance score above 0.75, while in the genomics dataset four markers have a high importance score, and most of the top 20 genes have an importance score aroung 0.7 or higher.
For clarity, we will focus on the first two latent components:
plot_top_features(
diablo_output,
mo_data = mo_set_complete,
label_cols = list(
"rnaseq" = "Name",
"metabolome" = "name"
),
latent_dimensions = paste("Component", 1:2)
)
For component 1, we can see three markers from the genomics dataset and three genes with a high importance score. From the metabolomics dataset, less than 20 features were retained.
12.5.4 Extracting top features
It is also possible to extract information about the top contributing features as a table with the function get_top_features(). The function offers two options:
Return for each dataset the top N features contributing to the latent dimensions of interest. N is specified through the
n_featuresargument;Return for each dataset features whose importance score is at least equal to a certain threshold. The threshold is specified through the
min_importanceargument.
The first option is preferred when many features contribute to a given factor. The second option has been implemented for cases where some latent dimensions are driven by a small number of features (so for example in the top 10 contributing features, only 3 features would actually have a high importance score). By default, the function returns the top contributing features from all datasets for all latent dimensions, but we can focus on some datasets and/or latent dimensions of interest by passing their names to the datasets and latent_dimensions arguments of the function. Additionally, we can pass to the function a MultiDataSet object (through the mo_data argument); the function will extract information about the features from the features metadata.
We illustrate the two possible options below. For the second option, we add features information to the table:
get_top_features(spls_output, n_features = 3) |>
head()
#> # A tibble: 6 × 5
#> feature_id dataset latent_dimension weight importance
#> <chr> <fct> <fct> <dbl> <dbl>
#> 1 ENSBTAG00000022715 rnaseq Component 1 -0.318 1
#> 2 ENSBTAG00000050618 rnaseq Component 1 -0.302 0.950
#> 3 ENSBTAG00000018223 rnaseq Component 1 -0.291 0.916
#> 4 HMDB00094 metabolome Component 1 0.451 1
#> 5 HMDB00169 metabolome Component 1 -0.327 0.726
#> 6 HMDB00159 metabolome Component 1 -0.312 0.692
get_top_features(
spls_output,
min_importance = 0.8,
mo_data = mo_set_complete
) |>
head()
#> # A tibble: 6 × 31
#> feature_id dataset latent_dimension weight importance chromosome p_value
#> <chr> <fct> <fct> <dbl> <dbl> <chr> <dbl>
#> 1 ENSBTAG0000002… rnaseq Component 1 -0.318 1 26 1.79e-33
#> 2 ENSBTAG0000005… rnaseq Component 1 -0.302 0.950 26 2.54e-33
#> 3 ENSBTAG0000001… rnaseq Component 1 -0.291 0.916 16 7.14e-31
#> 4 ENSBTAG0000004… rnaseq Component 1 -0.266 0.837 15 4.14e-36
#> 5 ENSBTAG0000003… rnaseq Component 1 -0.264 0.831 7 1.22e-35
#> 6 HMDB00094 metabo… Component 1 0.451 1 <NA> 2.48e-41
#> # ℹ 24 more variables: fdr <dbl>, start <int>, end <int>, width <int>,
#> # strand <fct>, Name <chr>, description <chr>, log_fc <dbl>, log_cpm <dbl>,
#> # f <dbl>, de_signif <chr>, de_status <chr>, hmdb_id <chr>, name <chr>,
#> # chemical_formula <chr>, monisotopic_molecular_weight <dbl>,
#> # cas_registry_number <chr>, smiles <chr>, inchikey <chr>, kegg_id <chr>,
#> # direct_parent <chr>, super_class <chr>, t_value <dbl>, padj <dbl>get_top_features(so2pls_output, n_features = 3) |>
head()
#> # A tibble: 6 × 5
#> feature_id dataset latent_dimension weight importance
#> <chr> <fct> <fct> <dbl> <dbl>
#> 1 ENSBTAG00000048835 rnaseq joint component 1 -0.268 1
#> 2 ENSBTAG00000049569 rnaseq joint component 1 -0.266 0.993
#> 3 ENSBTAG00000039037 rnaseq joint component 1 -0.248 0.925
#> 4 HMDB00094 metabolome joint component 1 0.751 1
#> 5 HMDB00042 metabolome joint component 1 0.335 0.445
#> 6 HMDB00357 metabolome joint component 1 0.319 0.424
get_top_features(
so2pls_output,
min_importance = 0.8,
mo_data = mo_set_complete
) |>
head()
#> # A tibble: 6 × 31
#> feature_id dataset latent_dimension weight importance chromosome p_value
#> <chr> <fct> <fct> <dbl> <dbl> <chr> <dbl>
#> 1 ENSBTAG0000004… rnaseq joint component… -0.268 1 24 5.43e-27
#> 2 ENSBTAG0000004… rnaseq joint component… -0.266 0.993 NKLS02001… 1.14e-27
#> 3 ENSBTAG0000003… rnaseq joint component… -0.248 0.925 NKLS02001… 1.45e-28
#> 4 ENSBTAG0000003… rnaseq joint component… -0.235 0.876 7 1.22e-35
#> 5 ENSBTAG0000005… rnaseq joint component… -0.227 0.847 NKLS02000… 1.31e-22
#> 6 ENSBTAG0000002… rnaseq joint component… -0.216 0.805 26 1.79e-33
#> # ℹ 24 more variables: fdr <dbl>, start <int>, end <int>, width <int>,
#> # strand <fct>, Name <chr>, description <chr>, log_fc <dbl>, log_cpm <dbl>,
#> # f <dbl>, de_signif <chr>, de_status <chr>, hmdb_id <chr>, name <chr>,
#> # chemical_formula <chr>, monisotopic_molecular_weight <dbl>,
#> # cas_registry_number <chr>, smiles <chr>, inchikey <chr>, kegg_id <chr>,
#> # direct_parent <chr>, super_class <chr>, t_value <dbl>, padj <dbl>get_top_features(mofa_output, n_features = 3) |>
head()
#> # A tibble: 6 × 5
#> feature_id dataset latent_dimension weight importance
#> <chr> <fct> <fct> <dbl> <dbl>
#> 1 BovineHD0100032240 snps Factor 1 -0.00268 1
#> 2 Hapmap55381-rs29025399 snps Factor 1 -0.00268 1
#> 3 BovineHD0100032254 snps Factor 1 0.00249 0.927
#> 4 ENSBTAG00000049569 rnaseq Factor 1 1.89 1
#> 5 ENSBTAG00000048835 rnaseq Factor 1 1.88 0.999
#> 6 ENSBTAG00000039037 rnaseq Factor 1 1.81 0.961
get_top_features(
mofa_output,
min_importance = 0.8,
mo_data = mo_set_complete
) |>
head()
#> # A tibble: 6 × 40
#> feature_id dataset latent_dimension weight importance chromosome position
#> <chr> <fct> <fct> <dbl> <dbl> <chr> <dbl>
#> 1 BovineHD0100… snps Factor 1 -0.00268 1 1 1.13e8
#> 2 Hapmap55381-… snps Factor 1 -0.00268 1 1 1.13e8
#> 3 BovineHD0100… snps Factor 1 0.00249 0.927 1 1.13e8
#> 4 BovineHD2700… snps Factor 1 -0.00225 0.839 27 3.57e7
#> 5 ENSBTAG00000… rnaseq Factor 1 1.89 1 NKLS02001… NA
#> 6 ENSBTAG00000… rnaseq Factor 1 1.88 0.999 24 NA
#> # ℹ 33 more variables: gen_train_score <dbl>, ref <chr>, alt <chr>,
#> # ilmn_strand <chr>, customer_strand <chr>, norm_id <dbl>, qtl_type <chr>,
#> # qtl_effect <dbl>, p_value <dbl>, fdr <dbl>, start <int>, end <int>,
#> # width <int>, strand <fct>, Name <chr>, description <chr>, log_fc <dbl>,
#> # log_cpm <dbl>, f <dbl>, de_signif <chr>, de_status <chr>, hmdb_id <chr>,
#> # name <chr>, chemical_formula <chr>, monisotopic_molecular_weight <dbl>,
#> # cas_registry_number <chr>, smiles <chr>, inchikey <chr>, kegg_id <chr>, …get_top_features(diablo_output, n_features = 3) |>
head()
#> # A tibble: 6 × 5
#> feature_id dataset latent_dimension weight importance
#> <chr> <fct> <fct> <dbl> <dbl>
#> 1 ARS-BFGL-NGS-27468 snps Component 1 -0.632 1
#> 2 BovineHD2300010006 snps Component 1 0.509 0.807
#> 3 BovineHD0300000351 snps Component 1 0.337 0.534
#> 4 ENSBTAG00000022715 rnaseq Component 1 0.465 1
#> 5 ENSBTAG00000050618 rnaseq Component 1 0.378 0.813
#> 6 ENSBTAG00000049888 rnaseq Component 1 0.305 0.655
get_top_features(
diablo_output,
min_importance = 0.8,
mo_data = mo_set_complete
) |>
head()
#> # A tibble: 6 × 40
#> feature_id dataset latent_dimension weight importance chromosome position
#> <chr> <fct> <fct> <dbl> <dbl> <chr> <dbl>
#> 1 ARS-BFGL-NGS-2… snps Component 1 -0.632 1 1 1.34e8
#> 2 BovineHD230001… snps Component 1 0.509 0.807 23 3.43e7
#> 3 ENSBTAG0000002… rnaseq Component 1 0.465 1 26 NA
#> 4 ENSBTAG0000005… rnaseq Component 1 0.378 0.813 26 NA
#> 5 HMDB00094 metabo… Component 1 -0.716 1 <NA> NA
#> 6 BovineHD070000… snps Component 2 -0.543 1 7 2.11e5
#> # ℹ 33 more variables: gen_train_score <dbl>, ref <chr>, alt <chr>,
#> # ilmn_strand <chr>, customer_strand <chr>, norm_id <dbl>, qtl_type <chr>,
#> # qtl_effect <dbl>, p_value <dbl>, fdr <dbl>, start <int>, end <int>,
#> # width <int>, strand <fct>, Name <chr>, description <chr>, log_fc <dbl>,
#> # log_cpm <dbl>, f <dbl>, de_signif <chr>, de_status <chr>, hmdb_id <chr>,
#> # name <chr>, chemical_formula <chr>, monisotopic_molecular_weight <dbl>,
#> # cas_registry_number <chr>, smiles <chr>, inchikey <chr>, kegg_id <chr>, …12.5.5 Extracting the selected features
For methods that perform feature selection (i.e. sPLS, sO2PLS and DIABLO), we can extract information about the selected features with the get_selected_features() function. As with get_top_features(), if a MultiDataSet object is passed to the function, it will extract information about the features from the features metadata tables and return them. If applied to the results of a method that does not perform feature selection, all features will be returned.
selected_features <- get_selected_features(spls_output)
dim(selected_features)
#> [1] 120 5
head(selected_features)
#> # A tibble: 6 × 5
#> feature_id dataset latent_dimension weight importance
#> <chr> <fct> <fct> <dbl> <dbl>
#> 1 ENSBTAG00000022715 rnaseq Component 1 -0.318 1
#> 2 ENSBTAG00000050618 rnaseq Component 1 -0.302 0.950
#> 3 ENSBTAG00000018223 rnaseq Component 1 -0.291 0.916
#> 4 ENSBTAG00000049888 rnaseq Component 1 -0.266 0.837
#> 5 ENSBTAG00000031647 rnaseq Component 1 -0.264 0.831
#> 6 ENSBTAG00000046158 rnaseq Component 1 -0.248 0.781selected_features <- get_selected_features(so2pls_output)
dim(selected_features)
#> [1] 1394 5
head(selected_features)
#> # A tibble: 6 × 5
#> feature_id dataset latent_dimension weight importance
#> <chr> <fct> <fct> <dbl> <dbl>
#> 1 ENSBTAG00000048835 rnaseq joint component 1 -0.268 1
#> 2 ENSBTAG00000049569 rnaseq joint component 1 -0.266 0.993
#> 3 ENSBTAG00000039037 rnaseq joint component 1 -0.248 0.925
#> 4 ENSBTAG00000031647 rnaseq joint component 1 -0.235 0.876
#> 5 ENSBTAG00000052012 rnaseq joint component 1 -0.227 0.847
#> 6 ENSBTAG00000022715 rnaseq joint component 1 -0.216 0.805selected_features <- get_selected_features(mofa_output)
dim(selected_features)
#> [1] 30735 5
head(selected_features)
#> # A tibble: 6 × 5
#> feature_id dataset latent_dimension weight importance
#> <chr> <fct> <fct> <dbl> <dbl>
#> 1 BovineHD0100032240 snps Factor 1 -0.00268 1
#> 2 Hapmap55381-rs29025399 snps Factor 1 -0.00268 1
#> 3 BovineHD0100032254 snps Factor 1 0.00249 0.927
#> 4 BovineHD2700010051 snps Factor 1 -0.00225 0.839
#> 5 Hapmap44084-BTA-48504 snps Factor 1 -0.00194 0.724
#> 6 BovineHD2100007115 snps Factor 1 0.00190 0.709selected_features <- get_selected_features(diablo_output)
dim(selected_features)
#> [1] 200 5
head(selected_features)
#> # A tibble: 6 × 5
#> feature_id dataset latent_dimension weight importance
#> <chr> <fct> <fct> <dbl> <dbl>
#> 1 ARS-BFGL-NGS-27468 snps Component 1 -0.632 1
#> 2 BovineHD2300010006 snps Component 1 0.509 0.807
#> 3 BovineHD0300000351 snps Component 1 0.337 0.534
#> 4 BovineHD1900011146 snps Component 1 -0.210 0.333
#> 5 BovineHD0900026231 snps Component 1 -0.191 0.302
#> 6 BovineHD1100030384 snps Component 1 -0.181 0.28612.5.6 Features measurements
Finally, we can investigate further some features of interest using the functions shown in Section 4.5. For the example, we will extract the top 3 contributing features from each dataset for the first latent dimension generated with the methods:
spls_top_features <- get_top_features(
spls_output,
n_features = 3,
latent_dimensions = "Component 1"
) |>
pull(feature_id)
spls_top_features
#> [1] "ENSBTAG00000022715" "ENSBTAG00000050618" "ENSBTAG00000018223"
#> [4] "HMDB00094" "HMDB00169" "HMDB00159"so2pls_top_features <- get_top_features(
so2pls_output,
n_features = 3,
latent_dimensions = "joint component 1"
) |>
pull(feature_id)
so2pls_top_features
#> [1] "ENSBTAG00000048835" "ENSBTAG00000049569" "ENSBTAG00000039037"
#> [4] "HMDB00094" "HMDB00042" "HMDB00357"mofa_top_features <- get_top_features(
mofa_output,
n_features = 3,
latent_dimensions = "Factor 1"
) |>
pull(feature_id)
mofa_top_features
#> [1] "BovineHD0100032240" "Hapmap55381-rs29025399" "BovineHD0100032254"
#> [4] "ENSBTAG00000049569" "ENSBTAG00000048835" "ENSBTAG00000039037"
#> [7] "HMDB00094" "HMDB00042" "HMDB00357"diablo_top_features <- get_top_features(
diablo_output,
n_features = 3,
latent_dimensions = "Component 1"
) |>
pull(feature_id)
diablo_top_features
#> [1] "ARS-BFGL-NGS-27468" "BovineHD2300010006" "BovineHD0300000351"
#> [4] "ENSBTAG00000022715" "ENSBTAG00000050618" "ENSBTAG00000049888"
#> [7] "HMDB00094" "HMDB00042" "HMDB00169"We can use the plot_data_heatmap() function to visualise the measurements for the features across the samples as a heatmap. By passing a MultiDataSet object to the function, we can add information about the samples and the features around the heatmap, and change the label used for the features.
plot_data_heatmap(
mo_set_complete,
spls_top_features,
center = TRUE,
scale = TRUE,
show_column_names = FALSE,
only_common_samples = TRUE,
samples_info = c("status", "day_on_feed"),
features_info = c("de_status"),
colours_list = list(
"status" = c("Control" = "gold", "BRD" = "lightblue"),
"day_on_feed" = colorRamp2(c(5, 70), c("white", "pink3")),
"de_status" = c("downregulated" = "deepskyblue",
"upregulated" = "chartreuse3")
),
label_cols = list(
"rnaseq" = "Name",
"metabolome" = "name"
),
truncate = 20
)
plot_data_heatmap(
mo_set_complete,
so2pls_top_features,
center = TRUE,
scale = TRUE,
show_column_names = FALSE,
only_common_samples = TRUE,
samples_info = c("status", "day_on_feed"),
features_info = c("de_status"),
colours_list = list(
"status" = c("Control" = "gold", "BRD" = "lightblue"),
"day_on_feed" = colorRamp2(c(5, 70), c("white", "pink3")),
"de_status" = c("downregulated" = "deepskyblue",
"upregulated" = "chartreuse3")
),
label_cols = list(
"rnaseq" = "Name",
"metabolome" = "name"
),
truncate = 20
)
plot_data_heatmap(
mo_set_complete,
mofa_top_features,
center = TRUE,
scale = TRUE,
show_column_names = FALSE,
only_common_samples = TRUE,
samples_info = c("status", "day_on_feed"),
features_info = c("de_status"),
colours_list = list(
"status" = c("Control" = "gold", "BRD" = "lightblue"),
"day_on_feed" = colorRamp2(c(5, 70), c("white", "pink3")),
"de_status" = c("downregulated" = "deepskyblue",
"upregulated" = "chartreuse3")
),
label_cols = list(
"rnaseq" = "Name",
"metabolome" = "name"
),
truncate = 20
)
plot_data_heatmap(
mo_set_complete,
diablo_top_features,
center = TRUE,
scale = TRUE,
show_column_names = FALSE,
only_common_samples = TRUE,
samples_info = c("status", "day_on_feed"),
features_info = c("de_status"),
colours_list = list(
"status" = c("Control" = "gold", "BRD" = "lightblue"),
"day_on_feed" = colorRamp2(c(5, 70), c("white", "pink3")),
"de_status" = c("downregulated" = "deepskyblue",
"upregulated" = "chartreuse3")
),
label_cols = list(
"rnaseq" = "Name",
"metabolome" = "name"
),
truncate = 20
)
Alternatively, we can display the features’ measurements against some samples covariate, with the plot_data_covariate() function:
plot_data_covariate(
mo_set_complete,
"status",
spls_top_features,
only_common_samples = TRUE,
label_cols = list(
"rnaseq" = "Name",
"metabolome" = "name"
)
)
plot_data_covariate(
mo_set_complete,
"status",
so2pls_top_features,
only_common_samples = TRUE,
label_cols = list(
"rnaseq" = "Name",
"metabolome" = "name"
)
)
plot_data_covariate(
mo_set_complete,
"status",
mofa_top_features,
only_common_samples = TRUE,
label_cols = list(
"rnaseq" = "Name",
"metabolome" = "name"
)
)
plot_data_covariate(
mo_set_complete,
"status",
diablo_top_features,
only_common_samples = TRUE,
label_cols = list(
"rnaseq" = "Name",
"metabolome" = "name"
)
)
12.6 Recap – targets list
For convenience, here is the list of targets that we created in this section:
Targets list for results interpretation
In the targets script, before the targets list, we will add the following code to create the colour palettes used in some of the plots:
## In the preamble of the _targets.R script
## Choosing colour palettes for the samples covariates
palette_status <- scale_colour_brewer(palette = "Set1")
palette_feedlot <- scale_colour_brewer(palette = "Set2")
palette_geno_comp <- scale_colour_brewer(palette = "Dark2")
palette_rnaseq_batch <- scale_colour_brewer(palette = "Paired")list(
## Generating standardised output
tar_target(
spls_output,
get_output(spls_final_run)
),
## Percentage of variance explained
tar_target(
spls_plot_variance_explained,
plot_variance_explained(spls_output)
),
## Samples score matrix plot
tar_target(
spls_samples_scores_plot,
plot_samples_score(
spls_output,
mo_data = mo_set_complete,
colour_upper = "status",
scale_colour_upper = palette_status,
shape_upper = "gender",
colour_lower = "feedlot",
scale_colour_lower = palette_feedlot
) +
theme(legend.box = "vertical")
),
## Distribution of features weight
tar_target(
spls_features_weight_distribution,
plot_features_weight_distr(spls_output)
),
## Plot of top contributing features
tar_target(
spls_top_features_plot,
plot_top_features(
spls_output,
mo_data = mo_set_complete,
label_cols = list(
"rnaseq" = "Name",
"metabolome" = "name"
)
)
),
## Table of selected features
tar_target(
spls_selected_features,
get_selected_features(spls_output)
)
)list(
## Generating standardised output
tar_target(
so2pls_output,
get_output(so2pls_final_run)
),
## Percentage of variance explained
tar_target(
so2pls_plot_variance_explained,
plot_variance_explained(so2pls_output, ncol = 1) +
theme(axis.text.x = element_text(size = 9, angle = 30, hjust = 1))
),
## Samples score matrix plot
tar_target(
so2pls_samples_scores_plot,
plot_samples_score(
so2pls_output,
latent_dimensions = "joint component 1",
mo_data = mo_set_complete,
colour_upper = "status",
scale_colour_upper = palette_status,
shape_upper = "gender"
)
),
## Distribution of features weight
tar_target(
so2pls_features_weight_distribution,
plot_features_weight_distr(so2pls_output) +
plot_layout(ncol = 2)
),
## Plot of top contributing features
tar_target(
so2pls_top_features_plot,
plot_top_features(
so2pls_output,
mo_data = mo_set_complete,
label_cols = list(
"rnaseq" = "Name",
"metabolome" = "name"
)
)
),
## Table of selected features
tar_target(
so2pls_selected_features,
get_selected_features(
so2pls_output,
latent_dimensions = "joint component 1"
)
)
)list(
## Generating standardised output
tar_target(
mofa_output,
get_output(mofa_trained)
),
## Percentage of variance explained
tar_target(
mofa_plot_variance_explained,
plot_variance_explained(mofa_output, ncol = 1) +
theme(axis.text.x = element_text(size = 9, angle = 30, hjust = 1))
),
## Samples score matrix plot
tar_target(
mofa_samples_scores_plot,
plot_samples_score(
mofa_output,
latent_dimensions = paste("Factor", 1:4),
mo_data = mo_set_complete,
colour_upper = "status",
scale_colour_upper = palette_status,
shape_upper = "gender",
colour_lower = "geno_comp_cluster",
scale_colour_lower = palette_geno_comp
) +
theme(legend.box = "vertical")
),
## Distribution of features weight
tar_target(
mofa_features_weight_distribution,
plot_features_weight_distr(
mofa_output,
latent_dimensions = paste("Factor", 1:4)
) +
plot_layout(ncol = 1)
),
## Plot of top contributing features
tar_target(
mofa_top_features_plot,
plot_top_features(
mofa_output,
mo_data = mo_set_complete,
label_cols = list(
"rnaseq" = "Name",
"metabolome" = "name"
),
latent_dimensions = paste("Factor", 1:2)
)
),
## Table of top contributing features
tar_target(
mofa_top_features,
get_top_features(
mofa_output,
min_importance = 0.8,
mo_data = mo_set_complete
)
)
)list(
## Generating standardised output
tar_target(
diablo_output,
get_output(diablo_final_run)
),
## Percentage of variance explained
tar_target(
diablo_plot_variance_explained,
plot_variance_explained(diablo_output, ncol = 2)
),
## Samples score matrix plot
tar_target(
diablo_samples_scores_plot,
plot_samples_score(
diablo_output,
mo_data = mo_set_complete,
colour_upper = "status",
scale_colour_upper = palette_status,
shape_upper = "gender",
colour_lower = "rnaseq_batch",
scale_colour_lower = palette_rnaseq_batch
) +
theme(legend.box = "vertical")
),
## Distribution of features weight
tar_target(
diablo_features_weight_distribution,
plot_features_weight_distr(
diablo_output
) +
plot_layout(ncol = 1)
),
## Plot of top contributing features
tar_target(
diablo_top_features_plot,
plot_top_features(
diablo_output,
mo_data = mo_set_complete,
label_cols = list(
"rnaseq" = "Name",
"metabolome" = "name"
),
latent_dimensions = paste("Component", 1:2)
)
),
## Table of top contributing features
tar_target(
diablo_selected_features,
get_selected_features(diablo_output)
)
)