GEO数据挖掘与重分析¶
一句话概述¶
从NCBI GEO数据库下载公共转录组数据，使用GEOquery/GEO2R工具进行标准化、差异表达分析和多数据集meta-analysis，实现文献数据的二次挖掘与验证。
核心知识点总览¶
| 知识点 | 关键内容 | 重要程度 |
|---|---|---|
| GEO数据库结构 | GSE/GPL/GSM层级关系 | ⭐⭐⭐⭐⭐ |
| GEOquery下载 | R包自动下载与解析 | ⭐⭐⭐⭐⭐ |
| 芯片数据标准化 | RMA/MAS5/quantile归一化 | ⭐⭐⭐⭐ |
| 差异表达分析 | limma对GEO数据的应用 | ⭐⭐⭐⭐⭐ |
| 探针注释 | 探针ID到基因Symbol映射 | ⭐⭐⭐⭐ |
| Meta-analysis | 多数据集联合分析/效应量整合 | ⭐⭐⭐⭐ |
| RNA-seq GEO数据 | SRA下载+重新定量 | ⭐⭐⭐ |
| 批次效应校正 | ComBat/SVA跨数据集整合 | ⭐⭐⭐ |
各步骤详解¶
第一步：GEO数据库检索与数据结构¶
白话解释： GEO（Gene Expression Omnibus）是全球最大的基因表达数据仓库。每个"数据集"（GSE）包含多个"样本"（GSM），使用特定"平台"（GPL，如Affymetrix HG-U133A芯片）。数据挖掘的第一步是找到与你研究相关的数据集。
技术细节：
# === GEO数据检索与下载 ===
library(GEOquery)
# 下载GSE数据集
gse <- getGEO("GSE68465", GSEMatrix = TRUE, getGPL = TRUE)
# 返回ExpressionSet对象列表（可能多个平台）
# 查看数据基本信息
eset <- gse[[1]]
show(eset)
dim(exprs(eset)) # 基因数 × 样本数
head(pData(eset)) # 临床/表型数据
head(fData(eset)) # 探针注释
# 提取表达矩阵
expr_matrix <- exprs(eset)
# 提取临床信息
clinical_info <- pData(eset)
# 通常需要整理：从杂乱的characteristics列中提取有用信息
clinical_info$group <- ifelse(grepl("tumor", clinical_info$`disease state:ch1`),
"Tumor", "Normal")
# === 批量下载多个GSE ===
gse_ids <- c("GSE68465", "GSE30219", "GSE31210")
gse_list <- lapply(gse_ids, function(id) {
getGEO(id, GSEMatrix = TRUE, getGPL = FALSE)[[1]]
})
第二步：数据标准化与质控¶
白话解释： GEO下载的数据可能已经标准化过（看experimentData），也可能是原始数据。需要确认数据状态，必要时重新标准化。常见问题包括：表达值范围异常、样本间分布不一致、异常离群样本等。
技术细节：
# === 质控与标准化 ===
library(limma)
library(affy)
# 1. 检查数据是否已标准化
summary(expr_matrix[1:5, 1:5])
# 如果值范围在2-15左右，可能已经log2转换
# 如果值很大（>100），可能需要log2转换
# 判断是否需要log转换
needs_log2 <- max(expr_matrix, na.rm = TRUE) > 30
if (needs_log2) {
expr_matrix <- log2(expr_matrix + 1)
}
# 2. Quantile归一化（如果样本间分布不一致）
library(preprocessCore)
expr_norm <- normalize.quantiles(expr_matrix)
dimnames(expr_norm) <- dimnames(expr_matrix)
# 3. 质控可视化
# 箱线图
boxplot(expr_norm[, 1:20], las = 2, main = "Sample Distribution")
# PCA检查批次效应和离群样本
pca <- prcomp(t(expr_norm), scale. = TRUE)
plot(pca$x[, 1:2], col = as.numeric(factor(clinical_info$group)),
pch = 19, main = "PCA")
# 4. 去除异常样本
# 基于PCA或样本间相关性识别离群值
sample_cor <- cor(expr_norm)
mean_cor <- colMeans(sample_cor)
outliers <- names(mean_cor[mean_cor < mean(mean_cor) - 2*sd(mean_cor)])
if (length(outliers) > 0) {
expr_norm <- expr_norm[, !colnames(expr_norm) %in% outliers]
clinical_info <- clinical_info[!rownames(clinical_info) %in% outliers, ]
}
第三步：探针ID到基因注释映射¶
白话解释： 芯片数据中每行是一个"探针"（probe），不是一个基因。需要将探针ID映射到基因Symbol。一个基因可能对应多个探针，此时通常取表达最高或最变异的那个探针代表该基因。
技术细节：
# === 探针注释 ===
# 方法1：从fData中直接获取
probe_anno <- fData(eset)
# 常见的基因列名：Gene.symbol, GENE_SYMBOL, gene_assignment
gene_col <- grep("symbol|gene", colnames(probe_anno), ignore.case = TRUE, value = TRUE)[1]
probe_to_gene <- probe_anno[, c("ID", gene_col)]
colnames(probe_to_gene) <- c("probe_id", "gene_symbol")
# 方法2：使用平台特异性注释包
# BiocManager::install("hgu133plus2.db")
library(hgu133plus2.db)
probe_ids <- rownames(expr_norm)
gene_symbols <- mapIds(hgu133plus2.db, keys = probe_ids,
column = "SYMBOL", keytype = "PROBEID",
multiVals = "first")
# 方法3：使用biomaRt
library(biomaRt)
ensembl <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
anno <- getBM(attributes = c("affy_hg_u133_plus_2", "hgnc_symbol"),
filters = "affy_hg_u133_plus_2",
values = probe_ids, mart = ensembl)
# === 多探针取一：最大均值法 ===
# 对多个探针映射到同一基因的情况，取平均表达最高的探针
expr_df <- data.frame(gene = gene_symbols, expr_norm, check.names = FALSE)
expr_df <- expr_df[!is.na(expr_df$gene) & expr_df$gene != "", ]
# 按基因取最大均值探针
expr_df$mean_expr <- rowMeans(expr_df[, -1])
expr_unique <- expr_df %>%
group_by(gene) %>%
slice_max(mean_expr, n = 1) %>%
ungroup()
# 整理为矩阵
rownames(expr_unique) <- expr_unique$gene
expr_gene <- as.matrix(expr_unique[, -c(1, ncol(expr_unique))])
第四步：差异表达分析¶
白话解释： 用limma做差异分析——跟原始研究用的方法类似，但你可以用自己定义的分组来发现原文没报告的差异基因。GEO2R是NCBI提供的在线版本，适合快速查看；自己写代码则更灵活。
技术细节：
# === limma差异分析 ===
library(limma)
# 确保分组信息正确
group <- factor(clinical_info$group, levels = c("Normal", "Tumor"))
# 设计矩阵
design <- model.matrix(~ 0 + group)
colnames(design) <- levels(group)
# 拟合
fit <- lmFit(expr_gene, design)
contrast <- makeContrasts(Tumor - Normal, levels = design)
fit2 <- contrasts.fit(fit, contrast)
fit3 <- eBayes(fit2)
# 结果
de_results <- topTable(fit3, number = Inf, sort.by = "P")
de_results$gene <- rownames(de_results)
# 筛选
sig_genes <- de_results[de_results$adj.P.Val < 0.05 & abs(de_results$logFC) > 1, ]
cat(sprintf("DEGs: %d (Up: %d, Down: %d)\n",
nrow(sig_genes), sum(sig_genes$logFC > 0), sum(sig_genes$logFC < 0)))
# 火山图
library(EnhancedVolcano)
EnhancedVolcano(de_results, lab = de_results$gene,
x = 'logFC', y = 'adj.P.Val',
title = "GSE68465 - Tumor vs Normal")
第五步：Meta-analysis——多数据集联合分析¶
白话解释： 单个数据集可能样本量不够、平台偏差大。Meta-analysis将多个独立GEO数据集的结果联合起来，通过效应量合并得到更可靠的差异基因。好比多项临床试验的荟萃分析，一篇文章的结论不够，多篇放在一起更有说服力。
技术细节：
# === Meta-analysis ===
library(metafor)
library(MetaIntegrator)
# 方法1：基于效应量的meta-analysis
# 对每个数据集分别做差异分析，获取效应量(logFC)和标准误
# 假设已有多个数据集的差异分析结果
study_results <- list(
GSE68465 = de_results_1, # 含logFC, SE, P列
GSE30219 = de_results_2,
GSE31210 = de_results_3
)
# 对每个基因做随机效应meta-analysis
meta_results <- data.frame()
common_genes <- Reduce(intersect, lapply(study_results, rownames))
for (gene in common_genes) {
effects <- sapply(study_results, function(x) x[gene, "logFC"])
ses <- sapply(study_results, function(x) x[gene, "SE"])
# 如果没有SE，用logFC/t-statistic估计
if (any(is.na(ses))) next
res <- rma(yi = effects, sei = ses, method = "REML")
meta_results <- rbind(meta_results, data.frame(
gene = gene,
meta_logFC = res$beta[1],
meta_se = res$se,
meta_pval = res$pval,
meta_I2 = res$I2, # 异质性
n_studies = length(effects)
))
}
meta_results$meta_fdr <- p.adjust(meta_results$meta_pval, method = "BH")
meta_sig <- meta_results[meta_results$meta_fdr < 0.05 & abs(meta_results$meta_logFC) > 0.5, ]
# 方法2：使用MetaIntegrator包
# 标准化的meta-analysis框架
metaObject <- list(
originalData = list(
GSE68465 = list(expr = expr1, pheno = pheno1, class = class1, keys = keys1),
GSE30219 = list(expr = expr2, pheno = pheno2, class = class2, keys = keys2)
)
)
metaResult <- runMetaAnalysis(metaObject)
第六步：RNA-seq GEO数据处理¶
白话解释： 越来越多的GEO数据是RNA-seq。这些通常存储在SRA中，需要下载FASTQ后重新定量分析。或者如果GEO中提供了count矩阵，可以直接使用DESeq2分析。
技术细节：
# === 从SRA下载RNA-seq数据 ===
# 使用SRA Toolkit
# 查找GSE对应的SRA project
# 如GSE157103对应PRJNA660427
# 批量下载
prefetch --max-size 50G SRR12345678
fastq-dump --split-files --gzip SRR12345678
# 或使用更快的fasterq-dump
fasterq-dump --split-files SRR12345678
gzip SRR12345678_1.fastq SRR12345678_2.fastq
# === 如果GEO提供了count矩阵 ===
# 从supplementary files下载
# 处理GEO中的RNA-seq count数据
library(DESeq2)
# 下载supplementary file
gse_supp <- getGEOSuppFiles("GSE157103")
counts <- read.csv("GSE157103_raw_counts.csv", row.names = 1)
# DESeq2分析
dds <- DESeqDataSetFromMatrix(countData = counts,
colData = sample_info,
design = ~ condition)
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "COVID", "Healthy"))
第七步：跨数据集整合与批次效应校正¶
白话解释： 当想合并多个GEO数据集的原始数据进行联合分析时，需要处理批次效应——不同实验室、不同时间、不同平台产生的数据有系统性差异。ComBat是最常用的批次校正方法。
技术细节：
# === 批次效应校正 ===
library(sva)
# 合并表达矩阵（只保留共同基因）
common_genes <- Reduce(intersect, list(rownames(expr1), rownames(expr2), rownames(expr3)))
merged_expr <- cbind(expr1[common_genes, ], expr2[common_genes, ], expr3[common_genes, ])
# 批次信息
batch <- c(rep(1, ncol(expr1)), rep(2, ncol(expr2)), rep(3, ncol(expr3)))
# 保留的生物学变量
group <- c(group1, group2, group3)
mod <- model.matrix(~ group)
# ComBat校正
combat_expr <- ComBat(dat = merged_expr, batch = batch, mod = mod, par.prior = TRUE)
# 验证：PCA应不再按批次聚类
pca_before <- prcomp(t(merged_expr))
pca_after <- prcomp(t(combat_expr))
par(mfrow = c(1, 2))
plot(pca_before$x[, 1:2], col = batch, main = "Before ComBat")
plot(pca_after$x[, 1:2], col = batch, main = "After ComBat")
实战命令速查¶
# GEO快速分析
library(GEOquery); library(limma)
gse <- getGEO("GSE_ID")[[1]]
expr <- exprs(gse); pheno <- pData(gse)
design <- model.matrix(~0+group)
fit <- eBayes(contrasts.fit(lmFit(expr, design), makeContrasts(A-B, levels=design)))
topTable(fit, number=Inf)
面试常问点¶
Q1: GEO数据挖掘的典型研究思路是什么？¶
A: (1) 确定研究问题→搜索相关GSE数据集；(2) 下载并质控数据；(3) 差异分析寻找候选基因；(4) 多数据集meta-analysis验证一致性；(5) 功能富集（GO/KEGG）理解生物学意义；(6) 生存分析评估临床价值；(7) 构建预后模型或基因签名。这种"纯干"研究成本低、速度快，适合初步探索和验证假设。
Q2: 如何处理GEO中不同平台的数据整合？¶
A: (1) 基因层面整合——将不同平台的探针映射到基因Symbol，只保留共同基因；(2) 效应量meta-analysis——各平台独立分析后在统计量层面整合（不直接合并原始数据）；(3) 如果必须合并原始数据，使用ComBat校正平台效应，但需确保生物学变量不与平台完全混淆（confounded）。
Q3: 为什么meta-analysis比单个大数据集分析更可靠？¶
A: (1) 增加样本量提高统计效力；(2) 跨平台/跨队列的一致信号更不可能是技术偶然；(3) 异质性分析（I²统计量）能识别不可靠的结果；(4) 减少单个研究的特异偏差（如样本选择偏差）。但前提是各研究的入排标准和生物学条件足够相似。
Q4: GEO数据挖掘文章的常见审稿问题？¶
A: (1) "结果缺乏实验验证"——至少需要qPCR/WB验证关键基因；(2) "样本临床信息不完整"——GEO经常缺失重要协变量；(3) "批次效应未处理"——跨数据集分析必须说明如何控制；(4) "多重检验校正不足"；(5) "结果在独立队列中未验证"——training+validation设计更有说服力。
Q5: 如何从GEO中选择高质量数据集？¶
A: 评估标准：(1) 样本量足够（每组至少15-20例）；(2) 平台注释完善（有对应的Bioconductor注释包）；(3) 临床信息详细（有分组、生存数据等）；(4) 原始数据可用（有raw data或supplementary count matrix）；(5) 发表在高质量期刊（实验可靠性）；(6) 研究设计清晰，与你的问题匹配。
易错点¶
1. 使用已标准化数据再次标准化¶
GEO中的"Series Matrix"数据通常已经做过某种标准化。再次做quantile normalization可能过度处理。应先确认数据处理状态。
2. 探针注释过时¶
GEO中平台注释可能多年未更新。建议使用最新的Bioconductor注释包或biomaRt重新映射。
3. 临床信息解析错误¶
GEO的phenoData格式混乱，同一变量在不同GSE中列名和编码方式不同。务必仔细检查分组是否正确。
4. 忽略数据中的技术混杂¶
芯片的扫描日期、RNA提取方法、样本保存时间等可能与生物学分组相关（如所有tumor样本在同一天处理）。这种confounding难以统计校正。
5. Meta-analysis中混入低质量研究¶
质量很差的数据集加入meta-analysis不会增加信息，反而增加噪声。应做敏感性分析（逐一去除各研究看结果是否稳定）。
补充知识¶
相关数据库与工具¶
- ArrayExpress：EBI的表达数据仓库（与GEO互补）
- TCGA/GTEx：癌症/正常组织大型数据集
- GEO2R：NCBI在线分析工具
- shinyGEO：交互式GEO分析Web应用
- networkanalyst：网络分析+GEO整合平台
引用推荐¶
- GEOquery: Davis & Meltzer, Bioinformatics, 2007
- limma: Ritchie et al., Nucleic Acids Research, 2015
- ComBat: Johnson et al., Biostatistics, 2007
- MetaIntegrator: Haynes et al., Pacific Symposium on Biocomputing, 2017（注意：该包已于2024年从CRAN存档，可从archive获取）