简介这份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而不是gastriccancerGEO的检索是基于分词匹配的连写反而可能漏掉标题里带空格的Series。检索出来之后不要急着全下。报告里做了一步过滤去掉没有重复试验的样品。所谓重复试验就是同一实验条件下有多个生物学重复。只有单样本的Series做差异分析时组内方差无法估计limma跑出来的P值不可信。常见做法是在GEO页面看Series的Sample列表每个分组至少要有3个样本才纳入。2.2 六套数据的平台与样本分布报告最终纳入的六套数据如下表。这张表本身就是选数据的决策依据值得逐列看Series平台正常样本肿瘤样本参考文献GSE64041GPL62445050Makowska et al, 2016GSE67764GPL170773030Yang et al, 2015GSE62232GPL5701081—GSE19665GPL5701010Nagae et al, 2010GSE46408GPL413366—GSE25097GPL10687243243Tung et al, 2011选平台时有几个判断GPL570Affymetrix 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 - NULLmapIds()的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后怎么整合取交集太粗暴——要求基因在六个数据集中都显著会漏掉只在部分数据集中显著的基因。报告用的是RobustRankAggregRRA它的零假设是每个基因在每个实验中的排名是随机的。如果一个基因在所有实验中都排在前列它随机出现的概率就极低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纵坐标是基因名红色logFC0绿色logFC0。# 差异基因热图 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检验过滤条件Pvalue0.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功能富集P0.05为筛选条件。DAVID的输入是基因列表背景选人类全基因组输出三个本体BP生物过程、CC细胞组分、MF分子功能。GO IDTermCountPValueGenesGO:0004784superoxide dismutase activity30.001692SOD1, NQO1, SOD3GO:0042629mast cell granule40.002400LAT, AKR1B1, PLA2G3, SRGNGO:1901687glutathione derivative biosynthetic process40.002628GSTM1, GSTM2, ESD, GSTP1看GO结果时Count是落在该Term上的差异基因数PValue是超几何检验的显著性。Count太小的Term比如只有1个基因即使P值显著生物学意义也有限。报告里superoxide dismutase activity有3个基因mast cell granule有4个这些是有解读价值的。5.2 PATHWAY富集KEGG通路与疾病关联PATHWAY富集用的是KEGG数据库报告里P0.05的通路包括Drug metabolism、Hepatocellular carcinoma、Systemic lupus erythematosus等。这里有个细节值得注意Hepatocellular carcinoma通路出现在胃癌分析里不是错误——KEGG的疾病通路往往共享基因SHC1、FZD7、GSTM1这些基因在多个癌症通路里都有角色。Pathway IDPathwayPValueMembershsa00983Drug metabolism0.019729TYMP, GSTM1, GSTP1, GSTM2hsa05225Hepatocellular carcinoma0.022852SHC1, FZD7, GSTM1, GSTM2, NQO1, GSTP1hsa05322Systemic lupus erythematosus0.029355CD28, H2AFY2, HIST1H3G, HIST3H2A, CD40解读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的先取log2mapIds()的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联合分析都会强制走一遍这个三层过滤而不是直接把差异基因列表丢给下游。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/22 13:15:55 · 阅读 1.2k 多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而不是gastriccancerGEO的检索是基于分词匹配的连写反而可能漏掉标题里带空格的Series。检索出来之后不要急着全下。报告里做了一步过滤去掉没有重复试验的样品。所谓重复试验就是同一实验条件下有多个生物学重复。只有单样本的Series做差异分析时组内方差无法估计limma跑出来的P值不可信。常见做法是在GEO页面看Series的Sample列表每个分组至少要有3个样本才纳入。2.2 六套数据的平台与样本分布报告最终纳入的六套数据如下表。这张表本身就是选数据的决策依据值得逐列看Series平台正常样本肿瘤样本参考文献GSE64041GPL62445050Makowska et al, 2016GSE67764GPL170773030Yang et al, 2015GSE62232GPL5701081—GSE19665GPL5701010Nagae et al, 2010GSE46408GPL413366—GSE25097GPL10687243243Tung et al, 2011选平台时有几个判断GPL570Affymetrix 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 - NULLmapIds()的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后怎么整合取交集太粗暴——要求基因在六个数据集中都显著会漏掉只在部分数据集中显著的基因。报告用的是RobustRankAggregRRA它的零假设是每个基因在每个实验中的排名是随机的。如果一个基因在所有实验中都排在前列它随机出现的概率就极低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纵坐标是基因名红色logFC0绿色logFC0。# 差异基因热图 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检验过滤条件Pvalue0.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功能富集P0.05为筛选条件。DAVID的输入是基因列表背景选人类全基因组输出三个本体BP生物过程、CC细胞组分、MF分子功能。GO IDTermCountPValueGenesGO:0004784superoxide dismutase activity30.001692SOD1, NQO1, SOD3GO:0042629mast cell granule40.002400LAT, AKR1B1, PLA2G3, SRGNGO:1901687glutathione derivative biosynthetic process40.002628GSTM1, GSTM2, ESD, GSTP1看GO结果时Count是落在该Term上的差异基因数PValue是超几何检验的显著性。Count太小的Term比如只有1个基因即使P值显著生物学意义也有限。报告里superoxide dismutase activity有3个基因mast cell granule有4个这些是有解读价值的。5.2 PATHWAY富集KEGG通路与疾病关联PATHWAY富集用的是KEGG数据库报告里P0.05的通路包括Drug metabolism、Hepatocellular carcinoma、Systemic lupus erythematosus等。这里有个细节值得注意Hepatocellular carcinoma通路出现在胃癌分析里不是错误——KEGG的疾病通路往往共享基因SHC1、FZD7、GSTM1这些基因在多个癌症通路里都有角色。Pathway IDPathwayPValueMembershsa00983Drug metabolism0.019729TYMP, GSTM1, GSTP1, GSTM2hsa05225Hepatocellular carcinoma0.022852SHC1, FZD7, GSTM1, GSTM2, NQO1, GSTP1hsa05322Systemic lupus erythematosus0.029355CD28, H2AFY2, HIST1H3G, HIST3H2A, CD40解读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的先取log2mapIds()的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联合分析都会强制走一遍这个三层过滤而不是直接把差异基因列表丢给下游。希望帮到你。本文还有配套的精品资源点击获取