Non-negative Integration of Gene Expression and Chromatin Accessibility of 10k PBMCs
2023-07-18
Contents
1 Introduction
This vignette demonstrates how to load a trained MuVI model using the MOFA2 R package, and perform a typical downstream analysis for multi-omics data. Briefly, MuVI is a multi-view latent variable model with domain-informed structured sparsity, that integrates (noisy) feature sets e.g. gene set annotations. In this analysis we rely on a pretrained MuVI model according to this python notebook on RNA-seq and ATAC-seq, in order to infer a joint latent space that is informed a priori by cell type gene set annotations.
library(data.table)
library(tibble)
library(ggplot2)
library(rhdf5)
# GSEA
library(msigdbr)
# make sure that `MOFA2` is imported last, to avoid collisions with functions from other packages
library(MOFA2)2 Model Loading
Load the MuVI model with the training data. Essentially, we call the underlying load_model from MOFA2 but rename the latent factors according to the prior information provided when training MuVI. For instance, instead of having Factor1, we have the actual label of the prior gene set informing this particular factor, e.g. Naive T Cell.
load_muvi <- function(file) {
factors_names <- h5read(file, "/factors/group_0")
model <- load_model(file,
sort_factors = FALSE,
remove_inactive_factors = FALSE)
factors_names(model) <- factors_names
return(model)
}
muvi <- load_muvi("10x_multiome_muvi.hdf5")
muvi## Trained MOFA with the following characteristics:
## Number of views: 2
## Views names: rna atac
## Number of features (per view): 3026 14896
## Number of groups: 1
## Groups names: group_0
## Number of samples (per group): 9534
## Number of factors: 243 Variance Decomposition
Variance decomposition provides two main insights.
First, the overall variance explained serves as a check for the quality of the data fit. In this case, the decomposition performed by MuVI explains slightly less than 50% of variance in the RNA and around 30% in the ATAC view.
plot_variance_explained(muvi, plot_total = TRUE)[[2]]
Second, it quantifies the percentage of variance explained by each factor in each view, thereby identifying shared and private factors. For instance, we see that the factor informed by the Naive T Cell gene set explains a large portion of the variance, and is also shared across RNA and ATAC.
plot_variance_explained(muvi)
4 Investigating Factors
Since the latent factors in MuVI are pre-labelled by design, we may directly inspect how they differentiate the cell population. In this case, we see that the latent factor informed by the Naive T Cell gene set assigns higher scores to the 4 groups of T cells. In addition, the Follicular B Cell factor helps distinguish B cells from the rest of the cell population. The same can be said about Nk Cells regarding natural killer cells and Dendritic Cell regarding mDCs and pDCs.
relevant_factors <-
c("Naive T Cell", "Follicular B Cell", "Nk Cells", "Dendritic Cell")
plot_factor(muvi,
factors = relevant_factors,
group_by = "rna.celltype",
color_by = "rna.celltype") +
theme(axis.text.x = element_text(
color = "black",
angle = 90,
vjust = 0.5,
hjust = 1
))
4.1 Factor Loadings
The factor loadings quantify the contribution of each feature (gene) to each latent factor, e.g. pathway.
For the rest of the analysis we focus on the Follicular B Cell factor/pathway.
We first plot the distribution of all factor loadings in the RNA. Notice we only have non-negative weights, as the model has been constrained during training. Among others we notice several B cell markers such as BANK1, IGHM, MS4A1, PAX5 and so on.
plot_weights(
muvi,
view = "rna",
factors = "Follicular B Cell",
nfeatures = 20,
text_size = 4
)
We get a clearer picture of the relevant features for this particular factor by plotting only the top loadings.
plot_top_weights(
muvi,
view = "rna",
factors = "Follicular B Cell",
sign = "positive",
nfeatures = 20,
)
4.2 Associating Factors and Features
Using MOFA2 we may also assess the correlation of the factor scores with respect to the most relevant features.
plot_data_scatter(
muvi,
view = "rna",
factor = "Follicular B Cell",
features = 6,
color_by = "rna.celltype",
add_lm = T,
dot_size = 1
)
4.3 Prior Knowledge Transfer
A powerful aspect of MuVI is pathway information transfer. In particular, we have only informed the RNA, but expect the ATAC view to be implicitly informed via shared factors. For instance, the corresponding gene symbol CD74 and BLK are also deemed highly relevant for the Follicular B Cell pathway, even in the uninformed ATAC view. Potentially, we may even extract a Follicular B Cell ATAC signature from the top loadings.
plot_top_weights(
muvi,
view = "atac",
factors = "Follicular B Cell",
sign = "positive",
nfeatures = 20,
)
4.4 Factor-specific Patterns
As opposed to scatterplots which focus mostly on two dimensions, we may render heatmaps to visualise the overall patterns of each factor. In this case, we observe that B cells are clearly identifiable by significantly larger Follicular B Cell scores.
plot_data_heatmap(
muvi,
view = "rna",
factor = "Follicular B Cell",
features = 25,
show_rownames = T,
show_colnames = F,
cluster_rows = T,
cluster_cols = F,
annotation_samples = "rna.celltype"
)
As we already notice a very strong correlation of the Follicular B Cell with BANK1 we may visualise a stripplot of the factor scores (y-axis) grouped by each cell type (x-axis) and coloured by the BANK1 expression.
plot_factor(muvi,
factors = "Follicular B Cell",
group_by = "rna.celltype",
color_by = "BANK1") +
theme(axis.text.x = element_text(
color = "black",
angle = 30,
vjust = 1,
hjust = 1
))
5 Non-linear Dimensionality Reduction
Similar to the principal components of a PCA decomposition, MuVI infers a compact latent space governed by domain knowledge in terms of feature (gene) sets. Therefore, we may apply clustering algorithms that learn non-linear manifolds such as t-SNE or UMAP. Here, we apply UMAP and further compress the latent space into two dimensions for better visualisation.
factors <- 1:20
muvi <- run_umap(muvi,
factors = factors,
n_neighbors = 15,
min_dist = 0.30)
plot_dimred(
muvi,
method = "UMAP",
color_by = "rna.celltype",
label = TRUE,
stroke = 0.05,
dot_size = 1,
legend = FALSE
)
We can then go over each latent factor, labelled by their corresponding gene set, and colour each sample by the inferred factor score. When comparing to the previous UMAP plot, we clearly see how each latent factor can be strongly identified with their corresponding cell type.
for (i in relevant_factors) {
p <- plot_dimred(
muvi,
method = "UMAP",
color_by = i,
stroke = 0.05,
dot_size = 1
) + ggtitle(i)
print(p)
}
6 Gene Set Enrichment Analysis
Although MuVI accommodates gene set information a priori, it is always advisable to perform a GSEA procedure after training in order to test whether the inferred latent factors (posterior gene sets) still resemble the prior gene sets. In cases where the prior information is not suitable or is severely noisy, the model will attempt to refine the prior gene set annotations given enough evidence from the data.
6.1 Significance to Prior Information
We perform GSEA on the prior gene set collection from the C8 category in MSigDB. The original gene set names contain the prefix HAY_BONE_MARROW_... but have been prettified for the plots, that is HAY_BONE_MARROW_NAIVE_T_CELL is renamed to Naive T Cell, and so on.
As the test results indicate, not only can we ensure that each inferred latent factor significantly matches the prior gene set, we also see that the prior gene set is typically ranked first as the most significant posterior gene set candidate.
msigdb.matrix <- msigdbr(species = "Homo sapiens",
category = "C8") %>% as.data.table %>% .[, id := 1] %>%
dcast(gs_name ~ gene_symbol, value.var = "id", fill = 0) %>%
remove_rownames %>%
column_to_rownames(var = "gs_name")
# GSEA on positive weights only because we have no negative weights
gsea.positive <- run_enrichment(muvi,
feature.sets = msigdb.matrix == 1,
view = "rna",
sign = "positive")
for (i in relevant_factors) {
p <-
plot_enrichment(gsea.positive, factor = i, max.pathways = 15) + ggtitle(i)
print(p)
}
We may also inspect the top loadings that contributed to the test results.
plot_enrichment_detailed(
gsea.positive,
factor = "Follicular B Cell",
max.genes = 10,
max.pathways = 5
)
6.2 Matching to Other Collections
We perform GSEA on another gene set collection, e.g. GO terms, in order to examine whether we can recover similar themes regarding each inferred factor.
Again, the test results point towards the same underlying biological processes, even when comparing to a different gene set collection.
msigdb.matrix <- msigdbr(species = "Homo sapiens",
category = "C5",
subcategory = "BP") %>% as.data.table %>% .[, id := 1] %>%
dcast(gs_name ~ gene_symbol, value.var = "id", fill = 0) %>%
remove_rownames %>%
column_to_rownames(var = "gs_name")
# GSEA on positive weights only because we have no negative weights
gsea.positive <- run_enrichment(muvi,
feature.sets = msigdb.matrix == 1,
view = "rna",
sign = "positive")
for (i in relevant_factors) {
p <-
plot_enrichment(gsea.positive, factor = i, max.pathways = 15) + ggtitle(i)
print(p)
}
plot_enrichment_detailed(
gsea.positive,
factor = "Follicular B Cell",
max.genes = 10,
max.pathways = 5
)
Session Info
sessionInfo()## R version 4.3.1 (2023-06-16)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Ventura 13.4.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Europe/Berlin
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] MOFA2_1.10.0 msigdbr_7.5.1 rhdf5_2.44.0 ggplot2_3.4.2
## [5] tibble_3.2.1 data.table_1.14.8 BiocStyle_2.28.0
##
## loaded via a namespace (and not attached):
## [1] tidyselect_1.2.0 dplyr_1.1.2 farver_2.1.1
## [4] filelock_1.0.2 fastmap_1.1.1 digest_0.6.33
## [7] lifecycle_1.0.3 magrittr_2.0.3 compiler_4.3.1
## [10] rlang_1.1.1 sass_0.4.6 tools_4.3.1
## [13] utf8_1.2.3 yaml_2.3.7 corrplot_0.92
## [16] knitr_1.43 ggsignif_0.6.4 S4Arrays_1.0.4
## [19] labeling_0.4.2 reticulate_1.30 DelayedArray_0.26.6
## [22] plyr_1.8.8 RColorBrewer_1.1-3 abind_1.4-5
## [25] HDF5Array_1.28.1 Rtsne_0.16 babelgene_22.9
## [28] withr_2.5.0 purrr_1.0.1 BiocGenerics_0.46.0
## [31] grid_4.3.1 stats4_4.3.1 fansi_1.0.4
## [34] ggpubr_0.6.0 colorspace_2.1-0 Rhdf5lib_1.22.0
## [37] scales_1.2.1 cli_3.6.1 rmarkdown_2.23
## [40] crayon_1.5.2 generics_0.1.3 rstudioapi_0.15.0
## [43] reshape2_1.4.4 cachem_1.0.8 stringr_1.5.0
## [46] splines_4.3.1 parallel_4.3.1 BiocManager_1.30.21
## [49] matrixStats_1.0.0 basilisk_1.11.2 vctrs_0.6.3
## [52] Matrix_1.6-0 jsonlite_1.8.7 carData_3.0-5
## [55] dir.expiry_1.8.0 bookdown_0.34 car_3.1-2
## [58] IRanges_2.34.1 S4Vectors_0.38.1 ggrepel_0.9.3
## [61] rstatix_0.7.2 irlba_2.3.5.1 magick_2.7.4
## [64] tidyr_1.3.0 jquerylib_0.1.4 glue_1.6.2
## [67] codetools_0.2-19 cowplot_1.1.1 uwot_0.1.16
## [70] RcppAnnoy_0.0.21 stringi_1.7.12 gtable_0.3.3
## [73] munsell_0.5.0 pillar_1.9.0 basilisk.utils_1.12.1
## [76] htmltools_0.5.5 rhdf5filters_1.12.1 R6_2.5.1
## [79] evaluate_0.21 lattice_0.21-8 highr_0.10
## [82] png_0.1-8 backports_1.4.1 pheatmap_1.0.12
## [85] broom_1.0.5 bslib_0.5.0 Rcpp_1.0.11
## [88] nlme_3.1-162 mgcv_1.9-0 xfun_0.39
## [91] MatrixGenerics_1.12.2 forcats_1.0.0 pkgconfig_2.0.3