# 第 09 章 功能富集分析

## 本章目标

学完本章，你将能够：

1. 理解富集分析的思想与统计学原理（超几何分布 / 费舍尔精确检验）；
2. 用 clusterProfiler 做 GO 富集，并分清 BP、CC、MF 三个层面；
3. 用 clusterProfiler 做 KEGG 通路富集，并了解 pathview 通路图；
4. 画出气泡图 / 柱状图并正确解读；
5. 理解 ORA 与 GSEA 的区别，能跑通 GSEA 并解读富集分数（ES）曲线。

## 正文

### 1. 富集分析是什么：从“一堆基因”到“一条通路”

第 07 章我们拿到了几百上千个差异基因，但“这些基因在功能上意味着什么？”——单个基因名看不出名堂。**功能富集分析**（functional enrichment analysis）回答的正是这个问题：给定一堆基因，看它们是否显著地集中在某些已知的功能类别或通路上。

统计思想用“袋子抽球”类比：想象一个袋子里有 N 个球，代表基因组里所有被测到的基因，其中 M 个是红球，代表“某条通路里的基因”。你随机抽了 n 个球（我们的差异基因），发现里面有 k 个红球：k 这么大是纯属巧合，还是说明这条通路真的和我们的基因列表有关？

这正是**超几何分布**（hypergeometric distribution）描述的“不放回抽球”概率问题；等价地，也可以用一张 2×2 列联表做**费舍尔精确检验**（Fisher's exact test）算出“纯属巧合”的 p 值。富集分析就是逐条通路（或 GO 条目）做这个检验，然后同样用 BH 校正控制多重检验。这一类方法统称 ORA（over-representation analysis，过度代表分析）。

两个容易踩的坑：

- **背景基因集（universe）**：默认是所有被测到的基因，而不是全基因组。背景选错了，结果会系统性偏差——用表达矩阵覆盖的基因当背景最稳妥；
- **输入列表的质量**：富集分析只在输入列表可靠时才有意义，所以务必使用第 07 章严格筛选出的差异基因。

### 2. GO 富集：三个层面

GO（Gene Ontology，基因本体）把基因功能分成三个层面：

- **BP**（biological process，生物学过程）：如炎症反应、细胞程序性死亡；
- **CC**（cellular component，细胞组分）：如线粒体、质膜；
- **MF**（molecular function，分子功能）：如半胱氨酸型内肽酶活性。

clusterProfiler 是 Bioconductor 上最主流的富集工具。标准流程如下：

```r
# BiocManager::install("clusterProfiler")
# BiocManager::install("org.Hs.eg.db")
library(clusterProfiler)
library(org.Hs.eg.db)

deg_entrez <- bitr(rownames(deg), fromType = "SYMBOL",
                   toType = "ENTREZID", OrgDb = org.Hs.eg.db)$ENTREZID

ego <- enrichGO(gene = deg_entrez,
                OrgDb = org.Hs.eg.db,
                keyType = "ENTREZID",      # 输入基因 ID 的类型
                ont = "BP",                # 可选 BP / CC / MF / ALL
                pAdjustMethod = "BH",
                pvalueCutoff = 0.05,       # 条目自身 p 值阈值
                qvalueCutoff = 0.2,        # q 值阈值（更严格）
                readable = TRUE)           # 结果中显示基因符号而非 Entrez ID

head(as.data.frame(ego))
```

关键参数解释：

- `bitr()`：**ID 转换**。enrichGO 默认接受 Entrez ID，而差异基因通常是基因符号（SYMBOL），需先转换（`fromType` / `toType` 声明来源与目标类型）；
- `keyType`：告诉 enrichGO 输入基因的 ID 类型；
- `pvalueCutoff` 与 `qvalueCutoff`：q 值是 p 值经多重检验校正后的版本（同第 07 章 adj.P 的思路）。clusterProfiler 默认 `qvalueCutoff = 0.2`，较宽松，富集检验通常只要求“信号存在”；
- `readable = TRUE`：把结果中的 Entrez ID 换回基因符号，方便阅读。

