MOFA与MOFA+多组学因子分析原理

可视化演示

展示多组学数据整合流程、贝叶斯因子分解模型、ARD稀疏权重约束热力图及MOFA与MOFA+版本功能对比，含星座网络背景和发光动画。

多组学背景

想象乐团里的小提琴、管乐与打击乐同时发声。如果只录下小提琴轨，你会漏掉管乐与打击乐之间的呼应；只有同时倾听全团，才能理解整首作品。多组学数据正是如此——同一批肿瘤病人的组织样本，既测了转录组（RNA 表达）、又测了甲基化谱与蛋白组。每种数据都像乐团的一个声部，虽然量纲、噪声和特征维度都不同，但它们毕竟来自同一批组织，必然受同一套生物学程序驱动。

于是，声部之间出现了“协同变异”：癌基因激活时，往往同时在 RNA 层面看到表达上调、在蛋白层面看到蛋白丰度增加、在甲基化层面看到调控区域甲基化异常。这种跨组学联动就是分析的核心对象。传统做法是把每个组学单独建模，相当于只听单一声部，会错过重要的协变信号；而把所有组学直接拼接成一个超长矩阵，又会让量纲差异最大的那个组学主导模型。MOFA 的处理方式是：让所有组学共享一个低维的“因子得分矩阵”，再给每个组学配备独立的“载荷矩阵”，从而既保留共性又兼顾个性。

$$X_{g} = W_{g} Z^{\top} + E_{g}, \quad g = 1, 2, \dots, G$$

1. X_g 数据矩阵：第 g 个组学的观测数据，n 行代表样本、p_g 列代表特征（如基因或甲基化位点），是乐评人听到的“声部原声”。

$$X_{g} \in \mathbb{R}^{n \times p_{g}}$$

2. Z 因子得分矩阵：n × K 的潜变量矩阵，所有组学共享，代表每个样本在 K 条“旋律线”上的坐标，是隐藏在主旋律中的样本差异。

$$Z \in \mathbb{R}^{n \times K}$$

3. W_g 载荷矩阵：K × p_g 的载荷矩阵，决定第 g 个组学如何由 K 条旋律线性叠加而成，相当于“每个声部对每条旋律的演奏音量”。

$$W_{g} \in \mathbb{R}^{K \times p_{g}}$$

4. E_g 噪声残差：捕获各组学独立的测量噪声与生物学特异波动，是乐谱之外的即兴杂音。

$$E_{g} \in \mathbb{R}^{n \times p_{g}}$$

回到乐团：MOFA 不是把三个声部的录音简单叠在一起，而是先承认“大家共同演绎了主旋律”，再允许每个声部有自己的独奏片段。这样，后续的因子解释与分子筛选才有了清晰的对象。

MOFA模型理论

乐评人该怎么拆解一场演出？他不会逐秒记录每个音符，而是抓住几条反复出现的旋律线索。MOFA 的建模思路如出一辙：把庞大的多组学矩阵压缩成 K 个因子，每个样本在这 K 个因子上都有自己的“坐标”，称为因子得分。某个样本在因子得分上接近，就说明它在对应的生物学程序上状态相似。

因子得分矩阵 Z 是全体组学共享的，因为样本是同一样的；但每个组学对因子的响应方式不同，于是 MOFA 为第 g 个组学设置独立的载荷矩阵 W_g。载荷值衡量“因子 k 对特征 j 的影响强度”。例如调控代谢的共享因子，在转录组的某些代谢酶基因上载荷很高，在甲基化组的启动子区域载荷可能相对温和，这正体现了不同分子层对同一程序的响应差异。

$$x_{ij} = \sum_{k=1}^{K} z_{ik} w_{kj} + \varepsilon_{ij}$$

1. x_ij 观测值：第 i 个样本第 j 个特征的表达量，例如病人 i 的基因 j 的 RNA 丰度。

$$x_{ij} \in \mathbb{R}$$

2. z_ik 因子得分：样本 i 在因子 k 上的投影，描述“这个样本离主旋律 k 有多近”。

$$z_{ik} \in \mathbb{R}$$

3. w_kj 载荷系数：因子 k 对特征 j 的贡献强度，相当于旋律线 k 在乐器 j 上的音量大小。

$$w_{kj} \in \mathbb{R}$$

