3. 批次校正基本概念：校正目标、数据类型适用范围
说实话，我刚接触批次校正时也被各种术语搞晕过。校正目标、数据类型适用范围，听起来像两个独立话题，其实它们是同一枚硬币的两面。你只有搞清楚“为什么要校”和“什么样的数据能校”，才能真正选对方法。
3.1 批次效应的本质与校正目标
先说说批次效应从哪来的。你去测一批样本，今天跑一板、明天跑一板，试剂批号不一样、操作手不一样、仪器状态不一样——这些都会引入系统性偏差。我曾在一次跨中心合作项目里见过：同一份细胞样本，在两个实验室测出来的表达谱，聚类后居然分成了两群。不是生物差异，纯粹是技术噪声。
校正的核心目标就三个：
- 消除技术变异——让相同生物学条件下的样本，表达谱尽可能一致。
- 保留生物变异——千万别把真实的组间差异也校掉了，那就本末倒置了。
- 提高可比性——使得来自不同批次、不同平台甚至不同实验室的数据可以放在一起分析。
嗯，这里要特别强调第二点。你想想看，如果一个校正方法把正常组和疾病组的差异都抹平了，那还要它干嘛？所以我个人习惯是，在正式校正前先画一张原始数据的PCA图，看看批次的分离方向和生物学分组的分离方向是否正交。如果它们混在一起，那校正的难度就大很多。
3.2 你要校正的是什么数据？
不同类型的数据，校正策略完全不同。我曾经踩过一个坑：把bulk RNA-seq的方法直接套在scRNA-seq上，结果细胞亚群结构全乱了。为什么呢？因为它们的特征维度、稀疏程度、噪声模型都不在一个频道上。
| 数据类型 | 典型特征 | 主要噪声来源 | 校正常见限制 |
|---|---|---|---|
| Bulk RNA-seq | 基因×样本的矩阵，每样本一个表达值 | 文库大小、GC含量、测序深度 | 需要足够的样本量；不能处理单细胞级别的异质性 |
| scRNA-seq | 基因×细胞的稀疏矩阵（dropout率高） | 捕获效率、扩增偏倚、批次特异性缺失 | 零膨胀模型、细胞数量差异大；部分方法要求至少2个批次 |
| Microarray | 探针×样本的连续信号值 | 芯片间归一化、染料效应、扫描仪偏移 | 探针集不兼容时需先交叉映射；无法解决平台间的系统偏差 |
你看表格里列出的限制，其实就对应了校正方法的设计初衷。Bulk RNA-seq 的校正思路偏向于全局缩放或分位数对齐，比如ComBat、limma等。而scRNA-seq需要处理大量缺失值，所以像Harmony、Seurat的CCA、以及Scanorama这些方法，都采用了不同的策略来对齐细胞嵌入空间，而不是直接对齐表达值。
我曾经犯过一个错：用Microarray的RMA归一化流程去处理bulk RNA-seq的count数据。结果发现低表达基因的分布严重变形。后来才搞清楚，不同平台的信号分布假设完全不一样——Microarray的强度值是log2变换后近似正态，而RNA-seq的count数据是负二项分布。校正前一定要先搞清楚输入数据的分布形态。
3.3 数据类型的适用范围：一张图讲清楚
为了让你一眼看清楚，我画了一张框架图。它展示了批次校正的基本流程以及不同数据类型对应的推荐方法范围。
这张图清楚地展示了：不是所有方法都能通吃。Microarray的历史最久，方法也最成熟，但千万别把它那一套直接搬到单细胞数据上。我碰到过有人在单细胞数据上用标准化后的ComBat——结果细胞群被强行拉平，本来稀有细胞亚群信号直接没了。
3.4 什么时候该做校正？
说白了，校正不是默认动作。我个人的经验是：先问自己三个问题。
- 批次和生物学分组是否完全混杂？ 比如批次1全是疾病组，批次2全是正常组——这种设计你几乎没法校正，因为技术效应和生物效应完全纠缠在一起。我曾经花了两周硬跑各种校正方法，结果全是徒劳。最后只能重新设计实验。
- 是否需要整合多个来源的数据？ 比如你要把GEO里三个不同研究的数据合并分析，那校正就是必选项，否则你会看到“研究来源”成为最大的聚类因子。
- 你的下游分析对批次敏感吗？ 如果是做差异表达，而且每个批次内都包含所有条件，那简单的线性模型加上批次作为协变量可能就够用了。但如果你是要做聚类、轨迹推断，那非要仔细校正不可。
一个小技巧：在跑任何校正方法之前，先画一个分面箱线图，x轴是样本，y轴是某个管家基因的表达，分面按批次。如果箱线图的中位数在不同批次间差了一倍以上，那说明批次的效应很强，需要认真对待。如果只差10%~20%，也许简单的分位数归一化就能解决。
3.5 校正前的数据准备：统一基因/探针命名
这一点常被忽略。不同批次的数据，可能用了不同版本的基因注释（比如Ensembl v75 vs v100），或者不同平台用了不同的探针ID。我见过最夸张的一次，同一个基因在一份数据里叫"BRCA1"，另一份叫"ENSG00000012048"。如果不先做ID映射，校正的结果就是一团乱麻。
一般流程是这样的：
# 假设你有一个表达矩阵 expr1 (基因×样本1)
# 另一个表达矩阵 expr2 (基因×样本2)
# 使用 biomaRt 或自定义映射表统一基因名
# 一个简单的示例：将Ensembl ID转为HGNC symbol
library(biomaRt)
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
gene_map <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
filters = "ensembl_gene_id",
values = rownames(expr1),
mart = ensembl)
# 保留有对应symbol的基因，去重复
# 然后合并表达矩阵
common_genes <- intersect(gene_map$hgnc_symbol, rownames(expr2))
expr1_common <- expr1[match(common_genes, gene_map$hgnc_symbol), ]
expr2_common <- expr2[common_genes, ]
merged_expr <- cbind(expr1_common, expr2_common)注意哦，这个步骤非常容易出错，尤其是当基因名有多个别名时。我会建议你用NCBI的官方symbol或者Ensembl的稳定ID作为基准，避免使用像"1-Dec"这种易混淆的符号。
3.6 一个新手常见的误区
很多人以为校正就是一步到位，把数据扔进函数就完事了。其实不是这样的。你校正完之后，必须做验证。我习惯的验证流程包括：
- 画校正前后的PCA图，检查批次是否还聚在一起。
- 计算批次间的Kolmogorov–Smirnov统计量（或者更简单的相关系数），看看校正后批次分布是否更接近。
- 保留一些已知不会受批次影响的“内参”基因（比如管家基因），检查它们在各个样本中的表达变异是否降低。
特别注意：如果你用的是scRNA-seq数据，校正之后一定要检查细胞亚群的比例是否改变。我有一个惨痛教训——用某个流行方法校正后，本来应该有5%的稀有细胞群，校正后只剩0.5%了。后来发现那个方法对低表达基因的偏差校正过度，把稀有群的信号稀释掉了。
好了，这一章我们聊了校正的目标和不同数据类型的脾气。记住一句话：没有银弹。每种方法都有它成立的假设，只有当你手上的数据符合这些假设时，校正才能帮你，而不是害你。下一回我们聊聊具体的方法论——ComBat系列，那才是真正动手的起点。