MOFA 1.2.0
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: 50
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.01
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
## Warning: Factor 4 is strongly correlated with the total expression for each sample in view_3
## Such (strong) factors usually appear when count-based assays are not properly normalised by library size.
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: 5
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.
sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.10-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.10-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.2.1 MOFAdata_1.1.0
## [3] MOFA_1.2.0 MultiAssayExperiment_1.12.0
## [5] SummarizedExperiment_1.16.0 DelayedArray_0.12.0
## [7] BiocParallel_1.20.0 matrixStats_0.55.0
## [9] Biobase_2.46.0 GenomicRanges_1.38.0
## [11] GenomeInfoDb_1.22.0 IRanges_2.20.0
## [13] S4Vectors_0.24.0 BiocGenerics_0.32.0
## [15] BiocStyle_2.14.0
##
## loaded via a namespace (and not attached):
## [1] ggrepel_0.8.1 Rcpp_1.0.2 lattice_0.20-38
## [4] assertthat_0.2.1 digest_0.6.22 foreach_1.4.7
## [7] R6_2.4.0 plyr_1.8.4 evaluate_0.14
## [10] pillar_1.4.2 zlibbioc_1.32.0 rlang_0.4.1
## [13] lazyeval_0.2.2 Matrix_1.2-17 reticulate_1.13
## [16] rmarkdown_1.16 labeling_0.3 stringr_1.4.0
## [19] pheatmap_1.0.12 RCurl_1.95-4.12 munsell_0.5.0
## [22] compiler_3.6.1 vipor_0.4.5 xfun_0.10
## [25] pkgconfig_2.0.3 ggbeeswarm_0.6.0 htmltools_0.4.0
## [28] tidyselect_0.2.5 tibble_2.1.3 GenomeInfoDbData_1.2.2
## [31] bookdown_0.14 codetools_0.2-16 reshape_0.8.8
## [34] withr_2.1.2 crayon_1.3.4 dplyr_0.8.3
## [37] bitops_1.0-6 grid_3.6.1 GGally_1.4.0
## [40] jsonlite_1.6 gtable_0.3.0 magrittr_1.5
## [43] scales_1.0.0 stringi_1.4.3 XVector_0.26.0
## [46] reshape2_1.4.3 doParallel_1.0.15 cowplot_1.0.0
## [49] Rhdf5lib_1.8.0 RColorBrewer_1.1-2 iterators_1.0.12
## [52] tools_3.6.1 glue_1.3.1 beeswarm_0.2.3
## [55] purrr_0.3.3 yaml_2.2.0 rhdf5_2.30.0
## [58] colorspace_1.4-1 BiocManager_1.30.9 corrplot_0.84
## [61] knitr_1.25