企业官网建设流程全解析

上篇把探针注释和表达矩阵的读取讲完了，不少朋友私信说，最头疼的不是“怎么把数据下载下来”，而是“数据拿到手之后，后面这一串到底怎么处理才敢放心交给下游分析”。做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数据里不同平台、不同数据类型适用的方法差异很大，我把几种常见的列成一张表：

对于微阵列芯片，我更推荐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$E

v$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数据时经常遇到的问题速查表，每一个都是真实踩过的坑，希望能帮大家少查几次搜索引擎。

6.1 踩坑实录：某次分组标签错位导致的结果反转

有一次我处理某个疾病数据集，差异分析跑完发现上调基因全是免疫相关通路，本来还挺高兴，结果在手动验证阶段发现，有两个tumor样本其实被错误标成了normal。原因是样本集里存在一个“normal_tumor”之类的中间命名，被grepl错误拾取。差一点就整篇分析推倒重来，还好当时做了topTable后的手动抽查。

所以我现在凡是遇到标题字段里包含否定前缀或者复合词的，一律先手动打印所有唯一标题再确认规则。分组向量生成之后，用table(group, title)交叉验证，真的能救一命。

6.2 踩坑实录：ComBat过度校正

另一次处理跨平台合并数据时，我直接对全部基因做了ComBat批次校正。校正完一看PCA，批次确实消失了，但连正常组和疾病组也完全重叠了。后来查文献才发现，因为两个平台的样本恰好分属两个分组，批次和分组完全混杂，ComBat把组间真实差异当作批次效应消掉了。后来我改用把批次作为协变量放进limma设计矩阵的办法，问题就解决了。

这个教训给我的启示是：批次效应的处理绝不是“跑一个ComBat”那么简单，必须想清楚你自己的数据里批次和分组的关系，再决定用哪种策略。

6.3 给初学者的建议：建立自己的预处理checklist

如果你刚开始接触GEO数据，我建议从第一个数据集开始，就养成跟着checklist走的习惯。我自己在项目里常用的流程是：

1. 检查表达矩阵维度、缺失值比例、列名和行名格式
2. 做样本相关性热图和PCA，识别离群样本
3. 根据平台选择归一化方法
4. 判断是否存在批次效应，若存在则明确与分组的关系
5. 提取分组信息并用match确保与表达矩阵列名对齐
6. 构建设计矩阵并检查因子水平顺序
7. 保存所有中间文件和RData
8. 用limma/edgeR/DESeq2跑差异分析，先抽查50个基因再做可视化

这个流程不复杂，但能保证每一步都有据可查。尤其是保存中间文件，很多人会忽略，等发现结果异常时就没有对照可查，只能从头再跑一遍，时间就这么浪费了。

做GEO数据预处理，真正拉开差距的往往不是那些炫酷的可视化技巧，而是每一步逻辑是否严密、顺序是否合理。数据清洗没有捷径，该查的缺失值、该做的归一化、该处理的批次效应，每一项都值得用文档记录下来。希望这篇下篇能把你的GEO分析流程补完，后续做差异分析、富集分析的时候，能少一点“数据重跑”的无奈，多一点“结果可解释”的从容。