企业官网建设流程全解析

简介：一份关于多个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 个数据集，我把关键信息整理成了一张表，方便对照理解：

注意一个细节：选的数据集立足于不同芯片平台，这增加了数据异质性。异质性这个东西对差异表达分析是个挑战，但 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 芯片用affy；GPL17077 这类全转录组芯片用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 文件跑了一晚上才发现后面代码报错，整个人会非常崩溃。从那以后我每次跑多数据集联合分析，都强制先跑通全流程小样再做正式分析，这步检查能省下大半天调错时间。希望帮到你。

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