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  RNA

Print number of features per view

data[,length(unique(feature)),by="view"]
##         view  V1
## 1:       RNA 887
## 2: motif_met 500
## 3: motif_acc 500

Load 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                             26
table(sample_metadata$lineage10x_2)
## 
##         Ectoderm         Endoderm         Epiblast         Mesoderm 
##              188              157              550              481 
## Primitive_Streak 
##              142

For 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 <- 10

Prepare 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: 7

5 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.003516151

Scatterplots 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