1.1 Option 1: base R approach

If using the base R approach, you simply need to provide a list of matrices where features are rows and samples are columns. Importantly, the samples need to be aligned. Missing values/assays should be filled with NAs.

data("CLL_data")
MOFAobject <- createMOFAobject(CLL_data)

## Creating MOFA object from list of matrices,
##  please make sure that samples are columns and features are rows...

MOFAobject

## Untrained MOFA model with the following characteristics: 
##   Number of views: 4 
##   View names: Drugs Methylation mRNA Mutations 
##   Number of features per view: 310 4248 5000 69
##   Number of samples: 200

1.2 Option 2: Bioconductor approach

If using the Bioconductor approach, you need to provide or create a MultiAssayExperiment object and then use it to build the MOFA object. For example, starting from a list of matrices where features are rows and samples are columns, this can be easily constructed as follows:

# Load data
# import list with mRNA, Methylation, Drug Response and Mutation data. 
data("CLL_data") 

# check dimensionalities, samples are columns, features are rows
lapply(CLL_data, dim) 

## $Drugs
## [1] 310 200
## 
## $Methylation
## [1] 4248  200
## 
## $mRNA
## [1] 5000  200
## 
## $Mutations
## [1]  69 200

# Load sample metadata: Sex and Diagnosis
data("CLL_covariates")
head(CLL_covariates)

##      Gender Diagnosis
## H045      m       CLL
## H109      m       CLL
## H024      m       CLL
## H056      m       CLL
## H079      m       CLL
## H164      f       CLL

# Create MultiAssayExperiment object 
mae_CLL <- MultiAssayExperiment(
  experiments = CLL_data, 
  colData = CLL_covariates
)

# Build the MOFA object
MOFAobject <- createMOFAobject(mae_CLL)
MOFAobject

## Untrained MOFA model with the following characteristics: 
##   Number of views: 4 
##   View names: Drugs Methylation mRNA Mutations 
##   Number of features per view: 310 4248 5000 69
##   Number of samples: 200

1.3 Overview of training data

The function plotDataOverview can be used to obtain an overview of the data stored in the object for training. For each sample it shows in which views data are available. The rows are the different views and columns are samples. Missing values are indicated by a grey bar.

plotDataOverview(MOFAobject)

2.1 Define options

DataOptions <- getDefaultDataOptions()
DataOptions 

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

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

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

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

2.2 Prepare MOFA

prepareMOFA internally performs a set of sanity checks and fills the DataOptions, TrainOptions and ModelOptions slots of the MOFAobject

MOFAobject <- prepareMOFA(
  MOFAobject, 
  DataOptions = DataOptions,
  ModelOptions = ModelOptions,
  TrainOptions = TrainOptions
)

## Checking data options...

## Checking training options...

## Checking model options...

Optionally, we can choose to regress out some (technical) covariates before training, using a simple linear model. For example, here we can choose to remove the effect of sex. Ideally, all undesired sources of variation should be removed a priori from the model. The reason ebing that, if strong technical factors exist, the model will “focus” on capturing the variability driven by the technical factors, and small sources of biological variability could be missed.
(Note: uncomment and running the function below will lead to a slight modification of the results)

# MOFAobject <- regressCovariates(
#   object = MOFAobject,
#   views = c("Drugs","Methylation","mRNA"),
#   covariates = MOFAobject@InputData$Gender
# )

2.3 Run MOFA

Now we are ready to train the MOFAobject, which is done with the function runMOFA. This step can take some time (around 15 min with default parameters). For illustration we provide an existing trained MOFAobject.
IMPORTANT NOTE: The software has evolved since the original publication and the results are not 100% reproducible with the last versions. Yet, the output should be very similar (if not improved) to the pre-existent model.

MOFAobject <- runMOFA(MOFAobject)

# Loading an existing trained model
filepath <- system.file("extdata", "CLL_model.hdf5", package = "MOFAdata")

MOFAobject <- loadModel(filepath, MOFAobject)

## [1] "Could not load the feature intercepts. Your features will be centered to zero mean"

MOFAobject