4. 残差 ε_ij：模型未能解释的部分，来自噪声或独立于因子的生物波动。

$$\varepsilon_{ij} \sim \mathcal{N}(0, \sigma_{g}^{2})$$

理解拆解之后，你会明白 MOFA 的优势：因子是“共享”的，所以不同组学的信息被拉到同一坐标系下比较；载荷是“独立”的，所以每个分子层保留了自身的响应特色。这就是“共享潜变量 + 独立载荷”的巧妙之处。

模型参数设定

乐评人下笔前要确定写几条旋律线索：太多，点评变得支离破碎；太少，又遗漏重要主题。模型设定阶段要回答同样的问题——先验怎么选、方差怎么分、因子数 K 取几。MOFA 采用概率框架，假设观测值由因子信号叠加高斯噪声生成。对载荷 W_g 施加零均值高斯先验，相当于在写乐评前“默认每个乐器与旋律关系不大”，除非数据给出强证据，否则权重不会轻易离开零。

方差分解是理解因子解释力的关键。整体变异同时包含“因子信号”和“噪声”两部分。若某个因子能解释某组学的大部分方差，说明它是该组学的重要驱动者；若所有因子加起来只能解释很小比例，说明这个组学与共享程序关系微弱。因子数 K 通常通过“方差解释率曲线”确定：绘制 K 从 1 到最大值时总解释方差的变化曲线，寻找拐点——曲线由陡变缓的位置，就是主旋律与杂音的边界。

$$\text{Var}(X_{g}) = \sum_{k=1}^{K} \text{Var}(w_{g k} z_{k}) + \text{Var}(E_{g})$$

1. 方差分解含义：左侧是第 g 个组学的整体变异，右侧第一项是所有因子贡献之和，第二项是噪声残留。

$$\text{Var}(X_{g}) = \underbrace{\text{Var}(W_{g} Z^{\top})}_{\text{因子信号}} + \underbrace{\text{Var}(E_{g})}_{\text{噪声}}$$

2. 单个因子项：因子 k 的方差贡献由载荷 w_gk 的平方与因子得分 z_k 的方差共同决定，两者越大，解释力越强。

$$\text{Var}(w_{g k} z_{k}) = w_{g k}^{2} \cdot \text{Var}(z_{k})$$

3. 解释率 VAF：因子解释方差占总方差的比例，是评估因子重要性的标准尺。

$$\text{VAF}_{g} = \frac{\sum_{k} w_{g k}^{2} \text{Var}(z_{k})}{\text{Var}(X_{g})}$$

好的模型设定就像一份详尽的乐评底稿：既不自作聪明地预设太多旋律，也不遗漏真正的主调。因子数 K 不是越大越好，而是让每条旋律都有明确意义、能稳定复现。

因子解释方法

终于到了乐评的核心段落：哪些旋律是全团的共同主题？哪些是某声部的独奏华彩？在 MOFA 中，这两种角色分别由共享因子与特异因子扮演。共享因子在多个组学上都有可观的方差贡献，例如“细胞增殖程序”在转录组与蛋白组上同步激活；特异因子则只在一个组学中显现，例如某个组织特有的甲基化修饰模式，在 RNA 层面几乎看不到对应信号。

识别方法很直观：对每个因子，分别计算它在各组学的平均载荷强度。如果多个组学的平均载荷都明显高于本底，则判为共享因子；如果只有一个组学强烈响应，其余组学接近零，则是特异因子。随后，取出每个因子载荷绝对值最高的那些特征做功能富集——比如对转录组因子取高权重基因做 GO/KEGG 通路分析，就能给这条旋律贴上“炎症应答”或“氧化磷酸化”等生物学标签。

$$\bar{w}_{g} = \frac{1}{p_{g}} \sum_{j=1}^{p_{g}} |w_{g j}|$$

1. 绝对值载荷：取绝对值是为了避免正负载荷相互抵消，只看实质贡献。

$$|w_{g j}| = \sqrt{w_{g j}^{2}}$$

2. 组内平均：把第 g 个组学所有特征上载荷的绝对值求平均，得到因子在该组学的整体活跃度。

$$\bar{w}_{g} = \frac{1}{p_{g}} \sum_{j=1}^{p_{g}} |w_{g j}|$$

