Demo 1: Multi-OMICs Factor Analysis (MOFA) and Unsupervised Data Learning
Nikolay Oskolkov (modified and updated by Quentin Read)
2025-06-04
Abstract
 Here we will discuss unsupervised manifold learning and introduce Multi-OMICs Factor Analysis (MOFA) as an example of unsupervised data exploration. Further, we will perform integrative OMICs analysis of chromatin accessibility, DNA methylation and transcriptome profiling data from the scNMT-seq (single-cell nucleosome, methylation and transcription sequencing) using MOFA.
Multi-OMICs Factor Analysis (MOFA)
See the original version of this notebook from the 2021 Sweden course for background information. (Schedule page with all course materials from 2021 Sweden course)
Prepare scNMT Data Set for MOFA
In this section we will apply MOFA to the single cell multi-OMICs data set scNMT which profiles chromatine accessebility (scATACseq), DNA methylation (scBSseq) and gene expression (scRNAseq) information from the same biological cells belonging to two types: differentiating mouse embryonic stem cells (ESCs) and embrionic bodies (EBs).
Load R packages.
library(mixOmics)## Loading required package: MASS## Loading required package: lattice## Loading required package: ggplot2##
## Loaded mixOmics 6.32.0
## Thank you for using mixOmics!
## Tutorials: http://mixomics.org
## Bookdown vignette: https://mixomicsteam.github.io/Bookdown
## Questions, issues: Follow the prompts at http://mixomics.org/contact-us
## Cite us: citation('mixOmics')library(MOFA2)##
## Attaching package: 'MOFA2'## The following object is masked from 'package:stats':
##
## predictWe will start with reading and having a look at the individual OMICs data sets:
Read single-cell RNAseq data.
Note: If you are running this on the demo server the data are located in a directory called
dataand you can run the first of the following two code chunks. If you are running this on your own machine, you may import the data directly from a Web URL, as in the second of the two code chunks.
scRNAseq<-read.delim("data/scRNAseq.txt", header = TRUE, check.names = FALSE, row.names = 1, sep="\t")
ens2genesymbol<-read.delim("data/ENSEMBLE_TO_GENE_SYMBOL_MOUSE.txt")scRNAseq<-read.delim("https://usda-ree-ars.github.io/SEAStatsData/omics/scRNAseq.txt", header = TRUE, check.names = FALSE, row.names = 1, sep="\t")
ens2genesymbol<-read.delim("https://usda-ree-ars.github.io/SEAStatsData/omics/ENSEMBLE_TO_GENE_SYMBOL_MOUSE.txt")Use a lookup table to match the IDs of the dataset with gene names. Display the first 5 rows and columns of the dataset.
ens2genesymbol<-ens2genesymbol[match(colnames(scRNAseq),as.character(ens2genesymbol$ensembl_gene_id)),]
colnames(scRNAseq)<-ens2genesymbol$external_gene_name
scRNAseq<-as.data.frame(t(scRNAseq))
scRNAseq[1:5,1:5]## ESC_A02 ESC_A03 ESC_A04 ESC_A05 ESC_A06
## Mrpl15 557.4201 130.8536 305.7356 417.1656 573.6504
## Lypla1 410.5958 138.8325 186.3237 106.0895 662.7008
## Tcea1 437.9119 356.3899 276.9121 161.8315 213.3068
## Atp6v1h 345.7199 500.5417 195.5884 133.0614 325.1376
## Oprk1 0.0000 0.0000 0.0000 0.0000 0.0000Read two additional omics datasets: single-cell BSseq (bisulfite) and
ATACseq. As above you may import them from the local folder
data if you are on the demo server, or from the Web
location if you are not.
scBSseq<-read.delim("data/scBSseq.txt", header = TRUE, check.names = FALSE, row.names = 1, sep="\t")
scATACseq<-read.delim("data/scATACseq.txt", header = TRUE, check.names = FALSE, row.names = 1, sep="\t")scBSseq<-read.delim("https://usda-ree-ars.github.io/SEAStatsData/omics/scBSseq.txt", header = TRUE, check.names = FALSE, row.names = 1, sep="\t")
scATACseq<-read.delim("https://usda-ree-ars.github.io/SEAStatsData/omics/scATACseq.txt", header = TRUE, check.names = FALSE, row.names = 1, sep="\t")Transpose the BSseq and ATACseq datasets so they are in the same format as the RNAseq. Display the first 5 rows and columns of each dataset.
scBSseq<-as.data.frame(t(scBSseq))
scBSseq[1:5,1:5]## ESC_A02 ESC_A03 ESC_A04 ESC_A05 ESC_A06
## 1_101606896 0.00000 0 0 0 0
## 1_101935054 88.88889 100 100 100 30
## 1_101935061 88.88889 100 100 100 90
## 1_102002184 100.00000 100 0 0 0
## 1_102238204 100.00000 100 100 100 100scATACseq<-as.data.frame(t(scATACseq))
scATACseq[1:5,1:5]## ESC_A02 ESC_A03 ESC_A04 ESC_A05 ESC_A06
## 1_100392139 100 100 100 100 67
## 1_100392151 100 100 100 50 67
## 1_100668590 50 86 44 45 28
## 1_100738967 83 88 83 83 88
## 1_100994324 0 0 11 75 0For the scRNAseq OMIC layer we will only select highly expressed genes (mean >= 1) in order to remove noisy features that might contaminate the further downstream analysis. We will also perform log-transform of the data which can be seen as a mild normalization. Look at the first 5 rows and columns of the updated dataset, and check its size.
scRNAseq<-scRNAseq[rowMeans(scRNAseq)>=1,]
scRNAseq<-log10(scRNAseq + 1)
scRNAseq[1:5,1:5]## ESC_A02 ESC_A03 ESC_A04 ESC_A05 ESC_A06
## Mrpl15 2.746961 2.120092 2.486764 2.621348 2.759404
## Lypla1 2.614471 2.145608 2.272593 2.029747 2.821972
## Tcea1 2.642377 2.553142 2.443907 2.211738 2.331036
## Atp6v1h 2.539979 2.700307 2.293558 2.127304 2.513401
## Oprk1 0.000000 0.000000 0.000000 0.000000 0.000000dim(scRNAseq)## [1] 12145 113Since scBSseq and scATACseq OMIC layers should be almost binary
(methylated vs unmethylated for scBSseq and open vs. close for
scATACseq) data sets, a good way to get rid of redundancy in the scBSseq
and scATACseq data (and thus perform feature pre-selection and reduce
dimensions) is to select only features with high variability. The
function nearZeroVar() from the mixOmics
package is used for this.
my_nearZeroVar<-nearZeroVar(as.data.frame(t(scBSseq)))
head(my_nearZeroVar$Metrics)## freqRatio percentUnique
## 1_101606896 112.00 1.7699115
## 1_105525627 35.00 5.3097345
## 1_105539120 26.25 5.3097345
## 1_10605572 0.00 0.8849558
## 1_107485472 0.00 0.8849558
## 1_109269874 20.00 7.9646018dim(my_nearZeroVar$Metrics)## [1] 3290 2scBSseq <- scBSseq[-which(rownames(scBSseq)%in%rownames(my_nearZeroVar$Metrics)),]
scBSseq[1:5,1:5]## ESC_A02 ESC_A03 ESC_A04 ESC_A05 ESC_A06
## 1_101935054 88.88889 100 100 100 30
## 1_101935061 88.88889 100 100 100 90
## 1_102002184 100.00000 100 0 0 0
## 1_102238204 100.00000 100 100 100 100
## 1_102255678 100.00000 100 100 100 0dim(scBSseq)## [1] 5285 113my_nearZeroVar<-nearZeroVar(as.data.frame(t(scATACseq)), uniqueCut = 1)
head(my_nearZeroVar$Metrics)## freqRatio percentUnique
## 1_10401731 0 0.8849558
## 1_128314944 0 0.8849558
## 1_13628215 0 0.8849558
## 1_178804829 0 0.8849558
## 1_183332775 0 0.8849558
## 1_183944178 0 0.8849558dim(my_nearZeroVar$Metrics)## [1] 91 2scATACseq <- scATACseq[-which(rownames(scATACseq)%in%rownames(my_nearZeroVar$Metrics)),]
scATACseq[1:5,1:5]## ESC_A02 ESC_A03 ESC_A04 ESC_A05 ESC_A06
## 1_100392139 100 100 100 100 67
## 1_100392151 100 100 100 50 67
## 1_100668590 50 86 44 45 28
## 1_100738967 83 88 83 83 88
## 1_100994324 0 0 11 75 0dim(scATACseq)## [1] 11709 113Let us now have a look at the histograms of individual OMICs layers in order to decide what distribution they follow and how we should model these distributions with MOFA:
hist(rowMeans(scRNAseq),breaks = 100, main = "scRNAseq")
hist(rowMeans(scBSseq), breaks = 100, main = "scBSseq")
hist(rowMeans(scATACseq), breaks = 100, main = "scATACseq")
We conclude that while scRNAseq data looks fairly Gaussian (or at least exponential), we should probably model the scBSseq and scATACseq data following Bernoulli distribution. They look quite bimodal, indicating the binary nature of the data, i.e. methylated vs. unmethylated for scBSseq and open vs. close for scATACseq. To make the scBSseq and scATACseq data purely Bernoulli-like, we will make them binary by encoding values below 50 as 0 and above 50 as 1. Since binary data typically have vey low variation compared to continuous data, we need to remove low-variance features in this case again:
scBSseq<-ifelse(scBSseq<50,0,1)
scBSseq[1:5,1:5]## ESC_A02 ESC_A03 ESC_A04 ESC_A05 ESC_A06
## 1_101935054 1 1 1 1 0
## 1_101935061 1 1 1 1 1
## 1_102002184 1 1 0 0 0
## 1_102238204 1 1 1 1 1
## 1_102255678 1 1 1 1 0my_nearZeroVar<-nearZeroVar(as.data.frame(t(scBSseq)))
head(my_nearZeroVar$Metrics)## freqRatio percentUnique
## 1_110547534 21.60 1.7699115
## 1_133172191 27.25 1.7699115
## 1_139117169 27.25 1.7699115
## 1_146093078 21.60 1.7699115
## 1_194260338 112.00 1.7699115
## 1_24611678 0.00 0.8849558dim(my_nearZeroVar$Metrics)## [1] 761 2scBSseq <- as.data.frame(scBSseq[-which(rownames(scBSseq)%in%rownames(my_nearZeroVar$Metrics)),])
scBSseq[1:5,1:5]## ESC_A02 ESC_A03 ESC_A04 ESC_A05 ESC_A06
## 1_101935054 1 1 1 1 0
## 1_101935061 1 1 1 1 1
## 1_102002184 1 1 0 0 0
## 1_102238204 1 1 1 1 1
## 1_102255678 1 1 1 1 0dim(scBSseq)## [1] 4524 113scATACseq<-ifelse(scATACseq<50,0,1)
my_nearZeroVar<-nearZeroVar(as.data.frame(t(scATACseq)))
head(my_nearZeroVar$Metrics)## freqRatio percentUnique
## 1_100392139 112.00 1.769912
## 1_100392151 27.25 1.769912
## 1_100738967 55.50 1.769912
## 1_101897864 55.50 1.769912
## 1_102283335 55.50 1.769912
## 1_102563492 112.00 1.769912dim(my_nearZeroVar$Metrics)## [1] 2238 2scATACseq <- as.data.frame(scATACseq[-which(rownames(scATACseq)%in%rownames(my_nearZeroVar$Metrics)),])
scATACseq[1:5,1:5]## ESC_A02 ESC_A03 ESC_A04 ESC_A05 ESC_A06
## 1_100668590 1 1 0 0 0
## 1_100994324 0 0 0 1 0
## 1_100994328 0 0 0 0 0
## 1_100994336 0 0 0 0 0
## 1_100994343 0 0 0 0 1dim(scATACseq)## [1] 9471 113Now the data cleaning step is finished and the OMICs are prepared to be integrated in an unsupervised way with MOFA.
Run MOFA on scNMT Data Set
Let us continue by creating a MOFA object. Combine the three Omics datasets into a list and look at the size of each matrix. MOFA allows for a handy overview of the data by displaying dimensions of each OMIC. Missing data are shown in gray but there are no missing values in this example dataset.
omics<-list(scRNAseq=as.matrix(scRNAseq),scBSseq=as.matrix(scBSseq),scATACseq=as.matrix(scATACseq))
lapply(omics,dim)## $scRNAseq
## [1] 12145 113
##
## $scBSseq
## [1] 4524 113
##
## $scATACseq
## [1] 9471 113MOFAobject <- create_mofa_from_matrix(omics)
plot_data_overview(MOFAobject, )
Now we will set some model training parameters such as distribution types corresponding to each OMIC, number of iterations, number of factors to be computed, etc.
#DEFINE DATA OPTIONS
DataOptions <- get_default_data_options(MOFAobject)
DataOptions## $scale_views
## [1] FALSE
##
## $scale_groups
## [1] FALSE
##
## $center_groups
## [1] TRUE
##
## $use_float32
## [1] TRUE
##
## $views
## [1] "scRNAseq" "scBSseq" "scATACseq"
##
## $groups
## [1] "group1"#DEFINE MODEL OPTIONS
ModelOptions <- get_default_model_options(MOFAobject)
mydistr <- c("gaussian","bernoulli","bernoulli")
names(mydistr) <- c("scRNAseq","scBSseq","scATACseq")
ModelOptions$likelihoods <- mydistr
ModelOptions$num_factors <- 20
ModelOptions## $likelihoods
## scRNAseq scBSseq scATACseq
## "gaussian" "bernoulli" "bernoulli"
##
## $num_factors
## [1] 20
##
## $spikeslab_factors
## [1] FALSE
##
## $spikeslab_weights
## [1] FALSE
##
## $ard_factors
## [1] FALSE
##
## $ard_weights
## [1] TRUE#DEFINE TRAIN OPTIONS
TrainOptions <- get_default_training_options(MOFAobject)
TrainOptions$seed <- 2018
# Automatically drop factors that explain less than 3% of variance in all omics
TrainOptions$drop_factor_threshold <- 0.03
TrainOptions$maxiter <- 1000
TrainOptions$verbose <- TRUE
TrainOptions## $maxiter
## [1] 1000
##
## $convergence_mode
## [1] "fast"
##
## $drop_factor_threshold
## [1] 0.03
##
## $verbose
## [1] TRUE
##
## $startELBO
## [1] 1
##
## $freqELBO
## [1] 5
##
## $stochastic
## [1] FALSE
##
## $gpu_mode
## [1] FALSE
##
## $gpu_device
## NULL
##
## $seed
## [1] 2018
##
## $outfile
## NULL
##
## $weight_views
## [1] FALSE
##
## $save_interrupted
## [1] FALSEFinally, we are ready to run MOFA.
NOTE: The first time you run MOFA on your machine, it will install a copy of Python locally. This is needed because the MOFA model is actually implemented in Python (see MOFA2 page for details). Installing Python will take 5-10 minutes. This process will only happen the first time you run MOFA. After Python is installed, MOFA will run on the example dataset. This will take about 5 minutes. If you want to skip installing and running MOFA, you may load the pre-fit MOFA object using the
readRDS()command below.
MOFAobject <- prepare_mofa(MOFAobject, data_options = DataOptions, model_options = ModelOptions, training_options = TrainOptions)
MOFAobject <- run_mofa(MOFAobject, use_basilisk = TRUE)As mentioned above, you may skip the code chunk directly above this
and load the pre-fit MOFA object. You may either import this from the
local folder data if you are on the demo server, or from
the Web location if you are not.
MOFAobject <- readRDS("data/MOFA2_singlecell.rds")MOFAobject <- readRDS(url("https://usda-ree-ars.github.io/SEAStatsData/omics/MOFA2_singlecell.rds"))View the contents of the MOFA object.
MOFAobject## Trained MOFA with the following characteristics:
## Number of views: 3
## Views names: scRNAseq scBSseq scATACseq
## Number of features (per view): 12145 4524 9471
## Number of groups: 1
## Groups names: group1
## Number of samples (per group): 113
## Number of factors: 3Let us check the covariance of the OMICs layers. The plots generated
by plot_variance_explained() are probably the most
important outputs of MOFA.
#ANALYZE RESULTS OF MOFA INTEGRATION
head(get_variance_explained(MOFAobject)$r2_total[[1]])## scRNAseq scBSseq scATACseq
## 13.339221 5.562496 4.713708head(get_variance_explained(MOFAobject)$r2_per_factor[[1]])## scRNAseq scBSseq scATACseq
## Factor1 6.923169 5.49092293 4.49206829
## Factor2 3.559488 0.03377795 0.16059279
## Factor3 3.318769 0.03517270 0.06973147get_variance_explained(MOFAobject)## $r2_total
## $r2_total$group1
## scRNAseq scBSseq scATACseq
## 13.339221 5.562496 4.713708
##
##
## $r2_per_factor
## $r2_per_factor$group1
## scRNAseq scBSseq scATACseq
## Factor1 6.923169 5.49092293 4.49206829
## Factor2 3.559488 0.03377795 0.16059279
## Factor3 3.318769 0.03517270 0.06973147MOFAobject@cache$variance_explained## $r2_total
## $r2_total$group1
## scRNAseq scBSseq scATACseq
## 13.339221 5.562496 4.713708
##
##
## $r2_per_factor
## $r2_per_factor$group1
## scRNAseq scBSseq scATACseq
## Factor1 6.923169 5.49092293 4.49206829
## Factor2 3.559488 0.03377795 0.16059279
## Factor3 3.318769 0.03517270 0.06973147plot_list <- plot_variance_explained(MOFAobject, x='view', y='factor', plot_total = T)
plot_list[[2]] #total variance
plot_list[[1]] #variance by factor
We can see that scRNAseq provides the largest variation (13%) in the total integrative OMICs scNMT data set, and scBSseq and scATACseq contribute around 5% of variation. We also observe that MOFA selected 3 Latent Factors (LFs); scRNAseq contributes to all of them while scBSseq and scATACseq contribute only to the first LF. Now let us interpret the LFs by displaying the features associated with each LF:
NumFactors<-dim(get_factors(MOFAobject)$group1)[2]
plot_weights(MOFAobject, view = "scRNAseq", factor = 1)
plot_top_weights(object = MOFAobject, view = "scRNAseq", factor = 1, nfeatures = 10)
plot_data_heatmap(object = MOFAobject, view = "scRNAseq", factor = "Factor1", features = 10, transpose = F, show_colnames=TRUE, show_rownames=TRUE,cluster_cols = TRUE)
plot_weights(MOFAobject, view = "scBSseq", factor = 1)
plot_top_weights(object = MOFAobject, view = "scBSseq", factor = 1, nfeatures = 10)
plot_data_heatmap(object = MOFAobject, view = "scBSseq", factor = "Factor1", features = 10, transpose = F, show_colnames=TRUE, show_rownames=TRUE,cluster_cols = TRUE)
plot_weights(MOFAobject, view = "scATACseq", factor = 1)
plot_top_weights(object = MOFAobject, view = "scATACseq", factor = 1, nfeatures = 10)
plot_data_heatmap(object = MOFAobject, view = "scATACseq", factor = "Factor1", features = 10, transpose = F, show_colnames=TRUE, show_rownames=TRUE,cluster_cols = TRUE)
Above we ranked by their importance and visualized features from each OMIC that contribute the most to the Factor 1 where all three OMICs (scRNAseq, scBSseq and scATACseq) are strongly correlated. Therefore one can expect those features to have something in common, i.e. they should also be correlated or involved in the same biological pathway.
Let us also display the cells in the low-dimensional latent space. for this purpose we will input MOFA factors into tSNE and UMAP, which are machine-learning dimension reduction and visualization methods. Display the cells as a consensus plot for scRNAseq, scBSseq and scATACseq:
controls<-c("EB_P1D12","EB_P1E12","EB_P1F12","EB_P1G12","EB_P2B12","EB_P2D12","EB_P2E12","EB_P2F12","EB_P2G12",
"ESC_B12","ESC_C12","ESC_D12")
colnames(scRNAseq)<-ifelse(colnames(scRNAseq)%in%controls,"CONTROL_CELL",colnames(scRNAseq))
cell_types<-matrix(unlist(strsplit(colnames(scRNAseq),"_")),ncol=2,byrow=TRUE)[,1]
sample_metadata <- data.frame(sample = samples_names(MOFAobject)[[1]],condition = cell_types)
samples_metadata(MOFAobject) <- sample_metadata
head(samples_metadata(MOFAobject), n=3)## sample condition group
## 1 ESC_A02 ESC group1
## 2 ESC_A03 ESC group1
## 3 ESC_A04 ESC group1set.seed(12345)
MOFAobject <- run_umap(MOFAobject, n_neighbors = round(sqrt(length(cell_types)),0), min_dist = 1)
MOFAobject <- run_tsne(MOFAobject, perplexity = round(sqrt(length(cell_types)),0))
plot_dimred(MOFAobject, method = "TSNE", dot_size = 3, color_by = "condition")
plot_dimred(MOFAobject, method = "UMAP", dot_size = 3, color_by = "condition")
We conclude that the MOFA unsupervised integrative OMICs approach can distinguish between all three types of cells: ESCs, EBs and Controls. Now let us see what MOFA factors contribute to different clusters:
plot_dimred(MOFAobject, method = "UMAP", dot_size = 3, color_by = "Factor1")
plot_dimred(MOFAobject, method = "UMAP", dot_size = 3, color_by = "Factor2")
plot_dimred(MOFAobject, method = "UMAP", dot_size = 3, color_by = "Factor3")
We can see that the ESC cluster is almost entirely driven by the MOFA Factor 1 meaning that all three OMICs (scRNAseq, scBSseq and scATACseq) contribute more or less equally into this cluster. In contrast the Control cluster is driven by the factor 2, which is mainly due to scRNAseq contribution. Finally, Factor 3 contributes more or less equally to all the three clusters, implying that scRNAseq impacts all the three clusters.
Session Info
sessionInfo()## R version 4.5.0 (2025-04-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 22631)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] MOFA2_1.18.0 mixOmics_6.32.0 ggplot2_3.5.2 lattice_0.22-6
## [5] MASS_7.3-65
##
## loaded via a namespace (and not attached):
## [1] gtable_0.3.6 dir.expiry_1.16.0 ellipse_0.5.0
## [4] xfun_0.52 bslib_0.9.0 corrplot_0.95
## [7] ggrepel_0.9.6 rhdf5_2.52.0 rhdf5filters_1.20.0
## [10] vctrs_0.6.5 tools_4.5.0 generics_0.1.4
## [13] stats4_4.5.0 parallel_4.5.0 tibble_3.2.1
## [16] rARPACK_0.11-0 pkgconfig_2.0.3 pheatmap_1.0.12
## [19] Matrix_1.7-3 RColorBrewer_1.1-3 S4Vectors_0.46.0
## [22] lifecycle_1.0.4 compiler_4.5.0 farver_2.1.2
## [25] stringr_1.5.1 codetools_0.2-20 htmltools_0.5.8.1
## [28] sass_0.4.10 yaml_2.3.10 pillar_1.10.2
## [31] crayon_1.5.3 jquerylib_0.1.4 tidyr_1.3.1
## [34] uwot_0.2.3 BiocParallel_1.42.0 cachem_1.1.0
## [37] DelayedArray_0.34.1 abind_1.4-8 RSpectra_0.16-2
## [40] basilisk_1.20.0 tidyselect_1.2.1 digest_0.6.37
## [43] Rtsne_0.17 stringi_1.8.7 dplyr_1.1.4
## [46] reshape2_1.4.4 purrr_1.0.4 labeling_0.4.3
## [49] forcats_1.0.0 cowplot_1.1.3 fastmap_1.2.0
## [52] grid_4.5.0 SparseArray_1.8.0 cli_3.6.5
## [55] magrittr_2.0.3 S4Arrays_1.8.0 h5mread_1.0.0
## [58] corpcor_1.6.10 withr_3.0.2 filelock_1.0.3
## [61] scales_1.4.0 XVector_0.48.0 rmarkdown_2.29
## [64] matrixStats_1.5.0 igraph_2.1.4 reticulate_1.42.0
## [67] gridExtra_2.3 png_0.1-8 HDF5Array_1.36.0
## [70] evaluate_1.0.3 knitr_1.50 IRanges_2.42.0
## [73] basilisk.utils_1.20.0 RcppAnnoy_0.0.22 rlang_1.1.6
## [76] Rcpp_1.0.14 glue_1.8.0 BiocGenerics_0.54.0
## [79] rstudioapi_0.17.1 jsonlite_2.0.0 Rhdf5lib_1.30.0
## [82] R6_2.6.1 plyr_1.8.9 MatrixGenerics_1.20.0