内容

本插图展示了如何使用MOFA，包括初始化，训练和下游分析。为了说明，我们使用MOFA出版物中使用的CLL数据。

库(MultiAssayExperiment)库(MOFA)库(MOFAdata)

1步骤1:加载数据并创建MOFA对象

有两种方式输入数据到MOFA:

1．1选项1:基本R方法

如果使用基本R方法，您只需要提供一个矩阵列表，其中特征是行，样本是列。重要的是，样本需要对齐。缺失值/检测应用NAs填充。

data("CLL_data") MOFAobject <- createMOFAobject(CLL_data)
##从矩阵列表中创建MOFA对象##请确保样本存储在列中，特征存储在行中…
MOFAobject
##未经训练的MOFA模型，具有以下特征:##视图数量:4 ##视图名称:药物甲基化mRNA突变##每个视图的特征数量:310 4248 5000 69 ##样本数量:200

1．2方案2:生物导体方法

如果使用Bioconductor方法，则需要提供或创建MultiAssayExperiment对象，然后使用它来构建MOFA对象。例如，从一个矩阵列表开始，其中特征是行，样本是列，这可以很容易地构造如下:

导入列表，包含mRNA、甲基化、药物反应和突变数据。数据("CLL_data") #检查维度，样本是列，特征是行
##[1] 310 200 ## ## $甲基化## [1]4248 200 ## ## $mRNA ##[1] 5000 200 ## ## $突变## [1]69 200
#加载样本元数据:性别和诊断数据(" cll_covariables ")头(cll_covariables)
性别诊断## H045 m CLL ## H109 m CLL ## H024 m CLL ## H056 m CLL ## H079 m CLL ## H164 f CLL
创建MultiAssayExperiment对象mae_CLL <- MultiAssayExperiment(experiments = CLL_data, colData = cll_covariables) #创建MOFA对象MOFAobject <- createMOFAobject(mae_CLL) MOFAobject
##未经训练的MOFA模型，具有以下特征:##视图数量:4 ##视图名称:药物甲基化mRNA突变##每个视图的特征数量:310 4248 5000 69 ##样本数量:200

1.3培训数据概述

这个函数plotDataOverview可用于获取存储在用于训练的对象中的数据的概览。对于每个示例，它显示哪些视图数据是可用的。行是不同的视图，列是示例。缺失的值用灰色条表示。

plotDataOverview (MOFAobject)

2步骤2:拟合MOFA模型

下一步是拟合模型。这部分管道是用Python实现的，所以首先要确保安装了相应的包(如果有问题，请参阅安装说明并阅读FAQ)。

2.1定义选项

2.1.1定义数据选项

用户需要定义的最重要的选项是:

  • scaleViews:逻辑指示是否缩放视图以具有单位方差。只要不同数据集的规模不是太高，这是不需要的。默认是假．

  • removeIncompleteSamples:逻辑指示是否删除在所有组学中未分析的样本。该模型可以处理缺失的分析，因此该选项不是必需的。默认是假．

DataOptions <- getDefaultDataOptions() DataOptions
## ## $removeIncompleteSamples ## [1] FALSE

2.1.2定义模型选项