3. 跨组学比较：比较 \bar{w}_1, ……, \bar{w}_G，若多个值都高则是共享因子，若只有一个高则是特异因子。

$$\bar{w}_{1} \approx \bar{w}_{2} \gg \bar{w}_{3} \quad \Rightarrow \quad \text{共享}$$

给因子命名之后，抽象的数字矩阵就变成了有生物学意义的叙事：因子 1 是“跨组学的增殖主旋律”，因子 3 是“甲基化独有的沉默华彩”。这样的解读，正是多组学整合相对单组学分析最迷人的地方。

组学权重解析

知道旋律还不够，乐评人还要点评每个声部的贡献度。组学载荷 W_g 就是这份“贡献清单”。特征 j 在因子 k 上的权重 w_kj 绝对值越大，说明该特征越直接地参与这条因子驱动的生物程序。因此，载荷权重是筛选关键分子的第一线索：对转录组数据，高权重基因常常就是该程序的核心效应分子；对甲基化数据，高权重位点往往定位在关键的调控区域。

跨组学比较权重还有额外收获。同一因子在不同组学上的载荷分布，能揭示“调控层级”的线索：如果 DNA 甲基化位点载荷高、但对应基因的转录组载荷低，说明该程序主要在表观层面调控；如果转录组和蛋白组载荷同步升高，则提示从 mRNA 到蛋白质的传递效率很高。这种跨层级的权重对比，是单组学分析永远看不到的风景。

$$\text{VAF}_{g k} = \frac{\sum_{j} w_{g k j}^{2} \sum_{i} z_{i k}^{2}}{\sum_{j} \sum_{i} x_{g i j}^{2}}$$

1. 载荷能量：因子 k 在组学 g 上的所有载荷平方和，反映该因子在特征空间的总作用强度。

$$\sum_{j} w_{g k j}^{2}$$

2. 得分能量：因子 k 在所有样本上的得分平方和，反映该因子在样本间的区分度。

$$\sum_{i} z_{i k}^{2}$$

3. 总体变异：组学 g 的全体数据平方和，作为归一化基准。

$$\sum_{j} \sum_{i} x_{g i j}^{2}$$

4. VAF_gk：两者相乘除以总变异，得到因子 k 对组学 g 的方差解释率，是分子筛选的“重要性排行榜”。

权重的意义不在数字本身，而在于它把“有没有关系”量化成了“关系有多强”。生物学家根据这份榜单，既能在海量基因里圈定少数关键分子用于后续验证，也能从跨组学载荷差异中读出深层的调控机制。

稀疏约束机制

一场演出如果所有乐器始终齐奏，听众反而什么都记不住。优秀的编曲会让每条旋律只由少数关键乐器承担，形成干净、鲜明的音色。稀疏约束正是这一思想的数学实现：它要求每个因子只由少数特征驱动，其余特征的载荷被强制压向零，从而滤除不相关的噪声特征，让因子更聚焦、更可解释。

MOFA 实现稀疏性的关键机制是 ARD（自动相关性确定）先验。每个权重 w_gkj 被赋予零均值高斯先验，但精度参数 α_gkj 由 Gamma 超先验控制。当数据不支持某个特征与因子存在关联时，α 会自适应地变得很大，反向方差 α^{-1} 趋近于零，将权重“拉”向零；反之，真正重要的特征其 α 保持在较小水平，权重得以保留。整个过程由变分推断自动完成，用户只需设定超参数，模型就在数据中自行学会“哪些乐器该闭嘴”。

$$w_{g k j} \sim \mathcal{N}(0, \alpha_{g k j}^{-1}), \quad \alpha_{g k j} \sim \text{Gamma}(a_{0}, b_{0})$$

1. 零均值高斯先验：权重默认围绕 0 分布，绝对值越大的权重需要越强的数据证据。

$$p(w_{g k j} \mid \alpha_{g k j}) = \mathcal{N}(0, \alpha_{g k j}^{-1})$$

2. 精度参数 α：α 越大，对应高斯分布越窄，权重被越强地压缩到 0 附近。

$$\alpha_{g k j} \to \infty \;\Rightarrow\; w_{g k j} \to 0$$

3. Gamma 超先验：α 本身也服从 Gamma 分布，模型在训练中从数据更新 α，而非人工固定。

