Vignette illustrating the use of MOFA on simulated data
9 February 2021
Package
MOFA 1.6.2
Contents
1 Simulate an example data set
To illustrate the MOFA workflow we simulate a small example data set with 3 different views. makeExampleData generates an untrained MOFAobject containing the simulated data. If you work on your own data use createMOFAobject to create the untrained MOFA object (see our vignettes for CLL and scMT data).
library(MOFA)By default the function makeExampleData produces a small data set containing 3 views with 100 features each and 50 samples. These parameters can be varied using the arguments n_views, n_features and n_samples.
set.seed(1234)
data <- makeExampleData()
MOFAobject <- createMOFAobject(data)## Creating MOFA object from a list of matrices
## please make sure that samples are stored in the columns and features in the rows...MOFAobject## Untrained MOFA model with the following characteristics:
## Number of views: 3
## View names: view_1 view_2 view_3
## Number of features per view: 100 100 100
## Number of samples: 502 Prepare MOFA: Set the training and model options
Once the untrained MOFAobject was created, we can specify details on data processing, model specifications and training options such as the number of factors, the likelihood models etc. Default option can be obtained using the functions getDefaultTrainOptions, getDefaultModelOptions and getDefaultDataOptions. We describe details on these option in our vignettes for CLL and scMT data.
Using prepareMOFA the model is set up for training.
TrainOptions <- getDefaultTrainOptions()
ModelOptions <- getDefaultModelOptions(MOFAobject)
DataOptions <- getDefaultDataOptions()
TrainOptions$DropFactorThreshold <- 0.013 Run MOFA
Once the MOFAobject is set up we can use runMOFA to train the model. As depending on the random initilization the results might differ, we recommend to use runMOFA multiple times (e.g. ten times, here we use a smaller number for illustration as the model training can take some time). As a next step we will compare the different fits and select the best model for downstream analyses.
n_inits <- 3
MOFAlist <- lapply(seq_len(n_inits), function(it) {
TrainOptions$seed <- 2018 + it
MOFAobject <- prepareMOFA(
MOFAobject,
DataOptions = DataOptions,
ModelOptions = ModelOptions,
TrainOptions = TrainOptions
)
runMOFA(MOFAobject)
})## Checking data options...## Checking training options...## Checking model options...## [1] "No output file provided, using a temporary file..."## Warning: Factor 2 is strongly correlated with the total expression for each sample in view_1## Such (strong) factors usually appear when count-based assays are not properly normalised by library size.## Checking data options...## Checking training options...## Checking model options...## [1] "No output file provided, using a temporary file..."## Warning: Factor 2 is strongly correlated with the total expression for each sample in view_1## Such (strong) factors usually appear when count-based assays are not properly normalised by library size.## Checking data options...## Checking training options...## Checking model options...## [1] "No output file provided, using a temporary file..."## Warning: Factor 2 is strongly correlated with the total expression for each sample in view_1## Such (strong) factors usually appear when count-based assays are not properly normalised by library size.4 Compare different random inits and select the best model
Having a list of trained models we can use compareModels to get an overview of how many factors were inferred in each run and what the optimized ELBO value is (a model with larger ELBO is preferred).
compareModels(MOFAlist)
With compareFactors we can get an overview of how robust the factors are between different model instances.
compareFactors(MOFAlist)
For down-stream analyses we recommned to choose the model with the best ELBO value as is done by selectModel.
MOFAobject <- selectModel(MOFAlist, plotit = FALSE)
MOFAobject## Trained MOFA model with the following characteristics:
## Number of views: 3
## View names: view_1 view_2 view_3
## Number of features per view: 100 100 100
## Number of samples: 50
## Number of factors: 55 Downstream analyses
On the trained MOFAobject we can now start looking into the inferred factors, its weights etc. Here the data was generated using five factors, whose activity patterns we can recover using plotVarianceExplained.
plotVarianceExplained(MOFAobject)
For details on downstream analyses please have a look at the vignettes on the CLL data and scMT data.
6 SessionInfo
sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.12-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.12-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] ggplot2_3.3.3 MOFAdata_1.6.0
## [3] MOFA_1.6.2 MultiAssayExperiment_1.16.0
## [5] SummarizedExperiment_1.20.0 Biobase_2.50.0
## [7] GenomicRanges_1.42.0 GenomeInfoDb_1.26.3
## [9] IRanges_2.24.1 S4Vectors_0.28.1
## [11] BiocGenerics_0.36.0 MatrixGenerics_1.2.1
## [13] matrixStats_0.58.0 BiocStyle_2.18.1
##
## loaded via a namespace (and not attached):
## [1] splines_4.0.3 jsonlite_1.7.2 foreach_1.5.1
## [4] assertthat_0.2.1 BiocManager_1.30.10 highr_0.8
## [7] GenomeInfoDbData_1.2.4 vipor_0.4.5 yaml_2.2.1
## [10] ggrepel_0.9.1 corrplot_0.84 pillar_1.4.7
## [13] lattice_0.20-41 glue_1.4.2 reticulate_1.18
## [16] digest_0.6.27 RColorBrewer_1.1-2 XVector_0.30.0
## [19] colorspace_2.0-0 cowplot_1.1.1 htmltools_0.5.1.1
## [22] Matrix_1.3-2 plyr_1.8.6 pkgconfig_2.0.3
## [25] pheatmap_1.0.12 magick_2.6.0 bookdown_0.21
## [28] zlibbioc_1.36.0 purrr_0.3.4 scales_1.1.1
## [31] tibble_3.0.6 mgcv_1.8-33 generics_0.1.0
## [34] farver_2.0.3 ellipsis_0.3.1 withr_2.4.1
## [37] magrittr_2.0.1 crayon_1.4.1 evaluate_0.14
## [40] GGally_2.1.0 nlme_3.1-152 doParallel_1.0.16
## [43] beeswarm_0.2.3 tools_4.0.3 lifecycle_0.2.0
## [46] stringr_1.4.0 Rhdf5lib_1.12.1 munsell_0.5.0
## [49] DelayedArray_0.16.1 compiler_4.0.3 rlang_0.4.10
## [52] rhdf5_2.34.0 grid_4.0.3 RCurl_1.98-1.2
## [55] iterators_1.0.13 rhdf5filters_1.2.0 rappdirs_0.3.3
## [58] bitops_1.0-6 labeling_0.4.2 rmarkdown_2.6
## [61] gtable_0.3.0 codetools_0.2-18 reshape_0.8.8
## [64] DBI_1.1.1 reshape2_1.4.4 R6_2.5.0
## [67] knitr_1.31 dplyr_1.0.4 stringi_1.5.3
## [70] ggbeeswarm_0.6.0 Rcpp_1.0.6 vctrs_0.3.6
## [73] tidyselect_1.1.0 xfun_0.20