MOFA+: integration of a time-course single-cell RNA-seq dataset
2020-09-03
Contents
1 Description
This vignette demonstrates the multi-group integration framework of MOFA+ on a single data modality.
We consider a data set of scRNA-seq experiments where 16,152 cells were isolated from a total of 8 mouse embryos from developmental stages E6.5, E6.75, E7.0 and E7.25 (two embryos per stage), spanning post-implantation and early gastrulation.
Despite differences in developmental time, embryos are expected to contain similar subpopulations of cells. Hence, MOFA+ should detect the existence of biological sources of variation that are shared between groups.
The data set we use here is a simplified subset version of the original data set, which can be visualised and downloaded from here.
2 How does the multi-group inference work?
Intuitively, this extension breaks the assumption of independent samples and allows inference across multiple groups, where groups are predefined sets of samples (i.e. different conditions, batches, cohorts, etc.).
Importantly, the model is not focused on capturing the differential changes between the groups (as for example when doing differential expression). The aim of the multi-group framework is to find out the sources of variability within each group and to charactersie which ones are shared between the different groups from those that are exclusive to a single group.
This is a rather advanced option that we disencourage if this is the first time that you are using MOFA. For more questions, please read the FAQ section.
3 Load libraries and data
3.1 Libraries
Make sure that MOFA2 is imported last, to avoid collisions with functions from other packages
library(Seurat)
library(ggplot2)
library(MOFA2)3.2 Data
The RNA expression data has been processed using Seurat. It has already been normalised and subset to the top 5,000 most variable genes (after regressing out the group effect).
load(url("ftp://ftp.ebi.ac.uk/pub/databases/mofa/scrna_gastrulation/gastrulation10x_seurat.RData"))
seurat## An object of class Seurat
## 5000 features across 16152 samples within 1 assay
## Active assay: RNA (5000 features, 5000 variable features)Define cell type colors for the visualisations
colors <- c(
"Epiblast" = "grey70",
"Primitive Streak" = "sandybrown",
"Mesoderm" = "violetred",
"ExE endoderm" = "#548B54",
"ExE ectoderm" = "black"
)4 Create MOFA object
Group cells according to the embryo and the stage they comne from
seurat@meta.data$stage_sample <- paste(seurat@meta.data$stage,seurat@meta.data$sample, sep="_")
unique(seurat@meta.data$stage_sample)## [1] "E6.5_1" "E6.5_5" "E7.0_10" "E7.0_14" "E7.25_26" "E7.25_27"MOFAobject <- create_mofa(seurat, groups = "stage_sample")4.1 Plot data overview
plot_data_overview(MOFAobject)
4.2 Define MOFA options
# Default data options
data_opts <- get_default_data_options(MOFAobject)
# Default model options
model_opts <- get_default_model_options(MOFAobject)
# Training options
train_opts <- get_default_training_options(MOFAobject)
train_opts$convergence_mode <- "fast"
train_opts$seed <- 424.3 Prepare the MOFA object
MOFAobject <- prepare_mofa(
object = MOFAobject,
data_options = data_opts,
model_options = model_opts,
training_options = train_opts
)5 Train the MOFA model
This step can take quite some time (~1h with standard CPU inference), so we provide a pre-trained model in the next chunk of code.
Note: if you train the model from scratch, the results will not be 100% reproducible with the pre-trained model (but hopefully similar enough!).
MOFAobject <- run_mofa(MOFAobject)5.1 Load pre-computed model
MOFA models are saved in hdf5 format and can be loaded into R with the function load_model. In this case, however, we provide the trained model as an RData file
# MOFAobject <- load_model(outfile)
load(url("ftp://ftp.ebi.ac.uk/pub/databases/mofa/scrna_gastrulation/gastrulation10x_mofa.RData"))6 Overview of the trained MOFA model
The MOFAobject consists of multiple slots where relevant data and information is stored. For descriptions, you can read the documentation by ?MOFA
slotNames(MOFAobject)## [1] "data" "intercepts" "imputed_data"
## [4] "samples_metadata" "features_metadata" "expectations"
## [7] "training_stats" "training_options" "stochastic_options"
## [10] "data_options" "model_options" "dimensions"
## [13] "on_disk" "dim_red" "cache"
## [16] "status"dimensions contains the dimensions of the model (K for factors, M for views, G for groups, N for samples (per group), D for features (per view)):
MOFAobject@dimensions## $M
## [1] 1
##
## $G
## [1] 6
##
## $N
## E6.5 (1) E6.5 (2) E7.0 (1) E7.0 (2) E7.25 (1) E7.25 (2)
## 337 1106 2152 1135 5537 5885
##
## $D
## RNA
## 5000
##
## $K
## [1] 10data contains the centered (and optionally scaled) input data matrices. It is a nested list of matrices where the first index refers to the view, and the second index to the group. In this case we have one view (RNA) and six groups
names(MOFAobject@data)## [1] "RNA"names(MOFAobject@data[[1]])## [1] "E6.5 (1)" "E6.5 (2)" "E7.0 (1)" "E7.0 (2)" "E7.25 (1)" "E7.25 (2)"samples_metadata contains the cell metadata that I have previously added to the MOFA object. Columns are:
- stage: developmental stage.
- lineage: cell type annotation.
- group: the grouping used for MOFA. In this case each embryo corresponds to a separate group.
head(MOFAobject@samples_metadata)## sample stage lineage group
## 1 cell_1 E6.5 Epiblast E6.5 (1)
## 2 cell_2 E6.5 Primitive Streak E6.5 (1)
## 3 cell_5 E6.5 ExE ectoderm E6.5 (1)
## 4 cell_6 E6.5 Epiblast E6.5 (1)
## 5 cell_8 E6.5 Epiblast E6.5 (1)
## 6 cell_9 E6.5 Epiblast E6.5 (1)6.1 Overview of training data
The function plot_data_overview can be used to obtain an overview of the input data.
It shows how many views (rows) and how many groups (columns) exist, what are their corresponding dimensionalities and how many missing information they have (grey bars).
In this case we have one view (RNA expression, a total of 5,000 genes) and 6 groups that correspond to different embryos at different stages of development, for a total of 16,152 cells.
plot_data_overview(MOFAobject, colors = c("RNA"="darkgreen"))
7 Plot variance explained
Quantifying the variance explained per factor across groups and views is probably the most important plot that MOFA+ generates. It summarises the (latent) signal from a complex heterogeneous data set in a single figure.
7.1 Plot variance explained per factor across groups
plot_variance_explained(MOFAobject, x="group", y="factor")
There is a lot of information contained in this plot. Factor 1 and Factor 2 explain a lot of variance across multiple groups. In contrast, Factor 4 increases its activity from E6.5 to E7.5, indicating that it captures a source of variation that emerges at E6.5.
We can also plot the total variance explained per group (with all factors) by adding the argument plot_total = TRUE. Notably, only 10 factors are sufficient to capture between 35% and 55% of the transcriptional variance per embryo
7.2 Plot total variance explained per group
plot_variance_explained(
MOFAobject,
x = "group",
y = "factor",
plot_total = T
)[[2]]
7.3 Plot variance explained for individual features
We can also inspect the variance explained by the MOFA factors for individual features. A high \(R^2\) implies that the MOFA factors captures most of the variation for this feature, whereas small values means that the variation for this feature is not explained by the model (i.e. it is considered as noise):
Variance explained by all factors, in each group:
features <- c("Rbp4","Ttr","Spink1","Mesp1","E130311K13Rik","Hey1")
plot_variance_explained_per_feature(
MOFAobject,
view = "RNA",
features = features
)
8 Characterisation of Factor 1
8.1 Visualisation of factor values
Each factor ordinates cells along a one-dimensional axis that is centered at zero. Samples with different signs indicate opposite phenotypes, with higher absolute value indicating a stronger phenotype. For example, if the \(k\)-th factor captures the variability associated with cell cycle, we could expect cells in Mitosis to be at one end of the factor (irrespective of the sign, only the relative positioning being of importance). In contrast, cells in G1 phase are expected to be at the other end of the factor. Cells with intermediate phenotype, or with no clear phenotype (i.e. no cell cycle genes profiled), are expected to be located around zero.
Let’s plot Factor 1 values and we color cells by lineage assignment. Clearly, this factors captures the emergence of ExE endoderm.
plot_factor(MOFAobject,
factor = 1,
color_by = "lineage" # lineage is a column in MOFAobject@samples.metadata
) + scale_fill_manual(values=colors)
Here are other ways of representing the same plot:
p <- plot_factor(MOFAobject,
factor = 1,
color_by = "lineage",
dodge = T, # dodge points with different colors
legend = F, # remove legend
add_violin = T # add violin plots
)
p <- p +
scale_fill_manual(values=colors) +
pOne can also change the default groups in the x-axis by some manually defined grouping structure. For example a column in the sample metadata
plot_factor(MOFAobject,
factor = 1,
color_by = "lineage",
group_by = "stage",
legend = F
) + scale_fill_manual(values=colors)
Combinations of factors can be plotted with plot_factors:
plot_factors(MOFAobject,
factors = c(1,4),
color_by = "lineage"
) + scale_fill_manual(values=colors)
8.2 Visualisation of weights
The weights provide a score for each gene on each factor. Genes with no association with the factor are expected to have values close to zero, whereas genes with strong association with the factor are expected to have large absolute values. The sign of the weight indicates the direction of the effect: a positive weight indicates that the feature is more active in the cells with positive factor values, and viceversa.
Let’s plot the distribution of weights for Factor 1.
plot_weights(MOFAobject,
view = "RNA",
factor = 1,
nfeatures = 10, # Top number of features to highlight
scale = T # Scale weights from -1 to 1
)
If you are not interested in the directionality of the effect, you can take the absolute value of the weights (abs=TRUE). You can also highlight some genes of interest using the argument manual:
plot_weights(MOFAobject,
view = "RNA",
factor = 1,
nfeatures = 5,
manual = list(c("Snai1","Mesp1","Phlda2"), c("Rhox5","Elf5")),
color_manual = c("darkgreen","red"),
scale = T,
abs = T
)
If you are not interested in the full distribution, but just on the top weights, you can do:
plot_top_weights(MOFAobject,
view = "RNA",
factor = 1,
nfeatures = 10,
scale = T,
abs = T
)
8.3 Visualisation of gene expression patterns given by the Factors
We expect negative weights For Factor 1 to be marker genes of ExE Endoderm. If we plot Factor 1 colouring cells by mRNA expresion of the genes with the largest weights:
genes <- c("Ttr","Apom")
for (i in genes) {
p <- plot_factor(MOFAobject,
factor = 1,
color_by = i
) + scale_colour_gradientn(colours = terrain.colors(10))
print(p)
}
The weights are useful to get an idea of which are top genes that drive the factors. However, to get an idea of how well Factors are associated with genomic features we can generate a scatter plot of the Factor values against mRNA expression for the genes with the highest weights:
Top genes with positive weights
p <- plot_data_scatter(MOFAobject,
view = "RNA",
factor = 1,
features = 6, # Number of features to show
sign = "positive", # select top 6 features with positive weights
color_by = "lineage", # color cells by lineage
add_lm = T, # add linear regression estimates
lm_per_group = F,
dot_size = 2
)
p <- p +
scale_fill_manual(values=colors) +
theme(legend.position = "none")
print(p)
Top genes with negative weights
p <- plot_data_scatter(MOFAobject,
view = "RNA",
factor = 1,
features = 6, # Number of features to show
sign = "negative", # select top 6 features with positive weights
color_by = "lineage", # color cells by lineage
add_lm = T, # add linear regression estimates
lm_per_group = F,
dot_size = 2
)
p <- p +
scale_fill_manual(values=colors) +
theme(legend.position = "none")
print(p)
9 Non-linear dimensionality reduction from the MOFA factors
The latent space inferred by MOFA can be employed as input to other single-cell algorithms that learn non-linear manifolds such as UMAP or t-SNE. This can be very useful to identify cellular populations and reconstruct complex pseudotime trajectories.
Run t-SNE
set.seed(42)
MOFAobject <- run_tsne(MOFAobject)In this data set, we see that the combination of MOFA factors have enough information to discriminate all cell types:
plot_dimred(MOFAobject,
method = "TSNE",
color_by = "lineage"
) + scale_fill_manual(values=colors)
10 sessionInfo
sessionInfo()## R version 4.0.0 (2020-04-24)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] MOFA2_0.99.0 ggplot2_3.3.2 Seurat_3.2.0 BiocStyle_2.16.0
##
## loaded via a namespace (and not attached):
## [1] Rtsne_0.15 ggbeeswarm_0.6.0 colorspace_1.4-1
## [4] ggsignif_0.6.0 deldir_0.1-25 rio_0.5.16
## [7] ellipsis_0.3.1 ggridges_0.5.2 spatstat.data_1.4-3
## [10] ggpubr_0.4.0 farver_2.0.3 leiden_0.3.3
## [13] listenv_0.8.0 ggrepel_0.8.2 codetools_0.2-16
## [16] splines_4.0.0 knitr_1.29 polyclip_1.10-0
## [19] jsonlite_1.7.0 broom_0.5.6 ica_1.0-2
## [22] cluster_2.1.0 png_0.1-7 pheatmap_1.0.12
## [25] uwot_0.1.8 shiny_1.5.0 HDF5Array_1.16.1
## [28] sctransform_0.2.1 BiocManager_1.30.10 compiler_4.0.0
## [31] httr_1.4.1 backports_1.1.8 Matrix_1.2-18
## [34] fastmap_1.0.1 lazyeval_0.2.2 later_1.1.0.1
## [37] htmltools_0.5.0 tools_4.0.0 rsvd_1.0.3
## [40] igraph_1.2.5 gtable_0.3.0 glue_1.4.1
## [43] RANN_2.6.1 reshape2_1.4.4 dplyr_1.0.0
## [46] Rcpp_1.0.4.6 spatstat_1.64-1 carData_3.0-4
## [49] cellranger_1.1.0 vctrs_0.3.1 ape_5.4
## [52] nlme_3.1-148 lmtest_0.9-37 xfun_0.15
## [55] stringr_1.4.0 globals_0.12.5 openxlsx_4.1.5
## [58] mime_0.9 miniUI_0.1.1.1 lifecycle_0.2.0
## [61] irlba_2.3.3 rstatix_0.6.0 goftest_1.2-2
## [64] future_1.17.0 MASS_7.3-51.6 zoo_1.8-8
## [67] scales_1.1.1 hms_0.5.3 promises_1.1.1
## [70] spatstat.utils_1.17-0 parallel_4.0.0 rhdf5_2.32.1
## [73] RColorBrewer_1.1-2 curl_4.3 yaml_2.2.1
## [76] reticulate_1.16 pbapply_1.4-2 gridExtra_2.3
## [79] ggrastr_0.1.9 rpart_4.1-15 stringi_1.4.6
## [82] S4Vectors_0.26.1 corrplot_0.84 BiocGenerics_0.34.0
## [85] zip_2.0.4 rlang_0.4.6 pkgconfig_2.0.3
## [88] matrixStats_0.56.0 evaluate_0.14 lattice_0.20-41
## [91] Rhdf5lib_1.10.0 ROCR_1.0-11 purrr_0.3.4
## [94] tensor_1.5 labeling_0.3 patchwork_1.0.1
## [97] htmlwidgets_1.5.1 cowplot_1.0.0 tidyselect_1.1.0
## [100] RcppAnnoy_0.0.16 plyr_1.8.6 magrittr_1.5
## [103] bookdown_0.20 R6_2.4.1 magick_2.4.0
## [106] IRanges_2.22.2 generics_0.0.2 DelayedArray_0.14.0
## [109] haven_2.3.1 foreign_0.8-80 pillar_1.4.4
## [112] withr_2.2.0 mgcv_1.8-31 fitdistrplus_1.1-1
## [115] survival_3.2-3 abind_1.4-5 tibble_3.0.1
## [118] future.apply_1.5.0 car_3.0-8 crayon_1.3.4
## [121] KernSmooth_2.23-17 plotly_4.9.2.1 rmarkdown_2.3
## [124] readxl_1.3.1 grid_4.0.0 data.table_1.12.8
## [127] forcats_0.5.0 digest_0.6.25 xtable_1.8-4
## [130] tidyr_1.1.0 httpuv_1.5.4 stats4_4.0.0
## [133] munsell_0.5.0 beeswarm_0.2.3 viridisLite_0.3.0
## [136] vipor_0.4.5