多组学整合与批次效应

引言

多组学整合（multi-omics integration）试图回答一个单组学无法回答的问题：「同一个生物系统在不同分子层面上是如何协同变化的」。转录组告诉你「哪些基因被转录」，蛋白组告诉你「实际有多少蛋白」，代谢组告诉你「最终产生了哪些小分子」，表观组告诉你「为什么这些基因被打开或关闭」。把它们放在一起，才能拼出从基因型到表型的完整链条。

但「放在一起」这件事在工程上极其困难。三个根本障碍贯穿始终：第一，量纲与分布完全不同——转录组是计数（负二项分布）、蛋白组是强度（对数正态）、代谢组是峰面积（高度偏态且缺失严重），任何直接的数值运算都缺乏统计基础。第二，样本量与维度严重不匹配——多组学实验通常只有几十个样本，而每个组学的特征是数万维，「维数灾难」让大多数多元统计方法失效。第三，批次效应在组学间不成比例——转录组与蛋白组的样本往往在不同时间、不同平台处理，批次结构不同，校正时需要区别对待。

最容易被低估的是批次效应。在多组学项目里，它往往不是「干扰项」而是「主导项」：一个未经校正的整合分析，聚类结果可能完全由「哪个批次」决定，而研究者却把它解读为「生物学亚型」。这不是统计细节，而是决定结论真伪的核心工程问题。

本文按「整合目标 → 共同挑战 → 批次效应 → 校正方法 → 整合策略 → 降维聚类 → 因子分析 → 网络与因果 → 验证与复现」的顺序展开。代码以 R（limma、sva、MOFA2、mixOmics）与 Python（scanpy、scikit-learn）为主。读完你应该能设计一个多组学项目，并判断整合结果是「真信号」还是「批次伪影」。

目录

1. 多组学整合的目标与三种范式
2. 组学数据的共同挑战：高维、小样本、异质
3. 批次效应：来源、检测与危害
4. 批次校正方法：ComBat、limma、RUV
5. 早期、中期与晚期整合策略取舍
6. 降维与聚类：PCA、NMF 与相似度网络
7. 多组学因子分析：MOFA 与 DIABLO
8. 调控网络与因果推断
9. 整合结果的验证与可复现性

1. 多组学整合的目标与三种范式

整合的目标可以归为三类，方法选择取决于目标：

「同一样本的多层描述」（样本为中心的整合）：把每个样本的多个组学特征拼成一个「多层特征向量」，用于分型、聚类、分类。典型问题：「这些肿瘤样本能分成几个亚型？每个亚型在多组学层面有什么共同特征？」代表方法：MOFA、iCluster、相似度网络融合（SNF）。

「跨组学的关联发现」（特征为中心的整合）：找「哪些转录本与哪些蛋白/代谢物协同变化」，进而推断调控关系。典型问题：「哪些 miRNA 调控哪些 mRNA？」「哪些代谢物与哪些基因表达相关？」代表方法：DIABLO、典型相关分析（CCA）、WGCNA 的多组学扩展。

「从基因型到表型的因果链」（因果推断）：用遗传变异作为工具变量，推断「基因表达 → 蛋白 → 代谢物 → 疾病」的因果方向。典型问题：「这个 eQTL 是通过影响基因表达来影响疾病的吗？」代表方法：孟德尔随机化（MR）、中介分析、eQTL 共定位。

三类目标的统计要求差异很大：第一类需要「降维 + 整合」，第二类需要「多变量关联检验」，第三类需要「因果假设检验」。先明确目标再选方法，是避免「用了复杂方法却回答了错误问题」的前提。

2. 组学数据的共同挑战：高维、小样本、异质

四个挑战决定了多组学整合的方法学形态：

高维小样本（p » n）。一个典型的多组学实验：n = 50 个样本，每个样本有 20000 个基因、5000 个蛋白、300 个代谢物。总特征数远超样本数，导致：协方差矩阵不可逆（无法直接做多元回归）、任何模型都容易过拟合、多重检验的校正极其严苛（2.5 万个特征 × 3 个组学 = 7.5 万次检验）。

应对手段是先降维再整合：用 PCA/因子分析把每个组学压到 10-50 个维度，再做跨组学关联。代价是「降维可能丢失与主要方差方向无关但生物学重要的信号」——一个低方差但关键的调控通路可能在前几个主成分里完全看不到。

异质分布。不同组学的数据分布形态不同，需要各自做标准化（transformation）：

