富集分析：GO / KEGG / GSEA

跟完这篇，您能把一份差异分析结果变成「这批基因富集在哪些通路」的结论：分清 ORA 与 GSEA 各自何时用，用 clusterProfiler 跑出 GO / KEGG / GSEA，读懂气泡图和富集曲线，并避开最常见的过度解读。

整篇的顺序是固定的四步，照做即可：

先选路线——ORA 还是 GSEA（决定要不要先卡差异基因）；
备好两样输入——差异分析结果 + 匹配的注释包；
跑富集——enrichGO / enrichKEGG（ORA）或 gseGO / gseKEGG（GSEA）；
读图、去冗余、克制地讲结论。

先分清：ORA 还是 GSEA（要不要先定差异基因）

富集分析回答同一个问题——「我这批基因，在哪些已知的功能通路上扎堆出现」。但有两条技术路线，区别就在于要不要先卡一份差异基因清单。

ORA（过表达分析，Over-Representation Analysis）：先按阈值挑出一份差异基因清单，再看每条通路里「被您挑中的基因」是不是显著多于随机预期。本质是超几何检验 / Fisher 精确检验。GO、KEGG 富集最常见的用法就是它。
GSEA（基因集富集分析，Gene Set Enrichment Analysis）：不挑差异基因，把全体基因按某个指标（如 log2FC）排好序，看某条通路的基因是整体偏向排名一头（上调端或下调端），还是均匀散布。

两者各有取舍：

	ORA	GSEA
输入	一份差异基因清单（先卡阈值）	全体基因 + 一个排序指标
检验对象	通路里差异基因是否显著偏多	通路基因是否集中在排序的一端
阈值依赖	强（换阈值结果会变）	无需卡差异阈值
方向信息	会丢（上下调混在一起会抵消）	保留（NES 正负即方向）
适合	已有明确的差异基因、想快速看功能	变化温和分散、担心卡阈值漏信息

一句话选型：手上已有一份差异基因、想直接看功能分布，用 ORA；担心「差一点点没过阈值」的一片小变化被漏掉，用 GSEA。 两个都跑、互相印证，也很常见。

开始前：准备数据和注释包

动手前把两样东西备齐：

一份差异分析结果：如 DESeq2 输出的表，至少含基因 ID、log2FoldChange、padj（GSEA 还会用到 stat）。
与物种、注释库匹配的注释包：人类常用 org.Hs.eg.db；KEGG 用物种代号，人类 hsa、小鼠 mmu。

有一个坑几乎人人踩：ID 类型要对上。enrichGO 可以直接吃 SYMBOL 或 ENSEMBL（用 keyType 指明就行），但 enrichKEGG / gseKEGG 默认要 ENTREZID。所以先看清自己手上是哪种 ID，需要时用 bitr() 转换：

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

# SYMBOL → ENTREZID，会丢一部分映射不上的基因，属正常现象
ids <- bitr(rownames(res), fromType = "SYMBOL",
            toType = "ENTREZID", OrgDb = org.Hs.eg.db)

# 留意丢了多大比例：映射率太低说明 ID 类型或物种可能选错了
nrow(ids) / nrow(res)

用 ORA 跑 GO 和 KEGG

ORA 的第一步是定差异基因。常用惯例是 padj < 0.05 且 |log2FC| > 1（约 2 倍变化）——这是惯例、不是铁律，样本量小或信号弱时可放宽到 |log2FC| > 0.585（1.5 倍），按您的数据和领域调。

GO 分三大类，报告时通常三类都看或按需取：

缩写	中文	回答的问题
BP	生物学过程	这些基因参与什么过程（如「炎症反应」）
CC	细胞组分	蛋白定位在哪（如「线粒体内膜」）
MF	分子功能	分子层面干什么（如「ATP 结合」）

用 clusterProfiler 跑，简例如下：

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

# 1. 从差异表里取出显著的上调 + 下调基因（以 SYMBOL 行名为例）
deg  <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
gene <- bitr(rownames(deg), fromType = "SYMBOL",
             toType = "ENTREZID", OrgDb = org.Hs.eg.db)$ENTREZID

# 2. GO 富集（ont = "BP" / "CC" / "MF" / "ALL"）
ego <- enrichGO(gene          = gene,
                OrgDb         = org.Hs.eg.db,
                keyType       = "ENTREZID",
                ont           = "BP",
                pAdjustMethod = "BH",     # 多重检验校正，BH 最常用
                pvalueCutoff  = 0.05,
                readable      = TRUE)     # 结果里把 ID 换回基因名

