跳到正文
百沐一下百沐一下
登录
本节目录

GSEA 基因集富集：原理与实操

跟完这篇，您能说清 GSEA 到底在算什么、为什么它不用先卡差异阈值，用 clusterProfiler 跑出 GO / KEGG / MSigDB 三种基因集库的 GSEA，看懂 NES 的正负和富集曲线，并绕开「喂错输入」这类最常见的错误。

GSEA 和 ORA 差在哪：为什么它不用先卡阈值

富集分析都在回答同一个问题——「我这批基因，在哪些已知功能通路上扎堆」。ORA（过表达分析）和 GSEA（基因集富集分析）是两条技术路线，最根本的区别就在要不要先挑一份差异基因清单。

  • ORA：先按阈值（如 padj < 0.05 且 |log2FC| > 1）切出一份差异基因，只拿这份清单去查通路。阈值以下的基因，一个都不参与。
  • GSEA：不挑差异基因。把全体基因按一个指标（如 log2FC）从大到小排好队，再看某条通路的成员是整体挤在队伍的一头（上调端或下调端），还是均匀散布。

关键差别在于 GSEA 用上了两样 ORA 丢掉的信息：阈值以下的微弱变化和每个基因的排名强弱。举个常见场景：某条通路有 30 个基因，全都温和上调，但没有一个的 |log2FC| 过得了 1。

ORA 看到的GSEA 看到的
这 30 个基因全被阈值挡在门外，等于不存在全都在，且都偏排序上端
结论这条通路没富集这条通路显著上调

所以 GSEA 特别适合变化温和、分散的体系（如很多临床样本），也天然免疫「差一点点没过阈值就被漏掉」的问题——它压根不设阈值这道门。代价是它需要一份完整的排序列表，而不是一份筛过的清单。ORA 与 GSEA 的整体选型见 富集分析：GO / KEGG / GSEA。

第一步：准备一个排好序的基因列表

GSEA 的输入不是清单，而是一个命名过、按指标降序排好的数值向量：名字是基因 ID，值是排序指标。这一步准备对了，后面基本不会错。

排序指标用什么？没有唯一正确答案，但它必须同时体现方向和强弱，且您全程只用同一套。常见三选：

排序指标含义特点
log2FoldChange变化幅度直观；但不含显著性，低表达噪声基因可能排得很靠前
stat（Wald 统计量）幅度与显著性方向兼有DESeq2 直接给出，常被推荐
sign(log2FC) × −log10(pvalue)p 值定强弱、FC 定方向两者兼顾，需自己算
# 以 DESeq2 结果 res 为例，用 stat 作排序指标（也可换成 log2FoldChange）
geneList <- res$stat
names(geneList) <- res$entrez            # 事先转好的 ENTREZID（见下节）

geneList <- na.omit(geneList)            # 去掉 NA
geneList <- sort(geneList, decreasing = TRUE)   # 必须降序

三条容易被忽略、却直接决定结果对不对的要点：

  1. 必须是全体基因，不是差异基因子集。喂进筛过的清单，GSEA 就退化了、失去意义。
  2. 必须降序排列，gseGO / gseKEGG 默认从上端开始走。
  3. 一个基因只能有一行。多探针映射到同一基因时先去重（常取绝对值最大或取平均），否则重复基因会干扰计分。

GSEA 是怎么算的：富集分数与 NES

理解了算法，读结果就不再靠猜。GSEA 沿着您排好的队伍从头走到尾，维护一个「跑分」：

  1. 每遇到一个属于目标通路的基因，跑分加一点（还会按这个基因的排序指标大小加权，排得越靠前加得越多）；
  2. 每遇到一个不属于该通路的基因，跑分减一点；
  3. 走完全程，跑分偏离 0 的最大幅度，就是这条通路的富集分数 ES（Enrichment Score）。

如果通路基因都挤在队首，跑分会一路冲高，ES 是个大正数；挤在队尾则 ES 是大负数；均匀散布则跑分在 0 附近晃、ES 接近 0。

两个由此衍生的关键量：

  • NES（标准化富集分数）：不同通路基因数不同，ES 不能直接比。GSEA 用置换（打乱标签重算）得到一个零分布，把 ES 除以该分布的均值幅度，得到 NES。NES 才是可以跨通路比较的量，也是结果里最该看的一列。
  • 显著性：同样靠置换检验给出 p 值，再对成百上千条通路做多重检验校正，得到 padj（FDR）。筛通路看 padj，不看原始 p 值。

还有一个常被用到的概念——leading edge（核心富集）子集：跑分到达峰值那一刻、已经出现过的那批通路基因，就是真正把这条通路「顶上去」的核心。讲机制时，落点应该是这批基因，而不是通路里的全部成员。

用 clusterProfiler 跑一遍

clusterProfiler 把上面的算法都封装好了，几行就能跑。注意 gseKEGG 默认要 ENTREZID，若手上是 SYMBOL / ENSEMBL，先用 bitr() 转好再命名向量。

library(clusterProfiler)
library(org.Hs.eg.db)

# GO 版（ont = "BP" / "CC" / "MF" / "ALL"）
gse <- gseGO(geneList     = geneList,
             OrgDb        = org.Hs.eg.db,
             keyType      = "ENTREZID",
             ont          = "BP",
             minGSSize    = 10,      # 通路太小结果不稳，默认 10
             maxGSSize    = 500,     # 通路太大信号被稀释，默认 500
             pvalueCutoff = 0.05)    # 按 padj 过滤，0.05 是常用惯例

