企业官网建设流程全解析

简介：一套面向生物信息学与医学科研人员的R语言GSEA绘图资源包，聚焦多组基因集富集分析的可视化输出，帮助用户快速生成符合SCI论文要求的富集曲线、通路气泡图与统计表格。压缩包共513个文件，容量41.63MB，主体包括293个PNG富集图片、110个XLS结果表单、102个HTML富集报告，以及R脚本、RNK排序文件、GMT基因集文件、CLS表型文件等输入与配置材料，几乎覆盖GSEA分析全流程所需的标准文件类型。已有423人学习使用，适用于需要复现多组GSEA富集结果、并将高分辨率图表用于论文发表的科研场景。资源内提供了可直接修改运行的R代码，基于GSEAPreranked、fgsea或clusterProfiler等常用库，结合gene_sets.gmt与cls表型文件可快速适配自身表达矩阵；同时PNG图片和HTML报告便于逐一核对KEGG等通路的富集方向、FDR值与核心基因，且内容涵盖MAPK、神经活性配体受体等经典通路，为数据解读和论文组图节省大量时间。

1. 多GSEA富集图：先搞明白你要画的到底是哪一张图

经常有同学拿着别人SCI论文里的“多GSEA富集图”来问：R 语言能不能一笔生成这种图，输入文件到底长什么样。答案是可以，但前提是先把“多GSEA”这个词拆清楚——它可能是同一个通路在多个对比组里的富集曲线拼图，也可能是多个显著通路按NES排序后的气泡图，还可能是一条通路在多个基因集数据库里的并排比较。三种图输入文件不一样，代码也不一样。这套资源的核心价值不在“会画图”，而在“输入文件怎么整理才能让GSEA结果有意义”。本文按我自己的跑通路径，把原理、输入文件、R代码、排查顺序完整过一遍，新手能直接照跑，熟手可以只看参数和踩坑部分。

2. GSEA 原理十行讲透：排序列表、富集分数与置换检验的关系

2.1 为什么 GSEA 不比“显著基因列表”，而比“全部基因的排序”

做GO/KEGG富集通常第一步是筛差异基因，然后拿差异基因去查通路，这叫 ORA（Over-Representation Analysis）。它的缺陷很明显：阈值一卡，排在 0.051 的基因就被丢掉了，而通路层面的协同变化往往就藏在这些“不显著但方向一致”的基因里。GSEA 的立场是：我不筛基因，我把所有基因按某个指标排个序，然后看某个基因集的成员在排序里是偏前还是偏后、是聚在一起还是散落各处。

具体到实现，排序指标最常用的是 log2FC 或 signal-to-noise。你把全基因组按这个值从大到小排好，然后从头到尾遍历一次，每碰到一个属于目标基因集的基因就加一点分数，碰不到就扣一点分数，加多少扣多少由表达量高低决定，遍历过程中记下最大偏离零点的地方，那个分数就是 ES（Enrichment Score）。为了消除基因集大小的影响，再除以一个归一化系数得到 NES，最后用置换检验算 p 值。

这里有个很多人忽略的点：GSEA 输出的是“整条通路的协同变化方向”。NES 为正表示这条通路整体在实验组里偏激活，为负就是整体抑制。所以多GSEA图里最常见的配色，就是红色表示激活、蓝色表示抑制，这不是审美选择，而是 NES 符号的直接映射。画图之前一定要确认自己的 geneList 排序方向是“实验组比对照组”，否则整张图的激活抑制方向全是反的。

2.2 多GSEA富集图的三种常见形态与选型

多GSEA图在SCI里常见三种形态，适用场景完全不同，做之前先选好，不然画到一半才发现图表达不了你的结论。第一种是经典富集图，也就是单条通路的绿色ES曲线加下面基因排序的条形热力区，用gseaplot2()一页拼多条通路，适合展示“某条通路在我关心的分组里显著激活”这种强结论。第二种是多通路气泡图：纵轴放通路名，横轴放分组，点大小表示基因集大小，颜色从蓝到红表示 NES 正负，适合做全局概览，一张图塞下几十条通路。第三种是同一通路在多个对比组里的 NES 条形图或热图，适合做“时间序列/多组织”的纵向比较。

选型的标准就一条：看你的核心结论是“某条通路在组间差异里排第一”还是“这批通路整体协调变化”。前者用经典富集图拼版，后者用气泡图。如果只是为了附录展示所有结果，气泡图信息密度最高，审稿人也最容易快速捕捉。三种形态的输入数据其实同源，都是 GSEA 跑完后的结果表，差别只在绘图入口。

