MEFISTO analysis of the scNMT-seq gastrulation data set
2021-04-22
Contents
1 Description
This vignette demonstrates how MEFISTO can be used to integrate single-cell multi-omics data sets. Here we used the scNMT-seq gastrulation data set, which consists of ~1,500 cells collected across four stages of mouse development, and can be downloaded from here. For roughly half of the cells we profiled RNA expression, DNA methylation and chromatin accessibility from the same cell, and for the other half only RNA expression was profiled. Thus, a significant challenge of this data set is the sparsity of the epigenetic readouts.
You may wonder, how can we use MEFISTO here if we don’t have continuous spatial or temporal covariates? Here we’ll use the RNA expression to derive a multivariate developmental trajectory, which we can then take as a proxy for real time, and use this manifold as a covariate in MEFISTO. This approach wil allow us to identify coordinated patterns of variation between the transcriptome and epigenome during development.
2 Load libraries
library(data.table)
library(purrr)
library(ggplot2)
library(ggpubr)
library(ggrepel)
# MOFA
library(MOFA2)3 Load data
Data is available for download here
The input data consists of three modalities: RNA expression quantified over genes; and DNA methylation and chromatin accessibility quantified over TF motifs. The script to do this step can be found here. Then, it is just a matter of doing some feature selection and preparing the input data for MEFISTO. The script for this step can be found here. For the purpose of this vignette we will load the multi-omics data already processed in a long data.table format
data <- fread("/Users/ricard/data/gastrulation/metaccrna/mefisto/vignette/scnmt_data.txt.gz")
head(data,n=3)## sample feature value view
## 1: E7.5_Plate1_C9 Sox17 0.5903 RNA
## 2: E7.5_Plate1_C9 Atp6v1h -0.3542 RNA
## 3: E7.5_Plate1_C9 Sulf1 -1.4312 RNAPrint number of features per view
data[,length(unique(feature)),by="view"]## view V1
## 1: RNA 887
## 2: motif_met 500
## 3: motif_acc 500Load cell metadata
sample_metadata <- fread("/Users/ricard/data/gastrulation/metaccrna/mefisto/vignette/scnmt_sample_metadata.txt")
colnames(sample_metadata)## [1] "sample" "embryo" "stage" "pass_rnaQC" "pass_metQC"
## [6] "pass_accQC" "lineage10x" "lineage10x_2" "UMAP1" "UMAP2"There are two cell type columns: lineage10x and lineage10x_2. The first one is obtained by mapping the scNMT-seq cells to the RNA gastrulation atlas. The problem with the most this annotation is that we have very few cells per cell type. In the second column these lineages are aggregated into bigger classes
table(sample_metadata$lineage10x)##
## Anterior_Primitive_Streak Caudal_epiblast
## 4 96
## Caudal_mesoderm Def._endoderm
## 5 44
## Epiblast ExE_mesoderm
## 548 1
## Gut Haematoendothelial_progenitors
## 38 5
## Intermediate_mesoderm Mesenchyme
## 23 27
## Mixed_mesoderm Nascent_mesoderm
## 121 222
## Notochord Paraxial_mesoderm
## 49 26
## Pharyngeal_mesoderm Primitive_Streak
## 17 42
## Rostral_neurectoderm Somitic_mesoderm
## 181 34
## Surface_ectoderm Visceral_endoderm
## 9 26table(sample_metadata$lineage10x_2)##
## Ectoderm Endoderm Epiblast Mesoderm
## 188 157 550 481
## Primitive_Streak
## 142For defining the multivariate pseudotime trajectory we applied UMAP on the RNA expression data. Due to the randomness of the algorithm, there is flexibility on how to calculate the latent manifold. For reproducibility purposes we included our precomputed UMAP coordinates in the samples metadata columns UMAP1 and UMAP2
ggscatter(sample_metadata, x="UMAP1", y="UMAP2", color="lineage10x_2") +
scale_color_manual(values=celltype.colors) +
theme_classic() +
ggplot_theme_NoAxes()
4 Create and train the MEFISTO object
Create the MEFISTO object
mefisto <- create_mofa_from_df(data)
mefisto## Untrained MOFA model with the following characteristics:
## Number of views: 3
## Views names: motif_acc motif_met RNA
## Number of features (per view): 500 500 887
## Number of groups: 1
## Groups names: single_group
## Number of samples (per group): 1518
## Add sample metadata to the model
samples_metadata(mefisto) <- sample_metadata %>%
.[sample%in%unlist(samples_names(mefisto))] %>%
setkey(sample) %>% .[unlist(samples_names(mefisto))]Add covariates
mefisto <- set_covariates(mefisto, c("UMAP1","UMAP2"))Define options
model_opts <- get_default_model_options(mefisto)
model_opts$num_factors <- 10Prepare MEFISTO object
mefisto <- prepare_mofa(
object = mefisto,
model_options = model_opts,
)Visualise data structure
plot_data_overview(mefisto)
Train the MEFISTO model (~30min in my MacBook)
mefisto <- run_mofa(mefisto)For reprodudicibility purposes, we also provide the precomputed MEFISTO model that we used for the manuscript:
mefisto <- readRDS("/Users/ricard/data/gastrulation/metaccrna/mefisto/vignette/scnmt_mefisto_model.rds")Select factors that explain at least 1% of variance in the RNA expression
r2 <- get_variance_explained(mefisto, views="RNA")[["r2_per_factor"]][[1]]
factors <- names(which(r2[,"RNA"]>1))
mefisto <- subset_factors(mefisto, factors)
mefisto## Trained MEFISTO with the following characteristics:
## Number of views: 3
## Views names: motif_acc motif_met RNA
## Number of features (per view): 250 250 946
## Number of groups: 1
## Groups names: single_group
## Number of samples (per group): 1518
## Number of covariates per sample: 2
## Number of factors: 75 Downstream analysis
5.1 Variance decomposition
Variance decomposition analysis. Plot the variance explained for each factor in each view
plot_variance_explained(mefisto, x="view", y="factor", max_r2 = 9)
5.2 Analyses of factors
Plot the smothness coefficient for each Factor. We note that Factors 1, 3 and 5 capture smooth patterns across the trajectory, whereas the remaining factors do not.
plot_smoothness(mefisto)
Visualise the Factor values overlayed onto the UMAP. Note that, as predicted from the smooth coefficients, Factors 1 and 3 capture smooth cell fate decision events, whereas Factor 2 is capturing some variation that is not associated with the developmental trajectory:
factors.dt <- get_factors(mefisto, factors=c(1,3,4))[[1]] %>%
as.data.table(keep.rownames = T) %>%
setnames("rn","sample")
to.plot <- sample_metadata[,c("sample","UMAP1","UMAP2")] %>%
merge(factors.dt, by="sample") %>%
melt(id.vars=c("UMAP1","UMAP2","sample"), variable.name="factor") %>%
.[,value:=value/max(abs(value)),by="factor"]
ggscatter(to.plot, x="UMAP1", y="UMAP2", fill="value", shape=21, stroke=0.15) +
facet_wrap(~factor) +
scale_fill_gradient2(low = "gray50", mid="gray90", high = "red") +
theme_classic() +
ggplot_theme_NoAxes()
6 Analyses of weights
Let’s fetch the weights for motif DNA methylation and the weights for motif chromatin accessibility, and remove the prefix to match the feature names. Then we can merge all into a single data.frame for plotting
w.met <- get_weights(mefisto, views="motif_met", factors=c(1,3), as.data.frame=T) %>%
as.data.table %>% .[,feature:=gsub("_met","",feature)] %>%
.[,value:=value/max(abs(value)),by=c("factor")]
w.acc <- get_weights(mefisto, views="motif_acc", factors=c(1,3), as.data.frame=T) %>%
as.data.table %>% .[,feature:=gsub("_acc","",feature)] %>%
.[,value:=value/max(abs(value)),by=c("factor")]
# Merge loadings
w.dt <- merge(
w.met[,c("feature","factor","value")],
w.acc[,c("feature","factor","value")],
by = c("feature","factor")
) %>% .[,feature:=strsplit(feature,"_") %>% map_chr(c(1))]
head(w.dt)## feature factor value.x value.y
## 1: ALX3 Factor1 -0.052179372 0.019081913
## 2: ALX3 Factor3 -0.022018014 0.024370113
## 3: ARGFX Factor1 -0.010555729 0.003068031
## 4: ARGFX Factor3 -0.001140567 0.001484293
## 5: ARNT..HIF1A Factor1 0.002651709 -0.161326074
## 6: ARNT..HIF1A Factor3 -0.040798583 -0.003516151Scatterplots of DNA methylation weights vs chromatin accessibility weights. Let’s highlight the top genes with the highest weights. The weights of these factors reflected the known negative relationship of DNA methylation and chromatin accessibility and identified key TFs associated with this process, including TBX6 and MSGN1 for the mesoderm fate and FOXA2 and HNF1 for the endoderm fate:
for (i in unique(w.dt$factor)) {
to.plot <- w.dt[factor==i]
to.label <- w.dt %>%
.[factor==i] %>%
.[,value:=abs(value.x)+abs(value.y)] %>% setorder(-value) %>% head(n=10)
p <- ggscatter(to.plot, x="value.x", y="value.y", size=1.5, add="reg.line", conf.int=TRUE) +
coord_cartesian(xlim=c(-1,1), ylim=c(-1,1)) +
scale_x_continuous(breaks=c(-1,0,1)) +
scale_y_continuous(breaks=c(-1,0,1)) +
geom_text_repel(data=to.label, aes(x=value.x, y=value.y, label=feature), size=3, max.overlaps=100) +
geom_vline(xintercept=0, linetype="dashed") +
geom_hline(yintercept=0, linetype="dashed") +
stat_cor(method = "pearson") +
labs(x="Methylation weights", y="Accessibility weights")
print(p)
}
7 Imputation
MEFISTO can exploit its generative model to perform data imputation and can use the posterior means of the Gaussian Process to obtain smooth estimates along the covariates.
Interpolate factor values (~1min)
mefisto <- interpolate_factors(mefisto, mefisto@covariates[[1]])
Z_interpol <- t(get_interpolated_factors(mefisto, only_mean = TRUE)[[1]]$mean)Imputation
imputed_data <- list()
for (i in views_names(mefisto)) {
imputed_data[[i]] <- tcrossprod(Z_interpol,get_weights(mefisto,i)[[1]]) %>% t
colnames(imputed_data[[i]]) <- colnames(mefisto@data[[i]][[1]])
}Plot DNA methylation and chromatin accessibility before and after imputation
for (view in c("motif_met","motif_acc")) {
# Define features to plot
# features.to.plot <- features_names(mefisto)[view]
features.to.plot <- list(paste0("MSGN1_412",gsub("motif","",view)))
names(features.to.plot) <- view
# Get original and imputed data
data <- get_data(mefisto, views=view, features=features.to.plot)[[1]][[1]]
data_imputed <- imputed_data[[view]][features.to.plot[[1]],,drop=F]
for (i in features.to.plot[[view]]) {
to.plot <- data.table(
sample = unlist(samples_names(mefisto)),
non_imputed = data[i,],
imputed = data_imputed[i,]
) %>%
setnames(c("sample",sprintf("%s (original)",i),sprintf(" %s (MEFISTO imputed)",i))) %>%
melt(id.vars="sample", value.name="value") %>%
merge(sample_metadata[,c("sample","UMAP1","UMAP2")], by="sample")
# Scale min/max values for visualisation
max.value <- max(to.plot[variable==sprintf(" %s (MEFISTO imputed)",i),value])
min.value <- min(to.plot[variable==sprintf(" %s (MEFISTO imputed)",i),value])
to.plot[value>max.value,value:=max.value]
to.plot[value<min.value,value:=min.value]
p <- ggscatter(to.plot, x="UMAP1", y="UMAP2", color="value") +
facet_wrap(~variable) +
theme_classic() +
ggplot_theme_NoAxes() +
theme(
legend.position = "none",
)
if (view=="motif_met") {
p <- p + scale_color_gradient2(low = "blue", mid="gray90", high = "red")
} else if (view=="motif_acc") {
p <- p + scale_color_gradient2(low = "yellow", mid="gray90", high = "purple")
}
print(p)
}
}
8 Gene set enrichment analysis
Sometimes factors cannot be easily characterised by simply inspecting the genes with the largest weight in each factor. Sometimes it is useful to combine information across genes and work instead with gene sets (i.e. biological pathways).
Here we applied GSEA on the RNA weights, and used this approach to link Factor 4 to cell cycle variation. For more details we refer the reader to our GSEA vignette
8.1 Load GSEA annotation
There are a large number of gene set annotations, and the right one to use will depend on your data set. Some generic and commonly used ones are MSigDB, Reactome and Gene Ontology.
We have manually processed some gene sets, which can be found in the MOFAdata package.
library(MOFAdata)
data("MSigDB_v6.0_C5_mouse") 8.2 Run enrichment analysis
enrichment.out <- run_enrichment(
object = mefisto,
view = "RNA",
feature.sets = MSigDB_v6.0_C5_mouse,
statistical.test = "parametric",
alpha = 0.01
)8.3 Visualise enrichment analysis
plot_enrichment(enrichment.out, factor=4, max.pathways = 15)
Session Info
sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## 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_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] MOFAdata_1.6.0 MOFA2_1.1.19 ggrepel_0.9.1 ggpubr_0.4.0
## [5] ggplot2_3.3.3 purrr_0.3.4 data.table_1.13.6 BiocStyle_2.18.1
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-152 matrixStats_0.58.0 filelock_1.0.2
## [4] RColorBrewer_1.1-2 tools_4.0.3 backports_1.2.1
## [7] R6_2.5.0 HDF5Array_1.18.1 uwot_0.1.10
## [10] DBI_1.1.1 BiocGenerics_0.36.0 mgcv_1.8-33
## [13] colorspace_2.0-0 rhdf5filters_1.2.0 withr_2.4.1
## [16] tidyselect_1.1.0 curl_4.3 compiler_4.0.3
## [19] basilisk.utils_1.2.2 DelayedArray_0.16.1 labeling_0.4.2
## [22] bookdown_0.21 scales_1.1.1 rappdirs_0.3.3
## [25] stringr_1.4.0 digest_0.6.27 foreign_0.8-81
## [28] rmarkdown_2.6 R.utils_2.10.1 basilisk_1.2.1
## [31] rio_0.5.16 pkgconfig_2.0.3 htmltools_0.5.1.1
## [34] MatrixGenerics_1.2.1 highr_0.8 rlang_0.4.10
## [37] readxl_1.3.1 generics_0.1.0 farver_2.0.3
## [40] jsonlite_1.7.2 dplyr_1.0.4 zip_2.1.1
## [43] car_3.0-10 R.oo_1.24.0 magrittr_2.0.1
## [46] Matrix_1.3-2 Rcpp_1.0.6 munsell_0.5.0
## [49] S4Vectors_0.28.1 Rhdf5lib_1.12.1 abind_1.4-5
## [52] reticulate_1.18 lifecycle_0.2.0 R.methodsS3_1.8.1
## [55] stringi_1.5.3 yaml_2.2.1 carData_3.0-4
## [58] rhdf5_2.34.0 Rtsne_0.15 plyr_1.8.6
## [61] grid_4.0.3 parallel_4.0.3 forcats_0.5.1
## [64] crayon_1.4.0 lattice_0.20-41 haven_2.3.1
## [67] cowplot_1.1.1 splines_4.0.3 hms_1.0.0
## [70] magick_2.6.0 knitr_1.31 pillar_1.4.7
## [73] ggsignif_0.6.0 reshape2_1.4.4 stats4_4.0.3
## [76] glue_1.4.2 evaluate_0.14 BiocManager_1.30.10
## [79] vctrs_0.3.6 cellranger_1.1.0 gtable_0.3.0
## [82] tidyr_1.1.2 assertthat_0.2.1 xfun_0.20
## [85] openxlsx_4.2.3 broom_0.7.4 rstatix_0.6.0
## [88] tibble_3.0.6 pheatmap_1.0.12 IRanges_2.24.1
## [91] corrplot_0.84 ellipsis_0.3.1