第三章 数据准备：打好ComBat校正的地基
说实话，很多人在跑ComBat之前翻车，十有八九都栽在数据准备这一步。我见过不下五次——师兄师弟拿着漂亮的表达矩阵来找我，结果校正完发现批量效应还在，仔细一查，啊，样本顺序和批次标签对不上。所以这一章咱们把准备工作做扎实，后面才能顺风顺水。
3.1 表达矩阵的格式要求
ComBat吃的输入，说白了就是一个标准的表达矩阵。行是基因（或探针），列是样本。数值必须是经过归一化的表达量，我一般习惯用log2(TPM+1)或log2(FPKM+1)。
# 一个正确的表达矩阵示例
# 行名: 基因ID, 列名: 样本名
Sample_1 Sample_2 Sample_3 Sample_4
GeneA 7.23 6.89 8.01 7.55
GeneB 4.12 4.55 3.98 4.33
GeneC 9.87 10.12 9.45 9.90三个要点你得记牢：
- 数值必须是连续的，不能有counts原始数据（除非你做了log转换）。
- 不能有负值，ComBat内部会做方差建模，负值会让模型跑飞。
- 不能有Inf或NaN，这个后面会专门讲缺失值处理。
我在项目中遇到过一位同事，直接拿原始counts塞进去，结果校正完所有样本的表达量变成了负几百。究其原因，counts的方差结构与log2尺度完全不同，ComBat假设方差齐性，原始counts显然不满足。所以请务必先做log归一化。
3.2 批次标签的准备工作
批次标签说白了就是一列向量，长度等于样本数，每个元素告诉模型这个样本来自哪一批。注意，批次标签的顺序必须与表达矩阵的列顺序严格一致。怎么强调都不过分——我曾经因为这个踩过大坑，代码跑完才发现标签错位，白白浪费两天。
# batch向量示例，假设5个样本来自两批
batch <- c(1, 1, 2, 2, 2)
# 或者用字符
batch <- c("Batch1", "Batch1", "Batch2", "Batch2", "Batch2")一个常见的误区是，有人把测序的lane号、日期、操作者都混在一起作为单一批次。我个人建议先做一个简单的PCA看看——如果批次标签能很好地把样本分开，说明这个标签值得用；如果分不开，反而要考虑是不是批次标签定义得太粗糙。
3.3 协变量的处理与编码
有时候你希望保留某些生物学变量的效应，比如年龄、性别、疾病状态。那你就得把这些变量作为协变量传给ComBat。这部分数据也需要整理好，格式上就是个data.frame，行数等于样本数，列是各个协变量。
# 协变量示例：年龄（连续） + 性别（分类）
covariate <- data.frame(
age = c(45, 52, 38, 61, 47),
sex = c("M", "F", "M", "F", "M")
)这里有两点值得注意：
- 分类变量最好做dummy编码，ComBat内部虽然能处理因子，但dummy编码能避免意外的排序问题。
- 小心混杂变量。如果批次和某个协变量高度相关（比如第1批全是男性，第2批全是女性），那校正时会很难分离真正的批次效应。这种情况我建议先检查一下设计矩阵的秩，必要时放弃彻底校正，或者改用更保守的方法。
你想想看，如果批次标签和协变量几乎一一对应，ComBat会把协变量的差异也当成批次纠正掉。结果就是你再也看不到性别差异了。所以我在每次校正前都会做一张热图，看看批次和协变量之间有没有明显关联。
3.4 缺失值处理
ComBat不允许表达矩阵中有任何缺失值。掉一个NaN进去，整个模型就会报错。所以必须提前处理。
常见的处理方法有三种：
| 方法 | 做法 | 优点 | 缺点 |
|---|---|---|---|
| 删除基因 | 去掉缺失值占比超过阈值的基因 | 简单粗暴，适用于缺失严重的基因 | 可能损失重要的低表达基因 |
| 填充均值/中位数 | 每个基因内用非缺失样本的均值填充 | 快速，不改变总体分布 | 会降低方差，可能引入假均匀性 |
| kNN填充 | 基于相似样本的加权平均填充 | 更精确，保留表达模式 | 计算量大，对小样本不稳定 |
我个人的习惯是：如果某个基因在超过20%的样本里都是缺失，那我直接删掉它。剩下的缺失，我用kNN（k=10）来填充。不过记得填充后要检查一下——我曾经遇到kNN把本来缺失的两个样本填成了完全一样的值，导致后续差异分析假阳性飙升。
最后，把你准备好的表达矩阵、批次向量、协变量（如果有）整理成RData或csv文件，给ComBat喂进去。
下面是整个数据准备流程的示意图，方便你对照检查自己的文件。
到这里，你的数据就已经准备好交给ComBat了。是不是感觉心里踏实多了？嗯，数据准备就像盖房子打地基，地基稳了，后面的校正才能干净利落。