$$p(\alpha_{g k j}) = \text{Gamma}(a_{0}, b_{0})$$

4. 稀疏结果：大量 α 收敛到较大值，对应权重为零；只有少数特征保留非零载荷，构成因子的“精简编制”。

稀疏约束如同指挥家果断地让部分乐器在特定段落休止。这样处理之后，每个因子的高权重基因数量大大减少，实验验证的目标更明确，因子之间的区分度也显著提升，这就是“少即是多”的建模智慧。

MOFA+扩展

原版 MOFA 处理的是“一个时间点、一群样本”的快照。但真实研究往往是纵向追踪——同一批病人在不同治疗阶段多次采样，或来自不同组织部位。这就像乐团开启多城巡演，每个城市（时间点/分组）的演出有全场共通的曲目风格，也会因地制宜地调整华彩段。MOFA+ 正是为这种复杂设计而生。

MOFA+ 的核心扩展是：保留共享的载荷矩阵 W_g，让每个时间点或分组拥有自己的因子得分 Z^{(t)}。由于载荷跨条件共享，因子的生物学含义在时间点之间保持一致，我们就能追踪同一因子如何随着治疗进程而激活或抑制。例如，某免疫浸润因子在基线、治疗中、缓解后三个时间点的得分逐渐下降，就提示免疫信号的动态变化与疗效相关。这种“载荷不变、得分变化”的结构，让纵向多组学分析真正可行。

$$X_{g}^{(t)} = W_{g} Z^{(t)\top} + E_{g}^{(t)}$$

1. 上标 (t)：代表第 t 个时间点或分组，例如治疗后第 3 个月的采样。

$$X_{g}^{(t)} \in \mathbb{R}^{n_{t} \times p_{g}}$$

2. 共享载荷 W_g：载荷跨所有时间点保持相同，保证因子在时间轴上“是同一个概念”。

$$W_{g} \in \mathbb{R}^{K \times p_{g}}$$

3. 时序得分 Z^{(t)}：每个时间点有自己的因子得分矩阵，展示样本在该时刻的因子状态。

$$Z^{(t)} \in \mathbb{R}^{n_{t} \times K}$$

4. 动态追踪：比较不同 t 的 z 值，即可描画因子的时间轨迹。

$$\Delta z_{k} = z_{k}^{(t_{2})} - z_{k}^{(t_{1})}$$

如果说 MOFA 是为单场音乐会写乐评，MOFA+ 就是为整轮巡演写长篇评论：既能识别每站共同的曲目风格，又能捕捉各站的即兴变化。这种灵活的建模方式，让多组学分析从静态描述走向动态机制研究。

设计权衡取舍

最后，任何乐评人都要面对现实限制：篇幅有限、精力有限，必须决定重点分析哪几条旋律、评价到什么深度。MOFA 的使用同样充满权衡。首先是因子数 K：K 太小会遗漏重要信号，K 太大则因子碎片化、难以解读、甚至过拟合噪声。常用做法是观察方差解释率曲线，在拐点附近取 K，并辅以交叉验证保证因子可重复。

其次是稀疏强度。稀疏约束过强，弱但真实的信号会被当作噪声抹掉；过弱则因子被大量零碎特征包裹，难以聚焦。实践中通常用若干候选超参数跑同类数据，比较因子的稳定性和富集显著性，再选择恰到好处的强度。最后是共享与特异的侧重：如果研究目标是跨组学共性机制，就把注意力放在共享因子上；如果目标是发现某个分子层的特异调控，则特异因子是主角。两者没有绝对优劣，取决于科学问题。

$$\text{VAF}(K) = 1 - \frac{\sum_{g} \left\| X_{g} - W_{g}^{(K)} Z^{(K)\top} \right\|_{F}^{2}}{\sum_{g} \left\| X_{g} \right\|_{F}^{2}}$$

1. 残差平方和：用 K 个因子重建数据后的误差，衡量未解释部分。

$$\left\| X_{g} - W_{g}^{(K)} Z^{(K)\top} \right\|_{F}^{2}$$

2. 总变异：所有组学的原始数据能量之和，作为对比基准。

$$\sum_{g} \left\| X_{g} \right\|_{F}^{2}$$

3. VAF(K)：K 个因子贡献的解释方差比例，随 K 增大而单调递增。

