介绍

MOFA是一种因子分析模型，它以完全无监督的方式为多元数据集的集成提供了一般框架。直观地说，MOFA可以被视为主成分分析（PCA）到多组学数据的通用和统计严格的推广。给定几个数据矩阵，在相同或重叠的样本集上测量多个'组学数据类型，MOFA根据（隐藏）因子推断可解释的低维数据表示。这些学习因素代表了数据模态变异的驱动源，从而有助于识别细胞状态或疾病亚组。

经过训练，模型输出可用于一系列下游分析，包括因子空间中样本的可视化，使用（基因集）富集分析的因子的自动注释，异常值的识别（例如由于样本交换）和对缺失值的估算。

简单来说MOFA支持使用不同组学的数据联合分析以发现样本间的异质性，它容许样本不同组学数据的缺失，MOFA将高维多组学数据集的异质性解构为一组捕捉全球变异来源的潜在因素。重要的是，这些因素可以在不同的组学中具有不同的活动模式。例如，批次效应可能影响RNA数据，但不影响甲基化数据。

步骤

MOFA的工作流程包括两个步骤：
（1）拟合步骤：用多组学数据训练模型，将异质性分解为少数潜在因子。
（2）下游分析：一旦推断出因子，就需要通过查看相应的权重，进行（基因集）富集分析，绘制因子，与已知协变量相关因子等来表征变异的技术或生物来源。此外，可以使用潜在因素来估算缺失值和预测临床结果。

1 数据导入

将不同组学数据组合成list，之后使用createMOFAobject函数创建MOFA对象，或者可以提前创建MultiAssayExperiment对象，然后用它来构建MOFA对象。

#使用内置CLL数据集，四种数据集的list，包括mRNA，甲基化，突变信息以及药物处理信息。
data("CLL_data")

class(CLL_data)
[1] "list"

MOFAobject <- createMOFAobject(CLL_data)

#或者使用MultiAssayExperiment创建MOFA对象

data("CLL_covariates")
head(CLL_covariates) #样本表型数据

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

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

对数据集进行总览

plotDataOverview(MOFAobject)

2 安装模型

对MOFA对象的参数进行一些调整，首先进行模型选项

DataOptions <- getDefaultDataOptions()
DataOptions 

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

ModelOptions <- getDefaultModelOptions(MOFAobject)
ModelOptions$numFactors <- 25
ModelOptions
#likelihood代表数据类型，高斯用于连续数据，伯努利用于二项数据，泊松用于计数数据，通常MOFA会自动识别数据类型
#numFactors代表因子数量（默认值是样本数量的0.5倍）
#sparsity代表稀疏性，默认TRUE

## $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
#包括最大迭代次数，数据容忍度
#DropFactorThreshold：超参数自动学习基于最小方差解释标准的因子数。DropFactorThreshold将删除解释所有视图中变化小于
#分数的因素。例如，值0.01意味着将丢弃解释所有视图中小于1％方差的因子。默认情况下，它为零，这意味着除非它们根本不
#解释任何差异，否则将保留所有因子。

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

3 prepareMOFA函数准备MOFA

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

4 运行MOFA

这里需要注意需要安装mofapy的python模块，使用reticulate指定python路径以及conda环境（注意python2环境）
MOFAobject <- runMOFA(MOFAobject)

模型分析

经过训练，我们可以对结果进行探索。在这里，作者提供了一个半自动化的管道来解开和表征所有已识别的变异来源。

第1部分：解开异质性

每组中每个因素解释的方差计算。这可能是MOFA生成的最重要的图，因为它在一个图中总结了数据集的整个异质性。在这里，我们可以看到在哪个组中，一个因素解释了变异，这可以通过调查这些权重来指导因素的进一步表征。

第2部分：个别因素的表征

• 检查具有最高负载的顶部特征：负载是特征重要性的度量，因此具有高负载的特征是驱动由该因子捕获的异质性的特征。
• 特征集富集分析（其中存在集合注释，例如mRNA视图的基因集）。
• 通过因子排序样本以揭示聚类和/或梯度：这与传统上使用主成分分析或t-SNE进行的类似。

还可以进行其他分析，包括估算缺失值和聚类样本。此外，这些因子可用于进一步分析，例如作为预测因子，例如用于预测临床结果或对患者进行分类，或控制不需要的变异来源。

#查看组间方差
r2 <- calculateVarianceExplained(MOFAobject)
r2$R2Total

#绘图展示
plotVarianceExplained(MOFAobject)

概览所有因子在不同组中的权重

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

我们观察到因子1和2具有大的非零权重。为了更详细地探索给定因子，我们可以使用plotWeights函数绘制单个因子的所有权重。例如，这里我们在Mutation数据中绘制因子1的所有权重。

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")
  
)

这里可以观察到IGHV和KLHL6具有较大的权重，随后使用plotTopWeights展示

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

同样我们还可以对不同view进行展示，只需要指定view即可。

#使用热图形式展示因子1中基因的表达模式
plotDataHeatmap(
  MOFAobject, 
  view = "mRNA", 
  factor = 1, 
  features = 20, 
  show_rownames = FALSE
)

随后进行注释

在这里，我们实现了一个runEnrichmentAnalysis从PCGSE包派生的特征集富集分析方法的函数。

该函数的输入是MOFA训练模型（MOFAmodel），执行特征集丰富的因子（字符向量），特征集（二进制矩阵）和关于如何执行分析的一组选项，另见文件 runEnrichmentAnalysis

#首先载入二进制基因注释矩阵
data("reactomeGS")

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

#画图查看因子富集情况
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))

这表明，因子4正在捕获与免疫反应相关的变异（可能是由于样品的T细胞污染），因子5与应激反应的差异有关，正如我们的论文中所讨论的。

可以使用该plotFactorScatter函数沿着感兴趣的因素可视化样本。我们可以使用模型中包含的特征（例如IGHV）来对样本进行着色。或者，外部协变量也可用于此目的。

#单组展示
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
)

#多组展示
plotFactorScatters(
  MOFAobject,
  factors = 1:3,
  color_by = "IGHV"
)

对于权重和因素的定制化的探索，你可以用“get”功能直接读取从模型中的变量：getWeights，getFactors及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)

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)

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

通过该predict功能，可以基于具有全部或部分因子的MOFA模型来预测完整视图。

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)

或者使用impute功能，可以根据MOFA模型估算所有缺失值。然后将插补数据存储在ImputedData slotMOFA对象中，并可通过该getImputedData功能访问。

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)

还可以使用clusterSamples函数根据它们在一些或所有潜在因子上的值对样本进行聚类。例如，可以使用该plotFactorScatters函数来可视化集群。

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
)