## Trained MOFA model with the following characteristics:
##   Number of views: 4 
##  View names: Drugs Methylation Mutations mRNA 
##   Number of features per view: 310 4248 69 5000 
##   Number of samples: 200 
##   Number of factors: 10

3.1 Part 1: Disentangling the heterogeneity: calculation of variance explained by each factor in each view

This is done by calculateVarianceExplained (to get the numerical values) and plotVarianceExplained (to get the plot). The resulting figure gives an overview of which factors are active in which view(s). If a factor is active in more than one view, this means that is capturing shared signal (co-variation) between features of different data modalities.
Here, for example Factor 1 is active in all data modalities, while Factor 4 is specific to mRNA.

# Calculate the variance explained (R2) per factor in each view 
r2 <- calculateVarianceExplained(MOFAobject)
r2$R2Total

##       Drugs Methylation   Mutations        mRNA 
##   0.4090497   0.2341916   0.2415335   0.3822602

# Variance explained by each factor in each view
head(r2$R2PerFactor)

##            Drugs  Methylation    Mutations       mRNA
## LF1 0.1507029688 1.732398e-01 1.537989e-01 0.07519230
## LF2 0.0350660758 5.844575e-03 8.195112e-02 0.04699171
## LF3 0.1123216630 3.789900e-05 2.242679e-05 0.01488067
## LF4 0.0001812843 1.925074e-04 1.934564e-07 0.09061719
## LF5 0.0605883335 1.837185e-05 3.587202e-06 0.02860510
## LF6 0.0340251701 1.181410e-05 3.004497e-06 0.04850234

# Plot it
plotVarianceExplained(MOFAobject)

3.2 Part 2: Characterisation of individual factors

plotWeightsHeatmap(
  MOFAobject, 
  view = "Mutations", 
  factors = 1:5,
  show_colnames = FALSE
)

plotWeights(
  MOFAobject, 
  view = "Mutations", 
  factor = 1, 
  nfeatures = 5
)

plotWeights(
  MOFAobject, 
  view = "Mutations", 
  factor = 1, 
  nfeatures = 5,
  manual = list(c("BRAF"),c("MED12")),
  color_manual = c("red","blue")
  
)

plotTopWeights(
  MOFAobject, 
  view="Mutations", 
  factor=1
)

plotTopWeights(
  MOFAobject, 
  view = "mRNA", 
  factor = 1
)

plotDataHeatmap(
  MOFAobject, 
  view = "mRNA", 
  factor = 1, 
  features = 20, 
  show_rownames = FALSE
)

# Load reactome annotations
data("reactomeGS") # binary matrix with feature sets in rows and features in columns

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

## Doing Feature Set Enrichment Analysis with the following options...

## View: mRNA

## Factors: LF1 LF2 LF3 LF4 LF5 LF6 LF7 LF8 LF9 LF10

## Number of feature sets: 501

## Local statistic: loading

## Transformation: abs.value

## Global statistic: mean.diff

## Statistical test: parametric

plotEnrichmentBars(gsea, alpha=0.01)

interestingFactors <- 4:5

fseaplots <- lapply(interestingFactors, function(factor) {
  plotEnrichment(
    MOFAobject,
    gsea,
    factor = factor,
    alpha = 0.01,
    max.pathways = 10 # The top number of pathways to display
  )
})

cowplot::plot_grid(fseaplots[[1]], fseaplots[[2]],
                   ncol = 1, labels = paste("Factor", interestingFactors))

3.3 Ordination of samples by factors to reveal clusters and gradients in the sample space

Samples can be visualized along factors of interest using the plotFactorScatter function. We can use features included in the model (such as IGHV or trisomy12) to color or shape the samples by. Alternatively, external covariates can also be used for this purpose.

plotFactorScatter(
  MOFAobject,
  factors = 1:2,
  color_by = "IGHV",      # color by the IGHV values that are part of the training data
  shape_by = "trisomy12"  # shape by the trisomy12 values that are part of the training data
)

Here we find again a clear separation of samples based on their IGHV status (color) along Factor 1 and by the absence or prescence of trisomy 12 (shape) along Factor 2 as indicated by the corresponding factor weights in the Mutations view.

An overview of pair-wise sctterplots for all or a subset of factors is produced by the plotFactorScatters function

