很多做生信的朋友拿到GEO数据，第一反应是赶紧跑差异分析、画火山图。结果跑出来的基因列表怎么看怎么不对劲，换一个差异分析工具结果就对不上，甚至换个样本组合结论就反转。这种问题十有八九不是下游算法的问题，而是预处理阶段埋下的雷。
这个系列的上篇我聊了GEO数据的下载、芯片探针注释、表达矩阵的读取这些基础工作，如果还没看可以翻一下。这篇下篇，重点讲拿到表达矩阵之后、正式做差异分析之前的那段路：样本质量到底靠不靠谱、数据要不要归一化、怎么处理批次效应，以及差异分析跑完怎么验证结果。这些步骤听起来不如画图热闹，但恰恰是决定后面所有结论能不能站住脚的关键节点。
这篇内容适合刚接触GEO数据挖掘、或者跑过流程但总觉得自己对结果没底的同学。我会尽量把每个步骤背后的判断逻辑讲清楚，不只是给代码。
1. 样本质量评估：先摸清这批数据的底细再动手
预处理的第一步不是急着归一化，而是先评估每个样本的质量。我见过不少人在这一步偷懒，结果到了差异分析阶段才发现某个样本其实是异常样本，整个分组比较被它带偏。这个顺序一定不能乱。
1.1 三个必须看的质控指标
芯片数据和测序数据的质控指标不太一样，但有几个共性的检查点：
- 表达值分布：单个样本的表达值分布如果明显偏离其他样本，比如整体表达量偏高或偏低，通常提示这个样本可能有问题。芯片数据可以看箱线图，测序数据可以看表达量密度分布图。
- 样本间相关性：计算所有样本两两之间的表达谱相关性。正常情况下，生物学重复之间的相关性应该明显高于不同组别之间的相关性。如果某个样本跟谁都不像，那它很可能是问题样本。
- PCA聚类：对所有样本做PCA降维，看它们在主成分空间里的分布是否合理。组内样本应当相对聚在一起，如果某个样本孤零零地漂在外面，要特别警惕。
质控这段代码在R里面写起来很快，核心就是用exp表达矩阵做几个基础计算和可视化。
# 样本相关性矩阵 cor_matrix <- cor(exp, method = "spearman") # 用pheatmap或ggplot2画相关性热图 pheatmap::pheatmap(cor_matrix, display_numbers = TRUE, number_format = "%.2f") # PCA计算 pca_res <- prcomp(t(exp), scale. = TRUE) pca_df <- data.frame(pca_res$x[, 1:2], group = group_info)有的朋友会问，prcomp里面的scale.要不要设成TRUE。我个人建议设成TRUE，因为GEO里很多矩阵是不同平台整合出来的，表达量量纲差异很大。用相关性矩阵做PCA的话，cor函数本身就相当于标准化了，但用prcomp时标准化选项还是加上更稳妥。
1.2 判定标准不是死板的，但异常样本必须处理
经常有人问我，相关性低于多少算不合格？PCA里离多远算离群？这个没有绝对标准,要看数据整体情况。
比如某组样本组内相关性普遍在0.85以上，有一个样本跟本组其他样本相关性只有0.6，跟另一组样本相关性却有0.75，那这个样本就很可疑——它可能样本标签标错了，也可能是原始数据质量差。这时候我会先查原始数据的处理记录，如果没法确认，就在下游分析中去掉它，而不是硬留着。
我不建议一上来就按照教科书说的"相关性低于0.8就剔除"。质性判断要看整体分布的。有些异质性很强的组织（比如肿瘤样本）组内相关性本来就不高，强行定一个阈值会把好样本误杀。
质控完了之后，记得把排名靠后、确实异常的样本剔掉，剩下的样本进入下一步。
2. 归一化到底是做什么，以及不同数据类型的正确姿势
归一化大概是整个预处理流程里被误解最多的环节。很多人只知道"要做归一化"，但不知道为什么做、什么时候做、不同数据该用哪套方法。搞清楚这些，你的预处理水平会明显上一个大台阶。
2.1 芯片数据：RMA和限值归一化的选择逻辑
芯片数据的表达值本质上是荧光信号强度经过背景校正之后的数值。不同芯片之间，甚至同一张芯片的不同区域之间，杂交效率和荧光检测效率都有系统性偏差。归一化的目的就是把这些技术偏差抹掉，尽量只保留生物学差异。
Affymetrix的芯片通常推荐做RMA标准化（Robust Multi-array Average），它同时完成背景校正、分位数归一化和探针汇总三个步骤。如果你用的是affy包或oligo包读取CEL文件，rma()一步就能完成。需要注意的是，RMA输出的数据已经是log2转换过的，之后不要再自己做log2。
Illumina和Agilent的芯片，GEO上通常提供的是已经处理过的基础表达矩阵，这类数据我会用limma包里面的normalizeBetweenArrays，方法选"quantile"或者"cyclicloess"。前者速度快、稳定；后者在样本量比较大的时候表现更好，能保留更多表达谱的分布信息。
# 芯片数据限值归一化 library(limma) exp_norm <- normalizeBetweenArrays(as.matrix(exp), method = "quantile")2.2 测序数据：标准化的本质在纠正文库大小
测序数据的归一化和芯片逻辑完全不同。RNA-seq得到的原始count数，受测序深度严重影响。一个样本测了3000万条reads，另一个只测了1500万条，前者几乎所有基因的count数都会翻倍，但这不代表基因表达量真的翻倍了。
所以测序数据的归一化核心是纠正文库大小这个技术偏差。常见做法包括：
- CPM/TPM：除以总reads数再乘以一百万，相当于把每个样本都换算到同一测序深度下。简单直观，适合样本间差异不大时使用。
- TMM：
edgeR包里的calcNormFactors用的就是这个方法，它基于假设大多数基因在两个样本间不差异表达，估算出一个缩放因子。比CPM更加稳健。 - DESeq2的median ratio：同样是基于"大多数基因不差异"的假设，用所有基因的表达量比值的中位数作为缩放因子。
代码上，如果后面准备用limma做差异分析，推荐用limma-voom流程。voom的作用是把count数据变成log2-CPM的形式，同时估算每个基因的方差异质性，给后续线性建模用。
# edgeR + voom 流程 library(edgeR) dge <- DGEList(counts = count_matrix, group = group_info) dge <- calcNormFactors(dge) design <- model.matrix(~ group_info) v <- voom(dge, design, plot = TRUE)2.3 一个常被忽略的问题：官方矩阵里的"已归一化"不等于能直接用
这是GEO数据挖掘里非常隐蔽的一个坑。很多GEO系列里提供的表达矩阵，文件说明里写着"normalized"或者"processed"，但这不代表你可以直接拿去做差异分析。
我遇到过的真实情况有这三种：
- 矩阵本身是log2转换后的比值值，但部分样本来自不同研究批次，批次效应赫然可见。
- 矩阵虽然做过最基础的归一化，但用的方法比较粗糙，其实是用MAS5处理的Affymetrix数据。MAS5和RMA这两种算法的数值分布差异非常大，混在一起做分析会出问题。
- 有的矩阵包含多个探针测得的数据，注释表里甚至有重复的基因名，需要先做探针汇总（collapse），再谈归一化。
所以拿到矩阵之后，先要搞清楚这个矩阵从哪来、经过了哪些处理。GEO里每个样本的Series Matrix File下载下来后，最好先看样本层面的!Sample_processed字段和系列层面的!Series_summary来理解数据背景。不能只相信文件名里的"normalized"字样。
3. 批次效应：什么时候必须处理，什么时候不该硬拆
批次效应是GEO数据挖掘绕不开的坎。GEO里的数据经常是不同实验室、不同时间、不同芯片批次生成的，这些技术差异会系统性地影响表达值。如果你想合并多个数据集做Meta分析，批次效应不处理，第一批PCA图就可以让你怀疑人生。
3.1 先判断你的场景是否需要批次矫正
不是所有场景都适合做批次矫正。这一点经常被忽略，我把它放在最前面说。
- 单个GEO系列，内部是同一个实验设计的：比如一个芯片平台、一次性杂交完成的。这种通常不需要做过度矫正，芯片本身的归一化就够了。如果强行做ComBat，有时候会把真实的生物学差异也拆掉一部分。
- 多个GEO系列合并：这是批次矫正的主战场。不同系列的样本在测序或芯片实验中本身就有平台、实验室、批次差异，必须处理。
- 需要验证模型：有些研究为了增加样本量，把多个数据集合并后重新分组做比较，这时批次矫正几乎是强制步骤。
有一个场景我不建议用ComBat：生物学差异恰好和批次完全混杂的时候。比如你的处理组全部来自某个数据集，对照组全部来自另一个数据集，这时候处理组和对照组之间的差异既包含生物学差异也包含批次差异，ComBat去不掉混杂的批次，反而可能把真正的生物学差异一并抹掉。遇到这种设计，正确的做法是只做样本级别的质控合并，再分析之前明确声明这个局限，或者通过数据集内部差异方向来做验证，而不是天真地用ComBat收拾残局。
3.2 ComBat与ComBat-seq的选择与实操
经典ComBat是面向芯片和连续型表达值的，假设表达值服从正态分布。RNA-seq的count数据分布是负二项或泊松型的，直接套经典ComBat会有偏差。这时候有两个选择：
sva包里的ComBat，适合处理芯片、log2-CPM或其他连续值。sva包里的ComBat_seq，针对RNA-seq的count数据，用负二项模型拟合，是目前处理测序批次效应的主流方法。
操作上，需要先构建batch向量和group向量。这里有个关键原则：矫正时务必把生物学分组信息传给group参数。
library(sva) # 连续型表达值（芯片/log2-CPM） exp_combat <- ComBat(dat = exp_norm, batch = batch_info, mod = model.matrix(~ group_info)) # RNA-seq count数据 dge_combat <- ComBat_seq(counts = count_matrix, batch = batch_info, group = group_info)为什么要传group？因为不传的话，ComBat会假设所有样本来自同一个分布，它在矫正批次时会把组间的真实差异也当成批次变异的一部分，导致矫正过度。传了分组，模型就会在有生物学分组信息的前提下，只拆掉"组内批次间"的偏差。
3.3 怎么验证批次矫正是否有效
有的人跑完ComBat心里没底，不知道矫正到什么程度才算到位。我的验证方法有三层：
- 矫正前的对照：先画矫正前的PCA，确认确实存在明显的批次分离。如果矫正前批次已经混在一起，那根本不必要矫正。
- 矫正后的PCA：看批次信息是否不再主导样本分布。同时看组间分离是否保留、甚至更清晰。
- 定量验证：计算矫正前后样本按批次分组时的组间差异（比如PERMANOVA的p值），确认批次解释的方差比例显著下降。
另外一个附加检查：矫正前后看目标基因的表达模式有没有异常反转。如果你关心的基因在原有数据里是有明确生物学逻辑的，矫正后这个逻辑还在，那批次矫正没有破坏核心信号。
还有一点要提醒：ComBat矫正后的表达矩阵，不再是原始的count整数了，它会输出浮点数。所以如果你后面要用DESeq2，需要谨慎。通常ComBat-seq矫正后的数据可以配合limma-voom使用，但如果你强制要用DESeq2的negative binomial模型做差异分析，用矫正后的非整数数值是不合适的。这点操作前一定想清楚。
4. 探针汇总与基因注释：预处理最琐碎也最容易翻车的环节
GEO里很多表达矩阵的每一行对应的是一个探针（probe），不是基因。做差异表达分析之前，必须把探针水平的数据映射到基因水平。这个环节很容易出错，而且错了很难发现，因为下游所有分析都会跟着错。
4.1 一个探针对应多个基因，一个基因对应多个探针
芯片平台都有这个特点。Affymetrix的3'端表达芯片，一个基因经常对应多个探针；而某些平台的一个探针可能同时比对到多个基因，这种多对多的关系网上填表的时候就很容易让你一脸懵。
处理原则是这样的：
- 一个基因多个探针：取平均值或中位数，或者选取方差最大的那个探针（系数变异最大、信息量最高）。我个人习惯选表达值方差最大的那个探针，因为信息量最高；如果差异分析用limma，也可以先不汇总，直接用
limma::avereps()按基因名折叠，原理一样但更省事。 - 一个探针对应多个基因：这个比较麻烦。优先查一下这个探针在平台注释里是否标了"cross-hybridization"（交叉杂交），如果标了，直接剔除，这类探针信号不可靠。剔除后如果只对应一个基因，那就归并；如果对应多个不同基因，一般也剔掉。
# 一个基因多个探针，取最大方差探针 probe_var <- apply(exp_matrix, 1, var, na.rm = TRUE) exp_matrix <- exp_matrix[order(probe_var, decreasing = TRUE), ] # 用重复基因名保留第一个（即方差最大的） exp_collapsed <- exp_matrix[!duplicated(gene_symbols), ]这里有个细节：order(probe_var, decreasing = TRUE)之后duplicated()保留的是第一次出现的行，也就是方差最大的那个探针，符合预期。
4.2 基因ID转换：SYMBOL、ENTREZID还是ENSEMBL
做富集分析时经常需要各种ID。大部分在线工具（比如某个富集分析网站或用R包做GSEA）都能自动帮你转换，但离线分析或者批量操作的时候，还是建议自己先把ID统一。
常用的转换工具是clusterProfiler包和org.Hs.eg.db注释包。比如人的数据：
library(clusterProfiler) library(org.Hs.eg.db) # 从SYMBOL转ENTREZID keytypes(org.Hs.eg.db) # 查看支持的ID类型 gene_entrez <- bitr(gene_symbols, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)转换过程中会有一些基因名匹配不到，这些基因后续做富集分析时会丢失。我通常建议保留原始SYMBOL和转换后的ENTREZID两列，后面随时可以反查。
另外一个常见坑：平台注释文件里面的基因名经常带引号，或者包含---这种填充字符。读进来之后先做一轮字符串清洗，不要直接拿去做匹配。
# 清洗注释列 gene_symbols <- gsub('"', '', gene_symbols) gene_symbols <- trimws(gene_symbols) gene_symbols[gene_symbols == "---"] <- NA4.3 我踩过的注释文件大坑
有一回我处理某平台的表达矩阵，注释文件下载下来是Excel格式。用readxl读进来之后发现基因名有几千个NA。一开始以为是文件损坏，后来仔细翻了原始平台说明，发现这个平台的注释表里基因名是放在第二列而非第一列的，而且分隔符有特殊情况。读进来之后必须按照平台说明里的列名认真核对，不能想当然认为第一列就是基因名。
这里也顺带说一句：很多朋友懒得下平台注释，直接用GEOquery包的getGEO函数里带的feature data。这个可以做，但要注意feature data里的注释信息有时不完整，特别是某些旧平台的注释是后来从外部映射上去的，缺失率可能很高。稳妥做法是去平台官网下载原始注释文件，和数据里的探针ID做一次全面匹配，这个注释表才是最终依据。
5. 差异表达分析实操：从limma流程到结果验证
预处理做完了，终于到了真正出结果的部分。但差异表达分析本身也有不少选择要做对，否则前面辛辛苦苦的预处理就白费了。
5.1 一次完整流程：探针汇总后的limma-voom
我用的最多的是limma流程，从探针汇总后的表达矩阵出发，一次性打通所有环节。下面是一个完整的参考代码框架：
library(limma) library(edgeR) # 1. 构建DGEList对象 dge <- DGEList(counts = count_matrix, group = group_info) # 2. 过滤低表达基因 keep <- filterByExpr(dge, group = group_info) dge <- dge[keep, , keep.lib.sizes = FALSE] # 3. 计算归一化因子 dge <- calcNormFactors(dge) # 4. 设计矩阵（这里是二分组比较，多分组也类似） design <- model.matrix(~ 0 + group_info) colnames(design) <- levels(factor(group_info)) # 5. voom转换 v <- voom(dge, design, plot = TRUE) # 6. 线性拟合 + 对比 fit <- lmFit(v, design) contrast_matrix <- makeContrasts( GroupA_vs_GroupB = GroupA - GroupB, levels = design ) fit2 <- contrasts.fit(fit, contrast_matrix) fit2 <- eBayes(fit2) # 7. 提取结果 results <- topTable(fit2, coef = "GroupA_vs_GroupB", number = Inf, sort.by = "none")这段代码几乎每个步骤都有讲究。filterByExpr是过滤低表达基因，因为测序数据里大量基因的count数可能全是0或者只有个位数，这些基因不仅没有统计功效，还会干扰多重检验校正。calcNormFactors是TMM归一化。voom返回的表达值已经是对数化的，直接用lmFit做线性建模。
5.2 多分组和配对设计怎么处理
GEO数据经常不止两组。比如某个疾病数据集里面有正常、早期、晚期三个分组，要做正常对晚期的比较。这种情况下，设计矩阵用~ 0 + group而不是~ group，是为了方便构建任意对比。makeContrasts里可以同时定义多个对比，比如正常vs早期、正常vs晚期。
如果数据是配对设计（同一个病人治疗前和治疗后），设计矩阵要包含个体信息：
design <- model.matrix(~ patient + treatment)这种设计矩阵会把个体间的差异量吸收掉，检验的只是治疗前后的效应，对配对数据来说是更准确的做法。如果没有包含patient这一项，就会被当成两组间完全随机比较，P值失去最合适的解释框架。
5.3 logFC和P值筛选阈值，不该只会看padj < 0.05
跑完topTable之后，常见的筛选条件是padj < 0.05和|logFC| > 1（即表达量变化大于两倍）。但实际处理中，这两个阈值的合理性要结合数据本身来判断。
- 如果是肿瘤组织对比正常组织，差异可能非常大（很多基因表达变化几十倍、上百倍），
logFC > 1偏保守，往往只保留噪音最小的少数基因。 - 如果是一个干预实验、药物处理对细胞系的影响，差异幅度可能没那么夸张，
logFC > 1可能把大量有生物学意义的基因筛掉了。
我个人习惯先看所有显著基因的数量，再结合预期生物学背景调整阈值。一般先用padj < 0.05筛选一遍，看数量是不是在合理范围（通常几百到几千）。如果明显太多或太少，再考虑放宽padj到0.1或者调整logFC阈值。不要死守默认值。
多重检验校正方法上，topTable默认做的是BH（Benjamini-Hochberg）校正，这也是绝大多数生信文章用的方法。如果你用的比较组数量非常多，建议考虑bonferroni，控制得更严格。
5.4 验证结果的三个层面
拿到显著基因列表后，我不会立刻忙着画热图，而是先做三层验证，确保结果不是随机噪音：
- 表达分布检查：随机挑几个显著基因，画出组间的表达分布（箱线图或小提琴图）。如果显著基因的分布几乎完全重叠，那这个显著性很可能是归一化假象或模型误设。这个检查能筛掉很多假阳性。
- 外部数据验证：看看这个基因在其他GEO数据集里是否也有类似差异方向。比如用同一个疾病的不同数据集做一个简单的meta方向验证。不需要多严格，方向一致就足以增加信心。
- 生物学合理性：把显著基因列表先做个简单的通路富集，看看富集到的通路是不是和这个疾病/组织的已知生物学机制吻合。如果富集出一堆风马牛不相及的通路，就该回头检查前面的注释和预处理了。
6. 预处理链条上最容易翻车的三个位置，以及我的处理习惯
最后把这轮内容里我认为最关键的提醒集中说一下。这个部分没有代码，但比代码更值钱，因为它来自实实在在踩过的坑。
6.1 样本顺序错位：表达矩阵列顺序和分组信息对不上
这是我见过最多、也最容易导致结论完全反转的问题。下载GEO矩阵后，表达矩阵里的样本顺序是平台按照自身顺序排的，而从GEOquery里提取的分组信息可能又是另一个顺序。如果直接按位置对应，A组样本全变成B组，差异结果会错到离谱。
我的处理习惯是：每次加载数据之后，第一步先核对样本ID一致性，用identical()函数或者all()做一次严格匹配。
stopifnot(all(rownames(pheno_data) == colnames(expr_matrix)))不通过就直接报错，宁可程序停下来也别带病运行。这句话看起来简单，但能救很多人一命。
6.2 多个探针对应同一个基因时，不同汇总方式结果差异很大
有人习惯"取均值"，有人习惯"取最大值"。这两种方式对同一个基因得到的表达量完全不同，进入下游差异分析之后，这个基因是否能达到显著阈值都可能不一样。
我自己的标准是：优先用avereps()（等价于取均值，但由limma内部处理更规范），如果基因的多个探针之间存在较大的表达差异，则改用最大方差探针。关键是要在论文方法部分写清楚，因为这会影响可重复性。报告里"平均值"和"最大方差探针"两行字，审稿人可能不在意，但对复现你结果的人来说影响巨大。
6.3 归一化的位置放错，整个流程白做
再强调一遍顺序问题，这是新手最容易搞混的：
芯片数据：读取 → 质控 → 归一化 → 探针汇总 → 差异分析。
测序数据：读取 → 质控 → 过滤低表达 → TMM归一化/voom → 探针汇总（如有） → 差异分析。
有人习惯先做探针汇总再做归一化，导致同一个基因的多个探针在归一化过程中被当作独立变量处理，归一化因子会被少量高表达探针主导，结果可能失真。
还有一点要留意：不要对GEO官方提供的连续型矩阵再做一次normalizeBetweenArrays。很多官方矩阵本身就做过归一化，重复做不仅没必要，还可能把原本正常的分布弄歪。判断依据是看矩阵里数值的分布。如果是明显log2化的连续值（范围可能在1到15之间），通常已经标准化过了，直接用即可。
写到这里，GEO数据预处理的主体流程就串起来了。这套流程不是什么地方都需要严格照做，但每一步为什么存在、什么时候可以跳过，我希望上面的解释已经让你心里有数。预处理没有激动人心的产出，它的价值全在于让后续的差异分析、可视化、功能注释建立在可靠的地基上。
我自己每次拿到新的GEO数据集，都会先重跑一遍质控和归一化的可视化了再决定往后走，哪怕这个数据已经看过很多遍。这个习惯帮我筛掉过不止一次因样本标签混淆带来的问题数据。也建议你把这些检查步骤固化到脚本里，每次换数据只改路径和分组文件，就能跑完一套完整质控报告，这会让你的数据挖掘过程稳定很多。
