多GEO数据联合分析实战：从差异基因到生存验证与hub基因筛选

简介：这份PDF面向生物信息学与数据挖掘方向的研究生、科研人员及生信初学者，聚焦多套GEO数据集联合分析的完整思路，帮助读者掌握从公共数据库挖掘胃癌差异基因并验证其临床价值的方法。资源为单个PDF文件，压缩包约3.8MB，内容涵盖GEO数据检索与过滤、CEL文件读取与log转换、多探针取均值、limma标准化与差异表达分析、RobustRankAggreg整合排名、pheatmap热图可视化，以及TCGA RNA-seq验证、生存分析和GO富集等模块。文中以胃癌为例，共识别出1210个差异基因（436个上调、774个下调），并给出差异基因列表、生存相关基因表及富集图表，读者可据此复现分析流程、理解多数据集整合策略，并借鉴其统计检验与可视化排错思路。目前已有129人学习，适合需要系统掌握GEO联合分析与TCGA验证的进阶学习者参考。

1. 从六套GEO芯片到一张生存曲线：这份生物信息分析报告能直接抄什么

手里攒着好几套胃癌表达谱，GSE64041、GSE67764、GSE62232、GSE19665、GSE46408、GSE25097，平台从GPL570到GPL6244各不相同，样本量从十几例到几百例，想合并分析却卡在批次效应和探针注释上——这是做肿瘤生信的人几乎都会撞上的墙。这份《数据挖掘与数据分析应用 数据挖掘之生物信息分析报告 多个GEO数据联合分析.pdf》就是一份把这条路走通的完整记录：从GEO检索、数据过滤、limma差异分析，到RobustRankAggreg跨数据集整合、TCGA验证、生存分析、GO与PATHWAY富集、蛋白互作网络，每一步都有参数、有结果表、有图。它适合两类人：刚接触GEO联合分析、需要一份能照着跑通全流程的新手；以及做过单数据集差异分析、想补上多数据集整合与生存验证这一环的熟手。下面我按自己复现时的顺序，把这份报告拆开讲。

2. GEO数据检索与下载：关键词、过滤条件与平台选择

2.1 检索式怎么写才不漏数据

GEO的检索语法比很多人想的要严格。报告里用的检索式是：

# GEO检索式：胃癌 + 人类物种限定 (gastric cancer) AND "Homo sapiens"[porgn:__txid9606]

这里有两个关键点。第一，[porgn:__txid9606]是NCBI Taxonomy的物种限定字段，9606就是人类，不加这个限定会混进小鼠、大鼠的胃癌模型数据，后面做注释时平台对不上，白跑一遍。第二，关键词用gastric cancer而不是gastriccancer，GEO的检索是基于分词匹配的，连写反而可能漏掉标题里带空格的Series。

检索出来之后不要急着全下。报告里做了一步过滤：去掉没有重复试验的样品。所谓重复试验，就是同一实验条件下有多个生物学重复。只有单样本的Series，做差异分析时组内方差无法估计，limma跑出来的P值不可信。常见做法是在GEO页面看Series的Sample列表，每个分组至少要有3个样本才纳入。

2.2 六套数据的平台与样本分布

报告最终纳入的六套数据如下表。这张表本身就是选数据的决策依据，值得逐列看：

选平台时有几个判断：GPL570（Affymetrix HG-U133 Plus 2.0）是最经典的芯片平台，注释全、探针多，GSE62232和GSE19665都用它，合并时注释可以直接复用。GPL6244和GPL17077是较新的平台，探针设计不同，需要各自的注释包。GPL4133和GPL10687属于Agilent和定制平台，样本量小，主要起补充验证作用。

提示：平台不一致时，不要试图把原始荧光值直接合并。正确做法是每个平台各自做标准化和差异分析，最后在基因水平上整合，这也是报告采用RobustRankAggreg的原因。

2.3 下载表达矩阵还是CEL文件

报告里区分了两种情况：有表达矩阵的直接下载矩阵，只有CEL文件的用affy包读取。这个判断很实际——GEO上不是每个Series都提供处理好的矩阵，有些只上传了原始CEL。

