引言：为什么需要多组学因子分析
随着高通量测序技术的发展，同一批样本往往同时测了多种组学数据：转录组(RNA-seq)、甲基化(DNA methylation)、蛋白质组(proteomics)、代谢组(metabolomics) 等。这些数据从不同层次刻画同一群样本，称之为多组学数据(multi-omics data)。
多组学数据的核心价值在于"信息互补"：转录组反映基因表达，甲基化反映表观调控，蛋白组反映最终功能分子。单独分析某一种组学往往只能看到"一角"，而整合起来才能还原完整的生物学图景。
但整合多组学面临两大难题：
- 维度灾难：每种组学都有成千上万个特征（基因、位点、蛋白），远超样本数。
- 组学间异质性：不同组学的数据分布、量纲、噪声结构完全不同，不能简单拼接后直接跑 PCA。
多组学因子分析(Multi-Omics Factor Analysis, MOFA) 就是为解决这两个问题而生的无监督方法。它把多种组学作为"多个视图(view)"，从中分解出共享的潜因子(latent factors)，既能降维、又能揭示组学间共享与组特异性的生物学信号。
MOFA 最初由 Argelaguet 等人提出（Nature Methods 2018），2020 年升级为 MOFA+，支持多组(multi-group)结构。R 中的实现是 MOFA2 包，官方文档见 https://biofam.github.io/MOFA2/。
MOFA 的核心原理：概率矩阵分解
1. 从单组学 PCA 到多组学 MOFA
先回顾 PCA。对单组学矩阵 $\mathbf{X}$（$N$ 个样本 × $D$ 个特征），PCA 把它近似分解为：
$$ \mathbf{X} \approx \mathbf{Z}\mathbf{W}^\top $$其中 $\mathbf{Z}$ 是 $N\times K$ 的因子得分(factor scores)（也叫潜变量），$\mathbf{W}$ 是 $D\times K$ 的载荷(loadings)，$K$ 是因子数。因子得分是每个样本在新坐标下的位置，载荷是每个特征对因子的贡献。
MOFA 把同样的思想推广到多组学：假设有 $M$ 个视图（组学），第 $m$ 个视图的观测矩阵 $\mathbf{X}_m$（$N$ 样本 × $D_m$ 特征）共享同一个因子得分矩阵 $\mathbf{Z}$，但每个视图有自己的载荷矩阵 $\mathbf{W}_m$：
$$ \mathbf{X}_m \approx \mathbf{Z}\mathbf{W}_m^\top + \boldsymbol{\varepsilon}_m $$关键点在于 $\mathbf{Z}$ 是共享的——所有组学都投影到同一组潜因子上。于是 MOFA 找到的因子，代表跨组学的共同变异模式。如果某个因子只在部分组学中有信号，则对应的载荷 $\mathbf{W}_m$ 在该组学中趋近于 0。
2. 概率模型：MOFA 的完整生成模型
MOFA 是概率模型（概率矩阵分解的贝叶斯版本），把观测看作是潜因子通过载荷映射，再叠加上噪声：
$$ \mathbf{X}_m = \mathbf{Z}\mathbf{W}_m^\top + \boldsymbol{\varepsilon}_m, \qquad \boldsymbol{\varepsilon}_m \sim N(0, \boldsymbol{\Sigma}_m) $$其中：
- 因子得分 $\mathbf{Z}$ 服从先验 $\mathbf{Z} \sim N(0, \mathbf{I})$；
- 载荷 $\mathbf{W}_m$ 上的先验是自动相关性判定(Automatic Relevance Determination, ARD)——每个因子在每个视图上有一个精度参数 $\alpha_{km}$。如果某因子对某视图不重要，$\alpha$ 会变大，把载荷收缩到 0，从而实现视图稀疏性（即因子只在相关组学中有载荷）。
3. ARD 先验与因子稀疏性
ARD 是 MOFA 的核心机制。每个载荷矩阵元素有先验：
$$ w_{dm} \sim N(0, \alpha_{km}^{-1}), \qquad \alpha_{km} \sim \text{Gamma}(a_0, b_0) $$$\alpha_{km}$ 是该因子-视图组合的精度参数，通过变分推断学习。若某因子在某视图不重要，$\alpha_{km}$ 会变得很大，强迫 $w_{dm}\to 0$。这正是 MOFA 能自动实现"因子只在部分组学中活跃“的原因。
4. MOFA+ 的多组扩展
MOFA+（MOFA2 支持）进一步允许**多组(multi-group)**结构：样本可以来自不同批次、不同队列或不同条件。此时因子得分变为 $\mathbf{Z}_{gk}$，即每个 group 有自己的因子值，但共享载荷结构：
$$ \mathbf{X}_{gm} = \mathbf{Z}_{gk}\mathbf{W}_{km}^\top + \boldsymbol{\varepsilon}_{gm} $$这意味着 MOFA 能区分共享因子(在所有组中都活跃)与组特异因子(只在某组中活跃)。这正是教学演示图里看到的：Factor2 只在 Disease 组有信号。
5. 似然与推断
给定各组学观测，模型的完全似然为：
$$ \log p(\mathbf{X} \mid \mathbf{Z}, \mathbf{W}) = \sum_m \log p(\mathbf{X}_m \mid \mathbf{Z}, \mathbf{W}_m) $$由于是贝叶斯模型，MOFA 通过变分推断(Variational Inference) 最大化证据下界(ELBO, Evidence Lower Bound) 来近似后验。ELBO 是训练收敛的判据，也用于模型比较（因子数选择）。
MOFA2 的后端实现是 Python 包 mofapy2，R 的 run_mofa() 通过 reticulate 调用它进行优化。
MOFA 与其它整合方法的区别
| 方法 | 核心思想 | 适用场景 | 特点 |
|---|---|---|---|
| MOFA | 概率矩阵分解 + ARD 稀疏 | 多种组学、可分组 | 因子可跨组学解释，支持多组，可做下游关联分析 |
| iCluster | 联合潜变量聚类 | 基因+拷贝数等 | 侧重于聚类发现亚型 |
| NMF | 非负矩阵分解 | 单一组学表达矩阵 | 非负约束，解释为共表达模块 |
| sMBPLS | 多块偏最小二乘 | 有监督关联 | 需要先验的响应变量 |
| JIVE | 联合+个体变异分解 | 两组学对齐样本 | 区分共享与特异成分，但不能分组 |
MOFA 的相对优势在于：无监督（无需标签）、组学间权重可解释（ARD 稀疏）、支持多组结构（MOFA+）、输出可直接用于下游统计（因子值与协变量关联）。
环境准备：安装 MOFA2 与 Python 后端
MOFA2 是 Bioconductor 包，安装时需要同时安装它的 Python 后端 mofapy2。
|
|
|
|
重要：
run_mofa()需要调用 Python 的mofapy2，因此必须让 R 的reticulate指向装有该模块的 Python 环境。有两种方式：
- 用
reticulate::use_condaenv("mofa", required = TRUE)指定 conda 环境；- 或在 R 启动前设置环境变量
Sys.setenv(RETICULATE_PYTHON = "/路径/python")。若已装 basilisk 且想自动管理环境，可调用
run_mofa(..., use_basilisk = TRUE)，但更推荐显式指定 Python 环境以便复现。
MOFA 完整操作流程
下面用 MOFA2 内置的 make_example_data() 合成一份三组学演示数据（RNA / 甲基化 / 蛋白），样本分成两组（Healthy / Disease）。合成数据的好处是可复现、无需下载，且真实结构已知，便于理解因子含义。真实数据的替换方式见文末。
步骤0：加载包并设定随机种子
|
|
步骤1：准备多组学数据
MOFA2 要求数据是一个 list，每个元素是一个组学视图的矩阵，行是特征(features)，列是样本(samples)。所有视图必须共享同一批样本且样本顺序一致。
我们用 make_example_data() 生成 3 个视图、2 组各 60 个样本、真实含 4 个潜在因子的数据：
|
|
make_example_data 的内部逻辑是：先生成真实的因子得分 $\mathbf{Z}$ 和载荷 $\mathbf{W}$，再用它们乘出无噪声信号并加观测噪声。这样我们知道因子"真相”，方便验证 MOFA 能否找回它们。
样本元数据：create_mofa() 需要样本分组信息。我们构造一个包含 group 和协变量 disease 的元数据框：
|
|
要点：元数据的行顺序必须与数据矩阵的列顺序完全一致。
make_example_data已经保证了这点，但用真实数据时务必自行核对。
步骤2：创建 MOFA 对象
|
|
这一步只是把数据组装成 MOFA 对象，尚未训练。
步骤3：配置模型参数
MOFA 有三组关键选项：data_options、model_options、training_options，每一组都有默认值。我们先取默认，再按需修改。
|
|
常用参数解读
数据选项 data_options：
center_groups：是否按组居中。不同组样本总表达水平常不同（批次效应），居中可避免它混入因子。scale_views：是否按视图缩放，使不同量纲的组学有可比性。
模型选项 model_options：
num_factors：初始因子数。建议设得比预期偏大（如 15），MOFA 训练后会自动移除不解释方差的因子。likelihoods：各视图的数据分布，默认"gaussian"。计数型数据（如 RNA 计数）可用"poisson"，二值可用"bernoulli"。spikeslab_weights：是否用 spike-and-slab 先验做特征层面的稀疏，适合特征数很多的场景。
训练选项 training_options：
seed：随机种子，保证可复现（MOFA 训练是随机初始化的）。convergence_mode：收敛标准，"slow"更严格但更慢。maxiter：最大迭代次数。
步骤4：准备模型
prepare_mofa() 把配置好的选项写入对象，并做最后校验。这一步之后对象会进入"可训练"状态。
|
|
步骤5：训练模型
|
|
训练过程会打印每次迭代的 ELBO 值，最终显示 Converged! 并给出模型概况。控制台会提示去除了多少个"不解释方差的因子"——这正是 ARD 稀疏机制在起作用。
结果解读
1. 模型概况
|
|
输出的要点：
- Number of views: 3（RNA / DNAme / Protein）
- Number of groups: 2（Healthy / Disease）
- Number of factors: 训练后保留的因子数（初始设了 8 个，不解释方差的会被剔除）
2. 因子得分与载荷的结构
MOFA2 的 get_factors() 返回每个样本在每个因子上的得分（$N\times K$ 矩阵），get_weights() 返回每个特征在每个因子上的载荷。由于是多组结构，get_factors() 返回一个按组分拆的列表。
|
|
3. 方差解释（最核心的输出）
MOFA 最重要的问题之一是"每个因子解释了多少方差"。用 get_variance_explained() 计算：
|
|
实际输出里 r2_total 显示每个组学在各组被解释的总方差（比如 RNA 在 Disease 组约 77%、Healthy 组约 56%），r2_per_factor 则给出每个因子对每个视图的方差贡献。你会发现 Factor1 主导了 Healthy 组所有视图的方差，而 Factor2 主导了 Disease 组的 RNA 和 Protein——这正是下方热图要讲的核心。需要提醒的是：不同因子解释的方差量级差异很大，个别因子（如这里的 Factor2/3/4）贡献较小，不代表它们无生物学意义，只是方差占比低。
用热图看更直观：
|
|
这张图正是 MOFA 的精髓所在（颜色越深代表该因子解释的方差比例越高）：
- Factor1（最下面一行）在 Healthy 组对 RNA、DNAme、Protein 都有较强信号（深紫），而在 Disease 组几乎不解释方差——这是一个在 Healthy 组占主导的因子；
- Factor2 只在 Disease 组对 RNA 和 Protein 有强信号（深蓝），在 Healthy 组几乎为 0——这是组特异因子（只在疾病组活跃），可能捕捉疾病相关的跨组学变异；
- Factor3 / Factor4 主要解释 DNAme 的方差（两个组的 DNAme 列都有信号），是视图特异因子，代表甲基化特有的变异维度。
注意：每个因子的"活性"是组特异且视图特异的——一个因子可能只在一组、且只在部分组学里活跃。能同时区分"跨组学共享"“组特异"“视图特异"信号，正是 MOFA+ 相对普通 PCA 的突出能力。这也提醒我们：解读因子时必须结合组与视图两个维度。
4. 解释因子：载荷的可视化
要回答"某个因子在生物学上代表什么”，需要看它在哪些特征上有大载荷。用 plot_weights() 展示因子在某个视图中的 Top 特征：
|
|
条形图按载荷大小排列，正负载荷分别表示该特征与因子正/负相关。
5. 因子与样本分组 / 协变量的关系
MOFA 找到的因子常与生物学状态相关。用 plot_factor() 按协变量着色，查看因子值是否区分组别：
|
|
可以看到 Factor1 在 Disease 组聚集在 0 附近，而 Healthy 组散布范围大——说明 Factor1 捕捉到两组间的结构差异。这是后续关联分析的基础。
6. 因子与协变量的定量关联
更严谨的做法是计算因子得分与协变量的回归。由于多组结构下 get_factors() 按组返回、列名带组前缀，这里演示合并两组后，用线性回归检验因子与疾病状态的关联（教学示例）：
|
|
若 p 值显著，说明该因子与疾病状态相关，可作为下游解读的线索。
进阶：模型诊断与选择
1. 选择合适的因子数
MOFA 靠 ARD 机制自动决定"有效因子”——初始给足因子数（如 15），训练后不解释方差的因子会自动被剔除。因此初始 num_factors 设大一些是安全的。也可比较不同因子数下的 ELBO：
|
|
实际项目中的经验做法：设置 num_factors = 15，用 convergence_mode = "slow" 跑一次，观察保留的因子数；若太多或太少再调整。
2. 检查是否有"技术性因子"
训练后 MOFA2 会给出 QC 警告，例如"某因子与总表达特征数强相关"。出现这类警告时，该因子可能捕捉的是文库大小 / 批次效应而非生物学信号，应结合载荷与协变量关联判断是否剔除。
3. 缺失值与其它似然
MOFA 天然支持缺失值（概率模型可从观测数据推断），这是它相比普通 PCA 的又一优势。若组学是计数数据（如 RNA 计数），把对应的 likelihoods 设为 "poisson"；二值数据设为 "bernoulli"。
MOFA 的优缺点
优点
- 真正的多组学整合：共享潜因子把不同组学放在同一坐标系，自动识别跨组学信号。
- 无监督、可解释：不需要标签即可降维，因子可通过载荷和协变量关联解释。
- 视图与组双重稀疏：ARD 先验自动决定"哪些因子在哪些视图/组中活跃"，避免人为设定。
- 支持多组结构：MOFA+ 能区分共享因子与组特异因子，适合多队列、多批次分析。
- 容忍缺失值：概率模型可处理部分组学缺失的样本。
- 下游分析友好：因子得分可直接用于 MANOVA、回归、聚类与可视化。
局限性
- 依赖 Python 后端：需要额外安装
mofapy2，环境配置是主要门槛。 - 计算开销：样本和特征很多时，变分推断耗时较长（可用
stochastic随机变分近似加速）。 - 对预处理敏感：中心化、缩放、归一化的选择会影响因子；原始 count 数据需先正确转换。
- 因子解释需要背景知识：MOFA 输出的是统计因子，其生物学语义需结合载荷与协变量关联人工解读。
- 高斯假设为主：默认
gaussian似然；非正态组学需显式指定其它似然。
注意事项
- 样本对齐：所有视图必须共享相同样本且列顺序一致，这是最常踩的坑。
- 预处理是前提：先做质量控制、批次校正（如
removeBatchEffect）与归一化，再喂给 MOFA，否则因子会被技术噪音主导。 - 设置随机种子：
train_opts$seed保证可复现，MOFA 训练是随机初始化的。 - 区分生物与技术因子：结合协变量关联和 QC 警告，识别并剔除技术性因子。
- 特征命名：给特征加视图前缀（如
RNA_、DNAme_），否则不同视图的特征同名会混淆。
从示例数据迁移到真实数据
把 data_list 换成你自己的多组学矩阵即可，格式要求是：
|
|
后续 create_mofa → prepare_mofa → run_mofa 完全一样。若数据是稀疏矩阵或过大，可用 create_mofa(..., extract_metadata = FALSE) 等方式控制内存。
参考文献与扩展
- Argelaguet, R., et al. (2018). Multi-Omics Factor Analysis—a framework for unsupervised integration of multi-omics data sets. Molecular Systems Biology, 14(6), e8124.
- Argelaguet, R., et al. (2020). MOFA+: a statistical framework for comprehensive integration of multi-modal single-cell data. Genome Biology, 21, 111.
- MOFA2 官方文档与教程：https://biofam.github.io/MOFA2/
- 相关扩展：
MOFA2的run_mofa底层是mofapy2；MEFISTO（时空因子分析）可通过mefisto_options配置。 - 集成到单细胞分析：MOFA+ 可与
Seurat/Scanpy的降维结果关联，用于多模态单细胞整合。
小结
多组学因子分析(MOFA)用概率矩阵分解 + ARD 稀疏先验，把多种组学整合进同一组共享潜因子，既能降维，又能揭示跨组学共享信号与组/视图特异性信号。本文用 MOFA2 包 + 合成三组学数据走通了完整流程：准备数据 → 创建对象 → 配置参数 → 训练 → 解读方差解释与载荷 → 因子与协变量关联 → 模型诊断。掌握这套流程，你就能在自己的多组学数据上挖掘跨组学的结构信号。