plotFactorScatters(
  MOFAobject,
  factors = 1:3,
  color_by = "IGHV"
)

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

A single factor can be visualised using the plotFactorBeeswarm function

plotFactorBeeswarm(
  MOFAobject,
  factors = 1,
  color_by = "IGHV"
)

3.4 Customized analysis

For customized exploration of weights and factors, you can directly fetch the variables from the model using ‘get’ functions: getWeights, getFactors and getTrainData:

MOFAweights <- getWeights(
  MOFAobject, 
  views = "all", 
  factors = "all", 
  as.data.frame = TRUE    # if TRUE, it outputs a long dataframe format. If FALSE, it outputs a wide matrix format
)
head(MOFAweights)

##   feature factor         value  view
## 1 D_001_1    LF1  0.0002667769 Drugs
## 2 D_001_2    LF1  0.0019258319 Drugs
## 3 D_001_3    LF1  0.0356598749 Drugs
## 4 D_001_4    LF1  0.0292638984 Drugs
## 5 D_001_5    LF1  0.0158091594 Drugs
## 6 D_002_1    LF1 -0.0161823261 Drugs

MOFAfactors <- getFactors(
  MOFAobject, 
  factors = c(1,2),
  as.data.frame = FALSE   # if TRUE, it outputs a long dataframe format. If FALSE, it outputs a wide matrix format
)
head(MOFAfactors)

##            LF1         LF2
## H045  1.277135  0.07556906
## H109  1.292226 -1.93887825
## H024 -1.438207  0.72267340
## H056 -1.102120  0.65452538
## H079  1.037262  0.51368436
## H164  1.563211  0.92332140

MOFAtrainData <- getTrainData(
  MOFAobject,
  as.data.frame = TRUE, 
  views = "Mutations"
)
head(MOFAtrainData)

##        view    feature sample value
## 1 Mutations gain2p25.3   H045     0
## 2 Mutations   gain3q26   H045     0
## 3 Mutations  del6p21.2   H045     0
## 4 Mutations    del6q21   H045     0
## 5 Mutations    del8p12   H045     0
## 6 Mutations   gain8q24   H045     1

4.1 Prediction of views

With the predict function, full views can be predicted based on the MOFA model with all or a subset of factors.

predictedDrugs <- predict(
  MOFAobject,
  view = "Drugs",
  factors = "all"
)[[1]]


# training data (incl. missing values)
drugData4Training <- getTrainData(MOFAobject, view="Drugs")[[1]]
pheatmap::pheatmap(drugData4Training[1:40,1:20],
         cluster_rows = FALSE, cluster_cols = FALSE,
         show_rownames = FALSE, show_colnames = FALSE)

# predicted data
pheatmap::pheatmap(predictedDrugs[1:40,1:20],
         cluster_rows = FALSE, cluster_cols = FALSE,
         show_rownames = FALSE, show_colnames = FALSE)

4.2 Imputation of missing observations

With the impute function all missing values are imputed based on the MOFA model. The imputed data is then stored in the ImputedData slot of the MOFAobject and can be accessed via the getImputedData function.

MOFAobject <- impute(MOFAobject)
imputedDrugs <- getImputedData(MOFAobject, view="Drugs")[[1]]

# training data (incl. missing values)
pheatmap::pheatmap(drugData4Training[1:40,1:20],
         cluster_rows = FALSE, cluster_cols = FALSE,
         show_rownames = FALSE, show_colnames = FALSE)

# imputed data
pheatmap::pheatmap(imputedDrugs[1:40,1:20],
         cluster_rows = FALSE, cluster_cols = FALSE,
         show_rownames = FALSE, show_colnames = FALSE)

4.3 Clustering of samples based on latent factors

Samples can be clustered according to their values on some or all latent factors using the clusterSamples function. Clusters can for example be visualised using the plotFactorScatters function

set.seed(1234)
clusters <- clusterSamples(
  MOFAobject, 
  k = 2,        # Number of clusters for the k-means function
  factors = 1   # factors to use for the clustering
)

plotFactorScatter(
  MOFAobject, 
  factors = 1:2, 
  color_by = clusters
)