Define data options

The most important options the user needs to define are:

• scaleViews: logical indicating whether to scale views to have unit variance. As long as the scale of the different data sets is not too high, this is not required. Default is FALSE.

• removeIncompleteSamples: logical indicating whether to remove samples that are not profiled in all omics. The model can cope with missing assays, so this option is not required. Default is FALSE.

DataOptions <- getDefaultDataOptions()
DataOptions 

## $scaleViews
## [1] FALSE
## 
## $removeIncompleteSamples
## [1] FALSE

Define model options

Next, we define model options. The most important are:

• numFactors: number of factors (default is 0.5 times the number of samples). By default, the model will only remove a factor if it explains exactly zero variance in the data. You can increase this threshold on minimum variance explained by setting TrainOptions$dropFactorThreshold to a value higher than zero.

• likelihoods: likelihood for each view. Usually we recommend gaussian for continuous data, bernoulli for binary data and poisson for count data. By default, the model tries to guess it from the data.

• sparsity: do you want to use sparsity? This makes the interpretation easier so it is recommended (Default is TRUE).

ModelOptions <- getDefaultModelOptions(MOFAobject)
ModelOptions$numFactors <- 25
ModelOptions

## $likelihood
##       Drugs Methylation        mRNA   Mutations 
##  "gaussian"  "gaussian"  "gaussian" "bernoulli" 
## 
## $numFactors
## [1] 25
## 
## $sparsity
## [1] TRUE

Define training options

Next, we define training options. The most important are:

• maxiter: maximum number of iterations. Ideally set it large enough and use the convergence criterion TrainOptions$tolerance.

• tolerance: convergence threshold based on change in the evidence lower bound. For an exploratory run you can use a value between 1.0 and 0.1, but for a “final” model we recommend a value of 0.01.

• DropFactorThreshold: hyperparameter to automatically learn the number of factors based on a minimum variance explained criteria. Factors explaining less than DropFactorThreshold fraction of variation in all views will be removed. For example, a value of 0.01 means that factors that explain less than 1% of variance in all views will be discarded. By default this it zero, meaning that all factors are kept unless they explain no variance at all.

TrainOptions <- getDefaultTrainOptions()
# Automatically drop factors that explain less than 2% of variance in all omics
TrainOptions$DropFactorThreshold <- 0.02
TrainOptions$seed <- 2017
TrainOptions

## $maxiter
## [1] 5000
## 
## $tolerance
## [1] 0.1
## 
## $DropFactorThreshold
## [1] 0.02
## 
## $verbose
## [1] FALSE
## 
## $seed
## [1] 2017

Inspection of top weighted features in the active views

To explore a given factor in more detail we can plot all weights for a single factor using the plotWeights function. For example, here we plot all weights from Factor 1 in the Mutation data. With nfeatures we can set how many features should be labelled (manual let’s you specify feautres manually to be labelled in the plot.)

plotWeights(
  
  MOFAobject, 

  view = "Mutations", 
  
  factor = 1, 
  
  nfeatures = 5

)

Features with large (absolute) weight on a given factor follow the pattern of covariation associated with the factor. In our case, this reveals a strong link to the IGHV status, hence recovering an important clinical marker in CLL (see our paper for details on the biological interpretation.) Note that the sign of the weight can only be interpreted relative to the signs of other weights and the factor values.

If you are only interested in looking at only the top features you can use the plotTopWeights function.

From the previous plots, we can clearly see that Factor 1 is associated to IGHV status. As the variance decomposition above told us that this factor is also relevant on all the other data modalities we can investigate its weights on other modalities, e.g. mRNA, to make connections of the IGHV-linked axes of variation to other molecular layers.

plotTopWeights(MOFAobject, "mRNA", 1)

Feature set enrichment analysis in the active views

Sometimes looking at the loadings of single features can be challenging, and often the combination of signal from functionally related sets of features (i.e. gene ontologies) is required.

To have an overview of the biological processes associated to the various factors run feature set enrichment analysis method (runEnrichmentAnalysis) derived from the PCGSE package.

The input of this function is a MOFA trained model (MOFAmodel), the factors for which to perform feature set enrichment (a character vector), the feature sets (a binary matrix) and a set of options regarding how the analysis should be performed, see also documentation of runEnrichmentAnalysis Let’s consider reactome annotations, which are contained in the package.

# Load reactome annotations
# binary matrix with feature sets in rows and feautres in columns

data("reactomeGS")

# perform enrichment analysis
fsea.out <- runEnrichmentAnalysis(
  
  MOFAobject,
  
  view = "mRNA",
  
  feature.sets = reactomeGS,
  
  alpha = 0.01

)

The next step is to visualise the results of the Gene Set Enrichment Analysis. There are two default plots:

1. General Overview: Barplot with number of enriched gene sets per view

plotEnrichmentBars(fsea.out, alpha=0.01)

From this we find enriched at a FDR of 1% gene sets on Factors 3-6 and 8. To look into which gene sets these are we can choose a factor of interest and visualize the most enriched gene sets as follows:

2. Factor-specific:

interestingFactors <- 4:5
for (factor in interestingFactors) {
  lineplot <- plotEnrichment(MOFAobject,
  fsea.out,
  factor = factor,
  alpha = 0.01,
  max.pathways = 10)
  
  print(lineplot)
}

This shows us that Factor 4 is capturing variation related to immune response (possibly due to T-cell contamination of the samples) and Factor 5 is related to differences in stress response, as discussed in our paper.