---
title: "Joint Dimensional Reduction methods for multi-omics data analysis with MOFA"
author: "Arnaud Gloaguen - Vincent Guillemot - Jimmy Vandel"
date: '`r format(Sys.Date(), "%B %d, %Y")`'
output: 
  rmdformats::material:
    code_folding: show
    use_bookdown: true
    thumbnails: false
    css: "ETBII.css"
header-includes:
  - \usepackage{subfig}
  - \usepackage{booktabs}
  - \usepackage{xcolor}
---
<style>
.textbox {
  border: 1px solid #ccc;
  padding: 10px;
  background-color: #f5f5f5;
}

.highlight {
  background-color: yellow;
  font-weight: bold;
}
</style>

```{r setup, eval = TRUE, echo = FALSE}

knitr::opts_knit$set(root.dir = ".")
knitr::opts_chunk$set(cache = TRUE)
options(knitr.table.format = "html")

```

<!-- ####################################################### -->

# Overview of the Tutorial


<div class="textbox">
  **General Statement**: 
  
In this tutorial, Principal Component Analysis (**PCA**), along with two multi-omics techniques : Regularized Generalized Canonical Correlation Analysis (**RGCCA**) and Multi-Omics Factor Analysis (**MOFA**), are going to be explored. These methods have in common that they all belong to the field of **multivariate statistics**, but also that, in different ways, they can be seen as **generalizations of PCA to more than one omic/matrix/block/table**. Indeed, similarly to PCA, they all aim at estimating a reduced/latent space that describes well the different omics. Some part of the literature refer to these methods as **Joint Dimension Reduction (JDR) methods** because their goal is to estimate this reduced space such that it is common to all the omic data used in the analysis (this is less true for RGCCA). Actually, applying PCA separatly on different omic data in not considered as a JDR method, but applying PCA on the concatenation of multiple omic data acquired on the same set of individuals is.
  
</div>

<br>

Bellow is figure **summarizing the key concept behind JDR methods**. This figure was taken from an interesting **benchmark** that compares, on different multi-omic data-sets, multiple JDR methods (among which RGCCA and MOFA).

```{r JDRprinciple, echo = FALSE, fig.align = "center", out.width = "100%", fig.cap = "Summarizing Joint Dimension Reduction (JDR) pinciple. A.: Multiomics are profiled from the same sample. Each omics corresponds to a different matrix Xi. jDR methods factorize the matrices $X^i$ into the product of a factor matrix $F$ and weight matrices $A^i$. These matrices can then be used to cluster samples and identify molecular processes. Figure and caption were taken from: Cantini, L., Zakeri, P., Hernandez, C. et al. [Benchmarking joint multi-omics dimensionality reduction approaches for the study of cancer](https://doi.org/10.1038/s41467-020-20430-7). Nat Commun 12, 124 (2021)."}

knitr::include_graphics("images/JDR_principle.PNG")

```