结果表的关键列：ID（GO 条目号）、Description（功能描述）、GeneRatio（差异基因中属于该条目的比例）、BgRatio（背景基因中属于该条目的比例）、pvalue、p.adjust、qvalue、geneID（参与基因）。

### 3. KEGG 富集：通路数据库

KEGG（Kyoto Encyclopedia of Genes and Genomes，京都基因与基因组百科全书）是**通路数据库**，它的条目是“通路图”而不是功能词条。KEGG 富集同样基于超几何检验，但基因要用 Entrez ID，并指定物种代码（人类为 `hsa`）：

```r
ekegg <- enrichKEGG(gene = deg_entrez,
                    organism = "hsa",      # 人类
                    keyType = "kegg",
                    pvalueCutoff = 0.05,
                    pAdjustMethod = "BH")
head(as.data.frame(ekegg))
```

如果结果为空（新手常见），先放宽 `pvalueCutoff` 看有没有信号，再检查基因 ID 是否转换正确。KEGG 通路图可以用 **pathview** 绘制：它把差异基因的 logFC 着色到通路图上（红 = 上调、绿 = 下调），直观展示基因在通路中的位置：

```r
# BiocManager::install("pathview")
library(pathview)
fc <- setNames(res$logFC, rownames(res))
fc_map <- bitr(names(fc), fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
fc_vec <- fc[fc_map$SYMBOL]; names(fc_vec) <- fc_map$ENTREZID   # 重复 ID 可去重，教程略
pathview(gene.data = fc_vec,
         pathway.id = "hsa04668",   # 示例：TNF 信号通路
         species = "hsa", out.suffix = "deg")
```

注意：pathview 需联网从 KEGG 下载通路图；KEGG 有使用条款（学术用途免费），论文中要正确引用。

### 4. 结果可视化：气泡图与柱状图

clusterProfiler 的结果对象自带绘图函数，最常用的是**气泡图**（dotplot）与**柱状图**（barplot）：

```r
dotplot(ego, showCategory = 15)   # 横轴 GeneRatio，点大小 = 基因数，颜色 = q 值
barplot(ego, showCategory = 15)   # 柱长 = 基因数/比例，颜色 = 校正 p 值
```

- **气泡图怎么读**：横轴 GeneRatio（差异基因中属于该条目的比例），纵轴是富集到的条目名，点越大包含的基因越多，颜色越红越显著。
- **柱状图怎么读**：柱长代表包含的基因数（或 GeneRatio），颜色代表显著性。

实际论文常把 GO 的 BP/CC/MF 三个层面并排展示，或用 `enrichplot::cnetplot()` 画“条目—基因”连接网络图。下图是富集结果示例（模拟数据绘制，仅示意样式）：

![图09 GO/KEGG 富集气泡图（示例，模拟数据）](../figures/fig09_go_kegg.png)

### 5. GSEA 进阶：用全部基因排序

ORA 有一个明显局限：它需要先设阈值筛出差异基因，阈值两端的基因被“一刀切”，且只回答“是否富集”，不回答“富集的方向”。**基因集富集分析**（gene set enrichment analysis, GSEA）换个思路：不筛基因，而是把**所有**基因按与表型的关联强度（如 log2FC）从高到低排序，再看某条通路里的基因是整体偏前（激活/上调）还是偏后（抑制/下调）。

两者的区别一句话：**ORA 用“取出来的差异基因”，GSEA 用“全部基因的排序”。**

GSEA 的核心统计量是**富集分数**（enrichment score, ES）：从排序顶端往下走，遇到通路基因加分、遇到非通路基因减分，累计曲线偏离 0 的最大值就是 ES。ES 为正表示通路基因集中在排序前段（整体上调）。之后用置换检验评估显著性，得到 NES（标准化 ES）和 FDR。