# 用affy读取CEL文件并做RMA标准化 library(affy) setwd("GSE64041_RAW") raw_data <- ReadAffy() # 读取当前目录下所有CEL文件 eset <- rma(raw_data) # RMA标准化：背景校正+分位数归一化+汇总 expr_matrix <- exprs(eset) # 提取表达矩阵 write.csv(expr_matrix, "GSE64041_expr.csv")

ReadAffy()会自动扫描工作目录下所有.CEL或.CEL.gz文件，不用手动列文件名。rma()是Affymetrix芯片的标准预处理流程，三步合一：背景校正去掉光学噪声，分位数归一化让各样本表达分布对齐，最后用中位数 polish 汇总探针集。如果你的数据已经是对数化过的矩阵，就跳过这一步，直接进入下一步的log处理。

3. 差异表达分析：limma标准化、探针合并与跨数据集整合

3.1 取log与多探针合并

芯片表达值如果没取log，分布会严重右偏，做t检验时方差大的基因会主导结果。报告的处理是：没有取log的值先取log2。判断方法很简单，看表达矩阵的最大值——如果超过100，基本就是原始尺度，需要取log；如果最大值在20以内，通常已经是对数尺度。

多探针合并是另一个容易翻车的点。同一个基因在芯片上往往有多个探针，报告的做法是取所有探针的平均值：

# 多探针取平均合并到基因水平 library(hgu133plus2.db) # GPL570对应的注释包 probe_ids <- rownames(expr_matrix) gene_symbols <- mapIds(hgu133plus2.db, keys = probe_ids, column = "SYMBOL", keytype = "PROBEID", multiVals = "first") # 多探针时先取第一个映射 expr_matrix <- expr_matrix[!is.na(gene_symbols), ] gene_symbols <- gene_symbols[!is.na(gene_symbols)] # 按基因名聚合，取平均值 expr_by_gene <- aggregate(expr_matrix, by = list(gene = gene_symbols), FUN = mean) rownames(expr_by_gene) <- expr_by_gene$gene expr_by_gene$gene <- NULL

mapIds()的multiVals参数控制一个探针映射到多个基因时的行为，这里用"first"取第一个。aggregate()按基因名分组求均值，得到每个基因一个表达值。注意聚合前要确保矩阵是数值型，如果读进来是data.frame带字符列，先转成matrix。

3.2 limma差异分析的标准流程

每个数据集单独做limma分析，结果保存在01_limma目录。limma的核心是设计矩阵和对比矩阵：

# limma差异表达分析 library(limma) group <- factor(c(rep("Normal", n_normal), rep("Tumor", n_tumor)), levels = c("Normal", "Tumor")) design <- model.matrix(~ 0 + group) # 无截距设计矩阵 colnames(design) <- levels(group) fit <- lmFit(expr_by_gene, design) # 线性拟合 contrast <- makeContrasts(Tumor - Normal, levels = design) fit2 <- contrasts.fit(fit, contrast) fit2 <- eBayes(fit2) # 经验贝叶斯平滑方差 deg <- topTable(fit2, number = Inf, adjust.method = "BH")

model.matrix(~ 0 + group)生成无截距的设计矩阵，每个分组一列，这样对比矩阵可以直接写Tumor - Normal。eBayes()是limma的精髓，它把每个基因的方差向全局方差收缩，样本量小时尤其重要，否则单基因方差估计噪声太大。topTable()的adjust.method = "BH"是Benjamini-Hochberg校正，控制FDR。

3.3 RobustRankAggreg：为什么不用取交集

六个数据集各自跑完limma后，怎么整合？取交集太粗暴——要求基因在六个数据集中都显著，会漏掉只在部分数据集中显著的基因。报告用的是RobustRankAggreg（RRA），它的零假设是：每个基因在每个实验中的排名是随机的。如果一个基因在所有实验中都排在前列，它随机出现的概率就极低，P值就小。

# RobustRankAggreg跨数据集整合 library(RobustRankAggreg) # 每个数据集的基因按logFC排序，提取基因名列表 gene_lists <- list() for (i in 1:6) { deg <- read.csv(paste0("01_limma/GSE", i, "_limma.csv")) deg <- deg[order(deg$logFC, decreasing = TRUE), ] gene_lists[[i]] <- deg$gene } # RRA分析 rra_result <- aggregateRanks(glist = gene_lists, N = 20000) # N为背景基因总数 rra_result$adjPval <- p.adjust(rra_result$Score, method = "bonferroni") sig_genes <- rra_result[rra_result$adjPval < 0.05, ]