关键原则：变换的目的是「让分布近似正态且方差稳定」，而不是「让不同组学可比」。后者是徒劳的——转录组与蛋白组的量纲本质上不可比，任何「拉齐量纲」的操作都会引入人为结构。跨组学比较应该在「变化方向与排序」层面做，而非原始数值层面。

缺失值模式不同。转录组的缺失极少（测序深度足够时），蛋白组的缺失常见（低丰度蛋白检不到），代谢组的缺失最多（20%-50% 常见）。缺失模式的差异会直接影响整合——若某个组学的缺失与「样本分组」相关，填补会引入偏差。

样本不匹配。并非所有样本都测了所有组学。部分重叠的样本结构（如 50 个样本做了转录组，其中 30 个做了蛋白组）是多组学研究的常态，需要专门的方法（如缺失组学下的因子分析）。

3. 批次效应：来源、检测与危害

批次效应（batch effect）指「由技术因素而非生物学因素导致的系统性差异」。它的来源包括：样本处理日期、操作人员、试剂批次、仪器、测序 lane、文库制备批次、甚至「样本在 96 孔板上的位置」。

为什么它在多组学里特别危险：单组学分析中，批次效应通常表现为「样本按批次聚集」，容易被察觉；而在多组学整合中，批次效应会伪装成「跨组学的一致性信号」——如果转录组和蛋白组样本在同一批处理，那么两个组学的批次效应会「协同」，整合分析会看到「多组学层面的一致变化」，从而误判为「真实的生物学亚型」。这是多组学文献中最常见的致命错误。

检测批次效应的手段：

# 1. 用 PCA 看主成分是否与批次相关（最直观）
pca <- prcomp(t(expr_matrix), scale. = TRUE)
plot(pca$x[,1:2], col = as.factor(metadata$batch), pch = 19)
# 若前几个主成分与批次高度相关 → 存在严重批次效应

# 2. 用 PVCA（Principal Variance Component Analysis）量化方差来源
library(pvca)
pvca_obj <- pvcaBatchAssess(assayData, metadata,
  threshold = 0.6, theProp = c("batch", "condition"))
# 输出各因素解释的方差比例：批次占 40% 就是严重问题

# 3. 用 kBET / LISI 等专门指标评估「批次混合程度」

判据的量化参考：批次解释的方差 < 5% 可忽略；5%-20% 需校正；> 20% 说明设计有严重问题（可能批次与生物学因素混杂），校正只能部分补救。

最危险的情形是「批次与分组完全混杂」（complete confounding）：如果所有对照组样本都在第一批处理、所有处理组都在第二批，那么「批次效应」与「处理效应」在数学上完全无法区分。这不是分析能解决的问题，只能重新设计实验。防患的做法是「随机化」——把不同组的样本均匀分配到各批次（如交替上样），让批次与分组正交。

4. 批次校正方法：ComBat、limma、RUV

三类方法，原理与适用场景不同：

ComBat（经验贝叶斯校正） 是最常用的方法，它把批次效应建模为「位置偏移 + 尺度缩放」，用经验贝叶斯估计每个基因在每个批次的校正参数，同时保留生物学协变量（在模型里显式指定）：

library(sva)
# mod：生物学协变量（如分组），必须显式传入，否则生物学信号会被当批次校正掉
mod <- model.matrix(~ condition, data = metadata)
expr_combat <- ComBat(dat = expr_matrix, batch = metadata$batch,
                      mod = mod, par.prior = TRUE, prior.plots = FALSE)

mod 参数是最关键的设置。若忘记传入生物学协变量，ComBat 会把「处理组与对照组的差异」当成批次效应一并抹掉——这是 ComBat 最经典的误用。ComBat 还有两个变体：ComBat-seq（针对 RNA-seq 计数，用负二项模型，不破坏整数性质）、ComBat 的非参数版本（par.prior = FALSE，适合小批次）。

limma 的 removeBatchEffect 把批次作为协变量放进线性模型，直接估计并移除批次贡献：

library(limma)
design <- model.matrix(~ condition + batch, data = metadata)
# 注意：这里把 batch 当协变量，条件效应会被保留
expr_corrected <- removeBatchEffect(expr_matrix, batch = metadata$batch,
                                    design = model.matrix(~ condition, metadata))

limma 的方法更透明（就是线性模型的系数移除），但不做「尺度校正」（只移除均值偏移，不移除方差差异）。ComBat 同时校正位置与尺度，效果通常更好，但对小批次（< 10 个样本）的估计不稳定。