# 3. KEGG 富集（organism 用物种代号，人类 = "hsa"；需要联网取 KEGG 库）
ekegg <- enrichKEGG(gene = gene, organism = "hsa", pvalueCutoff = 0.05)
# enrichKEGG 结果里是 ENTREZID，想看基因名再转一次：
ekegg <- setReadable(ekegg, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")

两个细节容易被忽略，却直接影响结果对不对：

背景基因集（universe）要对。ORA 默认拿注释库全体基因当背景。更稳妥的做法是把 universe 设成您实验里检测到 / 表达的全部基因（如过滤后进入 DESeq2 的那份基因列表），而不是全基因组，否则 p 值会偏乐观。
padj（校正后 p 值）才是判据，不是原始 p 值。通路成百上千，不校正会有大量假阳性；padj < 0.05 是常用惯例，按需可收紧到 0.01。

用 GSEA 跑（全体排序基因，不设阈值）

GSEA 不挑差异基因，它需要一个命名过、并按指标从大到小排好序的数值向量：名字是基因 ID，值是排序指标。

排序指标用什么？没有唯一正确答案，只要能同时体现「方向 + 强弱」、且全程用同一套即可。常见选择：

排序指标	含义	说明
`log2FoldChange`	变化倍数	最直观，但没考虑离散度
`stat`（Wald 统计量）	变化大小 + 显著性，带方向	DESeq2 直接输出，常用
`sign(log2FC) × -log10(pvalue)`	显著性带方向	也常见，突出显著性

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

# 1. 造排序向量：用 res$stat（或 log2FoldChange），名字用 ENTREZID
geneList <- res$stat
names(geneList) <- res$entrez               # 事先把行名转成 ENTREZID
geneList <- geneList[!is.na(names(geneList)) & !is.na(geneList)]
geneList <- geneList[!duplicated(names(geneList))]  # ID 必须唯一，先去重
geneList <- sort(geneList, decreasing = TRUE)       # 必须降序

# 2. GSEA（GO 版）
gse <- gseGO(geneList     = geneList,
             OrgDb        = org.Hs.eg.db,
             keyType      = "ENTREZID",
             ont          = "BP",
             pvalueCutoff = 0.05)

# 3. GSEA（KEGG 版）
gsekegg <- gseKEGG(geneList = geneList, organism = "hsa", pvalueCutoff = 0.05)

两点提醒：GSEA 要求基因 ID 唯一，重复的要先去重（上面 !duplicated 那步）；另外过小 / 过大的基因集会被默认过滤掉（由 minGSSize / maxGSSize 控制），这是为了统计稳健，一般无需改。

GSEA 结果里最关键的一列是 NES（标准化富集分数）：正值表示这条通路整体偏向排序上端（在您的对比里上调），负值偏向下端（下调）。这正是 ORA 丢掉、而 GSEA 保住的方向信息。

读懂气泡图和富集曲线

富集结果最常用三种图，clusterProfiler 一行就能出：

library(enrichplot)

dotplot(ego, showCategory = 15)   # 气泡图：一次看多条通路
barplot(ego, showCategory = 15)   # 条形图：更直观的排名

# GSEA 专用的富集曲线（跑分图）
gseaplot2(gse, geneSetID = 1, title = gse$Description[1])

气泡图 / 条形图上每个视觉元素代表什么，照下表读：

图上元素	代表什么	怎么读
纵轴	通路名称	一般按 `padj` 从显著到不显著排
横轴（GeneRatio）	命中比例	您清单里落在该通路的基因数 ÷ 清单里能注释上的总基因数
点的大小（Count）	命中基因个数	越大＝这条通路里被您挑中的基因越多
点的颜色	`padj`	越深＝越显著（图例会标明方向）

读图顺序建议：先看颜色（够不够显著）→ 再看大小和横轴（命中够不够多）→ 最后才落到通路名去讲生物学意义。 GSEA 的富集曲线则看两处：绿色跑分线的峰偏左（上调端）还是偏右（下调端），以及底部黑竖线（该通路的基因）是否在一端扎堆。

别过度解读：一堆通路里怎么挑对的讲

富集分析最容易出问题的不是跑，而是讲。一份结果动辄几十上百条显著通路，很容易挑出最符合预期的那条讲成故事。请守住几条底线：

GO 冗余要合并。GO 是有上下位关系的树，「免疫反应」和它一堆子术语会同时富集、彼此高度重叠。用 simplify(ego) 按语义相似度去冗余，别把同一件事重复报十遍。
别只挑符合假设的那条。先看 top 榜的整体图景（哪几类功能反复出现），再落到具体通路；与预期矛盾的、没入选的，也要如实纳入判断，而不是选择性失明。
上调、下调分开跑。ORA 把上下调基因混成一张清单，方向会互相抵消，通路层面看不出是激活还是抑制。要方向，就分开跑，或改用 GSEA 看 NES 正负。
显著 ≠ 重要。padj 小只说明「这种富集不像随机」，说明不了机制有多强；大基因集天然更容易显著。统计显著要和效应大小、生物学合理性一起看。
当心注释偏倚。研究越多的基因和通路，注释越全，也就越容易「中」。热门通路上榜，未必是您体系里最关键的。
记得 ID 映射会丢基因。bitr() 转换、注释库覆盖不全都会丢一部分基因，留意丢了多大比例，必要时在正文说明。

一句话：把富集当成「生成假设」的工具，不是「盖章结论」的工具。 它告诉您往哪几个方向去看，真正的机制还得靠后续验证。

在平台上跑一遍

在百沐一下不写代码也能做这一步：上传您的差异分析结果（或表达矩阵 + 分组表），用一句话说清目的，例如「对这批差异基因做 GO 和 KEGG 富集」或「用 stat 排序做 GSEA」，组学智能体会自动完成 ID 转换、背景设定、富集计算并出气泡图 / 富集曲线，还给可复现的 R 代码。用法见组学智能体使用指南。