aggregateRanks()的N参数是背景基因总数，人类一般填20000左右，这个值影响P值的计算，填错会系统性偏移。method = "bonferroni"是报告指定的校正方法，比BH更严格，适合这种多数据集整合场景。最终得到1210个差异基因，其中436个上调、774个下调。

3.4 热图：用pheatmap看一致性

差异基因的热图不只是好看，它能直接告诉你整合结果可不可信。报告用pheatmap画了最上调和最下调各20个基因的热图，横坐标是GEO ID，纵坐标是基因名，红色logFC>0，绿色logFC<0。

# 差异基因热图 library(pheatmap) top_up <- head(sig_genes[order(sig_genes$logFC, decreasing = TRUE), ], 20) top_down <- head(sig_genes[order(sig_genes$logFC), ], 20) top_genes <- c(top_up$gene, top_down$gene) # 从各数据集提取这些基因的logFC logfc_matrix <- extract_logfc(top_genes, datasets = 1:6) pheatmap(logfc_matrix, cluster_rows = FALSE, cluster_cols = FALSE, color = colorRampPalette(c("green", "black", "red"))(100), display_numbers = TRUE)

看热图时重点看两件事：上调基因是不是在所有数据集里都偏红，下调基因是不是都偏绿。如果某个基因在某个数据集里颜色反了，说明这个数据集的批次效应没处理好，或者样本分组有问题，需要回头检查。

4. TCGA验证与生存分析：从差异基因到临床相关基因

4.1 TCGA数据下载与矩阵合并

GEO得到的差异基因需要独立数据验证。报告从TCGA下载胃癌RNA-seq level 3数据，50个正常组织、374个肿瘤组织。TCGA的下载文件是每个样本一个FPKM文件，需要合并成矩阵：

# Perl合并TCGA单个样本FPKM文件 use strict; use warnings; my $dir = "TCGA_STAD_FPKM"; opendir(DH, $dir) or die $!; my @files = grep { /\.txt$/ } readdir(DH); closedir(DH); my %data; foreach my $file (@files) { open(FH, "$dir/$file") or die $!; while (<FH>) { chomp; next if /^gene_id/; my ($gene, $fpkm) = (split /\t/)[0, 1]; $data{$gene}{$file} = $fpkm; } close(FH); } # 输出合并矩阵 open(OUT, ">TCGA_STAD_matrix.txt"); print OUT "gene\t", join("\t", @files), "\n"; foreach my $gene (keys %data) { print OUT "$gene\t"; print OUT join("\t", map { $data{$gene}{$_} // "NA" } @files), "\n"; } close(OUT);

这段Perl的逻辑是：遍历目录下所有FPKM文件，按基因名聚合，每个基因一行，每个样本一列。// "NA"处理缺失值，避免某个样本缺基因时矩阵错位。合并后的矩阵用read.table()读进R，行名设为基因名。

4.2 Wilcoxon检验验证差异基因

TCGA数据是RNA-seq的FPKM，和GEO芯片的表达值尺度不同，不能直接比logFC。报告用的是Wilcoxon秩和检验，非参数方法，不假设正态分布，适合这种跨平台验证：

# Wilcoxon检验验证GEO差异基因 geo_deg <- read.csv("RRA_sig_genes.csv") tcga_expr <- read.table("TCGA_STAD_matrix.txt", header = TRUE, row.names = 1) tcga_expr <- log2(tcga_expr + 1) # FPKM取log normal_cols <- grep("Normal", colnames(tcga_expr)) tumor_cols <- grep("Tumor", colnames(tcga_expr)) results <- data.frame() for (gene in geo_deg$gene) { if (gene %in% rownames(tcga_expr)) { normal_val <- as.numeric(tcga_expr[gene, normal_cols]) tumor_val <- as.numeric(tcga_expr[gene, tumor_cols]) test <- wilcox.test(tumor_val, normal_val) results <- rbind(results, data.frame( gene = gene, pvalue = test$p.value, adjPvalue = p.adjust(test$p.value, method = "BH") )) } }

