简介一份关于多个GEO数据集联合分析的生物信息学报告PDF适合生物信息学初学者、医学数据分析人员及从事基因表达研究的科研工作者。报告以胃癌基因表达数据为例系统演示了从GEO数据库检索与筛选数据集、下载CEL文件与表达矩阵、对数转换、探针合并到使用limma包进行标准化和差异表达分析再通过RobustRankAggreg整合多实验排序结果最终识别出1210个差异基因的全过程。后续还展示了运用pheatmap绘制差异基因热图利用TCGA RNA-seq数据验证差异基因可靠性结合生存分析筛选与胃癌预后相关的272个基因以及DAVID工具进行GO功能富集分析等关键环节。资源包含1个PDF文件压缩包大小3.8MB图文并茂地呈现了数据表格、热图、生存曲线和富集结果可直接作为生信分析流程模板或课题参考资料。已有129人学习下载对希望掌握公共数据库挖掘与多数据集整合分析的读者具有较高参考价值。1. 多数据集联合分析一份能直接复现的胃癌生物信息报告做数据分析的人拿到 GEO 数据第一反应通常是跑个差异表达、画个热图就交差。但这份资源给的是另一条路把 6 套胃癌芯片数据联合起来做 RobustRankAggreg再用 TCGA 验证、接着做生存分析、GO/PATHWAY 富集、蛋白互作网络形成一套完整的数据挖掘闭环。它解决了单数据集样本量小、结论不可靠的典型痛点适合正在做生物信息课程设计、准备生信分析毕业课题、或者想系统走一遍多数据集联合分析流程的人。我看完最大的感受是这套流程每一步都有明确产出从 GEO 检索到最终生存相关基因列表全部可复现拿来改一改就能用到其它癌种上。下面按实际执行顺序把这套分析拆开讲透。2. 从 GEO 检索到数据质控六个数据集的筛选逻辑与预处理细节2.1 检索式怎么写限定物种和关键词组合GEO 检索最忌讳的是关键词太宽一搜出来几千个结果根本筛不动。这份资源里的检索式值得直接抄(gastric cancer) AND Homo sapiens[porgn:__txid9606]这个检索式里的门道在于[porgn:__txid9606]它把物种限定为人避免混入小鼠、大鼠的芯片数据。实际操作时我在 GEO 首页的搜索框里直接粘贴这个式子是能用的也可以用 GEO DataSets 的高级检索界面把 Organism 选成 Homo sapiens再把 Keywords 填 gastric cancer效果一样。筛完检索结果第一轮过滤标准是必须有重复试验样品——也就是每个分组至少要有生物学重复。这个过滤条件看着简单实际上能帮你剔除掉大量只有单样品或者没有对照组的低质量数据。那些只测了肿瘤组织没测正常组织的数据在差异表达分析里根本用不上。2.2 六个数据集的选择依据平台、样本量和文献支撑这份报告最终锁定了 6 个数据集我把关键信息整理成了一张表方便对照理解SeriesPlatformNormalTumor文献来源GSE64041GPL62445050Makowska et al, 2016GSE67764GPL170773030Yang et al, 2015GSE62232GPL5708181—GSE19665GPL5701010Nagae et al, 2010GSE46408GPL413366—GSE25097GPL10687243243Tung et al, 2011注意一个细节选的数据集立足于不同芯片平台这增加了数据异质性。异质性这个东西对差异表达分析是个挑战但 RobustRankAggreg 这种基于排序整合的方法恰恰能处理这种异质性因为每个数据集内部跑完 limma 后把基因按 logFC 排个序就行不直接比较表达值的大小。实际下载时有两种文件需要区分如果 GEO 页面提供了处理好的表达矩阵series matrix直接下载如果只有 CEL 原始文件就用 affy 包读取。两种文件对应完全不同的下游处理路径这个选错的话后面会绕很大一个弯子。2.3 表达矩阵的预处理取 log 和探针映射拿到表达矩阵之后第一步先看数据是否需要取 log。判断方法很简单如果矩阵里最大表达值在 1000 以上基本就是没取过 log 的芯片荧光信号强度值需要做 log2 转换。如果值域在 0-15 上下浮动说明已经取过 log 了直接往下走。expr_mat - read.table(GSE64041_series_matrix.txt, header TRUE, row.names 1, sep \t) # 判断是否需要取 log最大值远大于 50 说明是原始信号值 if (max(expr_mat, na.rm TRUE) 50) { expr_mat - log2(expr_mat 1) }这里用log2(x 1)而不是log2(x)是为了避免表达值为 0 时算出负无穷。加 1 是芯片数据分析里的常见做法虽然对高表达基因几乎没有影响但能保证低表达基因的数据仍然可用。探针到基因名的映射是另一个必踩的坑。一个基因往往对应多个探针如果直接把探针水平的数据拿去分析同一个基因会出现好几行后面 limma 跑完结果没法解释。处理逻辑很直接按基因名分组取每个组内所有探针表达值的平均值。# 假设 probe_map 是探针ID到基因Symbol的映射表 expr_mat$gene_symbol - probe_map$Symbol[match(rownames(expr_mat), probe_map$ProbeID)] expr_mat - na.omit(expr_mat) expr_agg - aggregate(. ~ gene_symbol, data expr_mat, FUN mean) rownames(expr_agg) - expr_agg$gene_symbol expr_agg - expr_agg[, -1]aggregate函数在这里的作用就是按基因分组求平均这样每个基因只剩一行表达值。几个数据集的表达矩阵都做完这步处理之后存成统一的格式为下一步 limma 分析做准备。这些代码在这份资源里都有对应的 R 脚本照着执行就能把原始数据变成可分析矩阵。3. limma 差异表达与 RobustRankAggreg 整合从单数据集到跨平台共识3.1 limma 分析的正确姿势分组矩阵和 contrast 设置limma 包做芯片差异表达分析核心逻辑是拟合线性模型再用经验贝叶斯方法修正方差。每个数据集单独跑别把数据混在一起因为芯片平台不同、批次不同、标准化基准不同直接合并会让平台效应淹没生物学差异。library(limma) # 假设 expr_agg 是处理好的表达矩阵列名是样本ID # group 向量标注每个样本属于 Normal 还是 Tumor group - factor(c(rep(Normal, 50), rep(Tumor, 50)), levels c(Normal, Tumor)) design - model.matrix(~ group) fit - lmFit(expr_agg, design) fit - eBayes(fit) result - topTable(fit, coef 2, number Inf, adjust.method BH)这里model.matrix(~ group)会生成两列设计矩阵第一列是截距第二列是 Tumor vs Normal 的对比。topTable里的coef 2就是取第二列系数对应肿瘤和正常的差异。每个数据集跑完之后保存一份结果到01_limma目录文件里至少要有 gene symbol、logFC、AveExpr、t 值、P 值、校正后 P 值这几列后面 RRA 分析要用 logFC 来排序。3.2 RobustRankAggreg 的原理和参数设定为什么排名整合比直接合并更稳RobustRankAggreg 的零假设是每个基因在每个实验中的排名是随机排列的。如果某个基因在所有实验里都稳定排在前面或后面它偏离随机程度的显著性就高是差异基因的可能性就大。这个思路避开了直接合并表达值时的平台校准难题。library(RobustRankAggreg) # rank_list 是一个列表每个元素是单个数据集中按 logFC 降序排列的基因名向量 # 上调基因取 logFC 降序排列下调基因取 logFC 升序排列 rank_list - list( GSE64041 rownames(result_GSE64041)[order(result_GSE64041$logFC, decreasing TRUE)], GSE67764 rownames(result_GSE67764)[order(result_GSE67764$logFC, decreasing TRUE)] ) # 六个数据集都放进去 aggregated - aggregateRanks(rank_list, method RRA)aggregateRanks的输出里有Score和P.value两列P 值经过 Bonferroni 校正后过滤条件用 0.05。这份资源里跑了出 1210 个差异基因436 个上调、774 个下调。这里要特别强调RobustRankAggreg 用的是每个数据集内部基因排名的相对顺序不要求不同芯片平台的表达值具有可比性这是它处理跨平台数据时最省心的原因。相比之下如果强行把所有数据集的表达值扔进一个 limma 模型芯片批次效应处理不好结果会很难看。# 取显著差异基因 sig_genes - aggregated[aggregated$P.value 0.05, ] # 区分上下调查看该基因在多数数据集中的 logFC 方向 up_genes - rownames(sig_genes)[sig_genes$Score 0] down_genes - rownames(sig_genes)[sig_genes$Score 0]3.3 pheatmap 热图的绘制细节只画头尾 20 个基因拿到 1210 个差异基因之后热图只挑了最上调和最下调的各 20 个基因来画。这个做法很聪明如果全画出来热图面积太大看不出表达模式只画极端变化的基因能一眼看出这些基因在几乎所有数据集里 logFC 方向一致。library(pheatmap) # logFC_matrix 是这 40 个基因在 6 个数据集中的 logFC 值矩阵 pheatmap(logFC_matrix, cluster_cols FALSE, color colorRampPalette(c(green, black, red))(100), display_numbers TRUE, number_format %.1f, fontsize_row 8)注意cluster_cols FALSE这里不需要对数据集聚类因为希望横轴保持数据集顺序直观对比同一个基因在不同数据集里的 logFC。颜色设置上绿-黑-红三段渐变能清楚区分负值、零值和正值。4. TCGA 验证与生存分析差异基因的临床价值闭环4.1 TCGA 数据下载和 FPKM 矩阵合并Perl 处理单样品文件的思路TCGA 下载下来的 RNA-seq 数据是一个样品一个 FPKM 文件几十个正常组织加几百个肿瘤组织就是几百个文件不可能手动一个个合并。资源里用了 Perl 来批量合并思路是把每个文件的文件名映射成样品名读文件、抽基因表达值、写进一个大矩阵。#!/usr/bin/perl use strict; use warnings; my files glob(*.txt); # 每个文件是一个样品的 FPKM 数据 my %matrix; foreach my $file (files) { open(my $fh, , $file) or die cannot open $file\n; my $sample $file; $sample ~ s/\.txt$//; # 去掉扩展名作为样品名 while (my $line $fh) { chomp $line; my ($gene, $fpkm) split /\t/, $line; $matrix{$gene}{$sample} $fpkm; } close($fh); } # 输出矩阵行是基因列是样品 my samples sort keys %{ $matrix{(keys %matrix)[0]} }; print gene\t, join(\t, samples), \n; foreach my $gene (sort keys %matrix) { print $gene; foreach my $sample (samples) { print \t, $matrix{$gene}{$sample} // NA; } print \n; }这个脚本的核心是嵌套哈希外层键是基因内层键是样品名值是对应的 FPKM。这样做一次遍历就能构建出完整矩阵不用反复横跳文件。合并前建议先确认每个文件的基因列表是一致的别出现不同样品基因版本不一样的情况那会导致矩阵出现大量 NA。4.2 Wilcoxon 检验验证差异基因非参数检验的选择理由TCGA 验证这一步的逻辑很简单GEO 里筛出来的差异基因在 TCGA 的独立队列里是不是也显著差异表达用的是 Wilcoxon 非参数检验而不是 t 检验因为 TCGA 的 FPKM 数据分布偏态严重不满足正态性假设非参数检验对分布没有要求更稳健。# expr_tcga 是合并好的 TCGA 表达矩阵group_tcga 标注 Normal/Tumor p_values - apply(expr_tcga[genes_to_validate, ], 1, function(x) { wilcox.test(x ~ group_tcga)$p.value }) adjusted_p - p.adjust(p_values, method BH)这份资源的验证结果里LYPD1、NQO1、MKI67、B3GNT5 这些基因在 TCGA 队列里仍然显著说明不是单个数据集偶然跑出来的假阳性。这一步的产出直接给差异基因增加了可信度。4.3 survival 包生存分析和 logrank 检验批量跑 1210 个基因的效率问题生存分析的输入有三个每个样品的生存时间、生存状态活着还是死亡、每个基因的表达值。分析目标是把 1210 个差异基因全部跑一遍 logrank 检验筛出其中与生存显著相关的基因P 0.01最终找到 272 个。library(survival) # surv_data 包含 time、status 和所有差异基因的表达值 results - data.frame(gene character(), pvalue numeric()) for (gene in rownames(expr_matrix)) { # 按表达值中位数把样本分成高表达组和低表达组 group - ifelse(expr_matrix[gene, ] median(expr_matrix[gene, ]), High, Low) fit - survdiff(Surv(time, status) ~ group, data surv_data) pval - 1 - pchisq(fit$chisq, df 1) results - rbind(results, data.frame(gene gene, pvalue pval)) }这里有个效率上的细节1210 个基因循环跑R 里这个操作要几分钟如果基因再多建议改用apply系列函数并行化。分组用中位数是默认做法当然也可以用上下三分位数但会损失一部分样本一般不推荐。# 筛选生存相关差异基因 survival_genes - results[results$pvalue 0.01, ] # 绘制单个基因的生存曲线 library(survminer) fit - survfit(Surv(time, status) ~ group, data surv_data) ggsurvplot(fit, pval TRUE, risk.table TRUE)TRIP13、TMEM106C、DTYMK、BIRC5、CDK1 这些基因的 P 值都跑到 1e-5 以下是典型的生存相关标志物候选。5. 常见问题与避坑实录六个数据集联合分析最容易翻车的五个地方5.1 探针注释不一致导致基因名对不齐现象合并多个数据集结果时发现同一基因在不同数据集里名字大小写不一致有的数据集里叫CD36另一个里变成了Cd36还有的带了版本号后缀直接导致 merge 之后基因数骤减。原因不同芯片平台的注释文件风格不同Affymetrix 的注释文件里基因 Symbol 通常是全大写但 Illumina 平台的注释可能首字母大写或全小写。解决在上游统一转换用toupper()把所有基因名转为大写再合并。同时在映射之前跑一遍去重检查看是否有同义名必要时用biomaRt包统一转换到 Entrez ID 再映射回 Symbol这个方案最稳。5.2 Series Matrix 文件读取后出现大段空值或乱码现象read.table读 series matrix 文件后表达矩阵中间出现 NULL 或者不均匀的行列数下游 limma 直接报错。原因GEO 的 series matrix 文件头部有大量以!开头的元信息行直接读会把元信息当成数据。解决先readLines找到!series_matrix_table_begin行号从下一行开始读或者用GEOquery包的getGEO()解析不需要手工处理表头。5.3 limma 结果的 logFC 方向在不同数据集里不一致现象同一个基因在 GSE64041 里是上调在 GSE67764 里变成下调RobustRankAggreg 整合后 P 值不显著基因被过滤掉。原因数据集的 grouping 定义不一致有的数据集里对照组在前有的实验组在前。芯片平台不同也会导致探针杂交信号的方向性差异。解决每个数据集跑 limma 之前先table(group)确认分组样本数和顺序再检查已知 marker 基因比如胃癌里的 CD36的 logFC 方向是否符合预期确认方向一致后再继续。5.4 TCGA 生存数据和表达数据的样品名对不上现象生存数据里的样品条形码是TCGA-BR-1234-01A表达矩阵里的列名可能是TCGA-BR-1234-01左对齐四个位置对不上merge 之后大量样本丢失。原因TCGA 不同数据文件输出的样品命名规范不完全一致有的带样本类型编码01A有的被截断。解决统一取前 12 位或前 15 位字符做匹配。一般建议取第 1 到第 15 位字符包含TCGA-XX-1234-01稳定且唯一标识样品。5.5 CEL 文件用错读取包导致表达量差异巨大现象同一套 CEL 文件用affy包的justRMA()和用oligo包的rma()跑出来的表达值分布差异肉眼可见同一个基因排序位置变化很大。原因不同芯片平台对应的 CEL 文件格式不同affy针对 Affymetrix 3 端表达芯片设计oligo更适合外显子芯片和基因芯片。混用会踩到格式解析问题。解决先查清楚 GSE 页面标注的 platform 类型。GPL570 这类老式 3 IVT 芯片用affyGPL17077 这类全转录组芯片用oligo。不确定的时候两个包都试跑一遍比较表达矩阵的箱线图分布选中间值在 6-9 之间、分布更紧凑的。6. GO 与 PATHWAY 富集、蛋白互作网络差异基因列表的生物学功能落地6.1 DAVID 做 GO 富集背景基因怎么选才不被喷GO 富集分析的统计原理是超几何分布核心参数有两个待分析基因列表和背景基因列表。背景基因选错富集结果全军覆没。常见错误是背景基因用了全基因组但实际上芯片只测了部分基因这会导致富集到一些假阳性通路。正确做法是用你自己的芯片或测序平台能检测到的所有基因作为背景。DAVID 操作上先把差异基因列表复制进去选择官方基因名种类背景选 Homo sapiens 全基因组P 值阈值设 0.05。这份资源的 GO 结果里superoxide dismutase activity、mast cell granule、glutathione derivative biosynthetic process 这几个 GO term 显著富集说明差异基因集中在氧化应激和谷胱甘肽代谢相关功能模块。这是胃癌发生机制里非常经典的通路。# 用 R 做 GO 富集的可替代方案clusterProfiler library(clusterProfiler) library(org.Hs.eg.db) # gene_list 是差异基因的 Entrez ID 向量 ego - enrichGO(gene gene_list, OrgDb org.Hs.eg.db, ont BP, pAdjustMethod BH, pvalueCutoff 0.05, qvalueCutoff 0.05, readable TRUE) barplot(ego, showCategory 15)6.2 PATHWAY 富集结果解读KEGG 富集和 GO 的结果互为印证PATHWAY 富集用的也是 DAVID 或者 clusterProfiler 里的 KEGG 通路数据库。胃癌差异基因最显著富集到 Drug metabolism - other enzymes、Hepatocellular carcinoma、Systemic lupus erythematosus、Cell adhesion molecules这几个通路里能同时看到 GSTM1、GSTP1、GSTM2、NQO1 这些谷胱甘肽转移酶家族的基因。和 GO 富集结果互相印证说明这些基因确实参与解毒代谢过程而不是随机散落的几个基因碰巧富集到一起。ekegg - enrichKEGG(gene gene_list, organism hsa, pvalueCutoff 0.05) dotplot(ekegg)这种 GO 和 KEGG 双维度验证的习惯值得保留如果 GO 富集到一组通路但 KEGG 结果完全不沾边先别急着写结论回头检查基因列表是不是有注释误差。6.3 String 构建蛋白互作网络核心节点的判定标准用 String 数据库跑蛋白互作网络输入的是 272 个生存相关差异基因输出一个 PPI 网络图。然后做一步关键转化用 R 读取 String 导出的 TSV 文件统计每个蛋白的邻接节点数目这个数就是 network degree度越高的蛋白在网络里越核心。string_data - read.table(string_interactions.tsv, header TRUE, sep \t) # 统计每个蛋白的互作伙伴数量 library(dplyr) degree - string_data %% group_by(node1) %% summarise(degree n()) %% arrange(desc(degree))这份资源跑完的结果显示 CFTR、SST、TIMP1 这些基因排在核心位置说明它们与大量生存相关基因有直接互作是网络中的关键枢纽节点更值得优先做实验验证。做这种网络分析有个值得坚持的习惯把互作置信度从默认的 0.4 调高到 0.7 再跑一次。低置信度的互作关系往往是文本挖掘出的共现关系没有直接实验证据保留下来只会把网络弄得很密核心节点反而不好识别。高置信度过滤后还能保留下来的互作关系基本都有实验证据支撑可信度完全不同。从实操角度讲这份资源对做复现的人最有用的是整套分析逻辑的顺序编排先 GEO 多数据集筛差异基因再 TCGA 独立队列验证接着用生存分析挑临床相关基因最后用富集分析和 PPI 网络解释生物学意义。每一步的产出都是下一步的输入链条非常清楚。真到自己拿着这套流程改到别的癌种上时有个经验值得分享数据下载完先做一个全流程小样测试拿 3 个基因跑通从 limma 到生存分析的完整代码确认代码没有路径错误和环境问题再放开跑全部基因。否则等你下载完几百个 CEL 文件跑了一晚上才发现后面代码报错整个人会非常崩溃。从那以后我每次跑多数据集联合分析都强制先跑通全流程小样再做正式分析这步检查能省下大半天调错时间。希望帮到你。本文还有配套的精品资源点击获取