$$\text{VAF}(K) \in [0, 1]$$

4. 拐点法则：当曲线从陡增变为平缓时取 K，即“再增加因子也只是锦上添花”。

$$\frac{\partial \text{VAF}}{\partial K} \to 0 \;\Rightarrow\; K^{\ast}$$

设计权衡的本质，是在解释力、简洁性与可重复性之间寻找平衡点。就像乐评人不会把每个装饰音都写进文章，好的 MOFA 分析也有节制——选择适当的 K、适度的稀疏与清晰的因子侧重，才能让结论既深刻又可信。

AI 应用场景

Python 演示

import numpy as np def soft_threshold(X, lam):    return np.sign(X) * np.maximum(np.abs(X) - lam, 0.0) def generate_multi_omics(n_samples=120, n_features=(300, 240, 180), K=3):    rng = np.random.default_rng(42)    Z_true = rng.normal(0, 1, size=(n_samples, K))    X_list = []    for g, nf in enumerate(n_features):        W_g = rng.normal(0, 1, size=(K, nf))        mask = rng.random((K, nf)) < 0.4        W_g[~mask] = 0.0        X_g = Z_true @ W_g + 0.3 * rng.normal(0, 1, size=(n_samples, nf))        X_list.append(X_g)    return X_list def sparse_factor_analysis(X_list, K=3, lam=0.2, max_iter=100):    rng = np.random.default_rng(7)    n_samples = X_list[0].shape[0]    Z = rng.normal(0, 1, size=(n_samples, K))    W_list = []    for X_g in X_list:        W_g = np.zeros((K, X_g.shape[1]))        W_list.append(W_g)    for it in range(max_iter):        for g, X_g in enumerate(X_list):            A = Z.T @ Z + 1e-6 * np.eye(K)            B = Z.T @ X_g            W_list[g] = soft_threshold(np.linalg.solve(A, B), lam)        S = np.zeros((K, K))        T = np.zeros((n_samples, K))        for g, X_g in enumerate(X_list):            S += W_list[g] @ W_list[g].T            T += X_g @ W_list[g].T        Z = T @ np.linalg.inv(S + 1e-6 * np.eye(K))    return Z, W_list def main():    X_list = generate_multi_omics()    K = 3    Z, W_list = sparse_factor_analysis(X_list, K=K)    print("因子得分矩阵形状:", Z.shape)    for g, W_g in enumerate(W_list):        non_zero = int(np.sum(np.abs(W_g) > 1e-6))        print("组学{} 权重矩阵: {}, 非零权重: {}".format(g + 1, W_g.shape, non_zero))    var_explained = []    for g, X_g in enumerate(X_list):        total_var = float(np.sum(X_g ** 2))        res = float(np.sum((X_g - Z @ W_list[g]) ** 2))        var_explained.append(round(1 - res / total_var, 3))    print("各组件方差解释率:", var_explained)    for k in range(K):        top_weights = [np.argsort(np.abs(W_list[g][k]))[-3:].tolist() for g in range(3)]        print("因子{}各组件高权重特征列:".format(k + 1), top_weights) if __name__ == "__main__":    main()

1. generate_multi_omics()模拟三个组学数据，每个组学由共享因子得分 Z_true 与稀疏权重矩阵相乘生成，内置 40% 的零权重模拟真实稀疏性。
2. soft_threshold(X, lam)软阈值算子，将绝对值小于 λ 的元素置零，是稀疏约束在优化中的直接体现，对应 ARD 的最终效果。
3. A = Z.T @ Z + 1e-6 * np.eye(K)给加性正则项，避免矩阵奇异无法求逆，稳定性技巧。
4. W_list[g] = soft_threshold(np.linalg.solve(A, B), lam)固定 Z 后解析更新各组学权重，再施加软阈值实现特征选择。
5. T += X_g @ W_list[g].T跨组学累积信息，让 Z 的更新同时接收所有组学的“投票”，保证是共享因子得分。
6. Z = T @ np.linalg.inv(S + 1e-6 * np.eye(K))固定所有 W 后重新估计因子得分，完成一轮交替迭代。
7. top_weights = ...对每个因子取各组学中权重绝对值最大的特征索引，模拟“高权重基因筛选”的下游分析。

常见误区

面试题