log2(FPKM + 1)的+1是为了避免log0。grep("Normal", colnames())从列名里识别分组，前提是列名里包含Normal/Tumor标识，TCGA的样本ID里有-01（肿瘤）和-11（正常）的后缀，也可以用这个来分。验证结果里LYPD1、NQO1、MKI67等基因的adjPvalue都在0.01以下，说明GEO整合结果在TCGA里站得住。

4.3 生存分析：survival包与log-rank检验

验证完差异，下一步是看这些基因和生存有没有关系。报告用survival包做生存分析，log-rank检验，过滤条件Pvalue<0.01，最终找到272个与胃癌生存相关的差异基因。

# 生存分析 library(survival) library(survminer) surv_data <- read.csv("TCGA_survival_data.csv") expr_data <- read.csv("TCGA_STAD_matrix.csv") merged <- merge(surv_data, expr_data, by = "sample") for (gene in sig_genes$gene) { if (gene %in% colnames(merged)) { # 按中位数分高低表达组 median_expr <- median(merged[[gene]], na.rm = TRUE) merged$group <- ifelse(merged[[gene]] > median_expr, "High", "Low") fit <- survfit(Surv(OS.time, OS) ~ group, data = merged) diff <- survdiff(Surv(OS.time, OS) ~ group, data = merged) pval <- 1 - pchisq(diff$chisq, length(diff$n) - 1) if (pval < 0.01) { ggsurvplot(fit, data = merged, pval = TRUE, palette = c("red", "blue"), title = gene, xlab = "Time (days)", ylab = "Overall Survival") } } }

Surv(OS.time, OS)构建生存对象，OS.time是生存时间，OS是生存状态（1为死亡，0为删失）。survfit()做Kaplan-Meier估计，survdiff()做log-rank检验。按中位数分高低表达组是最常用的分组方式，简单但有效。ggsurvplot()画出的曲线里，红色是高表达组，蓝色是低表达组，横坐标生存时间，纵坐标总生存率。

注意：生存分析里P值过滤到0.01比0.05更严格，因为这里做了几百次检验，多重比较问题严重。报告用0.01是在敏感性和特异性之间取的折中。

5. 功能富集与蛋白互作：GO、PATHWAY与String网络

5.1 GO富集：DAVID的使用与结果解读

GO富集回答的是“这些差异基因在功能上有什么共同点”。报告用DAVID做GO功能富集，P<0.05为筛选条件。DAVID的输入是基因列表，背景选人类全基因组，输出三个本体：BP（生物过程）、CC（细胞组分）、MF（分子功能）。

看GO结果时，Count是落在该Term上的差异基因数，PValue是超几何检验的显著性。Count太小的Term（比如只有1个基因）即使P值显著，生物学意义也有限。报告里superoxide dismutase activity有3个基因，mast cell granule有4个，这些是有解读价值的。

5.2 PATHWAY富集：KEGG通路与疾病关联

PATHWAY富集用的是KEGG数据库，报告里P<0.05的通路包括Drug metabolism、Hepatocellular carcinoma、Systemic lupus erythematosus等。这里有个细节值得注意：Hepatocellular carcinoma通路出现在胃癌分析里，不是错误——KEGG的疾病通路往往共享基因，SHC1、FZD7、GSTM1这些基因在多个癌症通路里都有角色。

解读PATHWAY结果时，不要只看P值最小的通路。要结合基因列表看：如果某个通路的成员基因里包含你关注的差异基因，这个通路就值得深挖。报告里Drug metabolism通路的四个基因全是GST家族，提示胃癌中谷胱甘肽代谢可能被扰动。

5.3 蛋白互作网络：String与邻接节点分析

蛋白互作网络用String构建，输入是生存相关基因。String的输出是互作关系对，报告用R绘制了邻接节点数目图——每个蛋白有多少个互作伙伴，邻接节点越多，说明它在网络里越核心。

# 蛋白互作网络邻接节点统计 library(igraph) edges <- read.table("string_interactions.tsv", header = TRUE, sep = "\t") g <- graph_from_data_frame(edges[, 1:2], directed = FALSE) degree_df <- data.frame( gene = names(degree(g)), degree = as.numeric(degree(g)) ) degree_df <- degree_df[order(degree_df$degree, decreasing = TRUE), ] # 画条形图 barplot(degree_df$degree[1:20], names.arg = degree_df$gene[1:20], las = 2, col = "steelblue", ylab = "Number of adjacent nodes", main = "Hub genes in PPI network")

