GEO数据预处理全流程：表达矩阵清洗、归一化与批次效应处理

上篇把探针注释和表达矩阵的读取讲完了不少朋友私信说最头疼的不是“怎么把数据下载下来”而是“数据拿到手之后后面这一串到底怎么处理才敢放心交给下游分析”。做GEO数据挖掘这几年我的体会是数据预处理的下半程才是真正的分水岭前面读进来的是毛坯后面洗出来的才是能进模型的成品。今天这篇就专门补上下篇的内容把表达矩阵清洗、归一化、批次效应处理、分组信息整理一直写到差异分析前的最后检查一条线走完。里面的代码我都用确定性项目跑过不是纸面教程是能直接拿过去用的那种。如果你是刚开始做GEO芯片或者转录组数据分析的这篇文章看完应该能少走很多弯路。1. 表达矩阵清洗先给自己定一套基线标准拿到上一步处理好的表达矩阵之后第一步不是急着画图而是先做样本和探针层面的清洗。这个环节看起来不起眼却是后面所有统计检验的地基。地基没打牢后面画出来的火山图再漂亮审稿人一句“你的样本量太少”或者“你的离群样本没有处理”就能把你所有工作打回原形。1.1 三个必须看的质控指标我在检查表达矩阵时从来都是按“样本-探针-整体分布”三个维度依次过一遍。样本维度的指标是平均表达量、内参基因表达、样本相关性。芯片数据一般会有阳性对照和阴性对照探针如果某个样本的阳性对照信号整体偏低基本可以判断这个样本的杂交质量不过关。如果是转录组数据则更看重比对率、基因检出数这些指标不过GEO里归档的转录组数据不一定都会给出完整质控报告所以我会优先看表达矩阵本身能不能通过后续清洗。探针维度的指标是缺失值比例和表达量方差。一个探针如果在超过70%的样本里都是缺失值它还留着的意义就不大。反过来一个探针在所有样本里表达量几乎不变方差接近0这种探针对后续差异分析也没什么贡献通常会在过滤那一步被筛掉。整体分布维度则是看表达值是否呈典型正态或长尾分布以及有没有离谱的低值或极值。你可以用summary函数快速扫一眼# 假设已经读入表达矩阵 expr行为探针列为样本 summary(expr)如果某个样本的中位数或者最大值明显偏离其他样本就要回到样本相关性那边去做进一步判断。1.2 缺失值处理能插补就插补别动不动就删样本GEO数据里出现缺失值太正常了尤其是芯片数据同一批次里总有那么几个探针在个别样本上亮不出来。处理缺失值有两个流派一个是删一个是补。我自己的原则是缺失比例高的探针直接删缺失比例低的探针用KNN插补。具体阈值没有硬性规定但我在项目里习惯用50%作为分界线。某探针在超过一半样本里都是缺失值说明这个探针的检测本身就不稳定插补出来的数字也没多大意义。KNN插补的原理其实很简单对某个探针的缺失值找到表达谱最相似的K个探针用它们的表达值加权平均来填充。实际操作里我会用R的impute包默认K取10就行跑起来也很快。library(impute) # 先去掉缺失比例过高的探针 miss_rate - apply(expr, 1, function(x) sum(is.na(x)) / length(x)) expr_filt - expr[miss_rate 0.5, ] # KNN插补剩余缺失值 expr_imp - impute.knn(as.matrix(expr_filt), k 10)$data这里有一个容易踩的坑KNN插补之前最好先把表达矩阵转成常规数值矩阵别带着数据框的rowname属性一起干。另外插补之后记得检查是否还有残留的NA如果还有很可能是K值太小或者数据本身太稀疏可以适当调整参数。1.3 离群样本识别PCA和相关性热图双管齐下清洗完缺失值之后还需要确认样本整体有没有“异类”。最直观的两个工具是主成分分析PCA和样本间相关性热图。我一般会先画相关性热图因为计算成本低且能直接看出哪个样本跟谁都不像。library(pheatmap) cor_mat - cor(expr_imp, method spearman) pheatmap(cor_mat, show_colnames TRUE, show_rownames TRUE)如果某个样本与其他样本的相关系数普遍低于0.8就要高度警惕。它可能是批次效应的产物也可能是杂交失败或者分组标签搞错了。再用PCA做一次确认pca_res - prcomp(t(expr_imp), scale. TRUE) plot(pca_res$x[, 1], pca_res$x[, 2], pch 19)PCA图里孤零零飘在角落的点大概率就是离群样本。不过不要一看到离群就删先去看一下这个样本的来源注释。如果只是技术原因导致的失败样本删掉是合理的如果是某个临床亚组的独特样本贸然删掉反而会失去重要信息。提示离群样本的处理一定要在分析记录里写清楚“为何保留/为何剔除”别怕麻烦这一步能让你在写论文方法部分时省很多事。2. 归一化把不同样本拉到同一把尺子上清洗完的矩阵只是“干净”了还不意味着“可比”。不同样本之间哪怕来自同一批次芯片扫描时的荧光强度也会受到很多非生物因素的影响比如杂交条件、扫描参数、试剂批次差异。如果不做归一化某些样本的绝对表达值偏高纯粹是技术原因后续所有比较都会失真。2.1 不归一化直接跑差异分析的后果我见过不少初学者直接拿原始表达量去跑差异分析然后发现某个基因上调了3倍兴奋得不行后来一查发现是那个样本本身整体信号就高基因并没有真的上调。这就是典型的没有归一化导致的假阳性。简单说归一化的目的是把所有样本的表达值分布调整到同一个基准上让样本之间残余的差异尽可能反映生物学差异而不是技术噪声。2.2 常见归一化方法该怎么选GEO数据里不同平台、不同数据类型适用的方法差异很大我把几种常见的列成一张表方法适用数据类型思路常用工具RMAAffymetrix芯片基于多芯片模型校正背景噪声和探针亲和度affy / oligoMAS5Affymetrix芯片单芯片归一化适合看绝对表达量affyQuantile各种芯片/表达矩阵让所有样本的表达值分位数分布一致preprocessCoreTMMRNA-seq基于经验分布计算缩放因子适合组间文库大小不一致edgeRCPM/TPMRNA-seq纯数值缩放适合样本间文库大小校正edgeR/genefilterVoomRNA-seq转微阵列思路先算cpm再估计均值-方差关系limma对于微阵列芯片我更推荐RMA或者Quantile原因是它们在校正样本间技术差异方面比较成熟。对于RNA-seq数据TMM配合后续的limma-voom流程是R社区里用得最广的组合。2.3 实战用limma的voom完成归一化很多人对limma的印象停留在差异分析其实它的voom函数就是专门为RNA-seq数据设计的一步到位的处理方案。它会先把原始计数转换为log2-CPM再根据均值-方差关系给每个基因一个权重。library(limma) library(edgeR) # 假设你有原始的count矩阵 count_mat行为基因列为样本 dge - DGEList(counts count_mat) dge - calcNormFactors(dge, method TMM) # 设计矩阵可以先不填voom必须在有样本信息的前提下运行 # 这里先用一个简单的分组向量来做演示 group - factor(c(ctrl, ctrl, treat, treat)) design - model.matrix(~ group) v - voom(dge, design, plot TRUE) expr_norm - v$Ev$E就是归一化后的log2-CPM表达矩阵。这个矩阵可以直接用于后续的PCA、热图、差异分析。如果你的GEO数据已经是过程值矩阵比如下载的是Series Matrix File并且数据本身来自芯片平台那就可以直接对表达矩阵做分位数归一化library(preprocessCore) expr_norm - normalize.quantiles(as.matrix(expr_imp)) colnames(expr_norm) - colnames(expr_imp) rownames(expr_norm) - rownames(expr_imp)归一化之后最好画一张箱线图来对比归一化前后的分布差异。归一化做得好的话不同样本的箱线图会高度一致中位数几乎落在同一条水平线上。boxplot(expr_imp[, 1:6], main Before normalization) boxplot(expr_norm[, 1:6], main After normalization)3. 批次效应真实世界绕不开的一关GEO数据有个很现实的问题一个数据集里的样本很可能不是同一次实验完成的。它们可能是某课题组在不同时间段跑的芯片也可能是从平台上下载后合并了多个子课题的数据。这种来源差异造成的系统性偏差就是批次效应。3.1 批次效应到底长什么样批次效应最常见的表现是不同批次的样本在PCA图上各自聚成一团哪怕这些样本的生物学分组是一样的。如果不处理你后续差异分析检出来的基因很可能有一部分不过是在反映批次的差异而不是疾病和正常样本的差异。有个经典场景某数据集总共30个样本其中15个正常、15个患病但15个患病样本全在批次2里。这个时候性别差异、年龄差异都可能混在批次差异里除非实验设计时做了平衡否则很难把真实的生物学效应和批次效应干净地剥离开。所以第一步永远是先做PCA并依据批次信息着色看看数据里到底存不存在明显的批次结构。3.2 什么时候校正什么时候不要校正我在项目里总结出一条经验看批次和分组的关系。如果批次和分组完全独立比如每个批次里都同时包含正常和患病样本那么校正批次效应是安全的。如果批次和分组几乎完全混杂比如批次1全是正常批次2全是患病那么用ComBat去校正就非常危险因为算法会把组间差异的一部分也当成批次效应消掉最后可能什么显著基因都跑不出来。对于混杂的情况比较稳妥的做法是不做全量批次校正而是把批次作为协变量放进差异分析模型里。3.3 用ComBat做批次校正的完整流程当批次效应真实存在且与分组独立时我优先使用sva包的ComBat函数。ComBat基于经验贝叶斯方法能够估计每个批次内的位置和尺度效应再把所有批次的数据调整到同一个参照水平。library(sva) # batch向量与样本一一对应比如 c(1,1,1,2,2,2) batch - c(1, 1, 1, 2, 2, 2) mod - model.matrix(~ group) expr_combat - ComBat(dat expr_norm, batch batch, mod mod)跑完之后建议再看一次PCA如果批次聚团现象明显减弱同时分组还能保持一定的分离那就说明校正效果不错。这里我要特别提醒的是ComBat之后不要再去跑一次分位数归一化否则等于重复校正容易把数据的真实方差压得太低。4. 分组信息整理最不起眼但翻车率最高的一步如果前面所有步骤都顺利最后却因为样本分组搞错导致结果全反那真是欲哭无泪。分组信息整理在GEO数据预处理里看起来最简单实际却是翻车率最高的一步——因为GEO的样本标题字段往往写得非常自由有人写“normal_1”有人写“Sample_N1”还有的干脆把所有信息塞在特征列里。4.1 从标题字段里提取分组信息假设我们从GEO下载了表达矩阵和对应的样本信息样本标题通常存在pData()对象里。常见的提取方式是根据字符串匹配来生成分组向量。library(GEOquery) gset - getGEO(GSE12345, GSEMatrix TRUE) pdata - pData(gset[[1]]) # 假设标题列里包含 normal 或 tumor 关键词 title - as.character(pdata$title) group - ifelse(grepl(tumor, title, ignore.case TRUE), tumor, normal)但这里有个隐患如果标题写的是“non_tumor”那么grepl就会把这个样本错误地归到tumor组。所以我每次写完都先table一下table(group, title)确认每一个标题和分组的对应关系符合预期。这一步多花30秒能省掉后面两小时返工。4.2 样本对齐顺序错了一切都白搭提取完分组之后最关键的操作是确认分组向量的顺序和表达矩阵的列顺序完全一致。这两个顺序如果没有对齐后面的差异分析就是把张三的标签贴到李四的样本上。最稳妥的方法是直接用样本名匹配而不是“第几个样本对应第几列”这种肉眼对齐。# 确认表达矩阵的列名与样本信息表完全一致 sample_order - colnames(expr_combat) pdata_ordered - pdata[match(sample_order, rownames(pdata)), ] # 再次确保分组和样本一一对应 stopifnot(all(rownames(pdata_ordered) sample_order)) group - factor(pdata_ordered$group, levels c(normal, tumor))stopifnot这行代码会在不满足条件时直接报错是防止后续分析在错误前提下继续运行的一道闸门。4.3 构建设计矩阵差异分析模型的入口分组信息对齐之后我们就可以构建设计矩阵了。设计矩阵是线性模型的基础告诉R你的表达数据要按哪些因素来解释。对于最简单的两种分组情况design - model.matrix(~ 0 group) colnames(design) - levels(group)这里用0 是为了去掉截距项得到两组各自的均值后续用contrasts来做两组比较很方便。如果你要包含年龄、性别、批次等协变量可以扩展公式design - model.matrix(~ batch age group)注意设计矩阵里分组变量必须放在最后而且不要用~ group这种带截距的形式做后续差异对比否则你得到的是“tumor相对于normal的倍数变化”如果分组水平顺序弄反了结果的正负号也会跟着反。5. 差异分析前的最后检查完成了上面所有步骤之后先别急着跑limma或者DESeq2多花几分钟做最后的检查。我一直认为差异分析的输出结果里最容易犯的错误不是代码写错而是数据处理流程留下了隐性问题等到结果出来了才发现解释不通。5.1 用箱线图、PCA和聚类树做最终确认一组理想的数据在预处理完成后应该满足三个条件样本间表达分布均匀、PCA上不同组别有可见分离趋势、聚类树里同一组的样本倾向于聚在一起。我会把这三张图都做出来存成一个PDF然后逐个检查pdf(final_check.pdf, width 10, height 8) # 箱线图 boxplot(expr_combat, col rainbow(ncol(expr_combat)), main Sample Boxplot) # PCA pca_res - prcomp(t(expr_combat), scale. TRUE) plot(pca_res$x[, 1], pca_res$x[, 2], col c(blue, red)[group], pch 19) legend(topright, legend levels(group), col c(blue, red), pch 19) # 聚类树 hc - hclust(dist(t(expr_combat)), method complete) plot(hc, labels sample_order) dev.off()如果聚类树里出现了明显跨分组的杂糅比如一个tumor样本和所有normal样本聚在一起先别急着怀疑数据回到分组表里查一下这个样本的临床注释很可能是样本标签写错了。5.2 选对差异分析工具limma、edgeR、DESeq2工具的选择应该基于数据类型和你前面的处理方式。如果是芯片数据或者经过voom转换的表达矩阵直接用limma。limma的优势在于稳健的方差估计和成熟的线性模型框架是目前微阵列和转录组数据差异分析的首选。如果是原始count数据但没有做voom转换可以选择edgeR的exactTest或者DESeq2的标准流程。DESeq2在处理小样本和低表达基因时表现更保守适合对假阳性率比较敏感的分析场景。这里有一个关键点不要混用流程。如果你已经用voom生成了log2-CPM矩阵就别再拿原始count去跑DESeq2否则等于把同样数据用两套逻辑分析结果反而互相打架。我用limma做两组比较时通常这样写fit - lmFit(expr_combat, design) contrast_matrix - makeContrasts(tumor - normal, levels design) fit2 - contrasts.fit(fit, contrast_matrix) fit3 - eBayes(fit2) results - topTable(fit3, coef 1, number Inf, adjust.method BH) head(results)coef 1对应第一个对比也就是tumor减去normal的log2倍数变化。正值表示在tumor组上调负值表示下调这个方向问题我在初学阶段也栽过跟头。5.3 保存工作区与结果导出完整的预处理流程跑完之后一定要把关键对象保存下来方便后面重新出图和复现。saveRDS(expr_combat, expr_combat.rds) saveRDS(group, group_factor.rds) saveRDS(design, design_matrix.rds) saveRDS(results, diff_results.rds)写结果文件时我会同时输出基因symbol、log2FC、P.Value和adj.P.Val方便后面做富集分析或者画火山图。write.csv(results, diff_analysis_results.csv)导出之后别急着关项目先随机抽几个基因在原矩阵里手动验证一下表达方向确认和table里展示的倍数变化方向一致。这一步听起来很原始但确实是抓错的最佳方式。6. 常见问题与排查技巧实录最后整理一个我在处理GEO数据时经常遇到的问题速查表每一个都是真实踩过的坑希望能帮大家少查几次搜索引擎。常见问题可能原因排查方向样本相关性热图出现明显暗块批次效应或样本标签错误检查批次变量按批次着色画图某个样本在PCA里孤立杂交失败、样品污染、分组异常查看原始质控和临床注释ComBat之后差异基因全部消失批次与分组混杂导致过度校正不要使用ComBat改为在设计矩阵中纳入批次协变量分组标签和表达矩阵列名不一致用位置而非样本名匹配用match函数按样本ID对齐火山图里全是正log2FC设计矩阵对比方向反了检查makeContrasts里两个组别的顺序表达值出现负数部分平台eg后表达值本身可为负确认平台注释负数可能是log2比例值并非异常重复样本合并未处理同一患者多时间点多部位取样确认是否要保留建议用平均表达值合并或建模处理6.1 踩坑实录某次分组标签错位导致的结果反转有一次我处理某个疾病数据集差异分析跑完发现上调基因全是免疫相关通路本来还挺高兴结果在手动验证阶段发现有两个tumor样本其实被错误标成了normal。原因是样本集里存在一个“normal_tumor”之类的中间命名被grepl错误拾取。差一点就整篇分析推倒重来还好当时做了topTable后的手动抽查。所以我现在凡是遇到标题字段里包含否定前缀或者复合词的一律先手动打印所有唯一标题再确认规则。分组向量生成之后用table(group, title)交叉验证真的能救一命。6.2 踩坑实录ComBat过度校正另一次处理跨平台合并数据时我直接对全部基因做了ComBat批次校正。校正完一看PCA批次确实消失了但连正常组和疾病组也完全重叠了。后来查文献才发现因为两个平台的样本恰好分属两个分组批次和分组完全混杂ComBat把组间真实差异当作批次效应消掉了。后来我改用把批次作为协变量放进limma设计矩阵的办法问题就解决了。这个教训给我的启示是批次效应的处理绝不是“跑一个ComBat”那么简单必须想清楚你自己的数据里批次和分组的关系再决定用哪种策略。6.3 给初学者的建议建立自己的预处理checklist如果你刚开始接触GEO数据我建议从第一个数据集开始就养成跟着checklist走的习惯。我自己在项目里常用的流程是检查表达矩阵维度、缺失值比例、列名和行名格式做样本相关性热图和PCA识别离群样本根据平台选择归一化方法判断是否存在批次效应若存在则明确与分组的关系提取分组信息并用match确保与表达矩阵列名对齐构建设计矩阵并检查因子水平顺序保存所有中间文件和RData用limma/edgeR/DESeq2跑差异分析先抽查50个基因再做可视化这个流程不复杂但能保证每一步都有据可查。尤其是保存中间文件很多人会忽略等发现结果异常时就没有对照可查只能从头再跑一遍时间就这么浪费了。做GEO数据预处理真正拉开差距的往往不是那些炫酷的可视化技巧而是每一步逻辑是否严密、顺序是否合理。数据清洗没有捷径该查的缺失值、该做的归一化、该处理的批次效应每一项都值得用文档记录下来。希望这篇下篇能把你的GEO分析流程补完后续做差异分析、富集分析的时候能少一点“数据重跑”的无奈多一点“结果可解释”的从容。