2.3 置换检验、多次运行与 set.seed 的关系

GSEA 的 p 值来自置换检验：把你基因名打乱，重新算一遍 ES，重复很多次得到一个“零分布”，再看你真实 ES 落在分布哪个位置。既然是随机置换，就存在每次运行结果不完全一致的问题。这不是 BUG，也不是代码错了，而是统计方法本身的性质。

正因为如此，在跑多组对比时，循环里每一组都必须在同一个set.seed()下运行，否则同一份数据你今天跑和明天跑，富集到通路的排序会略有浮动，虽然显著性结论一般不变，但“可复现性”这事在投稿时一旦被问到就很尴尬。我的习惯是固定set.seed(2024)这样每次跑出来完全一致。还有一个更稳妥的做法：置换次数足够多时结果就会稳定下来，但现在主流 R 包里置换次数已经内置固定，你能控制的只有随机种子，关于这一点后面的排查章节还会单独说。

3. 一张图要三份输入：表达谱、基因集与分组表的整理逻辑

3.1 从表达谱到 geneList：核心是“带名字的排序向量”

R 做 GSEA 的第一步不是读文件，而是构造一个命名向量geneList：每一个基因名对应一个排序指标值，向量按值降序排好。很多人卡在这里，因为不清楚该拿 counts 还是 TPM。我的经验是：手头有 DESeq2 或 edgeR 的差异结果时，直接用它的log2FoldChange列，这是最稳的；只有表达矩阵时，就用“实验组均值减对照组均值”作为排序值，也能跑，只是灵敏度略低。

# 读取表达矩阵：行=基因，列=样本 expr <- read.table("expr_matrix.txt", header = TRUE, row.names = 1) # 读取分组信息：两列，sample 和 condition group <- read.table("group.txt", header = TRUE) # 构造排序向量：以两组均值的差为排序指标 logFC <- rowMeans(expr[, group$sample[group$condition == "treat"]]) - rowMeans(expr[, group$sample[group$condition == "control"]]) deg <- data.frame(gene = rownames(expr), logFC = logFC) # 关键一步：去掉重复基因名，否则后续报错 deg <- deg[!duplicated(deg$gene), ] # 再丢掉缺失值，NA 会直接让 gsea 函数中断 deg <- deg[is.finite(deg$logFC), ] geneList <- deg$logFC names(geneList) <- deg$gene geneList <- sort(geneList, decreasing = TRUE) head(geneList)

这段代码的逻辑是先把表达矩阵按分组拆开，用rowMeans计算两组均值差作为排序指标，然后去重、去缺失、按降序排列。去重这一步最容易翻车：同一个基因名在矩阵里出现两次，names()就有重复，clusterProfiler 内部做匹配时会直接下标越界。参数上唯一需要你根据数据改的是rowMeans的分组条件，其余都是通用写法。如果用的是 DESeq2 结果，直接把logFC列替换成res$log2FoldChange，基因名列换成rownames(res)就行。

3.2 基因集格式与 ID 转换：symbol 和 ENTREZID 对不上就全白跑

GSEA 的基因集文件常见有两种：GMT 格式和 R 包内置的基因集。GMT 是纯文本，第一列通路名，第二列描述，第三列之后全是基因。但 R 里的gseKEGG()和gseGO()函数并不直接吃 GMT 文件，它们用的是 KEGG 和 GO 的最新注释，你只需要提供基因 ID。这里有个巨大的坑：gseKEGG()默认要求基因是 ENTREZID，而你手上十有八九是 gene symbol，不转换直接喂进去，返回结果全是空集。

library(clusterProfiler) library(org.Hs.eg.db) # 把 geneList 里的 symbol 转成 ENTREZID convert <- bitr(names(geneList), fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db) # 转换后可能丢基因，用保留下的基因重新做命名向量 geneList_entrez <- geneList[convert$SYMBOL] names(geneList_entrez) <- convert$ENTREZID # 排序向量要求值降序，必须再排一次 geneList_entrez <- sort(geneList_entrez, decreasing = TRUE)

