多组学数据堆在手里转录组、蛋白组、甲基化各做各的差异分析最后写讨论的时候发现三张表的结论对不上——这大概是每个做多组学的人都会经历的阶段。我最早接触多组学因子分析这套思路就是被这种“拼盘式分析”逼出来的每一层组学单独看都挺热闹合在一起却讲不出一个完整的故事。后来用上MOFAMulti-Omics Factor Analysis这套工具才慢慢摸到门道——它做的事情说穿了不复杂把多个组学层里共同变化的信号压缩成若干个潜在因子每个因子背后对应一组跨组学的特征组合再回头去解释这些因子和表型的关系。这篇文章我打算按“为什么这么设计、怎么准备数据、怎么跑、怎么解读、哪里会翻车”的顺序来写面向的是已经做过单组学差异分析、想往多组学整合走一步的读者。如果你手里正好有两层以上的组学数据样本能对上又不想一上来就陷进复杂的网络建模那这套方法的性价比相当高。下面所有代码和参数都是我在 R 环境里跑通的写法Python 侧的接口逻辑一致只是包名和调用方式不同我会在对应位置提一句。1. 先把问题想清楚为什么多组学整合总是卡在降维这一步1.1 拼盘式分析和共享因子分析差在哪我刚入门的时候做多组学的标准流程是转录组跑一遍差异表达蛋白组跑一遍差异丰度甲基化再跑一遍差异甲基化然后取交集画韦恩图。这个方法本身没错但它有个隐含假设——每一层组学都是独立的信息源交集就是“多重验证”。实际操作下来你会发现交集往往小得可怜而且交集里的基因未必是最重要的基因。原因是每层组学的噪声结构、动态范围、缺失机制都不一样用同一把尺子比如 P 0.05去卡等于强行假设它们可比。共享因子分析换了个思路不去追每个特征单独的显著性而是先假设存在若干个“潜在驱动因素”这些因素同时影响了多层组学里的一组特征。比如某个转录因子活性升高会导致它调控的 mRNA 上调、对应蛋白上调、启动子区甲基化下降——这三件事在数据里表现为同一方向上的一组协同变化。MOFA 干的就是把这个协同模式找出来用一个因子表示然后告诉你这个因子在每层组学里由哪些特征支撑。这个视角的转换带来两个直接好处。第一降噪单层组学里那种随机波动很难在多个层里同时出现同方向变化所以被自然过滤掉。第二可解释你不再需要对几十个差异基因逐个编故事而是对 5 到 15 个因子做解释每个因子背后是一组有共性的特征讨论部分的写作难度直接下降一个量级。需要说明的是这套方法不是万能的。它擅长的是无监督探索——找出数据里主要的变异轴然后你去赋予生物学意义。它不擅长的是因果推断和精确的差异检验。如果你老板要的是“A 组和 B 组之间哪些基因显著差异”那还是老老实实做差异分析MOFA 帮不上忙。这个边界一定要在动手之前想清楚否则跑完一圈回来发现回答不了原始问题时间就白费了。1.2 什么样的数据组合适合上这套方法不是所有组学拼在一起都有意义。我的经验是样本量是硬约束。MOFA 本质上是因子模型参数估计依赖样本间的协方差结构样本太少的时候因子会非常不稳定同一批数据换个随机种子结果就变样。业内的常见经验是每组学至少几十个样本起步如果只有十几个样本跑出来的因子基本没法信不如老老实实做单组学。组学层之间的生物学关联强度也很关键。转录组和蛋白组通常是最佳拍档因为存在明确的翻译关系共享因子比较容易解释。转录组加甲基化也可以但要注意甲基化和表达常呈负相关解读的时候方向别搞反。代谢组和微生物组这种“上下游”关系耦合往往比较弱跑出来的共享因子占比可能很低更多是各层自己的独立因子这时候要有心理准备——不是工具不行是数据本身关联就弱。还有一个容易忽略的点样本对齐。多层组学必须来自同一批样本或者至少是严格配对的设计。我见过有人把两组不同队列的转录组和蛋白组硬拼在一起靠基因名匹配样本结果因子结构完全被批次主导。这种情况先做批次校正也救不回来因为样本层面的对应关系已经错了。数据整理阶段花时间核对样本 ID 和顺序比后面调参有用得多。2. 模型骨架拆解它到底在算什么2.1 从 PCA 到因子分析一句话讲清区别很多人第一次看 MOFA 的输出会懵因为它的结果形式和 PCA 很像——都有样本得分、都有特征载荷。区别在于PCA 是对一个矩阵做正交分解追求用最少的维度重构原始数据因子分析允许每个特征有自己的噪声方差而且因子之间不要求正交更接近“潜变量生成模型”的思维。MOFA 的模型可以写成这样对第 m 层组学观测矩阵 X_m 近似等于 样本因子矩阵 Z 乘以该层的权重矩阵 W_m 的转置再加上噪声项。Z 是跨层共享的W_m 是每层独立的。这就解释了一个核心设计样本层面的因子得分是共享的特征层面的权重是分层的。换句话说“哪些样本像”这件事在多层组学里是统一的“为什么像”则由每层组学各自给出证据。用生活化的比喻假设你在评价一批餐厅Z 是每家餐厅在“性价比”“环境”“口味”三个维度上的得分这个是共通的W_m 则是不同点评平台相当于不同组学层对这些维度的具体打分项——大众点评可能更看重环境外卖平台更看重性价比。平台之间打分项不同但餐厅在这几个底层维度上的位置是一致的。实现层面MOFA 用的是贝叶斯变分推断而不是最大似然或最小二乘。这个选择决定了它的很多行为特征能处理缺失值、能给出不确定性、能通过先验做稀疏化。代价是计算量大样本和特征一多跑一次可能要几个小时甚至过夜。所以后面讲实操的时候我会重点说怎么在保证结果可信的前提下把时间压下来。2.2 稀疏先验为什么它能自动关掉无关组学这是 MOFA 最值得讲的设计。如果只用普通因子分析每个因子会在所有特征上都分到一点权重解读的时候你会得到一张密密麻麻的热图根本看不出重点。MOFA 给权重矩阵加了spike-and-slab 先验通俗说就是每个权重都有很大概率被压到接近零spike只有少数有强证据的权重才会被放出来slab。最终效果是每个因子只在少部分特征上有非零权重热图上呈现出清晰的块状结构。这个机制对多组学特别友好。假设你有一个因子主要由转录组驱动蛋白组只是跟着轻微变化那稀疏先验会把蛋白组那部分的权重压到零附近你一眼就能看出这个因子是“转录主导型”。如果没有稀疏化蛋白组会分到 5% 的方差解释你很难判断这 5% 是真实信号还是噪声。还有一个配套机制叫ARD自动相关性判定作用在因子层面。简单说如果某个因子对任何一层组学都没有实质贡献ARD 会把它的整体权重压向零相当于自动“弃用”这个因子。这意味着你不需要精确猜因子数——设一个偏大的数比如 15 或 20让模型自己决定保留几个。实际跑下来设 15 个最后往往剩 8 到 12 个活跃因子这个数字比拍脑袋定要靠谱。注意稀疏化不是越强越好。先验强度调太高因子会过度稀疏本来该有信号的特征被误杀因子解释的方差比例会掉得很明显。默认参数通常够用除非你的数据特征数特别大比如甲基化位点几十万才需要考虑调高稀疏强度。2.3 方差分解R² 是怎么算出来、又该怎么看MOFA 最常被引用的图就是方差解释条形图——每个因子在每层组学上解释了多少方差。计算方式不复杂用因子重构出的矩阵去逼近原始矩阵1 减去残差方差比上总方差就是 R²。公式写出来是 R² 1 - Var(X - ZW^T) / Var(X)按层、按因子分别算。看这张图有几个要点。第一总方差解释率反映数据里到底有多少结构可挖。如果所有因子加起来只解释了 20% 的方差说明你的数据里大部分是噪声或者与技术因素相关这时候硬解读因子风险很大。第二单因子的跨层分布是判断因子性质的依据。一个因子在转录组解释 15%、在蛋白组解释 8%那是跨层共享因子价值最高只在单层解释 20%、其他层接近零那是该层特异的因子解读时要谨慎因为它可能只是该层内部的批次效应。第三注意方差的绝对尺度。不同组学的特征数和噪声水平差异巨大转录组动辄两万个基因甲基化可能几十万个位点蛋白组只有几千个蛋白。特征多的层天然容易拿到更高的方差解释比例这不代表它在生物学上更重要。我一般会同时看“每层组学的方差被解释了多少”纵向比较和“每个因子在各层的分布”横向比较两个视角结合才不容易误判。表格对比一下三种典型情况的判读方式这个是我自己整理的经验不是文档里的标准说法现象可能原因处理建议所有因子总解释率低于 30%噪声大、样本异质性低、或预处理过度检查质控和批次考虑放宽特征筛选某层组学方差几乎全被单因子吃掉该层存在强技术因素如测序深度对该层单独做校正后重跑因子数远少于设定值ARD 判定多数因子无贡献正常现象看剩余因子的解释率因子在各层分布高度不均衡组学间关联弱或样本量差异大降低期望重点看主导层的独立因子3. 数据准备与质控八成的坑都埋在这一步3.1 矩阵对齐与缺失值处理的实际操作MOFA 的输入是一个列表每个元素是一个矩阵行是特征、列是样本这个方向千万别搞反。我第一次用的时候按“样本×特征”组织跑出来的因子全是乱的排查了半天才发现是维度问题。列名必须严格一致顺序也要一致——工具不会自动帮你按样本 ID 匹配它只认位置。缺失值的处理要分情况。MOFA 的贝叶斯框架对高斯似然下的缺失值用 NA 或 NaN 表示是支持的推断时会把这些位置边缘化掉所以零星缺失不用管。但如果某层组学缺失比例超过三成这个层的结果可信度就要打折扣了。另外二值型数据比如突变有无要用伯努利似然计数型数据比如微生物丰度用泊松似然这两种似然对缺失值的容忍度低得多缺了就是缺了得事先补或者删特征。我一般会做一轮“矩阵体检”检查每层的特征数、样本数、缺失率、是否有全零行、是否有方差为零的特征。方差为零的特征在所有样本里表达都一样必须删掉它们对因子没有贡献还会拖慢计算。全零行也要删那是技术失败的信号。这套检查写成脚本跑一遍两分钟的事能省掉后面几小时的困惑。check_matrix - function(mat, name) { cat(sprintf([%s] 特征数: %d, 样本数: %d\n, name, nrow(mat), ncol(mat))) cat(sprintf( 缺失率: %.2f%%\n, 100 * mean(is.na(mat)))) cat(sprintf( 零方差特征: %d\n, sum(apply(mat, 1, var, na.rm TRUE) 0, na.rm TRUE))) cat(sprintf( 全缺失特征: %d\n, sum(rowSums(is.na(mat)) ncol(mat)))) invisible(NULL) }3.2 归一化与批次校正做多了反而坏事这里我要说一个反直觉的观点MOFA 内部已经做了尺度标准化默认scale_views TRUE会把每层组学的总方差拉到同一量级避免特征多的层主导结果。所以外部的归一化不需要做到“方差完全一致”这种程度做过头反而会抹掉真实的生物学差异。我推荐的做法是每层组学各自做适合该平台的标准化然后交给 MOFA。转录组用 TPM 或 VST蛋白组用 log2 转换后中位数归一化甲基化用 M 值而不是 Beta 值Beta 值在两端有压缩会扭曲方差结构。做完这些就停手不要再做跨层的 Z-score 或者联合标准化——那样会把各层特有的结构压平共享因子反而找不出来。批次校正是另一个雷区。如果每个组学层的批次结构不一样比如转录组和蛋白组是不同时间跑的那批次效应会成为一个“伪因子”横跨所有层但毫无生物学意义。这种情况必须在各层内部单独做批次校正用 ComBat 或者把批次作为协变量回归掉然后再进 MOFA。但如果你把批次也放进 MOFA 的模型里作为协变量效果往往不如预处理来得干净。提示批次校正做完后建议用 PCA 快速看一眼每层组学的前两个主成分确认批次分离现象消失了再继续。这一步花五分钟能避免后面因子被批次主导却毫不知情的情况。3.3 特征筛选的尺度把握特征越多越好吗不是。计算复杂度和噪声都会随特征数上升。甲基化数据几十万个位点直接塞进去跑一天都不一定收敛而且大部分位点是无关的。我的经验阈值是每层组学控制在几千到一两万个特征。转录组取高变基因比如方差前 5000 到 10000蛋白组如果本来就只有几千个就全留甲基化按变异系数或者差异甲基化区域筛到几千个。筛选标准要事先定好、写进方法学不能跑完看结果不好再回头改——那属于数据窥探。筛选方式上我倾向用无监督的方差筛选而不是差异分析筛选。原因很简单如果先用分组信息筛特征再去做无监督的因子分析等于把答案提前泄露给模型因子的独立性就没了。这一点在审稿时经常被问值得提前注意。4. 完整实操流程从原始矩阵到可解释的因子4.1 环境搭建与对象构建R 环境下主流程是MOFA2这个包它依赖 Python 后端通过 basilisk 或者 reticulate 调用所以第一次安装会比较折腾尤其是国内网络环境下装 Python 依赖容易失败。我的建议是先把 Python 环境配好再装 R 包装完后跑一个最小示例确认能通。if (!requireNamespace(BiocManager, quietly TRUE)) install.packages(BiocManager) BiocManager::install(MOFA2) library(MOFA2) # 数据矩阵行是特征列是样本列顺序必须完全一致 data_list - list( RNA rna_mat, # 例如 log2(TPM1) Protein prot_mat, # log2 强度 Methylation meth_mat # M 值 ) # 核对样本顺序 stopifnot(identical(colnames(data_list$RNA), colnames(data_list$Protein))) stopifnot(identical(colnames(data_list$RNA), colnames(data_list$Methylation))) MOFAobject - create_mofa(data_list) print(MOFAobject)create_mofa之后可以调用plot_data_overview()看一眼数据覆盖情况它会画出每个样本在各层的缺失状态一眼能看出有没有样本在某层整列缺失。这一步很重要我曾经遇到过某个样本在蛋白组里完全没有数据跑出来的因子在这个样本上得分异常排查很久才发现是输入问题。4.2 关键参数设置与因子数选择参数分三组数据选项、模型选项、训练选项。我的常用配置如下每一项都有理由。# 数据选项 data_opts - get_default_data_options(MOFAobject) data_opts$scale_views - TRUE # 平衡各层总方差防止特征多的层主导 data_opts$center_groups - TRUE # 按组中心化 # 模型选项 model_opts - get_default_model_options(MOFAobject) model_opts$num_factors - 15 # 先设偏大让 ARD 自动裁剪 model_opts$likelihoods - c(RNA gaussian, Protein gaussian, Methylation gaussian) model_opts$spikeslab_weights - TRUE # 权重稀疏化 model_opts$ard_weights - TRUE # 因子层面自动相关性判定 model_opts$drop_factor_threshold - 0.01 # 解释率低于 1% 的因子丢弃 # 训练选项 train_opts - get_default_training_options(MOFAobject) train_opts$convergence_mode - slow # 慢速收敛结果更稳 train_opts$maxiter - 1000 train_opts$seed - 42 # 固定种子保证可复现 train_opts$num_threads - 8 train_opts$verbose - TRUE关于num_factors我的做法是分两轮。第一轮设 15 到 20看 ARD 砍掉多少同时观察方差解释的累积曲线是否已经饱和。如果 12 个因子之后曲线基本走平说明 15 够用如果还在上升说明结构没挖完可以加到 25 再试。第二轮用最终确定的因子数重跑一次保证结果稳定。这个两步法比一次性猜准要可靠。convergence_mode的选择要权衡时间。fast适合快速探索确认流程能跑通medium是日常默认slow在写论文出最终结果时用因为它对变分推断的收敛更严格因子排序和权重更稳定。我用fast和slow跑同一批数据因子结构大体一致但个别因子的方差解释比例能差 2 到 3 个百分点这个量级在解读边界因子时是有影响的。4.3 训练与结果提取MOFAobject - prepare_mofa(MOFAobject, data_options data_opts, model_options model_opts, training_options train_opts) # 训练模型会存成 hdf5 文件方便复现和分享 MOFAobject - run_mofa(MOFAobject, outfile mofa_model.hdf5, use_basilisk TRUE)训练时间取决于数据规模。三层组学、每层几千特征、一两百样本8 线程大概十几分钟到半小时。如果超过两小时还没收敛八成是特征数太多或者因子数设得过大回头砍特征比干等更划算。训练完成后提取结果的三件套# 1. 样本因子得分 (N x K) factors - get_factors(MOFAobject, factors all)[[1]] # 2. 特征权重按层列出 weights - get_weights(MOFAobject, views all, factors all) # 3. 方差解释 r2 - calculate_variance_explained(MOFAobject) plot_variance_explained(MOFAobject, x view, y factor)factors这个矩阵是后面所有下游分析的起点。行是样本列是因子每一列可以当成一个连续的“综合指标”去和表型做关联。这里要注意因子得分是无序的不能因为因子 1 排在最前面就认为它最重要——排序是按方差解释量来的但生物学意义要看权重构成。4.4 因子解读的完整套路解读一个因子我一般走四步。第一步看方差解释分布确定它是跨层共享还是单层特异。第二步看权重排序取每层权重绝对值前 30 到 50 的特征看它们有没有功能上的共性。第三步做富集分析把权重高的特征对应的基因拿去做通路富集这时候往往能看出主题。第四步关联表型把因子得分和临床变量、分组信息做相关或差异检验。# 查看某个因子在每层的 top 权重特征 plot_weights(MOFAobject, view RNA, factor 1, nfeatures 30) # 因子得分与表型的关联 library(ggplot2) df - data.frame(factor1 factors[, 1], group sample_metadata$group) ggplot(df, aes(x group, y factor1, fill group)) geom_boxplot() theme_classic() labs(y Factor 1 score)富集分析这一步有个技巧不要只拿权重最高的那几个特征去做容易得到过于宽泛的结果比如全是“细胞周期”。我的做法是取权重分布的正负两端各取前 100 个分别做富集看正负方向的功能是否有差异。很多因子其实是双向的——一端上调一端下调单看一端会丢掉一半信息。可视化方面plot_data_heatmap能展示因子相关的特征在样本上的表达模式plot_factor能看因子得分在样本上的分布两者结合是标准的“因子故事”配图组合。热图默认按因子得分排序样本能明显看到块状结构这正是稀疏先验在起作用。5. 常见问题排查那些手册里不会写的情况5.1 结果不稳定、换种子就变样这是新手最容易慌的情况。同一个数据集种子 42 和种子 123 跑出来因子顺序不同、个别因子结构也不一样。先别急着怀疑数据。原因通常有两个。第一数据里的信号本身不够强几个方差解释率相近的因子在竞争随机初始化不同就会收敛到不同局部解。解决办法是把num_factors调小一点让模型专注于最强的那几根轴或者把convergence_mode设成slow多迭代。第二样本量相对特征数偏小参数估计的方差大。这种情况再怎么调参也只是缓解根本上是数据规模的限制。判断是否可信的一个土办法用三个不同种子各跑一次看因子间的一致性。如果主要因子的权重相关性在 0.8 以上说明结构是稳的如果相关性只有 0.3 到 0.5那这些因子就不能当结论写。我一般会做这个稳定性检查写进补充材料审稿人问起来也有据可依。5.2 因子被某一层组学独占跑完发现前几个因子全是转录组在解释蛋白组和甲基化几乎为零。这种情况先检查scale_views是不是开着。如果开着还这样那可能是数据本身的问题——某一层的噪声水平明显高于其他层导致共享信号被淹没。处理思路有两种。一种是对该层做更强的降噪比如提高特征筛选阈值、更严格的低表达过滤、或者用平滑方法甲基化可以按区域聚合而不是单 CpG 位点。另一种是接受现实分开解读跨层共享因子往往只占少数大量因子本来就是层特异的这在生物学上完全合理。强行追求“所有因子都跨层”反而会误导自己。5.3 因子和任何表型都对不上这个也很常见尤其是探索性研究。我的建议是先分清是“因子没意义”还是“表型选错了”。检查方式是把因子得分和所有可用的元数据都做一遍关联包括年龄、性别、批次、测序深度这些技术变量。如果因子和技术变量强相关那它是技术伪影需要校正重跑。如果确实和任何变量都无关也别急着放弃。无监督方法找到的因子可能是数据里真实存在但你的研究设计没有覆盖到的变异轴比如细胞类型组成的变化、或者某个未知的亚群结构。这时候可以拿因子得分去做聚类看能不能分出新的样本亚群再用亚群的差异特征反推生物学意义。这条路我走过一次最后发现某个因子对应的是样本的细胞纯度差异虽然和研究假设无关但反过来解释了为什么某些样本的信号偏弱。为了便于查阅把我踩过的坑整理成速查表现象排查顺序典型解法训练不收敛或超时特征数 → 因子数 → 线程数削减特征至几千、因子数降到 10-15因子 1 解释率异常高40%检查该层是否含强技术因素层内批次校正后重跑大量因子被 ARD 丢弃检查特征筛选是否过度放宽筛选阈值保留更多变异特征因子得分与分组完全重合检查是否泄漏了标签信息确认特征筛选为无监督方式热图无明显块状结构检查稀疏化是否生效确认spikeslab_weights TRUE某样本因子得分极端检查该样本各层数据完整性缺失过多则剔除该样本5.4 计算资源不够怎么办三层组学、每层一万特征内存占用能到几十 GB。如果机器扛不住有几个实用手段。降特征是最有效的把每层压到三千到五千内存需求直接降一个量级而主要因子结构通常不会变。分组处理如果样本来自多个批次或条件可以用 MOFA 的 group 机制分组训练它会为每组估计独立的因子均值内存和时间的增长比整体训练温和。分步训练先只用两层组学跑出因子骨架再把第三层加进来做投影虽然严格性差一点但应急够用。我个人的配置参考16 GB 内存的机器三层组学各三千特征、150 样本8 线程、medium收敛模式大约二十分钟跑完峰值内存 6 到 8 GB。这个规模对大多数课题来说是够用的。6. 一点个人体会这套方法真正的价值在哪用了一段时间之后我对 MOFA 的定位有了比较清晰的认识。它最大的价值不是给出一个“标准答案”而是逼着你从整体结构的角度去看数据。以前做单组学习惯盯着 top 差异基因思路是“找最显著的那个点”现在会先问“这批样本主要的变异方向有几个每个方向由什么驱动”。这个思维方式的转变比工具本身更有价值。另一个体会是关于解读的克制。因子分析出来的东西是相关性结构不是因果。我看到过不少把因子直接命名成某个通路名字的做法风险很大——一个因子可能同时包含多个通路的特征硬套一个名字会限制后续思考。我自己的习惯是给因子起一个描述性的中性名字比如“因子 3增殖相关跨转录与蛋白”保留解释空间等下游验证做完再下结论。还有一点关于结果的呈现。如果要把这套分析写进论文或者分享给同行方差解释图、因子权重热图、因子-表型关联图这三张是基本配置。补充材料里放上种子稳定性检查、参数敏感性分析会让结果的可信度高很多。这些内容不占正文篇幅但能省掉很多解释成本。最后分享一个我在实际项目里总结的小技巧跑完 MOFA 之后把因子得分当成新的变量回过头去做单组学的差异分析用因子得分分组而不是用原始分组。这样做的好处是你能找到“驱动某个变异轴的具体基因”而不是“两组之间的所有差异基因”后者往往包含大量与该变异轴无关的噪声。这个方法我用过几次找出来的基因列表短而集中后续验证的成功率高不少。它不是什么标准流程但作为思路扩展我觉得挺值得一试。