graph_from_data_frame()把边列表转成图对象，degree()返回每个节点的度。报告里CFTR、SST、TIMP1的邻接节点数最多，位于网络核心。这些hub基因往往是调控的关键节点，后续做实验验证时优先选它们。

提示：String的互作关系有置信度评分（combined_score），默认阈值是0.4，做网络分析时建议调到0.7以上，只保留高置信度互作，否则网络太密，hub基因不突出。

6. 复现时最容易翻车的五个地方

现象一：limma跑出来的差异基因特别少，甚至为零。原因通常是分组因子设反了，或者设计矩阵的列名和对比矩阵对不上。makeContrasts(Tumor - Normal)里的Tumor和Normal必须和colnames(design)完全一致，大小写都不能差。 解决：跑完model.matrix()后先打印colnames(design)确认，再写对比矩阵。

现象二：多探针合并后基因数对不上，或者表达值异常大。原因是聚合前没检查矩阵是否已经取log，或者探针注释包里有多重映射没处理干净。 解决：聚合前先summary(expr_matrix)看数值范围，超过100的先取log2；mapIds()的multiVals参数明确设为"first"或"list"，不要用默认值。

现象三：RobustRankAggreg的P值全部接近1，没有显著基因。原因是aggregateRanks()的N参数填错了。N是背景基因总数，填成实际基因数（比如5000）会导致P值计算偏大。 解决：人类数据统一填20000，小鼠填20000左右也适用，这个值是近似值，影响的是P值的绝对大小，不影响排序。

现象四：TCGA验证时Wilcoxon检验全部不显著。原因是FPKM矩阵没有取log，或者正常/肿瘤样本的列名识别错了。 解决：log2(FPKM + 1)是必须的；用table(grepl("Normal", colnames(tcga_expr)))确认分组列数对不对，TCGA的样本ID里-11结尾是正常，-01结尾是肿瘤。

现象五：生存曲线画不出来，或者P值和预期相反。原因是生存时间和生存状态的数据类型不对，或者高低表达分组的中位数计算包含了NA。 解决：Surv()要求time是数值型、status是0/1数值型，字符型会报错；median(merged[[gene]], na.rm = TRUE)一定要加na.rm，否则有NA时中位数返回NA，分组全错。

7. 从差异基因到hub基因：一个筛选优先级的小技巧

跑完上面整套流程，你手里会有几个列表：1210个RRA差异基因、272个生存相关基因、若干GO和PATHWAY富集结果、一个蛋白互作网络。接下来最实际的问题是：如果只能选几个基因做后续实验验证，选哪个？

我自己的习惯是做一个三层过滤。第一层，取RRA差异基因和生存相关基因的交集，这批基因既差异表达又和生存相关，临床价值最高。第二层，在这个交集里看蛋白互作网络的邻接节点数，优先选degree排名前20%的hub基因。第三层，看GO/PATHWAY富集，如果hub基因落在某个显著富集的通路里，这个基因就值得优先做。

# 三层过滤筛选候选基因 rra_genes <- read.csv("RRA_sig_genes.csv")$gene surv_genes <- read.csv("04_TCGAsurvival/survival.xlsx")$gene degree_df <- read.csv("ppi_degree.csv") # 第一层：交集 candidate <- intersect(rra_genes, surv_genes) # 第二层：hub基因 degree_df <- degree_df[degree_df$gene %in% candidate, ] hub_cutoff <- quantile(degree_df$degree, 0.8) hub_genes <- degree_df$gene[degree_df$degree >= hub_cutoff] # 第三层：通路注释 pathway_genes <- unique(unlist(strsplit( read.csv("06.PATHWAY/PATHWAY.xlsx")$members_input_overlap, ";"))) final_genes <- intersect(hub_genes, pathway_genes)

这个三层过滤的逻辑是：差异表达保证基因在疾病和正常之间有区别，生存相关保证这个区别有临床后果，hub位置保证它在调控网络里不是边缘角色，通路注释保证它有已知的生物学背景。四者都满足的基因，做实验验证的成功率会高很多。

从那以后我每次跑完GEO联合分析，都会强制走一遍这个三层过滤，而不是直接把差异基因列表丢给下游。希望帮到你。

本文还有配套的精品资源，点击获取