The tutorial is organized as follow:

  - First a brief presentation of [Multi-Omics Factor Analysis (MOFA)](https://gitlab.com/ifb-elixirfr/training/etbii/-/blob/main/THEMATIC-SESSIONS/DimensionReduction/01_Theory/MOFA.pdf?ref_type=heads).
  - In section \@ref(tcga), we present the biological context data-set used as example in this tutorial.
  - In section \@ref(mofa), you will apply MOFA on the data-set.

<!-- ####################################################### -->

# Libraries 

## Load libraries

<div class="textbox">
**DIY**: Load required libraries : ....

```{r loadLibrariesMOFA, class.source = 'fold-hide', results = "hide", message = FALSE, warning = FALSE}

library("FactoMineR") #for PCA
library("factoextra") #for fviz_pca_ind
library("MOFA2")      #for MOFA
library("DT")
library("rmarkdown")
library("knitr")
library("rmdformats")
library("bookdown")

```
</div>
<br>

<!-- ####################################################### -->


# TCGA breast cancer dataset {#tcga}
This dataset includes measurements of mRNA, miRNA, and protein abundances taken from samples representing different molecular subtypes of breast cancer.

This dataset was obtained from [Mibiomics](https://urldefense.com/v3/__https://gitlab.univ-nantes.fr/combi-ls2n/mibiomics/-/tree/master/data/data_breast__;!!JFdNOqOXpB6UZW0!omJ3QEzdDCejWLznykgU9pq7Jor0wApVNg_dpRe8PkELgqIGzToeAJXNeBjDOMbliS16gycA8TT22nmLMEQi87ViEcvs9g3cuLAy0w$ )

<div class="textbox">
Download the dataset from the [IFB/SIB summer school's GitHub repository](https://urldefense.com/v3/__https://github.com/sib-swiss/summer-school-multiomics-data-analysis-and-integration/tree/master/general/datasets/BreastCancer_WGCNA__;!!JFdNOqOXpB6UZW0!omJ3QEzdDCejWLznykgU9pq7Jor0wApVNg_dpRe8PkELgqIGzToeAJXNeBjDOMbliS16gycA8TT22nmLMEQi87ViEcvs9g35Qd_UDQ$ )
</div>
<br>

For this dataset, we will use 3 modalities: 

* `mrna.csv`, contains mRNA expression data
* `mirna.csv`, contains miRNA expression data
* `protein.csv`, contains protein abundance data
 
We will also use the file `clinical.csv`, which contains samples metadata, including molecular subtypes.

<div class="textbox">

Again, in the upcoming sections, our primary focus will be on TCGA dataset. However, in the rest of the week, **the procedure we will discuss can also be applied to other datasets**. 

</div>
<br>

[TCGA](https://urldefense.com/v3/__https://www.cancer.gov/ccg/research/genome-sequencing/tcga__;!!JFdNOqOXpB6UZW0!omJ3QEzdDCejWLznykgU9pq7Jor0wApVNg_dpRe8PkELgqIGzToeAJXNeBjDOMbliS16gycA8TT22nmLMEQi87ViEcvs9g31Vnr1eA$ ).

```{r TCGAoverview, echo = FALSE, fig.align = "center", out.width = "100%", fig.cap = " [The Cancer Genome Atlas (TCGA) overview](https://www.cancer.gov/ccg/research/genome-sequencing/tcga)"}

knitr::include_graphics("images/blog_genomic_data_commons.png")
```


<div class="textbox">
Our objective is to build a multi-omic molecular signature discriminating TCGA patients
</div>
<br>

<!-- ####################################################### -->

# Input data {#loadata}

## In summary

For **MOFA**, data must be a list of **matrices** with specific shapes: 

- **variables** (e.g. genes, proteins ...) are displayed on the **rows** 
- **samples** (e.g. patients, organisms ...) are displayed on the **columns**

## Load your dataset

Load all the modalities of the data-set into R, except for the meta/clinical data.
You will probably need to read a csv/tsv file: look at the documentation for the function ```read.table```:

```{r doc_read.table}
?read.table
```

If your data is not stored into a csv/tsv file and you don't know how to load it, don't hesitate to ask us.

<div class="textbox">
**DIY**: Load all the modalities, except for the meta/clinical data, as a list of **matrices**.
```{r loadData, class.source = 'fold-hide'}

data_path  = "../data/TCGA/"
modalities    = c("mrna", "mirna", "protein")
data_files    = paste0(data_path, modalities,  ".csv")
blocks        = lapply(data_files, function(x) t(read.table(x, sep = ",", header = T, row.names = 1)))
names(blocks) = modalities
```
</div>
<br>


<div class="textbox">
**DIY**: Load separately the meta/clinical data.
```{r loadClinical, class.source = 'fold-hide'}
clinical_data = read.table(paste0(data_path,"clinical.csv"), sep = ",", header= T, row.names = 1)

# We rename the variables as they are too long
colnames(clinical_data) = c("age_at_diagnosis", "neoplasm_subdivision", "gender", "vital_status", "subtype")


```
</div>
<br>
  
<!-- ####################################################### -->

# Multi-Omics Factor Analysis (MOFA) {#mofa}

Multi-Omics Factor Analysis (MOFA) is a multi-omics data integration approach giving a computational framework for unsupervised discovery of the principal axes of biological and technical variation when multiple omics assays are applied to the same samples. Its extension, MOFA+, is designed for integration of single-cell multi-modal data.

![MOFA: model overview and downstream analyses](https://urldefense.com/v3/__https://i.imgur.com/4gG6vXG.jpg__;!!JFdNOqOXpB6UZW0!omJ3QEzdDCejWLznykgU9pq7Jor0wApVNg_dpRe8PkELgqIGzToeAJXNeBjDOMbliS16gycA8TT22nmLMEQi87ViEcvs9g0-pe0d1g$ )
Argelaguet et al., Mol Syst Biol, 2018.

*MOFA decomposes each omics into the product of a factor matrix and omics-specific weight matrices.*

If you are interested going further on the subject, here are an non exhaustive list of publications that might interest you: 

* Link to [MOFA article](https://urldefense.com/v3/__https://www.embopress.org/doi/full/10.15252/msb.20178124__;!!JFdNOqOXpB6UZW0!omJ3QEzdDCejWLznykgU9pq7Jor0wApVNg_dpRe8PkELgqIGzToeAJXNeBjDOMbliS16gycA8TT22nmLMEQi87ViEcvs9g3I70LpYQ$)

* Link to [MOFA+ article](https://urldefense.com/v3/__https://genomebiology.biomedcentral.com/articles/10.1186/s13059-020-02015-1__;!!JFdNOqOXpB6UZW0!omJ3QEzdDCejWLznykgU9pq7Jor0wApVNg_dpRe8PkELgqIGzToeAJXNeBjDOMbliS16gycA8TT22nmLMEQi87ViEcvs9g3tlD2erg$)

* Link to [MOFA tutorials](https://biofam.github.io/MOFA2/tutorials.html)

* Link to [MOFA FAQ](https://biofam.github.io/MOFA2/faq.html)

* Link to [Python code](https://github.com/bioFAM/mofapy2/tree/f4e4978c1fda7ed98ae2351324b23dd03119bddf)

* Link to [Documentation of MOFA2'Bioconductor package](https://bioconductor.org/packages/release/bioc/manuals/MOFA2/man/MOFA2.pdf)

* If you want to deep dive into MOFA's equations, good LUCK !! Here are some links that can help you. First, [Ricard Argelaguet's PhD manusccript](https://www.dropbox.com/scl/fi/doqrcgh477c5zgyyyvrx1/thesis.pdf?rlkey=wcj2mqh8dmltfsfrovjkamnh5&e=1&dl=0), the first author of MOFA, along with the [appendix of MOFA](https://pmc.ncbi.nlm.nih.gov/articles/instance/6010767/bin/MSB-14-e8124-s001.pdf) and the [appendix of MOFA+](https://pmc.ncbi.nlm.nih.gov/articles/instance/7212577/bin/13059_2020_2015_MOESM2_ESM.pdf). [Here]() is also a presentation of the paper [Multi-Omics Transfer Learning (MOTL)](https://doi.org/10.1186/s13059-025-03675-7) with a theoretical part that heavily rely on MOFA. 

We will perform the MOFA analysis and interpret the results using an object containing various data modalities.

## First steps

### Vocabulary

In this MOFA section, we will use the following vocabulary :

* view or modality : omic type that will compose a block of data in the analysis (eg transcriptomic, proteomic, genomic...)
* feature : variable of a specific view (eg transcipt, protein, locus...)
* factor values : compose the Z matrix (similar to components/scores in PCA)
* weight values : compose the W matrix (similar to wieghts/loadings in PCA)

### Input data

You should use a multiview matrix that contains all data cleaned and normalized. Each view of the matrix is a block of data corresponding to a modality, with features as lines and samples as columns. No concordance between lines across modalities is required, but columns should correspond to the same sample across modalities (you can still have some missing data). 

MOFA accept several formats for this multiview matrix: data.frame, list of matrices, multiAssayExperiments...

<div class="textbox">
**DIY**: write the code to create such a multiview matrix for your data.

```{r MOFA_data_loading, class.source = 'fold-hide', results = "hide", message = FALSE, warning = FALSE}
#we just have to transpose the corresponding blocks 
multiview_mat=list(mrna=t(blocks$mrna), mirna=t(blocks$mirna), protein=t(blocks$protein))
```
</div>
<br>


## Prepare MOFA

Before running MOFA, you should first create a MOFA object including the different data modalities and then specify options.

### Create MOFA object

You can create a MOFA object from an artificially generated dataset with ```make_example_data``` function or using a real dataset as the multiview matrix you just created. Each sub-matrix contains one modality with sample as columns.  


<div class="textbox">
**DIY**: use the ```create_mofa``` function to create a *MOFAobject* containing input data, then use the ```plot_data_overview``` function over this object to check available information for each modality.

```{r MOFAobject_creation, message=FALSE, warning=FALSE}
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                          Create the MOFA object                            ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##

MOFAobject <- create_mofa(multiview_mat)

#Plot overview of the input data, including the number of samples, views, features, and the missing assays.
plot_data_overview(MOFAobject) # everything matches here, but we could have kept non matching samples with MOFA
```

NB: by default no group is defined among samples, however you can specify some groups if necessary using "groups" option in the ```create_mofa``` function. Group information is useful to analyse single cell dataset by specifying groups of sample (*ie* cell) coming from the same biological sample or to take into account batch or experimental conditions.
</div>
<br>


### MOFA options {.tabset}

We have to define MOFA options before running it, they are grouped in 3 categories :

* Data options
* Model options
* Training options

You will have to create a list for each of these option categories.




#### Data options

**Main data options are:**

* scale_views: logical indicating whether to scale views to have the same unit variance. As long as the scale differences between the views is not too high, this is not required. Default is FALSE.
* scale_groups: logical indicating whether to scale groups to have the same unit variance. As long as the scale differences between the groups is not too high, this is not required. Default is FALSE.
* center_groups: logical indicating whether to center groups to have the same mean. As long as the mean differences between the groups is not too high, this is not required. Default is TRUE.

*Beware*,  scale_views and center_groups are only used for Gaussian likelihoods, see [code](https://github.com/bioFAM/mofapy2/blob/f4e4978c1fda7ed98ae2351324b23dd03119bddf/mofapy2/build_model/utils.py#L78). 

You can also change the name of views and groups through the 'data_options' field of the MOFA object.

<div class="textbox">
**DIY**: retrieve the list of default data options associated to your MOFAobject in a dedicated variable (eg named *data_options*) with the ```get_default_data_options``` function and modify any fields of this list if necessary.

```{r MOFA_dataoptions}
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                              MOFA data options                             ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##

data_options <- get_default_data_options(MOFAobject)
data_options$scale_views <- TRUE
data_options
```
</div>
<br>



#### Model options

**Main model options are:**

* likelihoods: data nature per view (possible options are “gaussian”, “poisson”, “bernoulli”). 
'gaussian' is for continuous data, 'bernoulli' for binary data and 'poisson' for discrete data. By default they are inferred automatically.
* num_factors: number of factors to use for the factorization

You can also specified through 'model_options' the use of sparsity priors on the weight W, and which type of prior to use (none, Spike-and-slab, Automatic Relevance Determination (ARD) or both). In the original MOFA approach, sparsity is only applied on the W matrix (both with ARD and spike-and-slab). In MOFA+ sparsity prior is also applied to the factor matrix Z when considering several groups. 

Corresponding options are : spikeslab_factors, spikeslab_weights, ard_factors and ard_weights.

<div class="textbox">
**DIY**: retrieve the list of default model options associated to your MOFAobject in a dedicated variable (eg named *model_options*) with the ```get_default_model_options``` function and check first the components of these list and modify if necessary the sparsity levels to consider in your model. Likelihoods and number of factors will be modified in the next steps.

```{r MOFA_get_modeloptions}
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                              MOFA model options                            ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##

model_options <- get_default_model_options(MOFAobject)
model_options$spikeslab_weights <- TRUE
model_options
```
</div>
<br>


**Likelihoods to use.**
<div class="textbox">

The model flexibly handles different data distribution, it is important to choose the likelihoods corresponding to the data (you can visually check data distribution for that). You can also let MOFA deciding for you, but it's better if you don't. [Here](https://github.com/bioFAM/mofapy2/blob/f4e4978c1fda7ed98ae2351324b23dd03119bddf/mofapy2/build_model/utils.py#L99) is the function used by MOFA to automatically decipher likelihoods.

**DIY**: check visually the distribution of each data modality using basic function such as ```hist```.

*Visual inspection of mRNA data*

```{r MOFA_countdistribmRNA}
## mRNA histogram
###################################

hist(MOFAobject@data$mrna$group1 , breaks = 100, main = "Histogram of mRNA count values stored in MOFA object")
```

*Visual inspection of miRNA data*

```{r MOFA_countdistribmiRNA}
## miRNA histogram
###################################

hist(MOFAobject@data$mirna$group1 , breaks = 100, main = "Histogram of miRNA values stored in MOFA object")
```
*Visual inspection of miRNA data ==> Do not forget that we are going to work with centered data*

```{r MOFA_countdistribmiRNA_centered}
## miRNA histogram
###################################

hist(scale(t(MOFAobject@data$mirna$group1), center = T, scale = F) , breaks = 100, main = "Histogram of centered miRNA values stored in MOFA object")
```


*Visual inspection of protein data*

```{r MOFA_countdistribProtein}
## protein histogram
###################################

hist(MOFAobject@data$protein$group1 , breaks = 100, main = "Histogram of protein values stored in MOFA object")
```


*Open question: Should we model all these data modalities with a gaussian distribution?*

**DIY**: if you decide it, change the corresponding likelihood for a given modality in the *model_options* variable.

To change this option, you can use *model_options$likelihoods* <- "gaussian" or "poisson" or "bernoulli"

```{r MOFA_changedistribMutations}
## as an example
#model_options$likelihoods["mrna"] <- "gaussian" or "poisson" or "bernoulli"
```

You can also transform your data to get a distribution closer than the expected distribution (for gaussian likelihood), it's better to anticipate this step before creating your MOFAobject, or you can modify the corresponding fields in ```MOFAobject@data```.

```{r MOFA_changedistribmRNA}
## as an example
#MOFAobject@data$mirna$group1<-log(MOFAobject@data$mirna$group1)
#hist(MOFAobject@data$mirna$group1 , breaks = 100, main = "Histogram of log miRNA count values stored in MOFA object")
```
</div>
<br>


**Number of factors**

<div class="textbox">

Defining the number of factors to extract is not an easy task. One can think about a Cross-validation procedure, i.e. setting this number to maximize a certain prediction performance. However, when you are doing unsupervised analysis this is no longer possible. Stability of the results might be a solution. You might change the initial condition several times and evaluate which factors are robust and which are not. We will try something is that direction later. Or you can rely on unsupervised analysis of each block and make assumptions on what will happen on your cohort. In that sense, Principal Component Analysis (PCA) is of great help. 
</div>
<br>


<div class="textbox">
**DIY**: Run a PCA with the `PCA` function of the `FactoMineR` package. You will have to set:

- `ncp`: the number of components.
- `scale.unit`: a boolean specifying if you want or not to scale (to unit variance) all your variables separately to make them comparable. You may skip this step if you previously performed a more expert normalization. Actually, scaling here could biais the extrapolation of results from PCA to MOFA as MOFA does not allow to perform such scaling. Though, [MOFA guidelines](https://biofam.github.io/MOFA2/faq.html) suggested that for *For count-based data such as RNA-seq or ATAC-seq we recommend size factor normalisation + variance stabilisation*. In order to be coherent between PCA and MOFA normalisation steps and making the hypothesis that an operation allowing to harmonize variance across features (as genes) has already been performed, this will be put to FALSE.
- `graph`: to avoid unwanted graphical output.
</div>
<br>


<div class="textbox">
  **Important Reminder**: Centering all your variables when performing PCA is mandatory, otherwize you are not extracting the components with the highest variance anymore. 
```{r mrnaPCA, class.source = 'fold-hide', fig.align = 'center'}

mrna_pca_results <- PCA(t(multiview_mat$mrna), scale.unit = FALSE, ncp = 25, graph = F) 

```
</div>
<br>

<div class="textbox">
  **Important Reminder**: Looking at ```mrna_pca_results$eig```, 
what is the number of components from which subsequent components represent less than 2% of the total variance ? How much of the total vairance do they represent ?
```{r mrnaPCA_cumulative_variance, class.source = 'fold-hide', fig.align = 'center'}

nb_mrna_comp              = max(which(mrna_pca_results$eig[, "percentage of variance"]>2))
mrna_ratio_total_variance = mrna_pca_results$eig[nb_mrna_comp, "cumulative percentage of variance"]

print(paste("Number of PCA components from which subsequent components represent less than 2% of the total variance:", nb_mrna_comp))
print(paste("This represent the following pourcentage of the total variance:", mrna_ratio_total_variance))


```
</div>
<br>

<div class="textbox">
  **Important Reminder**: Same question for the 2 other blocks
```{r other_omicsPCA_cumulative_variance, class.source = 'fold-hide', fig.align = 'center'}

#miRNA
mirna_pca_results <- PCA(t(multiview_mat$mirna), scale.unit = FALSE, ncp = 25, graph = F) 

nb_mirna_comp              = max(which(mirna_pca_results$eig[, "percentage of variance"]>2))
mirna_ratio_total_variance = mirna_pca_results$eig[nb_mirna_comp, "cumulative percentage of variance"]

print(paste("Number of PCA components from which subsequent components represent less than 2% of the total variance:", nb_mirna_comp))
print(paste("This represent the following pourcentage of the total variance:", mirna_ratio_total_variance))


#Protein
protein_pca_results <- PCA(t(multiview_mat$protein), scale.unit = FALSE, ncp = 25, graph = F) 

nb_protein_comp              = max(which(protein_pca_results$eig[, "percentage of variance"]>2))
protein_ratio_total_variance = protein_pca_results$eig[nb_mirna_comp, "cumulative percentage of variance"]

print(paste("Number of PCA components from which subsequent components represent less than 2% of the total variance:", nb_protein_comp))
print(paste("This represent the following pourcentage of the total variance:", protein_ratio_total_variance))

```
</div>
<br>

In the worst case scenario where no commonalities are found across omics, if we decide to extract 7+11+10 and we decide to drop every component explaining less than 2%, we get a result close to PCA, representing at least 46% of the total variance of each bloc. We hope factors will be dropped due to commonalities found accross omics.

<div class="textbox">
**DIY**: set the number of factors in the corresponding field in your *model_options* variable.

```{r MOFA_set_modeloptions}
## Model options
###################################

model_options$num_factors <- 28  
```
</div>
<br>

#### Training options

**Main training options are:**

* maxiter: maximum number of iterations (default is 1000). 

* convergence_mode: "fast" ($5\times10^{-4}$), "medium" ($5\times10^{-5}$), "slow" ($5\times10^{-6}$). For exploration, the fast mode is good enough.

* startELBO: initial iteration to compute the ELBO (the objective function used to assess convergence)

* freqELBO: frequency of computations of the ELBO (the objective function used to assess convergence)

* drop_factor_threshold: minimum fraction of variance explained criteria to drop factors while training (eg 0.01 for 1%). Default is -1 (no dropping of factors). One can change it to discard factors that explain a low percentage of variance in all views. Explained variance per factor and view is computed as follows:

$$R(\mathbf{z}_k, RNA) = 1 - \frac{\sum_j^{p_{RNA}}\sum_i^{n_{samples}} \left(\hat{y}_{i,j}^{rna} - z_{ik}\times w_{kj}^{RNA}\right)^2}{\sum_j^{p_{RNA}}\sum_i^{n_{samples}} \hat{y}_{i,j}^{rna^2}}$$
where $\mathbf{z}_k$ is the $k^{th}$ factor, $z_{ik}$ it $i^{th}$ element, $w_{kj}^{RNA}$ the weight associated with factor $k$, feature $j$ on the $RNA$ view. $\hat{y}_{i,j}^{RNA}$ is the pseudo-data associated with inividual $j$ on feature $j$ of the RNA view. If RNA's likelihood is Gaussian, pseudo-data egals original data, otherwise, it undergoes a dynamics transformation (all along the algorithm's iterations) to make it look like Gaussian data. Total explained variance for view RNA is computed as:

$$R(\mathbf{Z}, RNA) = 1 - \frac{\sum_j^{p_{RNA}}\sum_i^{n_{samples}} \left(\hat{y}_{i,j}^{rna} - \sum_k^Kz_{ik}\times w_{kj}^{RNA}\right)^2}{\sum_j^{p_{RNA}}\sum_i^{n_{samples}} \hat{y}_{i,j}^{rna^2}}$$


* stochastic: logical indicating whether to use stochastic variational inference (only required for very big data sets, default is FALSE).

* seed : numeric indicating the seed for reproducibility (default is 42).

<div class="textbox">
**DIY**: retrieve the list of default training options associated to your MOFAobject in a dedicated variable (*eg* named *train_options*) with the ```get_default_training_options``` function and modify this list if necessary. As mentioned earlier, let us drop factors if they represent less than 2% of the total variance accross all views.

```{r MOFA_trainingoptions}
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                          MOFA training options                             ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##

train_options <- get_default_training_options(MOFAobject)
train_options$convergence_mode <- "fast"
train_options$freqELBO <- 5
train_options$startELBO <- 1
train_options$seed <- 53
train_options$drop_factor_threshold <- 0.02

train_options
```
</div>
<br>


## Run MOFA

### Prepare the object and run MOFA

<div class="textbox">
**DIY**: use the ```prepare_mofa``` function to wrap all the previous variables including *MOFAobject*, *data_options*, *model_options* and *train_options.* Then run MOFA using the ```run_mofa``` function on the *MOFAobject*.  

```{r prepare_MOFA, warning=FALSE, message=FALSE}
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                          Prepare the MOFA object                           ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##

prepare_MOFAobject <- prepare_mofa(MOFAobject,
                                   data_options     = data_options,
                                   model_options    = model_options,
                                   training_options = train_options
)
```




```{r run_MOFA, warning=FALSE, message=FALSE}
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                                  Run MOFA                                  ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##

MOFAobject <- run_mofa(prepare_MOFAobject, "MOFA_models/MOFA2_TCGA.hdf5", use_basilisk = TRUE)

```

**NB:** You can also load a pre-computed MOFA model with the ```load_Model``` function from an hdf5 file.
</div>
<br>

## Results interpretation

### Get familiar with the MOFA object

#### Slots of MOFA object

Resulting *MOFAobject* is a S4 class used to store all relevant data to analyse a MOFA model.

<div class="textbox">
**DIY**: display available attributes using the ```slotNames``` function.

```{r MOFA_available_slots}
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                            MOFA object attributes                          ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##

slotNames(MOFAobject)
```

Object attributes accessible using '@' :

* data : input data (either a list of matrices or a MultiAssayExperiment)
* intercepts : feature intercepts
* imputed_data: imputed data (filled by running impute on the object, which is not necessary)
* samples_ and features_metadata (basically groups)
* expectations: expectations of the factors and weights ($Z and $W matrices)
* training_stats: statistics about iterative training process (time, elbo..)
* data_options: data processing options (given as parameters)
* model_options: model options (given as parameters)
* training_options: training options (given as parameters)
* stochastic_options: training options (given as parameters is stochastic option is TRUE in training_options)
* dimensions : dimensions for the number of views (M), groups (G), samples (N), features per view (D) and factors (K)
* on_disk : logical indicating whether data is loaded from disk.
* dim_red : non-linear dimensionality reduction manifolds.
* cache : variance explained by factors.
* status : auxiliary variable indicating whether the model has been trained.
* mefisto_options, covariates, covariates_warped and interpolated_Z are optional slot to store sample covariate only for training in MEFISTO
</div>
<br>

#### Functions to get infos from MOFA object

Some interesting functions to access these attributes: 

* ```factors_names()```: to get or set factor names
* ```features_names()```: to get or set feature names
* ```samples_names()```: to get or set sample names
* ```views_names()```: to get or set view names
* ```get_dimensions()```: to get dimensions (number of samples, features, etc.)
* ```get_factors()```: to get model factors
* ```get_weights()```: to get model weights

<div class="textbox">
**DIY**: try some of these functions on *MOFAobject*.

```{r MOFA_basic_function_tests}
## Function tests
###################################
views_names(MOFAobject)
get_dimensions(MOFAobject)
```
</div>
<br>

#### Add sample metadata to the model

<div class="textbox">
**DIY**: use the ```samples_metadata``` function to access and rewrite (or add) metadata on the *MOFAobject*.

```{r MOFA_addMetadata}
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                                 Metadata                                   ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##

## Add sample metadata to the model
temp_metadata <- cbind(rownames(clinical_data),clinical_data)
colnames(temp_metadata)[1]="sample" #one metadata column should contains IDs of samples matching with input data and labelled as "sample"
samples_metadata(MOFAobject) <- temp_metadata
rm(temp_metadata)
```
</div>
<br>

### Results

To interpret these results it is important to have information about the samples [Biological context (TCGA)].

#### General factors plots {.tabset}

First, we'll analyze all the factors at a high level to get an overview of the main trends in decomposition.

##### Factors correlation plot

<div class="textbox">
**DIY**: use the ```plot_factor_cor``` to display correlations between factors.

```{r MOFA_CorrelationPlot}
## Factor correlation plot
###################################

plot_factor_cor(MOFAobject)
```

Normally, factors should not be too correlated. If they are, you've probably chosen too many factors, creating redundancy between some of them.
</div>
<br>

##### Variance decomposition

Representing percentage of variance explained by each factor in each view is a key plot of MOFA and should always be done before inspecting factor values or weights.

<div class="textbox">
**DIY**: use the ```plot_variance_explained``` to display variance explained by each factor for each view.

```{r MOFA_varianceExplainedFactors, warning=FALSE}
## Variance explained values
###################################

plot_variance_explained(MOFAobject, max_r2=20)
#plot_variance_explained(MOFAobject, plot_total = T)[[1]]
```

**DIY**: you can get the numerical values from the cache attribute of the *MOFAobject*, in the ```variance_explained``` field.

```{r MOFA_varianceExplained, warning=FALSE}
## Variance explained values
###################################

MOFAobject@cache$variance_explained$r2_per_factor
```

As opposed to PCA, a factor does not necessarily explain less variance than the previous one. 

We can ask whether the model provides a good fit to the data. To check it we can plot the total variance explained (using all factors). The resulting values will depend on the nature of the data set, the number of samples, the number of factors, etc....

We can identify some rules about this amount: 

* Noisy datasets will result in a lower amounts of variance explained.
* Higher the number of factors --> higher the total variance explained.

The objective is not to reach 100%, even with a large number of factors. Explaining ~40% for a view thanks to a linear model is already good.

**DIY**: use the ```plot_variance_explained``` to display total variance explained for each view.

```{r MOFA_varianceHisto, warning=FALSE}
## Variance explained histogram
###################################

plot_variance_explained(MOFAobject, plot_total = T)[[2]]
```

The latent factors inferred by MOFA represent the underlying principal sources of variation among the samples.

Some factor can capture a source of variability that is present across all data modalities. Thus, its etiology is likely to be something very important for the studied phenomenon.
Whereas some factor can capture a very strong source of variation in only a specific modality
</div>
<br>

We can also access the variance explained per sample (only available with Gaussian data).

<div class="textbox">
**DIY**: if your data are exclusively Gaussian (regarding model_options parameters) use the ```calculate_variance_explained_per_sample``` function to display variance explained per sample.

```{r MOFA_varianceExplainedSample, warning=FALSE}
## Variance explained per sample
###################################

#calculate_variance_explained(MOFAobject)
datatable(calculate_variance_explained_per_sample(MOFAobject)$group1) #only with pure gaussian data
```
</div>
<br>

##### Summary plot of factors

We can look at associations (correlations) between MOFA factors and some metadata such as Gender, Age... these associations are a good starting point to interpret factors.

<div class="textbox">
**DIY**: use the ```correlate_factors_with_covariates``` function to search such associations with available metadata.

```{r MOFA_metadataCorrelations, fig.width = 10, warning=FALSE}
## Factor correlation with metadata
###################################

correlate_factors_with_covariates(MOFAobject, 
  covariates = names(clinical_data), 
  plot="log_pval"
)
```
</div>
<br>

ATTENTION: to compute such correlations, categorical are simply converted to numeric with ```as.numeric()```.

We can also stay at the sample level and display factor values for each sample (Z matrix) crossed with metadata.

<div class="textbox">
**DIY**: use the ```plot_factor``` function to display factors values and use metadata to color dots.

Using metadata variable significant in the previous correlation analysis for a given factor to color dots should give you a sharp separation for the same factor on this plot.

```{r MOFA_factorValues, fig.width = 10, warning=FALSE}
## Plot factor values
###################################

p <- plot_factor(MOFAobject, 
  factors = c(1:5),
  color_by = "subtype",
  dot_size = 3,        # change dot size
  dodge = T,           # dodge points with different colors
  legend = T,          # remove legend
#  add_violin = T,      # add violin plots (here we have too few points to make violin, let's do boxplots)
#  violin_alpha = 0.25  # transparency of violin plots
  add_boxplot = T
)

# The output of plot_factor is a ggplot2 object that we can edit
print(p)
```
</div>
<br>


##### Combinations of factors

We can study pairwise combination of factors values using the ```plot_factors``` (with an s) function to detect pairwise association between factors and a metadata variable (similar to PCA representation).

<div class="textbox">
**DIY**:try ```plot_factors``` on the *MOFAobject* still coloring dots using metadata.

```{r MOFA_factorCombination, fig.width = 10, fig.height = 9, warning = FALSE, message=FALSE}
## Plot combinations of factor values
###################################

plot_factors(MOFAobject, 
  factors = 1:5,
  color_by = "subtype"
)
```
</div>
<br>



#### Explore individual factors {.tabset}

After this general overview about all factors, you can explore in details a specific factor of interest.

##### Useful functions

**To explore loadings (W matrix) **

* plot_top_weights(): plot top weights for a given factor and view
* plot_weights(): plot all weights for a given factor and view
* plot_weights_heatmap(): plot a heatmap of the weights from multiple factors of a given view

**To explore data**

* plot_data_heatmap(): heatmap of the training data using only top features of a given factor
* plot_data_scatter(): scatterplot of the data using only top features of a given factor

The use of theses functions is illustrated in the next sections.

##### Table of weights/factors

The weights measure how strong each feature relates to each factor : 

* no association = 0
* strong association = high absolute values

The sign of the weight gives the direction of the effect (*eg* a positive weight indicates that the feature has higher levels in the sample with positive factor values).

<div class="textbox">
**DIY**: access the 'W' field of the ```expectation``` *MOFAobject* attribute and display weights associated to a given modality. Use the ```datatable``` function for a nicer visualization.

```{r MOFA_weightTable}
## Weight values
###################################

datatable(MOFAobject@expectations$W$mrna[order(MOFAobject@expectations$W$mrna[,1], decreasing = TRUE ),]) %>% formatRound(c(1:dim(MOFAobject@expectations$W$mrna)[2]),3)
```
</div>
<br>


```{r MOFA_factorToStudy}
## Factor to study
###################################

factorToStudy=1
```

In the next section we will focus on factor `r factorToStudy`, however feel free to choose another factor.

##### Exploration of factor `r factorToStudy`

**Factor values**

We can first display factor `r factorToStudy` values for each sample (Z matrix).

<div class="textbox">
**DIY**: use the ```plot_factor``` function to display factor `r factorToStudy` values and color dots according to this value.

You can also group the sample according to qualitative metadata variable to identify biological association.

```{r MOFA_factorValuePlot}
## Plot factor values
###################################

plot_factor(MOFAobject, 
  factors =  factorToStudy, 
  color_by = paste("Factor",factorToStudy,sep=""),
  group_by = "subtype"
) 
```
</div>
<br>

**Weight of associated features**

Then we can search for features associated to factor `r factorToStudy` in each modality to give a "signature" of the factor. 

<div class="textbox">
**DIY**: visualize weight distribution for factor `r factorToStudy` in the different modalities using ```plot_weights``` function. 

```{r MOFA_distribWeights}
## Plot weight values distribution
###################################

plot_weights(MOFAobject,
 view = "mrna",
 factor = factorToStudy,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```

**NB:** you can specify a given number of features to display in this distribution through 'nfeatures' option.
<br/>
<br/>

**DIY**: use the ```plot_top_weights``` function to focus on top weights (up and down) for factor `r factorToStudy` in each modality. 

The weights measure how strong each feature relates to each factor (no association = 0 ; strong association = high absolute values). The sign of the weight gives the direction of the effect (e.g. a positive weight indicates that the feature has higher levels in the sample with positive factor values).

```{r MOFA_topWeights}
## Plot top weight values
###################################

plot_top_weights(MOFAobject,
 view = "mrna",
 factor = factorToStudy,
 nfeatures = 10,     # Top number of features to highlight
 scale = T           # Scale weights from -1 to 1
)
```
</div>
<br>

**Raw values of associated features**

Instead of looking at weight values associated to each feature, we can also display raw values of features associated to factor `r factorToStudy`.

<div class="textbox">
**DIY**: use the ```plot_data_heatmap``` function to display raw values of the top25 features associated to factor `r factorToStudy` in each modality. You can also cluster features and/or sample based on these values.

```{r MOFA_featureHeatmap}
## Plot feature heatmap
###################################

plot_data_heatmap(MOFAobject, 
  view = "mrna",
  factor = factorToStudy,  
  features = 25,
  cluster_rows = TRUE, cluster_cols = FALSE,
  show_rownames = TRUE, show_colnames = TRUE,
  denoise=TRUE,
  scale = "row"
)
```
</div>
<br>

To focus on a limited number of features you can also display raw values as a scatter plot and cross information with metadata or factors. 

<div class="textbox">
**DIY**: use the ```plot_data_scatter``` function to plot raw values of top4 features having positive factor `r factorToStudy` values in each modality, and color samples based on metadata. 

```{r MOFA_posWeightScatter}
## Plot positive weight feature scatter
###################################

plot_data_scatter(MOFAobject, 
  view = "mrna",
  factor = factorToStudy,  
  features = 4,
  sign = "positive",
  color_by = "subtype",
) + labs(y="mRNA expression")
```
<br/>
<br/>

**DIY**: do the same for top4 features having negative factor `r factorToStudy` values in each modality, and color samples based on metadata. 

```{r MOFA_negWeightScatter}
## Plot negative weight feature scatter
###################################

plot_data_scatter(MOFAobject, 
  view = "mrna",
  factor = factorToStudy,  
  features = 4,
  sign = "negative",
  color_by = "subtype",
) + labs(y="mRNA expression")
```
</div>
<br>

### To go further ...

Before going any further, consider [Regularized Generalized Canonical Correlation Analysis (RGCCA)] if you haven't tried them yet.

#### Play with data pre-processing and MOFA options

You can decide to modify the pre-processing of the data and/or the model/training options to check their effect on the factors.

#### Cluster samples

MOFA is an unsupervised approach, as in PCA you can try to cluster sample based on their projection in the MOFA factor space.

<div class="textbox">
**DIY**: The function ```cluster_samples``` is based on [k-means](https://en.wikipedia.org/wiki/K-means_clustering#/media/File:K-means_convergence.gif). Use this function to function identify clusters of samples based on their projection all infered factors. Then use the ```plot_factors``` function to display the result of this clustering by coloring samples according to clusters. 

```{r MOFA_cluster, message = FALSE, warning = FALSE, class.source = 'fold-hide'}
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                                 Clustering                                 ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##

## cluster samples based on factor projection
clusters <- cluster_samples(MOFAobject, k=4, factors = "all")

## plot sample projection on factor space with cluster information
plot_factors(MOFAobject,
  factors=1:2,
  color_by=clusters$cluster)

```
</div>
<br>


<div class="textbox">
**DIY**: Try different value of ```k```, the number of clusters to extract. With the ```proxy::dist``` function, compute the distance from each sample to the barycenter of its clusters and compute the averaged [silhouette score](https://en.wikipedia.org/wiki/Silhouette_(clustering)) out of each clustering. Find the number of extracted cluster leading to the best silhouette score.

```{r MOFA_cluster_silhouette, message = FALSE, warning = FALSE, class.source = 'fold-hide'}

if (!require("proxy")){
  install.packages("proxy")
}

if (!require("Silhouette")){
  install.packages("Silhouette")
}

res_silhouette  = c()
tested_k_values = 2:10
for (k in tested_k_values){
    ## cluster samples based on factor projection
    clusters       = cluster_samples(MOFAobject, k=k, factors = "all")
    dist_matrix    = proxy::dist(MOFAobject@expectations$Z$group1, clusters$centers)
    sil            = Silhouette(dist_matrix)
    res_silhouette = c(res_silhouette, (mean(sil$sil_width)))
}
names(res_silhouette) = tested_k_values
print(res_silhouette)

```
</div>
<br>


<div class="textbox">
**DIY**: Re-lauch the clustering on the optimal number of clusters and compute both a contingency table between TCGA' subtypes and the inferred clusters and the [Rand Index](https://en.wikipedia.org/wiki/Rand_index) between the two. 

```{r MOFA_cluster_contingency_RI, message = FALSE, warning = FALSE, class.source = 'fold-hide'}

if (!require("aricode")){
  install.packages("aricode")
}

clusters = cluster_samples(MOFAobject, k=tested_k_values[which.max(res_silhouette)], factors = "all")

print("Contingency Table")
sortPairs(clinical_data$subtype, clusters$cluster, spMat = T)$spMat

print("Rand Index")
RI(clinical_data$subtype, clusters$cluster)

```
</div>
<br>


#### Evaluate stability of number of extracted factors


<div class="textbox">
**DIY**: MOFA is not randomly initialized as MOFA's latent space $\mathbf{Z}$ and weight matrix $\mathbf{W}$ are initialized based on PCA. So if you launch 2 subsequent MOFA without dropping factors, results will be exactly identicals. However, when you drop factors, if at a specific iteration, several factors are to be dropped, a single one is and picked randomly. Using this variability, apply multiple MOFA models with different seeds. Each time allow factors to be dropped.

```{r MOFA_nbFactors_stability, message = FALSE, warning = FALSE, results = "hide", class.source = 'fold-hide'}

list_MOFA_models          = list()
number_of_initializations = 10

for (i in 1:number_of_initializations){
  prepare_MOFAobject@training_options$seed = 53*i
  list_MOFA_models[[i]]                    = run_mofa(object = prepare_MOFAobject, use_basilisk = TRUE) 
}

```
</div>
<br>

<div class="textbox">
**DIY**: Use the function ```compare_factors``` to compare these models according to their factors.

```{r MOFA_nbFactors_stability_compare, message = FALSE, warning = FALSE, class.source = 'fold-hide', fig.width = 15, fig.height = 15}

compare_factors(list_MOFA_models)

```
</div>
<br>


#### Ontology enrichment on features

<div class="textbox">
**DIY**: You can run gene set enrichment analysis using group of genes having a high weight for a given factor to get indications about the biological meaning of factors.

```{r MOFA_enrichment, message = FALSE, warning = FALSE, class.source = 'fold-hide', fig.width = 20, fig.height = 10}
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##
##                                 Enrichment                                 ##
##––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––––##

if (!require("dplyr")){
  install.packages("dplyr")
}

if (!require("fgsea")){
  install.packages("fgsea")
}

annot_table  = read.csv("../data/TCGA/annotation_table_mrna.csv", row.names = 1)
# check for the first component

res = c()
for (comp in 1:dim(MOFAobject@expectations$W$mrna)[2]){
  a1           = MOFAobject@expectations$W$mrna[, comp]
  
  ## Check order
  all(annot_table$symbols_old == gsub("^mRNA_", "", names(a1)))
  names(a1)     = annot_table$entrez
  kegg_pathways = readRDS("../data/TCGA/kegg_pathways.rds")
  
  kegg_hsa_name <- kegg_pathways %>%
    dplyr::select(pathway_id, pathway_name) %>%
    dplyr::distinct()
  
  sum(names(a1) %in% kegg_pathways$entrez_gene)
  
  fgsea_res <- fgseaMultilevel(
    pathways = split(kegg_pathways$entrez_gene, kegg_pathways$pathway_id),
    stats = a1**2,
    minSize = 15,
    maxSize = 500,
    eps = 0
  )
  
  fgsea_res <- fgsea_res |>
    left_join(kegg_hsa_name, by = c("pathway" = "pathway_id")) %>%
    relocate(pathway_name, .after = 1)
  
  res = rbind(res, fgsea_res$pval)
}

res_pval_adjusted = matrix(p.adjust(c(res), method = "BH"), nrow = dim(res)[1], ncol = dim(res)[2], byrow = F)
colnames(res_pval_adjusted) = fgsea_res$pathway_name
rownames(res_pval_adjusted) = paste0("Factor_", 1:dim(res_pval_adjusted)[1])

res_pval_adjusted[which(res_pval_adjusted > 0.05, arr.ind = T)] = NA

pheatmap::pheatmap(-log10(res_pval_adjusted), cluster_cols = F, cluster_rows = F)


```
</div>
<br>

