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.
library(data.table)
library(purrr)
library(ggplot2)
library(ggpubr)
library(ggrepel)
# MOFA
library(MOFA2)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()