RUV（Remove Unwanted Variation） 走的是另一条路：用「负对照」估计不需要的变异。负对照有两类：RUVg 用「管家基因」（假设它们不该变化），RUVs 用「样本内重复」（技术重复的差异纯属噪声）。它不需要预先知道批次标签，适合「批次未知但存在系统性变异」的场景：

library(RUVSeq)
# RUVs：用技术重复样本估计不需要的变异因子
differences <- makeGroups(metadata$sample_group)
ruv_result <- RUVs(expr_matrix, cIdx = genes_of_interest,
                   k = 3, differences, isLog = TRUE)
expr_ruv <- ruv_result$normalizedCounts

k 参数（要移除的变异因子数）是 RUV 的核心权衡：k 太小校正不足，k 太大把真实生物学信号也移除。经验做法是「尝试 k = 1-5，看负对照的方差是否降到平台、而正对照（已知的差异基因）是否仍保留」。

Harmony 的思路值得单独说明：它在「降维后的嵌入空间」做校正，通过迭代聚类找到「同一细胞类型跨批次的对应点」并拉近。这一方法在 单细胞测序数据分析 里是标配，但它只适用于「有明确对应关系」的场景（同一细胞类型在不同批次都存在）。对于多组学整合，Harmony 的适用性有限——不同组学之间没有「同一特征的对应点」。

5. 早期、中期与晚期整合策略取舍

按「在哪个阶段合并数据」分类：

早期整合（early / concatenation）：把各组学的特征矩阵横向拼接成一个矩阵，再做单一模型。优点是最简单、能捕捉跨组学的交互；缺点是量纲不可比的问题被放大——若不做权重处理，特征数多的组学（转录组 20000 维）会完全主导结果（蛋白组 5000 维）。

中期整合（intermediate / 因子模型）：假设各组学由「共同潜在因子」驱动，从各组学中同时估计这些因子。代表方法：MOFA、iCluster、相似度网络融合。优点是有明确的统计模型、能处理部分缺失的组学；缺点是计算复杂、超参数多。

晚期整合（late / 结果层面）：各组学独立分析，最后比较结果（如「差异基因与差异蛋白的重叠」）。优点是简单、各组学方法成熟、可解释性强；缺点是丢失跨组学的协同信息——两个组学各自都不显著但协同显著的变化，晚期整合永远看不到。

实践建议：从晚期整合开始。先做各组学的独立分析，确认每个组学都有信号，再用中期方法寻找协同模式。直接上中期/早期整合而跳过单组学质控，是「用复杂方法掩盖数据问题」的典型。

6. 降维与聚类：PCA、NMF 与相似度网络

多组学整合的降维有几种思路：

各组学分别 PCA 后拼接：把每个组学的前 N 个主成分拼起来，做联合聚类。简单但主成分的「方差解释比例」在不同组学间不可比，需要标准化。

NMF（非负矩阵分解）：把非负的数据矩阵分解为「特征 × 模块」与「模块 × 样本」两个非负矩阵，天然给出「模块（可解释为通路/程序）」与「样本权重」。NMF 的模块比 PCA 的主成分更可解释（因为非负，模块是「叠加」而非「抵消」），是「多组学模块发现」的常用工具。

相似度网络融合（SNF）：为每个组学单独构建样本相似度网络，再用非线性方法融合成一个网络，最后聚类：

# SNF 的 Python 实现（snfpy）
import snf
from sklearn.metrics import pairwise_distances

# 每个组学一个特征矩阵 → 相似度矩阵
affinities = []
for omics_matrix in [rna_matrix, protein_matrix, methyl_matrix]:
    dist = pairwise_distances(omics_matrix, metric="euclidean")
    affinities.append(snf.make_affinity(dist, K=20, mu=0.5))

# 融合 → 谱聚类
fused = snf.snf(affinities, K=20)
labels = snf.spectral_clustering(fused, n_clusters=3)

SNF 的优势是不假设组学间有线性关系，且对「某个组学在某个样本上缺失」较鲁棒。它输出的聚类标签可以直接用于分型验证（如「这个亚型在三个组学层面都有一致特征」）。

一个必须做的诊断是「聚类是否由单一组学驱动」：分别用每个组学做聚类，看标签的一致性。若融合聚类的结果与「仅用转录组聚类」几乎一致，说明其他组学没有贡献新信息——这时应如实报告，而不是宣称「多组学整合发现了新亚型」。

7. 多组学因子分析：MOFA 与 DIABLO

MOFA（Multi-Omics Factor Analysis） 是中期整合的代表方法。它把每个组学建模为「若干潜在因子」的线性组合，用变分推断同时估计因子与载荷：