# KEGG 版（organism 用物种代号：人类 hsa、小鼠 mmu）
gsekegg <- gseKEGG(geneList = geneList, organism = "hsa", pvalueCutoff = 0.05)

minGSSize = 10 和 maxGSSize = 500 是 clusterProfiler 的默认值，也是通行做法：太小的通路几个基因就能偶然显著，太大的通路（几百上千基因）信号会被稀释。这两条线是惯例、按数据调，不是铁律。

基因集库怎么选

跑 GSEA 前得先定「拿哪套通路来比」。不同库覆盖面和风格差别很大，常用这几种：

基因集库特点在 R 里怎么取
GO（BP/CC/MF）覆盖最广、最细，但条目多、冗余高gseGO + org.Hs.eg.db
KEGG经典代谢 / 信号通路，图谱直观gseKEGG(organism = "hsa")
MSigDB Hallmark（H）50 条经过精心整理、去冗余的「主题」通路，信噪比高msigdbr + GSEA()
Reactome通路层级清晰、更新较勤ReactomePA::gsePathway

新手建议从 MSigDB 的 Hallmark（H 类）起步：它只有 50 条，是 Broad 研究所从大量重叠通路里归并出的高质量「主题」，结果好读、不容易被一堆近义的 GO 条目淹没。取用要靠 msigdbr 包把基因集整理成 clusterProfiler 认的两列格式：

library(msigdbr)
library(clusterProfiler)

# 取人类 Hallmark 基因集（H 类），整理成「通路名 - 基因」两列
h_sets    <- msigdbr(species = "Homo sapiens", category = "H")
term2gene <- h_sets[, c("gs_name", "entrez_gene")]

gse_h <- GSEA(geneList, TERM2GENE = term2gene, pvalueCutoff = 0.05)

MSigDB 还有 C2（含 Reactome / KEGG / WikiPathways 等整理通路）、C5（GO）、C7（免疫）等多个子集，把上面的 category 换成对应代号即可。选库没有标准答案，报告时写清用了哪个库、哪个版本，比选哪个更重要。

结果怎么读：NES、padj 与富集曲线

gseGO() 返回的结果表，逐列看懂它：

列含义怎么用
Description通路名称讲故事的落点，但要放到最后看
setSize该通路在您列表里命中的基因数太小的通路（十几个）结果别轻信
enrichmentScore原始富集分数 ES跨通路不可比，看 NES
NES标准化富集分数**正 = 偏排序上端（上调），负 = 偏下端（下调）；`
pvalue置换检验原始 p 值不直接当判据
p.adjust多重检验校正后筛通路就看这一列
core_enrichmentleading-edge 核心基因讲机制时落到这批基因

配套的图，clusterProfiler / enrichplot 一行就能出：

library(enrichplot)

# 富集曲线（跑分图）：细看某一条通路
gseaplot2(gse, geneSetID = 1, title = gse$Description[1])

# 点图：一次看多条，按激活 / 抑制分面
dotplot(gse, showCategory = 15, split = ".sign") + facet_grid(. ~ .sign)

富集曲线是 GSEA 最标志性的图，分三层看：

图上元素代表什么怎么读
顶部绿色曲线跑分（running ES）峰在左（上调端）= NES 正；峰在右（下调端）= NES 负
中间黑色竖线该通路基因在排序中的位置在某一端扎堆＝这条通路确实偏向那头
底部灰色柱全体基因的排序指标由正到负穿过 0，划开上调区与下调区

读图顺序建议：先看 padj（够不够显著）→ 再看 NES 的正负和大小（方向与强弱）→ 看曲线的峰和黑线是否真的在一端扎堆 → 最后才落到通路名去讲生物学。 关于富集图与统计表的通用读法，见 读懂分析结果。

一点关于显著性阈值的诚实说明：clusterProfiler 默认 pvalueCutoff = 0.05；而 Broad 的经典 GSEA 文档在样本量有限、只做假设生成时，建议放宽到 FDR < 0.25。两者都是惯例、按目的与数据调——您是要严格锁定还是探索性找方向，会决定这条线划在哪。

三个最常见的坑

  1. 喂了差异基因子集。这是最致命也最常见的错。GSEA 的前提就是全体基因排序，喂筛过的清单，算法直接失去意义。要跑 GSEA 就给完整列表，要用清单就改回 ORA。
  2. 排序指标只用了 p 值、丢了方向。-log10(pvalue) 不带符号，上调下调会混在一起，NES 方向失真。带方向的指标（stat、log2FC、或 sign(log2FC) × −log10(p)）才对。
  3. 把 NES 跨数据集硬比。NES 是相对该次分析自己的置换分布标准化的，换一份数据、换一套基因集库，绝对数值不宜直接比大小，比的是「在各自分析里谁更靠前、方向如何」。

再叠一句和所有富集分析共通的提醒：把 GSEA 当成生成假设的工具，不是盖章结论的工具。 显著 ≠ 机制强，热门通路更容易上榜；它告诉您往哪几个方向看，真正的机制还得靠后续验证。

在平台上跑一遍

在百沐一下不写代码也能做这一步：上传差异分析结果（含基因 ID、log2FC、p 值），用一句话说清目的，例如「用 stat 排序对全体基因做 GSEA，基因集库用 Hallmark」，组学智能体会自动完成 ID 转换、排序、GSEA 计算，并出富集曲线和点图，还给可复现的 R 代码。用法见 组学智能体使用指南。

相关