clusterProfiler 中跑 GSEA（以 KEGG 为例）：

```r
genelist <- sort(setNames(res$logFC, rownames(res)), decreasing = TRUE)  # 全部基因按 logFC 降序
gl_map <- bitr(names(genelist), fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
genelist <- genelist[gl_map$SYMBOL]
names(genelist) <- gl_map$ENTREZID

gsea_res <- gseKEGG(geneList = genelist,
                    organism = "hsa",
                    minGSSize = 10, maxGSSize = 500,
                    pvalueCutoff = 0.25,     # GSEA 常用较宽松的阈值
                    seed = TRUE)
head(as.data.frame(gsea_res))

gseaplot2(gsea_res, geneSetID = 1, title = "示例通路 GSEA 曲线")
ggsave("../figures/fig10_gsea.png", width = 7, height = 5, dpi = 300)
```

`gseaplot2` 的输出包含三部分：顶部是排序基因的位置标记（黑色竖线代表该通路基因）；中部是 ES 累积曲线（富集分数随排序位置的变化）；底部是排序度量值（如 logFC）。曲线在**左侧**冲高，说明通路基因整体偏前（上调）。**leading edge**（前沿子集）指曲线达到最大偏移之前的那段基因——“驱动富集”的核心基因，值得重点关注。

![图10 GSEA 富集曲线（示例，模拟数据）](../figures/fig10_gsea.png)

GO 版本的 GSEA 用法类似：`gseGO(geneList = genelist, OrgDb = org.Hs.eg.db, ont = "BP", ...)`。fgsea 包是 GSEA 的更快独立实现，进阶可选。GSEA 需要全部基因的排序向量，这正是第 07 章保留全基因结果表的原因。

### 6. 结果解读注意

- **读懂通路名的生物学含义**：比如富集到"cytokine-cytokine receptor interaction"（细胞因子—细胞因子受体相互作用）提示炎症信号活跃，与焦亡（促炎性细胞死亡）的课题主题吻合。不要只报编号。
- **富集 ≠ 因果**：某通路“富集”只说明通路里的基因在你的列表中出现得比随机多，不代表通路一定被激活或抑制；方向性问题要靠 GSEA 回答。
- **避免过度解读**：一个 GO 条目描述很长时，看它包含的具体基因（geneID 列）来核实；显著但只包含 1–2 个基因的条目通常不可靠。GO 条目有父子层级，可能造成“重复富集”，别当成独立发现。
- **分层报告**：GO 的 BP/CC/MF、KEGG、GSEA 回答不同问题，论文方法部分要交代各自的输入基因、数据库版本与阈值。

## 本章小结

- 富集分析用超几何分布 / 费舍尔精确检验回答“我的基因是否集中在某些通路”，记得做 BH 校正；背景基因集要选对。
- GO 分 BP（过程）/CC（组分）/MF（功能）三层；clusterProfiler 流程：`bitr` 转 ID → `enrichGO` → `dotplot`。
- KEGG 用 `enrichKEGG(gene, organism = "hsa")`，通路图可用 pathview 把 logFC 着色上去。
- ORA 用差异基因子集；GSEA 用全部基因的排序，输出 ES 曲线，能回答富集方向。
- 解读时看 geneID 核实、区分相关与因果、报告时明确输入与阈值。

## 思考题

1. ORA 里“背景基因集”选错会怎样？比如用全基因组做背景，而你的表达矩阵只覆盖 2 万个基因。
2. 为什么 GSEA 不需要先筛差异基因？它和 ORA 各适合什么场景？
3. 富集到"NF-kappa B signaling pathway"（q = 0.01，含 5 个基因）和富集到"Apoptosis"（q = 0.05，含 40 个基因），你会优先报告哪个？为什么？
4. 如果我们的焦亡基因只有 8 个，直接拿它们做富集分析合理吗？怎样做更稳妥？