# MOFA2 的 Python 接口
from mofapy2.run.entry_point import entry_point

ent = entry_point()
ent.set_data_options(scale_views=False)
ent.set_data_matrix(data=[rna, protein, metabolite],
                    views=["RNA", "Protein", "Metabolite"])
ent.set_model_options(factors=10, spikeslab_weights=True)
ent.set_train_options(iter=1000, convergence_mode="fast", seed=42)
ent.build(); ent.run()
ent.save("mofa_model.hdf5")

MOFA 输出的关键诊断：

• 方差解释比例（R² per view per factor）：哪个因子解释了哪个组学的多少方差。一个「多组学因子」应该在多个组学里都有解释力——若某因子只在转录组里有 R²，那它只是转录组的内部结构，不是多组学因子。
• 因子与样本协变量的关联：因子是否与分组、临床变量相关（用线性模型检验）。
• 载荷（weights）：每个组学中哪些特征驱动该因子（可用于通路富集）。

DIABLO（mixOmics 包）走的是「有监督整合」：给定分组标签，寻找「能最大程度区分组别」的跨组学特征组合。它的优势是直接服务于「分类/生物标志物发现」，且给出「各组学特征之间的相关性网络」：

library(mixOmics)
# 设计矩阵：指定各组学之间的关联强度（0-1）
design <- matrix(1, ncol = 3, nrow = 3,
                 dimnames = list(c("RNA","Protein","Metab"), c("RNA","Protein","Metab")))
diag(design) <- 0

diablo_res <- block.splsda(X = list(RNA=rna, Protein=prot, Metab=metab),
                           Y = factor(metadata$condition),
                           ncomp = 3, design = design, keepX = list(RNA=50, Protein=30, Metab=20))
# 交叉验证选最优组分与特征数
perf_res <- perf(diablo_res, validation = "Mfold", folds = 5, nrepeat = 10)

DIABLO 的最大风险是过拟合：在有监督框架下，只要特征数足够多，总能找到「完美区分组别」的组合。必须用严格的交叉验证评估性能（perf 函数），并报告「交叉验证的分类错误率」而非「训练集准确率」。若交叉验证错误率接近随机水平（如二分类 50%），说明没有真实的判别信号。

8. 调控网络与因果推断

关联分析只能给出「A 与 B 相关」，而生物学关心的是「A 调控 B」。从关联到因果需要额外的信息：

eQTL 与共定位（colocalization）。eQTL 分析找「遗传变异与基因表达的关联」，若某变异同时是 eQTL（影响表达）与 GWAS 信号（影响疾病），则提示「变异通过影响表达来影响疾病」。共定位分析（coloc 包）用贝叶斯方法评估「两个信号是否由同一个因果变异驱动」，比「简单地看是否重叠」严格得多：

library(coloc)
res <- coloc.abf(dataset1 = list(beta = eqtl_beta, varbeta = eqtl_varbeta,
                                 snp = snps, N = eqtl_n, type = "quant"),
                 dataset2 = list(beta = gwas_beta, varbeta = gwas_varbeta,
                                 snp = snps, N = gwas_n, type = "cc",
                                 s = case_fraction))
# 看 PP.H4（共享因果变异的后验概率），> 0.8 为强证据

孟德尔随机化（MR） 用遗传变异作为「工具变量」推断因果关系。核心假设有三条：工具变量与暴露强相关、与混杂因素独立、只通过暴露影响结局（排除多效性）。第三条假设最难满足——基因多效性（一个变异影响多个性状）在人类基因组中极其普遍，因此 MR 必须做多效性检验（如 MR-Egger、MR-PRESSO）。

中介分析（mediation analysis） 用于「检验 A → B → C 的间接效应」。在多组学语境下，可以检验「基因型 → 表达 → 蛋白 → 疾病」的链条。但中介分析对模型假设极其敏感（无未测混杂、正确的函数形式），用观测数据做中介分析得出的因果结论必须谨慎。

一个务实的原则：多组学网络推断的结果应定位为「假设生成」而非「结论」。共表达网络、调控网络给出的边是「统计关联」，需要实验验证（如敲除实验、ChIP 验证）才能升级为「调控关系」。

9. 整合结果的验证与可复现性

多组学整合的结论必须经过验证，且验证标准比单组学更高（因为引入了更多自由度）：

内部验证。用交叉验证评估整合模型的实际预测能力（不是拟合能力）。特别注意「批次与分组混杂」的假验证——若随机划分训练/测试集时批次结构不同，交叉验证会给出虚高的性能。

