GEO数据挖掘实战：探针转换、批次校正与差异分析

这篇文章最有意思的点是，作者真的把GEO的每个环节都跑了一遍，不是那种纸上谈兵的教程。我自己试过之后发现，最坑的其实不是下载数据，而是探针到基因的转换这一步——平台注释文件经常过时，一不小心就会丢大半基因。另外，那个ComBat批次校正虽然好用，但你要是样本量太小，它反而会引入假信号。这些细节文章里都提到了，挺实在的。

---

第一次用GEO，我差点被数据搞崩溃

说实话，我第一次打开NCBI的GEO主页时，整个人是懵的。满屏的GPL、GSM、GSE编号，跟密码似的。而且你搜一个关键词，出来几千个结果，根本不知道哪个数据集靠谱。后来踩了几个坑，才摸清楚门道。

GEO是个好东西，2000年就上线了，存了超过100亿个基因表达测量值，覆盖100多个物种，1500多个实验室往里灌数据。但问题也在这儿——数据太杂，质量参差不齐。有的数据集只有两个样本，有的批次效应明显到一眼就能看出来。

三层数据结构，其实没那么复杂

GEO的数据用三层模型组织，每个编号都有特定含义：

平台（GPL编号）

平台就是实验用的芯片或测序仪。比如GPL96对应的是Affymetrix的HG-U133A芯片。每个平台都有自己的探针注释文件，这个文件是关键，后面会提到它有多坑。

样本（GSM编号）

每个样本一个GSM编号，里面包含测量数据和元数据——生物来源、处理条件啥的。同一个平台可以对应几十上百个样本。

数据集（GSE编号）

一个GSE就是一项研究的所有样本集合。这是你下载数据的主要单位。比如GSE12345，里面可能包含正常组和疾病组的芯片数据。

数据存储有三种形式：SOFT格式（纯文本，包含完整元数据和表达矩阵）、Series Matrix（压缩的矩阵文件，行是基因，列是样本）、原始数据（CEL或FASTQ文件）。多数时候你只需要Series Matrix就够了，除非要做自己的标准化。

下载数据，我试了四种方式

网页检索

NCB上直接搜关键词或GSE编号，用过滤器缩小范围。这个最直观，但一次只能下载一个数据集，批量操作不行。

R语言（推荐）

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("GEOquery")
library(GEOquery)
gset <- getGEO("GSE12345", destdir = "./data/", GSEMatrix = TRUE)

GEOquery包是我最常用的，自动下载、解析、返回ExpressionSet对象，后面直接接limma或DESeq2就行。缺点是偶尔会因为网络问题中断，我习惯用getOption设置超时时间。

Python

pip install GEOparse
import GEOparse
gse = GEOparse.get_GEO("GSE12345", destdir="./")

Python生态用户可以用GEOparse，跟pandas、scikit-learn集成很方便。但它的文档比R包少，遇到格式不标准的GSE会报错。

FTP直接下载

适合批量下载原始数据。路径结构是ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE5nnn/GSE5290/soft/GSE5290.soft.gz。我一般只在需要CEL文件时用这个，因为10MB以下的微阵列数据直接FTP比API快。

质量控制，别信那些“直接跑差异分析”的教程

每次看到有人拿到数据就跑差异分析，我都想提醒一句：先做质控！不然结果全是噪声。

我通常从三个层面看：

样本水平：PCA和层次聚类。离群样本一眼就能看出来。有一次我跑了个GSE，PCA图上对照组和处理组完全混在一起，仔细一看，原来有两个样本标签贴反了。

探针水平：信号强度过滤。我会剔除在超过90%样本中表达值低于1的探针。这个阈值可以根据平台调整，但默认1是个不错的起点。

批次效应：如果样本来自不同批次，一定要看PCA上批次是否分离。分离了就要校正。R包arrayQualityMetrics能自动生成质控报告，省很多事。

数据预处理，探针转基因是第一大坑

探针→基因转换

这是最容易出问题的一步。一个基因可能对应多个探针，常见做法是取中位数或均值。但有个坑：平台注释文件有时候会过时，或者某些探针没有注释到基因。我处理过一个Affymetrix平台，探针总数5.4万个，能匹配到基因的只有2.8万个。剩下那些就丢了，你敢信？