bitr的参数逻辑很直观：fromType是现状，toType是目标，OrgDb按物种选，人类就是org.Hs.eg.db，小鼠是org.Mm.eg.db。转换过程会丢弃一部分没有对应关系的基因，这是正常的，别看到基因数变少就慌。真正需要警惕的是转换后names(geneList_entrez)出现重复，因为多个 symbol 可能映射到同一个 ENTREZID，此时务必按 3.1 的去重逻辑再处理一次。另外，不同物种的 KEGG 缩写不一样，人类是hsa，小鼠是mmu，大鼠是rno，写错会直接报错或者富集出别的物种通路。

3.3 分组表与多组对比：多GSEA图的数据结构设计

多GSEA和多组对比的数据组织方式，决定了你后期画图的顺手程度。我的做法是：一开始就不做“一个对比一个脚本”，而是把所有对比结果合成一张长表，每一行对应“某个分组对比 + 某条通路 + 它的 NES/p 值”，靠一列contrast区分来源。这样画气泡图、拼多图、筛选显著通路都在同一张表上操作。

# 以三次对比为例，构造一个所有结果的长表 group_info <- data.frame( contrast = c("treat_vs_ctrl", "drug_vs_ctrl", "drug_vs_treat"), treat = c("treat", "drug", "drug"), control = c("ctrl", "ctrl", "treat") )

这张表本身不参与 GSEA，它只是循环的驱动。真正跑分析时，每次从一个总表达矩阵里按当前对比的两个条件取子集，重新算均值差，生成新的geneList，跑完把结果加上contrast列存起来。这样做的好处是等所有对比跑完，你手里的结果表天然就是“多”GSEA 的数据基础，不需要后期做各种rbind大缝合，也不容易因为列名不一致报错。

4. 从单组到多组：R 跑 GSEA 的完整代码与关键参数

4.1 单次 GSEA 的最小完整代码：先跑通再扩展

无论最终要画多复杂的图，第一步一定是先把单次 GSEA 跑通。clusteProfiler 的gseKEGG()是最常用的入口，因为 KEGG 通路数据库覆盖面广、注释稳定。这里的输入就是 3.2 里构造的geneList_entrez，注意必须是命名向量、数值型、已按降序排好，三个条件缺一不可。

library(clusterProfiler) library(enrichplot) set.seed(2024) # 固定随机种子，保证结果可复现 kk <- gseKEGG( geneList = geneList_entrez, organism = "hsa", minGSSize = 10, maxGSSize = 500, pvalueCutoff = 0.05, seed = TRUE ) # 查看最核心的几列 result_df <- as.data.frame(kk@result) head(result_df[, c("ID", "Description", "setSize", "NES", "pvalue", "p.adjust")])

minGSSize和maxGSSize控制纳入分析的基因集大小范围。minGSSize=10表示少于 10 个基因的通路直接跳过，能有效过滤那些太小、经不起置换检验的碎片通路；maxGSSize=500则是为了过滤那些超大且宽泛的通路，比如“代谢通路”这种动辄几百上千基因的集合，富集结果很难解读。seed = TRUE和前面set.seed()是双保险，保证置换检验可复现。跑出来的kk是一个对象，真正的结果在kk@result里，后续画图直接从kk@result取数据即可。

4.2 多组对比循环：一个函数解决所有 Contrast

多GSEA图如果是一个个对比手动跑，不仅慢，而且容易在某个环节混入不一致的参数。我习惯写一个小的循环函数，把 3.3 里的group_info和总表达矩阵作为输入，一次输出所有对比的结果。