独立数据验证。用独立队列（不同医院、不同平台）验证发现的亚型或标志物。这是唯一能排除「批次伪影」的方法——如果一个亚型只在原队列中成立，在新队列里消失，那么它极可能是批次效应的产物。

可复现性的工程要求。多组学项目涉及的工具与参数比单组学更多，因此可复现性挑战更大：

多组学项目的可复现性清单
  1. 每个组学的原始数据与预处理参数（版本、软件、参数文件）
  2. 样本元数据：批次、处理日期、平台、协变量（必须完整）
  3. 随机种子：降维、聚类、交叉验证都涉及随机性
  4. 校正方法与其超参数（k 值、ComBat 的 mod 公式）
  5. 中间产物：每个组学校正前后的矩阵都保存
  6. 版本锁定：容器镜像 + 依赖清单

参见 生信可复现性工程 了解如何用容器与流程引擎实现这套要求。特别要强调**「保存校正前后的矩阵」**：批次校正的效果需要对照才能评估，只保留校正后的数据会让「校正是否过度」无法回溯。

过度校正的识别。批次校正的共同风险是「把生物学信号也校掉」。检测方法：把已知的生物学差异（如阳性对照基因、已知的疾病标志物）在校正前后做对比——若某个已知的强信号在校正后消失，说明校正过度。

权衡取舍

常见坑清单

1. 批次与分组完全混杂：校正无法分离两者；只能重新设计实验，分析层面无解。
2. ComBat 忘记传 mod：生物学差异被当批次抹掉；必须显式指定生物学协变量。
3. 早期整合不做权重：高维组学完全主导结果；先各组学降维或标准化到相同维度。
4. 用「多组学一致变化」证明真实信号：同批处理的组学批次效应会协同；必须先校正再比较。
5. DIABLO 只报训练集准确率：过拟合导致虚高性能；必须用交叉验证并报告错误率。
6. RUV 的 k 设太大：真实生物学信号被移除；用负对照方差平台 + 正对照保留双重判据。
7. 缺失值填补忽略 MNAR：低丰度蛋白的缺失被均值填补，信号被抹平；按缺失类型分别处理。
8. 聚类后不检查「是否单一组学驱动」：宣称多组学新发现实际只是转录组结构；分别聚类对比。
9. 中介分析忽略未测混杂：因果结论不可靠；用 MR 或敏感性分析检验稳健性。
10. 只保留校正后矩阵：无法回溯校正是否过度；校正前后都要保存。

小结

多组学整合的困难不在「用什么方法」，而在「数据能否支撑结论」。高维小样本让模型极易过拟合，异质分布让跨组学数值比较缺乏统计基础，而批次效应在多组学场景下会伪装成「跨组学一致性信号」——这三点决定了多组学分析必须以「严谨的实验设计与充分的质控」为前提，而非以「复杂的方法」为卖点。

工程上最该建立的三个认知：第一，批次是设计的产物，不是分析的产物。随机化上样、把组间比较放在同批次内，比任何校正算法都有效；校正只能补救，不能创造。第二，整合策略应从简单到复杂。先做各组学独立分析（晚期整合），确认每个组学都有信号，再用中期方法寻找协同模式。第三，验证标准随自由度上升。多组学引入了更多参数与选择，因此必须有独立队列验证——只在原队列成立的发现，很可能只是批次伪影。

下一步可以看 生信可复现性工程 了解如何把多组学的预处理、校正与整合步骤纳入可审计的流程，或 RNA-seq 转录组分析 了解单组学差异分析的基础——多组学整合的第一步，永远是先把每个组学单独做对。如果你的整合聚类结果「完美分开但生物学上讲不通」，第一件事就是检查批次与分组是否混杂。

常见问题

Q：样本量多少才够做多组学整合？
没有绝对数字，但有经验参考：无监督整合（MOFA、SNF）通常需要 ≥ 50 个样本才能稳定估计因子；有监督整合（DIABLO）每组至少 15-20 个样本，否则交叉验证的方差极大。样本量不足时，务实的选择是「减少组学数量、只做两两关联分析」，而不是硬上复杂的整合模型。

Q：不同组学的样本不完全重叠怎么办？
三种处理：一是只保留所有组学都测了的样本（样本量损失大）；二是用容忍缺失的方法（MOFA 支持部分缺失的组学）；三是分别建模后再对齐（晚期整合）。不要对缺失的组学做「样本层面的插补」——用其他样本的平均值填一个缺失样本的多组学特征，等于人为制造了一个「平均样本」，会严重扭曲聚类结果。