解决方案：用最新的注释包，比如hgu133a.db，定期去Bioconductor更新。实在不行就手动去NCBI官网下载最新版。

标准化

芯片数据我推荐RMA（affy::rma()）或分位数标准化（limma::normalizeBetweenArrays()）。RNA-seq数据用DESeq2的vst转换或者edgeR的TMM标准化。

注意：RMA只适用于Affymetrix芯片。Agilent的数据要用limma::normalizeWithinArrays()。

批次效应校正

ComBat（sva包里的）是目前最常用的。用法如下：

library(sva)
combat_data <- ComBat(dat = exprs(eset), batch = batch_info, mod = model_matrix)

但有个限制：如果你只有两个批次，每批次样本少于3个，ComBat的估计会不稳定，甚至可能放大噪声。这时候我会选择removeBatchEffect（limma包），或者在差异分析模型里直接把批次作为协变量。

差异表达分析，选对工具很关键

实验设计

先构建分组向量，明确对照组和处理组。对比矩阵用makeContrasts定义。

方法选择

我常用limma处理芯片数据，因为它对异常值不敏感，而且支持voom转换直接处理RNA-seq计数。跑完结果后，筛选标准：调整后p值<0.05，log2(FC)绝对值≥1。

功能富集分析，clusterProfiler是真爱

拿到差异基因后，下一步就是看它们参与了哪些通路。我用clusterProfiler包做GO和KEGG富集，代码很简单：

library(clusterProfiler)
library(org.Hs.eg.db)
ego <- enrichGO(gene = rownames(top_genes),
                OrgDb = org.Hs.eg.db,
                keyType = "SYMBOL",
                ont = "BP",
                pAdjustMethod = "fdr",
                pvalueCutoff = 0.05)

但注意：如果差异基因数量太少（比如只有几十个），富集结果往往不显著。这时候可以适当放宽p值或logFC的阈值，但一定要在论文里说明。

可视化，让数据自己说话

火山图展示差异基因的显著性和变化幅度。热图展示表达模式。富集分析的结果用气泡图或网络图。

我用ggplot2画火山图，pheatmap画热图。注意：热图最好按行标准化，否则不同基因的表达量范围不同，颜色会失真。

实用技巧和踩过的坑

数据选择策略

• 优先找跟自己研究物种、组织、处理条件匹配的数据集
• 芯片数据选Affymetrix或Agilent，注释完善
• 每组至少3个生物学重复，越多越好

常见问题

1. 批次效应过强：用ComBat校正，或者在模型中加入批次。如果批次效应大到完全覆盖了生物学差异，那就只保留同一批次的数据。
2. 差异基因太少：检查质控是否太严。logFC阈值从1.5降到1，但要小心假阳性。也可能是实验设计的问题，处理效应本身就不强。
3. 富集结果不理想：换数据库试试，GO不行就KEGG，再不行Reactome。另外检查探针→基因转换是否准确，这一步丢基因太多，富集肯定没戏。

工具推荐

我的建议：从简单数据集开始练手

别一上来就挑战多因素多批次的数据。先找一个单因素设计的芯片数据（比如正常vs疾病），跑通完整流程——下载、质控、预处理、差异分析、富集。回头再看这篇文章，你会发现每一步都是有道理的。

等熟练了，再尝试整合多组学数据，或者做meta分析。但记住：GEO数据再大也是别人做的实验，你合成的数据永远替代不了自己的实验验证。分析结果顶多帮你提出假说，别把它当成真理。

常见问题（FAQ）

GEO数据下载后，探针转基因怎么那么麻烦？

平台注释文件常过时，很多探针无法匹配到基因。建议使用最新Bioconductor注释包（如hgu133a.db）或手动从NCBI下载最新注释文件。

批次效应太强怎么办？ComBat校正好用吗？

ComBat常用但样本量小（如每批<3个）时可能放大噪声。此时可用limma的removeBatchEffect或在模型中将批次作为协变量。

GEO数据那么多，怎么选靠谱的数据集？

优先选研究物种、组织、处理条件匹配的，每组至少3个重复。Affymetrix或Agilent芯片注释更完善，避免选择样本量过小或批次效应明显的数据。