run_multi_gsea <- function(expr, group_df, group_info, organism = "hsa") { results <- list() for (i in seq_len(nrow(group_info))) { ct <- group_info$contrast[i] treat_cond <- group_info$treat[i] ctrl_cond <- group_info$control[i] # 按当前对比取样本 treat_samples <- group_df$sample[group_df$condition == treat_cond] ctrl_samples <- group_df$sample[group_df$condition == ctrl_cond] # 构造 geneList logfc <- rowMeans(expr[, treat_samples]) - rowMeans(expr[, ctrl_samples]) deg <- data.frame(gene = rownames(expr), logfc = logfc) deg <- deg[!duplicated(deg$gene) & is.finite(deg$logfc), ] gl <- deg$logfc names(gl) <- deg$gene gl <- sort(gl, decreasing = TRUE) # ID 转换 convert <- bitr(names(gl), fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db) gl_entrez <- gl[convert$SYMBOL] names(gl_entrez) <- convert$ENTREZID gl_entrez <- sort(gl_entrez, decreasing = TRUE) # 跑 GSEA，注意每个对比都在同一个随机种子下运行 set.seed(2024) kk <- gseKEGG(geneList = gl_entrez, organism = organism, minGSSize = 10, maxGSSize = 500, pvalueCutoff = 0.05) tmp <- as.data.frame(kk@result) tmp$contrast <- ct results[[ct]] <- tmp } do.call(rbind, results) } res_all <- run_multi_gsea(expr = expr, group_df = group, group_info = group_info, organism = "hsa")

这个函数里最容易忽略的是符号方向：treat_samples在rowMeans里写在前面，所以 NES 为正表示相对于对照组，当前处理组这条通路是激活的。如果你某个对比想表达的是“抑制”，那只要在画图时把contrast的顺序调整一下就行，不必反向重跑数据。另一个容易忽略的是半径：每次循环里重新生成gl后，bitr可能丢弃部分基因，后续 KEGG 富集的对象是转换后剩下的基因，所以如果你发现结果里setSize偏小，别怀疑代码，去查原始矩阵里基因名是不是混了版本号之类的脏字符。

4.3 画多GSEA图的两个出口：拼版富集图与多通路气泡图

跑完多组对比后，就到了真正“出图”的环节。如果核心结论是“某条通路在某个对比下显著激活”，那就用gseaplot2画经典富集图并横向拼版；如果核心结论是“多个通路的整体变化模式”，那就用 ggplot 画气泡图。两条路径都需要先把结果表过滤出显著通路。

# 路径1：经典富集图拼版（展示特定通路） kk_sig <- kk # 用单次跑出来的结果对象 p1 <- gseaplot2(kk_sig, geneSetID = c("hsa04110", "hsa04115", "hsa04218"), pvalue_table = TRUE, title = "Cell Cycle related pathways") ggsave("gsea_classic_panel.pdf", p1, width = 9, height = 7, dpi = 300)

# 路径2：多通路气泡图（展示所有对比的全局模式） library(ggplot2) sig <- res_all[res_all$p.adjust < 0.05, ] # 每个对比取NES绝对值最大的前15条通路 top <- do.call(rbind, lapply(unique(sig$contrast), function(ct) { sub <- sig[sig$contrast == ct, ] sub[order(abs(sub$NES), decreasing = TRUE), ][1:15, ] })) p2 <- ggplot(top, aes(x = contrast, y = Description, size = setSize, color = NES)) + geom_point(alpha = 0.8) + scale_color_gradient2(low = "#2166AC", mid = "#F7F7F7", high = "#B2182B", midpoint = 0) + scale_size(range = c(2, 8)) + labs(x = "Comparison", y = "Pathway", color = "NES", size = "Gene count") + theme_bw(base_size = 11) + theme(axis.text.y = element_text(size = 9)) ggsave("gsea_bubble_multi.pdf", p2, width = 8, height = 10, dpi = 300)

气泡图的参数里最需要调的是scale_size(range = c(2, 8))，如果通路setSize从 20 到 400 跨度很大，默认的点大小会让小通路几乎看不见，建议先看range(top$setSize)再决定范围。颜色方面midpoint = 0保证 NES 正负配色对称。alpha = 0.8是为了让重叠的点透出底下的信息，多个对比之间的通路名有重叠时非常有用。最后提醒一句：gseaplot2里的geneSetID填的是结果表的行号或 ID 号，不是通路名的字符串，填错会直接空白图。

5. 常见问题排查：五个让 GSEA 结果翻车的具体场景

5.1 现象：每次跑同样数据，结果不完全一样

这是我最常被问的问题。现象是同一份数据两次运行，p 值小数点后几位有变化，偶尔通路排序有微小浮动。原因就是 2.3 里说的置换检验的随机性，多个对比循环里没有固定set.seed()，或者每次单独跑脚本时种子不同。解决方法是：在跑 GSEA 的代码最开头统一加set.seed(2024)，如果用了多个包并行，还要注意并行核心数也会影响随机序列。固定种子后如果两次运行仍不一致，检查代码里是不是有排序函数用到了order()但有个别基因是NA。

5.2 现象：报错 “subscript out of bounds”

报错场景几乎全在gseKEGG()那一行。原因是geneList的names()里有重复，或者包含空字符串。重复基因名出现过一次，内部做排序索引时就会越界。解决：在构造geneList后跑一行any(duplicated(names(geneList)))检查，返回 TRUE 就回去看表达矩阵。矩阵第一列如果是基因名，但某些行是“---”或者空，去重删掉即可。另一个隐藏来源是bitr转换后多个 symbol 映射到同一 ENTREZID，同样要重新去重。

5.3 现象：富集结果全是“没有显著通路”

现象是跑完head(result_df)里 p 值全是 1，一条显著通路都没有，但直觉上数据有明显差异。原因最常见是geneList排序方向反了，把所有值取反后 GSEA 的富集方向完全颠倒，原本激活的通路变成抑制，再遇上通路本身双向都被扰动，p 值就被稀释了。另一个常见原因是pvalueCutoff卡太严，默认 0.05 对某些数据本身就苛刻。解决：先别改阈值，把pvalueCutoff临时设为 1，看整体 p 值分布，如果全部集中在 0.4 以上，说明数据本身没有通路级别的协调信号，问题在表达矩阵或分组标签上，而不在统计阈值。

5.4 现象：画出来的通路名是 “hsa04110” 而不是可读名称

gseKEGG的结果里ID列是hsa04110这种编号，Description列才是 “Cell cycle” 这类名称。很多人在gseaplot2里直接传geneSetID = c("hsa04110")时发现图是能出，但标题或图例显示编号，审稿人看着不舒服。解决：如果最终只展示几条关键通路，直接在result_df里用Description过滤，再用过滤出的行号传给geneSetID；如果画气泡图，直接用Description作 y 轴，不存在这个问题。还有个小坑：GO 分析的ID形如GO:0007049，画图前建议把冒号去掉或换成下划线，某些 PDF 设备对冒号渲染不友好。

5.5 现象：绘图正常但富集方向和文献完全相反

现象是同样的数据集、同样的通路，你的 NES 是正，文献里是负。原因不是代码，是排序指标的符号定义。你用treat - ctrl，文献用的是ctrl - treat，两个结果的 NES 必然互为相反数。解决：画图之前先用一个已知通路做 sanity check，比如找到你数据里一条公认激活的通路，确认它的 NES 符号符合实验方向。这个检查成本极低，却能在画图阶段避免整张图返工。另外注意 DESeq2 的log2FoldChange自带方向，不要自己在外面再乘一个负号。

6. 出图前的最后一步：SCI 版式、字体、配色与我的固定习惯

6.1 一页多图的版式控制：patchwork 优于手动拼图

多GSEA图做拼版时，我早期用gridExtra::grid.arrange()，后来换成patchwork，原因在于patchwork控制图例和间距更灵活。比如三张经典富集图竖向叠放，需要共享绘图区但保留各自的小标题，patchwork用一行代码就能对齐。横向拼多组对比时，注意保持每个子图的坐标尺度一致，否则通路名长度不同会导致图面积参差不齐。如果画的是气泡图长表，ggsave的高度按通路条数估算，10 条通路给 6 英寸，多 5 条加 1 英寸，这个经验值比固定高度靠谱。

6.2 字体与配色的固定参数：直接用能减少一轮修改

期刊投稿第一阶段最烦的不是分析，而是字体不统一被打回。我的固定做法是在theme()里统一base_family = "Arial"，PDF 输出后商用字体和 Times 混排的问题就不会出现。配色上，气泡图我习惯用RColorBrewer的RdYlBu反转版，NES 正负分明，色盲友好。如果期刊要求灰度打印，那图例就不能只靠颜色区分激活和抑制，必须在通路名后加“(+) / (−)”标记。这些细节我会提前写进一个theme_sci()函数里，全篇文章所有图都从这一个函数输出，保证字级、线宽、图例位置完全一致。

6.3 输出格式与最后检查：PDF 矢量图优先

多GSEA图如果字符密集，比如气泡图有十多个通路名，导出 300dpi 的 TIFF 容易发虚。我一般用ggsave(..., device = "pdf")输出矢量稿，然后交给期刊系统转位图。ggsave里设width和height用英寸，不用像素。提交前最后一轮检查我会做三件事：第一，通路名没有截断，特别是长名称如 “Hepatitis C” 这类带空格的，在 y 轴显示时改成换行；第二，图例字号不小于正文的 7pt；第三，所有子图的contrast顺序一致，比如统一 “treat_vs_ctrl → drug_vs_ctrl → drug_vs_treat”，别在最后拼版时调乱顺序。这套流程我反复用了很多个项目，翻过车的地方全在这几处，提前写在代码里就不用来回改。希望帮到你。

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