接下来，我们定义模型选项。最重要的是:

  • numFactors:因子数量(默认为样本数量的0.5倍)。默认情况下，模型只会删除一个因子，如果它恰好解释了数据中的零方差。您可以通过设置来增加最小方差的阈值TrainOptions dropFactorThreshold美元到一个大于零的值。

  • 可能:各视图的可能性。对于连续数据，我们通常推荐高斯分布，对于二进制数据推荐伯努利分布，对于计数数据推荐泊松分布。默认情况下，模型尝试从数据中猜测。

  • 稀疏:你想使用稀疏性吗?这使得解释更容易，所以推荐使用(默认为真正的）.

ModelOptions <- getDefaultModelOptions(MOFAobject
##药物甲基化mRNA突变##“高斯”“高斯”“高斯”“伯努利”## ## $numFactors ##[1] 25 ## ## $稀疏## [1]TRUE

2.1.3定义培训选项

接下来，我们定义培训选项。最重要的是:

  • 麦克斯特:最大迭代次数。理想情况下，设置它足够大，并使用收敛准则美元TrainOptions宽容．

  • 宽容:基于证据下界变化的收敛阈值。对于探索性运行，您可以使用1.0到0.1之间的值，但对于“最终”模型，我们建议使用0.01的值。

  • DropFactorThreshold:超参数，自动学习基于最小方差解释标准的因素数量。解释不足的因素DropFactorThreshold所有视图的部分变化将被删除。例如，值为0.01意味着在所有视图中解释方差小于1%的因素将被丢弃。默认情况下，它为零，这意味着所有的因素都被保留，除非它们解释不了任何方差。

TrainOptions <- getDefaultTrainOptions() #自动删除在所有组学中解释方差小于2%的因素
# # $麦克斯特# # # # # # 5000[1]公差0.1 # # # # # # [1]DropFactorThreshold # #美元[1]0.02 # # # # $ verbose # #[1]假# # # # $种子# # 2017年[1]

2．2准备遐

prepareMOFA在内部执行一组健全性检查并填充DataOptions，TrainOptions而且ModelOptions槽位号MOFAobject

MOFAobject <- prepareMOFA(MOFAobject, DataOptions = DataOptions, ModelOptions = ModelOptions, TrainOptions = TrainOptions)
##检查数据选项…
##检查培训选项…
##检查模型选项…

可选地，我们可以选择在训练之前使用简单的线性模型回归出一些(技术)协变量。例如，这里我们可以选择去除性别的影响。理想情况下，所有不希望出现的变化源都应该从模型中先验地去除。原因在于，如果存在强大的技术因素，该模型将“集中”于捕捉由技术因素驱动的变异性，而可能会错过生物变异性的小来源。
(注意:取消注释并运行下面的函数将导致对结果的轻微修改)

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

2．3运行遐

现在我们准备训练MOFAobject，这是用函数完成的runMOFA．这个步骤可能需要一些时间(默认参数下大约15分钟)。为了说明，我们提供了一个现有的训练有素的MOFAobject．
重要提示:该软件自最初的发布以来一直在发展，结果并不是100%可重复的最后版本。然而，输出应该与先前存在的模型非常相似(如果没有改进的话)。

MOFAobject <- runMOFA(MOFAobject)
#加载一个已训练好的模型文件路径<- system. #文件(“extdata”、“CLL_model。hdf5"， package = "MOFAdata") MOFAobject <- loadModel(filepath, MOFAobject)
##[1] "无法加载功能拦截。你的特征会集中在均值为零的地方"
MOFAobject
##训练的MOFA模型具有以下特征:##视图数量:4 ##视图名称:Drugs Methylation Mutations mRNA ##每个视图的特征数量:310 4248 69 5000 ##样本数量:200 ##因子数量:10

3.步骤3:分析一个训练好的MOFA模型

训练结束后，我们可以探索MOFA的结果。在这里，我们提供了一个半自动化的管道来分解和描述所有已识别的变化来源(因素)。

第1部分:理清异构性

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

第二部分:个体因素的特征

其他分析，包括缺失值的imputation和样本的聚类也可用。请参见下面对这些功能的简短说明。此外，这些因子可用于进一步分析，例如作为预测因子，例如用于预测临床结果或对患者进行分类，或用于控制不需要的变异源。说明这一点的小插图很快就会出现。

3．1第一部分:解开异质性:计算方差解释的每个因素在每个观点

这是由calculateVarianceExplained(以获得数值)和plotVarianceExplained(得到情节)。结果图给出了哪些因素在哪些视图中处于活动状态的概述。如果一个因子在多个视图中处于活动状态，这意味着它正在捕获不同数据形态特征之间的共享信号(共变)。
例如，在这里，因子1在所有数据模式中都是活跃的，而因子4是特定于mRNA的。

#计算每个视图中每个因素的方差解释(R2) R2 <- calculateVarianceExplained(MOFAobject) R2 $R2Total
药物甲基化突变mRNA ## # 0.4090497 0.2341916 0.2415335 0.3822602
#由每个视图头中的每个因素解释的方差(r2$R2PerFactor)
##药物甲基化突变mRNA ## LF1 0.1507029688 1.732398e-01 1.537989e-01 0.07519230 ## LF2 0.0350660758 5.844575e-03 8.19511216630 3.789900e-05 2.242679e-05 0.01488067 ## LF4 0.0001812843 1.925074e-04 1.934564e-07 0.09061719 ## LF5 0.0605883335 1.83718202e -05 3.004497e-06 0.04850234
# Plot it plotVarianceExplained(MOFAobject)

3.2第二部分:个体因素的特征

3.2.1之上活动视图中顶部加权特征的检查

要获得给定视图中所有因素权重的概览，可以使用plotWeightsHeatmap函数。例如，这里我们绘制了突变数据中所有因素的所有权重:

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

我们观察到因子1和因子2具有较大的非零权重。为了更详细地研究给定的因素，我们可以使用plotWeights函数。例如，这里我们绘制了Mutation数据中Factor 1的所有权重。与nfeatures我们可以设置应该标记多少个特征(手册让我们手动指定特征，并在图中标记。)

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

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

在给定因素上具有较大(绝对)权重的特征遵循与该因素相关的共变模式。在我们的病例中，这揭示了与IGHV状态的强烈联系，因此在CLL中恢复了一个重要的临床标记(有关生物学解释的详细信息，请参阅我们的论文)。请注意，权重的符号只能相对于其他权重的符号和因子值来解释。

如果您只对顶部功能感兴趣，则可以使用plotTopWeights函数。
例如，我们在这里绘制因子1上负载最大的突变。右边的标志表示装载方向(正/负)。

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

同样，在给定因素中具有较大权重的特征意味着它们遵循与该因素相关的共变模式。例如，这里的因子与肿瘤起源的b细胞有关，与IGHV状态的大权重一致(详情请参阅我们的手稿)。

从前面的图中，我们可以清楚地看到因子1与IGHV状态相关。正如上面的方差分解告诉我们的，这个因素也与所有其他数据模态相关，我们可以研究它在其他模态(如mRNA)上的权重，以建立ighv链接的变化轴与其他分子层的联系。

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

最后，与其查看“抽象”权重，不如观察原始数据中顶部特征的协调异质性。可以使用plotDataHeatmap函数。在该图中，样本(按行排列)根据它们在因子(这里是因子1)上的值排序。在这里，这显示了样本在沿因子的前20个加权基因中的基因表达的清晰模式。

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

3.2.2活动视图中的特征集丰富分析

有时查看单个特征的加载可能具有挑战性，并且通常需要来自功能相关的特征集(即基因本体)的信号组合。

这里我们实现了一个特征集富集分析方法函数(runEnrichmentAnalysis)源自PCGSE包．

该函数的输入是一个MOFA训练模型(MOFAmodel)，执行特征集丰富的因素(一个字符向量)，特征集(一个二进制矩阵)和一组关于应该如何执行分析的选项，请参见文档runEnrichmentAnalysis

函数说明此函数的使用reactome注释，它们包含在包中。根据您的数据，其他基因或特征集可能是有用的，您可以准备您的自定义特征集并使用feature.sets函数的参数。

#加载反应组注释数据(“reactomeGS”)#二进制矩阵，特征集在行中，特征集在列中#执行富集分析gsea <- runEnrichmentAnalysis(MOFAobject, view = "mRNA"，特征。set = reactomeGS, alpha = 0.01)
使用以下选项进行特征集丰富分析…
##视图:mRNA
因素:LF1 LF2 LF3 LF4 LF5 LF6 LF7 LF8 LF9 LF10
##功能集数量:501
##本地统计:加载中
##转换:abs.value
##全局统计:mean.diff
##统计检验:参数

下一步是将基因集富集分析的结果可视化。有几种方法:

  1. 绘制每个因子富集基因集的数量
plotEnrichmentBars (gseaα= 0.01)

从中我们发现因子3-6和8在FDR中富集了1%的基因集。为了研究这些是哪些基因集，我们可以选择一个感兴趣的因子，并将最丰富的基因集可视化如下:

  1. 为每个因素绘制最丰富的路径
interestingFactors <- 4:5 fseaplots <- lapply(interestingFactors, function(factor) {plotEnrichment(MOFAobject, gsea, factor = factor, alpha = 0.01, max. lapply(interestingFactors, function(factor) {plotEnrichment(MOFAobject, gsea, factor = factor, alpha = 0.01, max. lapply)cowplot::plot_grid(fseaplots[[1]]， fseaplots[[2]]， ncol = 1, labels = paste("Factor"， interestingFactors))

这表明因子4正在捕捉与免疫反应相关的变异(可能是由于t细胞污染了样本)，因子5与应激反应的差异有关，正如我们在论文中讨论的那样。

3．3通过因子对样本进行排序，揭示样本空间中的聚类和梯度

样例可以沿着感兴趣的因素进行可视化plotFactorScatter函数。我们可以使用模型中包含的特征(如IGHV或trisomy12)来着色或塑造样品。另外，外部协变量也可以用于此目的。

plotFactorScatter(MOFAobject, factors = 1:2, color_by = "IGHV"， #颜色由IGHV值组成的训练数据shape_by = "trisomy12" #形状由trisomy12值组成的训练数据)

在这里，我们再次发现，根据因子1的IGHV状态(颜色)和由突变视图中相应的因子权重指示的因子2的12三体(形状)的有无，样本之间有明显的分离。

对所有或部分因素的成对散点图进行概述plotFactorScatters函数

plotFactorScatters(MOFAobject, factors = 1:3, color_by = "IGHV")
##已注册S3方法被“GGally”覆盖:##方法从## +。gg ggplot2

单个因素可以使用plotFactorBeeswarm函数

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

3.4自定义分析

对于权重和因子的自定义探索，您可以使用' get '函数直接从模型中获取变量:getWeights，getFactors而且getTrainData：

MOFAweights <- getWeights(MOFAobject, views = "all"， factors = "all"， as.data.frame = TRUE #如果为TRUE，则输出一个长数据帧格式。如果为FALSE，则输出一个宽矩阵格式。
##特征因子值视图## 1 D_001_1 LF1 0.0002667769毒品## 2 D_001_2 LF1 0.0019258319毒品## 3 D_001_3 LF1 0.0356598749毒品## 4 D_001_4 LF1 0.0292638984毒品## 5 D_001_5 LF1 0.0158091594毒品## 6 D_002_1 LF1 -0.0161823261毒品
MOFAfactors <- getFactors(MOFAobject, factors = c(1,2)， as.data.frame = FALSE #如果为TRUE，则输出一个长数据帧格式。如果为FALSE，则输出一个宽矩阵格式。
## lf1 lf2 ## 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)
##查看特征样本值## 1突变gain2p25.3 H045 0 ## 2突变gain3q26 H045 0 ## 3突变del6p21.2 H045 0 ## 4突变del6q21 H045 0 ## 5突变del8p12 H045 0 ## 6突变gain8q24 H045 1

4进一步的功能

4．1视图预测

与预测函数中，基于MOFA模型的全部或子集因子可以预测完整视图。

predictedDrugs <- predict(MOFAobject, view="Drugs"， factors =" all")[[1]] #训练数据(包括缺失值)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)

# predictedDrugs[1:40,1:20]， cluster_rows = FALSE, cluster_cols = FALSE, show_rownames = FALSE, show_colnames = FALSE)

4.2缺失观测值的归因

与嫁祸于函数中所有缺失值均基于MOFA模型进行估算。然后将估算的数据存储在ImputedData槽，并可通过getImputedData函数。

MOFAobject <- impute(MOFAobject) imputedDrugs <- getImputedData(MOFAobject, view="Drugs")[[1]] #训练数据(包括缺失值)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.３基于潜在因素的样本聚类

样本可以根据它们在某些或所有潜在因素上的值进行聚类clusterSamples函数。类可以可视化集群plotFactorScatters函数

set.seed(1234) clusters <- clusterSamples(MOFAobject, k = 2， # k-means函数因子= 1 #用于聚类的因子的集群数量)plotFactorScatter(MOFAobject, factors = 1:2, color_by = clusters)

5SessionInfo

sessionInfo ()
## R版本4.0.3(2020-10-10)##平台:x86_64-pc-linux-gnu(64位)##运行在Ubuntu 18.04.5 LTS ## ##矩阵产品:默认## BLAS: /home/biocbuild/bbs-3.12-bioc/R/lib/libRblas。所以## LAPACK: /home/biocbuild/bbs-3.12-bioc/R/lib/libRlapack。所以## ## 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_phone = c# # [11] LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION=C ## ##附加的基本包:##[1]并行stats4统计图形grDevices utils数据集##[8]方法基础## ##其他附加包:[1] MultiAssayExperiment_1.16.0 SummarizedExperiment_1.20.0 ## [5] Biobase_2.50.0 GenomicRanges_1.42.0 ## [7] GenomeInfoDb_1.26.3 IRanges_2.24.1 ## [9] S4Vectors_0.28.1 BiocGenerics_0.36.0 ## [11] MatrixGenerics_1.2.1 matrixStats_0.58.0 ## [13] BiocStyle_2.18.1 ## ##通过命名空间加载(并且没有附加):## [4] BiocManager_1.30.10 highr_0.8 GenomeInfoDbData_1.2.4 ## [7] vipor_0.4.5 yaml_2.2.1 ggrepel_0.9.1 ## [10] corrplot_0.84 pillar_1.4.7 lattice_0.20-41 ## [13] glue_1.4.2 reticulate_1.18 digest_0.6.27 ## [19] cowplot_1.1.1 htmltools_0.5.1.1 Matrix_1.3-2 ## [22] plyr_1.8.6 pkgconfig_2.0.3 pheatmap_1.0.12 ## [25] magick_2.6.0 bookdown_0.21 zlibbioc_1.36.0 ## [28] purrr_0.3.4[40] beeswarm_0.2.3 tools_4.0.3 lifecycle_0.2.0 ## [43] string_1 .4.0 Rhdf5lib_1.12.1 munsell_0.5.0 ## [46] DelayedArray_0.16.1 compiler_4.0.3 rlang_0.4.10 ## [52] rhdf5filters_1.0 grid_4.0.3 RCurl_1.98-1.2 ## [52] iterators_1.0.13 rhdf5filters_1.2.0 bitops_1.0-6 ## [55] labeling_0.4.2 rmarkdown_2.6 gtable_0.3.0 ## [58]codetools_0.2-18 reshape_0.8.8 DBI_1.1.1 ## [61] reshape2_1.4.4 R6_2.5.0 knitr_1.31 ## [64] dplyr_1.0.4 stringi_1.5.3 ggbeeswarm_0.6.0 ## [67] Rcpp_1.0.6 vctrs_0.3.6 tidyselect_1.1.0 ## [70] xfun_0.20