第4讲：功能富集分析

从差异基因到生物学解释

本讲概述

学习目标： - 理解GO和KEGG数据库的结构 - 掌握ORA超几何检验原理 - 理解GSEA与ORA的区别 - 掌握clusterProfiler的使用

重点难点： - 🔑 重点：ORA原理、clusterProfiler使用、基因ID转换 - ⚠️ 难点：GSEA原理、KEGG的ID转换、结果生物学解释

1. 为什么需要功能富集分析？

1.1 问题背景

差异表达分析后，你可能得到： - 上调基因：1,500个 - 下调基因：1,200个

面对数千个基因，如何理解其生物学意义？

1.2 分析目标

差异基因列表      →      功能富集分析      →      生物学洞察
┌──────────┐           ┌──────────┐           ┌──────────┐
│ Gene A   │           │ GO:1234  │           │ 细胞增殖 │
│ Gene B   │    →      │ KEGG:567 │    →      │ 信号通路 │
│ Gene C   │           │ GO:5678  │           │ 代谢过程 │
│ ...      │           │ ...      │           │ ...      │
└──────────┘           └──────────┘           └──────────┘

目标	说明
降维	从基因列表到功能类别
解释	识别受影响的生物学过程
假设生成	发现新的研究方向

2. 功能注释数据库

2.1 Gene Ontology (GO)

三大分支：

分支	缩写	描述	示例
细胞组分	CC	基因产物所在位置	nucleus, membrane
分子功能	MF	基因产物的生化活性	kinase activity
生物过程	BP	参与的生物学过程	cell division

GO结构：有向无环图（DAG）

            生物过程 (Biological Process)
                   │
         ┌─────────┴─────────┐
         │                   │
    细胞增殖            细胞凋亡
         │
    ┌────┴────┐
    │         │
DNA复制    有丝分裂

GO注释特点： - 一个基因可有多个GO注释 - 注释遵循”真路径规则”（True Path Rule） - 父节点的注释包含所有子节点

2.2 KEGG通路

数据库	内容	应用
PATHWAY	代谢通路、信号通路	通路富集分析
GENES	基因序列信息	基因注释
DISEASE	疾病相关基因	疾病研究

常见KEGG通路： - 代谢通路：糖酵解、脂肪酸代谢 - 信号通路：MAPK、PI3K-Akt、Wnt - 疾病通路：Pathways in cancer - 免疫通路：T细胞受体、B细胞受体

2.3 其他数据库

数据库	特点	用途
Reactome	人工注释，详细	通路分析
MSigDB	基因集集合	GSEA
DO (Disease Ontology)	疾病本体	疾病富集

3. 基因ID转换 🔑

3.1 为什么需要ID转换？

不同数据库使用不同的基因标识符：

ID类型	示例	使用场景
SYMBOL	TP53, BRCA1	人类可读
ENTREZID	7157, 672	NCBI数据库
ENSEMBL	ENSG00000141510	Ensembl数据库
UNIPROT	P04637	蛋白数据库
KEGG	hsa:7157	KEGG通路

3.2 查看支持的ID类型

library(org.Hs.eg.db)

# 查看人类基因注释包支持的所有ID类型
keytypes(org.Hs.eg.db)

# 常用类型包括：
# - SYMBOL: 基因符号
# - ENTREZID: Entrez Gene ID
# - ENSEMBL: Ensembl Gene ID
# - UNIPROT: UniProt ID
# - REFSEQ: RefSeq ID

3.3 bitr进行ID转换

library(clusterProfiler)

# SYMBOL → ENTREZID
gene_df <- bitr(gene_symbols,
                fromType = "SYMBOL",
                toType = "ENTREZID",
                OrgDb = org.Hs.eg.db)

# 多目标转换：SYMBOL → ENTREZID + UNIPROT
gene_df <- bitr(gene_symbols,
                fromType = "SYMBOL",
                toType = c("ENTREZID", "UNIPROT"),
                OrgDb = org.Hs.eg.db)

⚠️ 注意：ID转换可能不完全

部分基因可能没有对应的目标ID（返回NA）
一个SYMBOL可能对应多个ENTREZID（需去重）
转换后建议检查转换成功率

4. 过表达分析（ORA）🔑

4.1 ORA原理

核心问题：某GO/通路中的差异基因是否比随机期望更多？

超几何检验：

背景基因：N个（如20,000，所有检测到的基因）
背景中属于某通路的基因：M个

差异基因：n个
差异基因中属于该通路的基因：k个

P(X=k) = C(M,k) × C(N-M, n-k) / C(N,n)

4.2 ORA分析步骤

1. 获取差异基因列表（|log2FC|>1, padj<0.05）
        ↓
2. 选择背景基因集（所有检测到的基因）
        ↓
3. 基因ID转换（如SYMBOL → ENTREZID）
        ↓
4. 对每个GO/通路进行超几何检验
        ↓
5. 多重检验校正（BH）
        ↓
6. 筛选显著富集的条目（p.adjust < 0.05）

4.3 富集比值

富集比值 = (k/n) / (M/N)

- k/n：差异基因中属于该通路的比例
- M/N：背景中属于该通路的比例

富集比值 > 1：该通路被富集
富集比值 = 2：差异基因在该通路中的比例是背景的2倍

4.4 clusterProfiler ORA 🔑

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

# 1. 基因ID转换（SYMBOL → ENTREZID）
gene_df <- bitr(deg_genes,
                fromType = "SYMBOL",
                toType = "ENTREZID",
                OrgDb = org.Hs.eg.db)

# 2. GO富集分析
ego <- enrichGO(
  gene = gene_df$ENTREZID,    # 基因列表
  OrgDb = org.Hs.eg.db,       # 物种注释包
  keyType = "ENTREZID",       # 输入基因ID类型
  ont = "BP",                 # BP/MF/CC/ALL
  pAdjustMethod = "BH",       # 多重检验校正方法
  pvalueCutoff = 0.05,        # p值阈值
  qvalueCutoff = 0.2          # q值阈值
)

# 3. 查看结果
head(ego)

# 4. 导出结果
df <- as.data.frame(ego)
write.csv(df, "GO_enrichment_results.csv")

4.5 ORA结果解读

列名	含义	说明
ID	GO/通路ID	唯一标识符，如GO:0007049
Description	描述	功能名称，如”cell cycle”
GeneRatio	k/n	差异基因中属于该通路的比率，如45/120
BgRatio	M/N	背景中属于该通路的比率，如500/20000
pvalue	原始p值	富集显著性
p.adjust	校正p值	BH校正后，用于筛选
qvalue	FDR控制	更严格的校正
geneID	基因列表	属于该通路的差异基因ID
Count	k	属于该通路的差异基因数

提示

选择富集阈值： - 宽松：p.adjust < 0.05（探索性分析，返回较多结果） - 标准：p.adjust < 0.05 且 qvalue < 0.2（常用） - 严格：qvalue < 0.05（发表级别，结果更可靠）

4.6 GO三大分支分别分析

# 生物过程 (BP) - 最常用
ego_bp <- enrichGO(gene = gene_list, OrgDb = org.Hs.eg.db, ont = "BP")

# 分子功能 (MF)
ego_mf <- enrichGO(gene = gene_list, OrgDb = org.Hs.eg.db, ont = "MF")

# 细胞组分 (CC)
ego_cc <- enrichGO(gene = gene_list, OrgDb = org.Hs.eg.db, ont = "CC")

# 全部分析 (ALL)
ego_all <- enrichGO(gene = gene_list, OrgDb = org.Hs.eg.db, ont = "ALL")

5. KEGG通路富集 🔑⚠️

5.1 KEGG分析的特殊性

注意：enrichKEGG() 支持的ID类型与GO不同： - 支持：kegg, ncbi-geneid, ncbi-proteinid, uniprot - 不支持：SYMBOL, ENSEMBL

5.2 KEGG分析流程

# 方法1：使用ENTREZID（推荐）
gene_df <- bitr(gene_symbols,
                fromType = "SYMBOL",
                toType = "ENTREZID",
                OrgDb = org.Hs.eg.db)

kk <- enrichKEGG(
  gene = gene_df$ENTREZID,
  organism = "hsa",         # hsa=human, mmu=mouse, rno=rat
  keyType = "kegg",         # 或"ncbi-geneid"
  pAdjustMethod = "BH",
  pvalueCutoff = 0.05
)

# 方法2：使用UNIPROT（某些情况更稳定）
gene_df <- bitr(gene_symbols,
                fromType = "SYMBOL",
                toType = "UNIPROT",
                OrgDb = org.Hs.eg.db)

kk <- enrichKEGG(
  gene = gene_df$UNIPROT,
  organism = "hsa",
  keyType = "uniprot",
  pAdjustMethod = "BH",
  pvalueCutoff = 0.05
)

5.3 常用物种代码

物种	代码	说明
人类	hsa	Homo sapiens
小鼠	mmu	Mus musculus
大鼠	rno	Rattus norvegicus
斑马鱼	dre	Danio rerio
果蝇	dme	Drosophila melanogaster

⚠️ KEGG分析常见问题

网络问题：KEGG数据库在国外，可能被防火墙阻挡
ID转换失败：部分基因可能没有KEGG注释
版本问题：KEGG定期更新，结果可能随时间变化

6. GSEA简介 ⚠️

6.1 ORA的局限性

需要预设阈值筛选差异基因（丢失阈值附近的基因）
丢失微小但协调变化的基因
忽略基因表达变化的方向

6.2 GSEA原理

Gene Set Enrichment Analysis：

根据所有基因排序（如log2FC从高到低）
评估基因集成员是否富集在排序列表顶部或底部
计算富集分数（Enrichment Score, ES）

排序的基因列表（按log2FC从高到低）：

log2FC:  +5    +3    +1    0    -1    -3    -5
         │     │     │    │     │     │     │
Gene:   A───B───C───D───E───F───G───H───I───J
         ↑     ↑                       ↑
       通路X成员                     通路Y成员

富集分数计算：
- 通路X富集在顶部（正相关）→ ES > 0
- 通路Y富集在底部（负相关）→ ES < 0

6.3 GSEA分析

# 准备排序基因列表（所有基因，按log2FC排序）
geneList <- res$log2FoldChange
names(geneList) <- rownames(res)
geneList <- sort(geneList, decreasing = TRUE)

# 基因ID转换
geneList_df <- bitr(names(geneList),
                    fromType = "SYMBOL",
                    toType = "ENTREZID",
                    OrgDb = org.Hs.eg.db)

# 创建命名向量（ENTREZID作为名称，log2FC作为值）
names(geneList) <- geneList_df$ENTREZID[match(names(geneList), geneList_df$SYMBOL)]
geneList <- geneList[!is.na(names(geneList))]

# GSEA分析
gsea <- gseGO(
  geneList = geneList,
  OrgDb = org.Hs.eg.db,
  ont = "BP",
  pAdjustMethod = "BH",
  pvalueCutoff = 0.05,
  minGSSize = 10,     # 基因集最小基因数
  maxGSSize = 500     # 基因集最大基因数
)

6.4 GSEA结果解读

列名	含义
ID	基因集ID
setSize	基因集中检测到的基因数
enrichmentScore	富集分数（ES）
NES	标准化富集分数（不同基因集可比）
pvalue	显著性
p.adjust	校正后p值
leading_edge	核心富集基因

NES解释： - NES > 0：基因集富集在排序列表顶部（上调基因） - NES < 0：基因集富集在排序列表底部（下调基因） - |NES| > 1：通常认为有意义

6.5 ORA vs GSEA

特性	ORA	GSEA
输入	差异基因列表	全部基因（排序）
阈值依赖	是	否
检测微弱信号	弱	强
计算量	小	大
适用场景	快速筛选	深入分析

7. 结果可视化

7.1 条形图（Barplot）

# 显示Top 20富集条目
barplot(ego, showCategory = 20, title = "GO BP Enrichment")

7.2 点图（Dotplot）🔑

# 按p.adjust排序显示Top 15
dotplot(ego, showCategory = 15, orderBy = "p.adjust")

# 按基因数排序
dotplot(ego, showCategory = 15, orderBy = "Count")

点图解读： - 点大小：富集基因数（Count） - 颜色：显著性（p.adjust） - X轴：GeneRatio

7.3 网络图（Cnetplot）🔑

# 富集条目与基因的关联网络
# foldChange参数显示基因的表达变化方向
cnetplot(ego,
         categorySize = "pvalue",
         foldChange = geneList,
         showCategory = 5)

用途： - 显示富集条目间的基因重叠 - 查看哪些基因参与多个通路 - 结合表达变化方向（红色=上调，蓝色=下调）

7.4 GSEA mountain plot

# 绘制指定基因集的GSEA图
gseaplot2(gsea, geneSetID = 1, title = gsea$Description[1])

# 同时绘制多个基因集
gseaplot2(gsea, geneSetID = 1:3)

图例解读： - 绿线：富集分数曲线 - 黑色竖线：基因集中基因在排序列表中的位置 - 热图：基因表达值分布

7.5 热图（Heatplot）

# 富集通路中基因的表达热图
heatplot(ego, foldChange = geneList, showCategory = 5)

8. 完整分析流程示例

# 完整功能富集分析流程
library(DESeq2)
library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)

# 1. 差异基因筛选（实验2已完成）
deg_up <- rownames(subset(res, padj < 0.05 & log2FoldChange > 1))
deg_down <- rownames(subset(res, padj < 0.05 & log2FoldChange < -1))

# 2. ID转换
deg_up_df <- bitr(deg_up, fromType = "SYMBOL",
                  toType = "ENTREZID", OrgDb = org.Hs.eg.db)

# 3. GO富集分析
ego_up <- enrichGO(gene = deg_up_df$ENTREZID,
                   OrgDb = org.Hs.eg.db,
                   ont = "BP",
                   pvalueCutoff = 0.05,
                   qvalueCutoff = 0.2)

# 4. 可视化
dotplot(ego_up, showCategory = 15)
barplot(ego_up, showCategory = 15)

# 5. KEGG富集分析
kk_up <- enrichKEGG(gene = deg_up_df$ENTREZID,
                    organism = "hsa",
                    pvalueCutoff = 0.05)

dotplot(kk_up, showCategory = 15)

# 6. 保存结果
write.csv(as.data.frame(ego_up), "GO_BP_upregulated.csv")
write.csv(as.data.frame(kk_up), "KEGG_upregulated.csv")

9. 结果解释与报告

9.1 撰写富集分析结果

示例描述：

对差异表达基因进行GO生物过程富集分析（clusterProfiler，
p.adjust < 0.05）。上调基因显著富集于"细胞周期调控"
（GeneRatio=45/120, p.adjust=2.3×10^-10）和"DNA复制"
（GeneRatio=32/120, p.adjust=5.6×10^-8）等增殖相关通路。
下调基因则富集于"细胞凋亡"（GeneRatio=28/98,
p.adjust=8.1×10^-6）和"免疫应答"（GeneRatio=35/98,
p.adjust=3.2×10^-5）。KEGG分析显示上调基因富集于
"Pathways in cancer"（p.adjust=1.5×10^-5）和"Cell cycle"
（p.adjust=8.7×10^-4）。

这些结果表明，ATC相较于PTC可能通过促进细胞增殖
和抑制凋亡来增强恶性表型。

9.2 注意事项 ⚠️

注意事项	说明
选择合适背景	使用所有检测到的基因作为背景，而非全基因组
关注GeneRatio	高富集比（>2）更可靠
结合表达方向	上调和下调基因分别富集分析
避免过度解读	富集≠因果，需实验验证
文献验证	关键通路需结合已有文献支持
多重比较	需进行BH校正，控制FDR

9.3 常见问题排查

问题	可能原因	解决方案
富集结果为空	ID转换失败	检查keyType是否正确
	阈值过严	放宽pvalueCutoff
结果过多	阈值过松	使用qvalue或更严格的p.adjust
通路名称不显示	网络问题	检查KEGG数据库连接
基因数过少	差异基因少	放宽差异基因筛选阈值

10. 本讲小结

要点	内容
功能数据库	GO（CC/MF/BP）、KEGG是主要资源
ID转换	使用`bitr()`进行SYMBOL↔︎ENTREZID↔︎UNIPROT转换
ORA原理	超几何检验评估基因集是否过度代表
GSEA优势	不依赖阈值，可检测协调的微小变化
可视化	条形图、点图、cnetplot、GSEA mountain plot
结果解读	关注显著性、GeneRatio、生物学意义

思考题

为什么富集分析需要设置背景基因集？

点击查看答案
背景基因集决定了”预期”的随机富集水平。使用所有检测到的基因作为背景，而非全基因组，因为未检测到的基因不应该参与计算。
ORA和GSEA各自适用于什么场景？

点击查看答案
ORA适用于快速筛选已有明确差异基因列表的场景；GSEA适用于不预设阈值、希望检测协调微小变化的深入分析场景。
富集分析p值<0.05就一定有意义吗？

点击查看答案
不一定。需要：1) 使用校正后的p值(p.adjust)；2) 考虑GeneRatio（富集比）；3) 结合生物学知识判断；4) 需要文献支持。
如何判断一个富集结果是否可靠？
点击查看答案
1. p.adjust < 0.05且qvalue < 0.2；2) GeneRatio > 2；3) 富集基因数Count >= 5；4) 与文献报道一致；5) 多个数据库(GO/KEGG)一致。

参考文献

Yu G, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284-287.
Subramanian A, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. PNAS. 2005;102(43):15545-15550.
Ashburner M, et al. Gene ontology: tool for the unification of biology. Nat Genet. 2000;25(1):25-29.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27-30.
Yu G, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation. 2021;2(3):100141.

--- title: "第4讲：功能富集分析" subtitle: "从差异基因到生物学解释" --- ## 本讲概述 **学习目标**： - 理解GO和KEGG数据库的结构 - 掌握ORA超几何检验原理 - 理解GSEA与ORA的区别 - 掌握clusterProfiler的使用 **重点难点**： - 🔑 **重点**：ORA原理、clusterProfiler使用、基因ID转换 - ⚠️ **难点**：GSEA原理、KEGG的ID转换、结果生物学解释 --- ## 1. 为什么需要功能富集分析？ ### 1.1 问题背景差异表达分析后，你可能得到： - 上调基因：1,500个 - 下调基因：1,200个 **面对数千个基因，如何理解其生物学意义？** ### 1.2 分析目标 ``` 差异基因列表 → 功能富集分析 → 生物学洞察 ┌──────────┐ ┌──────────┐ ┌──────────┐ │ Gene A │ │ GO:1234 │ │ 细胞增殖 │ │ Gene B │ → │ KEGG:567 │ → │ 信号通路 │ │ Gene C │ │ GO:5678 │ │ 代谢过程 │ │ ... │ │ ... │ │ ... │ └──────────┘ └──────────┘ └──────────┘ ``` | 目标 | 说明 | |------|------| | 降维 | 从基因列表到功能类别 | | 解释 | 识别受影响的生物学过程 | | 假设生成 | 发现新的研究方向 | --- ## 2. 功能注释数据库 ### 2.1 Gene Ontology (GO) **三大分支**： | 分支 | 缩写 | 描述 | 示例 | |------|------|------|------| | **细胞组分** | CC | 基因产物所在位置 | nucleus, membrane | | **分子功能** | MF | 基因产物的生化活性 | kinase activity | | **生物过程** | BP | 参与的生物学过程 | cell division | **GO结构**：有向无环图（DAG） ``` 生物过程 (Biological Process) │ ┌─────────┴─────────┐ │ │ 细胞增殖细胞凋亡 │ ┌────┴────┐ │ │ DNA复制有丝分裂 ``` **GO注释特点**： - 一个基因可有多个GO注释 - 注释遵循"真路径规则"（True Path Rule） - 父节点的注释包含所有子节点 ### 2.2 KEGG通路 | 数据库 | 内容 | 应用 | |--------|------|------| | **PATHWAY** | 代谢通路、信号通路 | 通路富集分析 | | **GENES** | 基因序列信息 | 基因注释 | | **DISEASE** | 疾病相关基因 | 疾病研究 | **常见KEGG通路**： - **代谢通路**：糖酵解、脂肪酸代谢 - **信号通路**：MAPK、PI3K-Akt、Wnt - **疾病通路**：Pathways in cancer - **免疫通路**：T细胞受体、B细胞受体 ### 2.3 其他数据库 | 数据库 | 特点 | 用途 | |--------|------|------| | **Reactome** | 人工注释，详细 | 通路分析 | | **MSigDB** | 基因集集合 | GSEA | | **DO (Disease Ontology)** | 疾病本体 | 疾病富集 | --- ## 3. 基因ID转换 🔑 ### 3.1 为什么需要ID转换？不同数据库使用不同的基因标识符： | ID类型 | 示例 | 使用场景 | |--------|------|----------| | **SYMBOL** | TP53, BRCA1 | 人类可读 | | **ENTREZID** | 7157, 672 | NCBI数据库 | | **ENSEMBL** | ENSG00000141510 | Ensembl数据库 | | **UNIPROT** | P04637 | 蛋白数据库 | | **KEGG** | hsa:7157 | KEGG通路 | ### 3.2 查看支持的ID类型 ```r library(org.Hs.eg.db) # 查看人类基因注释包支持的所有ID类型 keytypes(org.Hs.eg.db) # 常用类型包括： # - SYMBOL: 基因符号 # - ENTREZID: Entrez Gene ID # - ENSEMBL: Ensembl Gene ID # - UNIPROT: UniProt ID # - REFSEQ: RefSeq ID ``` ### 3.3 bitr进行ID转换 ```r library(clusterProfiler) # SYMBOL → ENTREZID gene_df <- bitr(gene_symbols, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db) # 多目标转换：SYMBOL → ENTREZID + UNIPROT gene_df <- bitr(gene_symbols, fromType = "SYMBOL", toType = c("ENTREZID", "UNIPROT"), OrgDb = org.Hs.eg.db) ``` ::: {.callout-warning} ### ⚠️ 注意：ID转换可能不完全 - 部分基因可能没有对应的目标ID（返回NA） - 一个SYMBOL可能对应多个ENTREZID（需去重） - 转换后建议检查转换成功率 ::: --- ## 4. 过表达分析（ORA）🔑 ### 4.1 ORA原理 **核心问题**：某GO/通路中的差异基因是否比随机期望更多？ **超几何检验**： ``` 背景基因：N个（如20,000，所有检测到的基因）背景中属于某通路的基因：M个差异基因：n个差异基因中属于该通路的基因：k个 P(X=k) = C(M,k) × C(N-M, n-k) / C(N,n) ``` ### 4.2 ORA分析步骤 ``` 1. 获取差异基因列表（|log2FC|>1, padj<0.05） ↓ 2. 选择背景基因集（所有检测到的基因） ↓ 3. 基因ID转换（如SYMBOL → ENTREZID） ↓ 4. 对每个GO/通路进行超几何检验 ↓ 5. 多重检验校正（BH） ↓ 6. 筛选显著富集的条目（p.adjust < 0.05） ``` ### 4.3 富集比值 ``` 富集比值 = (k/n) / (M/N) - k/n：差异基因中属于该通路的比例 - M/N：背景中属于该通路的比例富集比值 > 1：该通路被富集富集比值 = 2：差异基因在该通路中的比例是背景的2倍 ``` ### 4.4 clusterProfiler ORA 🔑 ```r library(clusterProfiler) library(org.Hs.eg.db) # 1. 基因ID转换（SYMBOL → ENTREZID） gene_df <- bitr(deg_genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db) # 2. GO富集分析 ego <- enrichGO( gene = gene_df$ENTREZID, # 基因列表 OrgDb = org.Hs.eg.db, # 物种注释包 keyType = "ENTREZID", # 输入基因ID类型 ont = "BP", # BP/MF/CC/ALL pAdjustMethod = "BH", # 多重检验校正方法 pvalueCutoff = 0.05, # p值阈值 qvalueCutoff = 0.2 # q值阈值 ) # 3. 查看结果 head(ego) # 4. 导出结果 df <- as.data.frame(ego) write.csv(df, "GO_enrichment_results.csv") ``` ### 4.5 ORA结果解读 | 列名 | 含义 | 说明 | |------|------|------| | ID | GO/通路ID | 唯一标识符，如GO:0007049 | | Description | 描述 | 功能名称，如"cell cycle" | | GeneRatio | k/n | 差异基因中属于该通路的比率，如45/120 | | BgRatio | M/N | 背景中属于该通路的比率，如500/20000 | | pvalue | 原始p值 | 富集显著性 | | **p.adjust** | 校正p值 | **BH校正后，用于筛选** | | qvalue | FDR控制 | 更严格的校正 | | geneID | 基因列表 | 属于该通路的差异基因ID | | Count | k | 属于该通路的差异基因数 | ::: {.callout-tip} **选择富集阈值**： - 宽松：p.adjust < 0.05（探索性分析，返回较多结果） - 标准：p.adjust < 0.05 且 qvalue < 0.2（常用） - 严格：qvalue < 0.05（发表级别，结果更可靠） ::: ### 4.6 GO三大分支分别分析 ```r # 生物过程 (BP) - 最常用 ego_bp <- enrichGO(gene = gene_list, OrgDb = org.Hs.eg.db, ont = "BP") # 分子功能 (MF) ego_mf <- enrichGO(gene = gene_list, OrgDb = org.Hs.eg.db, ont = "MF") # 细胞组分 (CC) ego_cc <- enrichGO(gene = gene_list, OrgDb = org.Hs.eg.db, ont = "CC") # 全部分析 (ALL) ego_all <- enrichGO(gene = gene_list, OrgDb = org.Hs.eg.db, ont = "ALL") ``` --- ## 5. KEGG通路富集 🔑⚠️ ### 5.1 KEGG分析的特殊性 **注意**：`enrichKEGG()` 支持的ID类型与GO不同： - 支持：kegg, ncbi-geneid, ncbi-proteinid, uniprot - **不支持**：SYMBOL, ENSEMBL ### 5.2 KEGG分析流程 ```r # 方法1：使用ENTREZID（推荐） gene_df <- bitr(gene_symbols, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db) kk <- enrichKEGG( gene = gene_df$ENTREZID, organism = "hsa", # hsa=human, mmu=mouse, rno=rat keyType = "kegg", # 或"ncbi-geneid" pAdjustMethod = "BH", pvalueCutoff = 0.05 ) # 方法2：使用UNIPROT（某些情况更稳定） gene_df <- bitr(gene_symbols, fromType = "SYMBOL", toType = "UNIPROT", OrgDb = org.Hs.eg.db) kk <- enrichKEGG( gene = gene_df$UNIPROT, organism = "hsa", keyType = "uniprot", pAdjustMethod = "BH", pvalueCutoff = 0.05 ) ``` ### 5.3 常用物种代码 | 物种 | 代码 | 说明 | |------|------|------| | 人类 | hsa | Homo sapiens | | 小鼠 | mmu | Mus musculus | | 大鼠 | rno | Rattus norvegicus | | 斑马鱼 | dre | Danio rerio | | 果蝇 | dme | Drosophila melanogaster | ::: {.callout-warning} ### ⚠️ KEGG分析常见问题 1. **网络问题**：KEGG数据库在国外，可能被防火墙阻挡 2. **ID转换失败**：部分基因可能没有KEGG注释 3. **版本问题**：KEGG定期更新，结果可能随时间变化 ::: --- ## 6. GSEA简介 ⚠️ ### 6.1 ORA的局限性 - 需要预设阈值筛选差异基因（丢失阈值附近的基因） - 丢失微小但协调变化的基因 - 忽略基因表达变化的方向 ### 6.2 GSEA原理 **Gene Set Enrichment Analysis**： 1. 根据所有基因排序（如log2FC从高到低） 2. 评估基因集成员是否富集在排序列表顶部或底部 3. 计算富集分数（Enrichment Score, ES） ``` 排序的基因列表（按log2FC从高到低）： log2FC: +5 +3 +1 0 -1 -3 -5 │ │ │ │ │ │ │ Gene: A───B───C───D───E───F───G───H───I───J ↑ ↑ ↑ 通路X成员通路Y成员富集分数计算： - 通路X富集在顶部（正相关）→ ES > 0 - 通路Y富集在底部（负相关）→ ES < 0 ``` ### 6.3 GSEA分析 ```r # 准备排序基因列表（所有基因，按log2FC排序） geneList <- res$log2FoldChange names(geneList) <- rownames(res) geneList <- sort(geneList, decreasing = TRUE) # 基因ID转换 geneList_df <- bitr(names(geneList), fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db) # 创建命名向量（ENTREZID作为名称，log2FC作为值） names(geneList) <- geneList_df$ENTREZID[match(names(geneList), geneList_df$SYMBOL)] geneList <- geneList[!is.na(names(geneList))] # GSEA分析 gsea <- gseGO( geneList = geneList, OrgDb = org.Hs.eg.db, ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, minGSSize = 10, # 基因集最小基因数 maxGSSize = 500 # 基因集最大基因数 ) ``` ### 6.4 GSEA结果解读 | 列名 | 含义 | |------|------| | ID | 基因集ID | | setSize | 基因集中检测到的基因数 | | enrichmentScore | 富集分数（ES） | | **NES** | 标准化富集分数（不同基因集可比） | | pvalue | 显著性 | | p.adjust | 校正后p值 | | leading_edge | 核心富集基因 | **NES解释**： - NES > 0：基因集富集在排序列表顶部（上调基因） - NES < 0：基因集富集在排序列表底部（下调基因） - |NES| > 1：通常认为有意义 ### 6.5 ORA vs GSEA | 特性 | ORA | GSEA | |------|-----|------| | 输入 | 差异基因列表 | 全部基因（排序） | | 阈值依赖 | **是** | **否** | | 检测微弱信号 | 弱 | **强** | | 计算量 | 小 | 大 | | 适用场景 | 快速筛选 | **深入分析** | --- ## 7. 结果可视化 ### 7.1 条形图（Barplot） ```r # 显示Top 20富集条目 barplot(ego, showCategory = 20, title = "GO BP Enrichment") ``` ### 7.2 点图（Dotplot）🔑 ```r # 按p.adjust排序显示Top 15 dotplot(ego, showCategory = 15, orderBy = "p.adjust") # 按基因数排序 dotplot(ego, showCategory = 15, orderBy = "Count") ``` **点图解读**： - **点大小**：富集基因数（Count） - **颜色**：显著性（p.adjust） - **X轴**：GeneRatio ### 7.3 网络图（Cnetplot）🔑 ```r # 富集条目与基因的关联网络 # foldChange参数显示基因的表达变化方向 cnetplot(ego, categorySize = "pvalue", foldChange = geneList, showCategory = 5) ``` **用途**： - 显示富集条目间的基因重叠 - 查看哪些基因参与多个通路 - 结合表达变化方向（红色=上调，蓝色=下调） ### 7.4 GSEA mountain plot ```r # 绘制指定基因集的GSEA图 gseaplot2(gsea, geneSetID = 1, title = gsea$Description[1]) # 同时绘制多个基因集 gseaplot2(gsea, geneSetID = 1:3) ``` **图例解读**： - 绿线：富集分数曲线 - 黑色竖线：基因集中基因在排序列表中的位置 - 热图：基因表达值分布 ### 7.5 热图（Heatplot） ```r # 富集通路中基因的表达热图 heatplot(ego, foldChange = geneList, showCategory = 5) ``` --- ## 8. 完整分析流程示例 ```r # 完整功能富集分析流程 library(DESeq2) library(clusterProfiler) library(org.Hs.eg.db) library(enrichplot) # 1. 差异基因筛选（实验2已完成） deg_up <- rownames(subset(res, padj < 0.05 & log2FoldChange > 1)) deg_down <- rownames(subset(res, padj < 0.05 & log2FoldChange < -1)) # 2. ID转换 deg_up_df <- bitr(deg_up, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db) # 3. GO富集分析 ego_up <- enrichGO(gene = deg_up_df$ENTREZID, OrgDb = org.Hs.eg.db, ont = "BP", pvalueCutoff = 0.05, qvalueCutoff = 0.2) # 4. 可视化 dotplot(ego_up, showCategory = 15) barplot(ego_up, showCategory = 15) # 5. KEGG富集分析 kk_up <- enrichKEGG(gene = deg_up_df$ENTREZID, organism = "hsa", pvalueCutoff = 0.05) dotplot(kk_up, showCategory = 15) # 6. 保存结果 write.csv(as.data.frame(ego_up), "GO_BP_upregulated.csv") write.csv(as.data.frame(kk_up), "KEGG_upregulated.csv") ``` --- ## 9. 结果解释与报告 ### 9.1 撰写富集分析结果 **示例描述**： ```markdown 对差异表达基因进行GO生物过程富集分析（clusterProfiler， p.adjust < 0.05）。上调基因显著富集于"细胞周期调控" （GeneRatio=45/120, p.adjust=2.3×10^-10）和"DNA复制" （GeneRatio=32/120, p.adjust=5.6×10^-8）等增殖相关通路。下调基因则富集于"细胞凋亡"（GeneRatio=28/98, p.adjust=8.1×10^-6）和"免疫应答"（GeneRatio=35/98, p.adjust=3.2×10^-5）。KEGG分析显示上调基因富集于 "Pathways in cancer"（p.adjust=1.5×10^-5）和"Cell cycle" （p.adjust=8.7×10^-4）。这些结果表明，ATC相较于PTC可能通过促进细胞增殖和抑制凋亡来增强恶性表型。 ``` ### 9.2 注意事项 ⚠️ | 注意事项 | 说明 | |----------|------| | **选择合适背景** | 使用所有检测到的基因作为背景，而非全基因组 | | **关注GeneRatio** | 高富集比（>2）更可靠 | | **结合表达方向** | 上调和下调基因分别富集分析 | | **避免过度解读** | 富集≠因果，需实验验证 | | **文献验证** | 关键通路需结合已有文献支持 | | **多重比较** | 需进行BH校正，控制FDR | ### 9.3 常见问题排查 | 问题 | 可能原因 | 解决方案 | |------|----------|----------| | 富集结果为空 | ID转换失败 | 检查keyType是否正确 | | | 阈值过严 | 放宽pvalueCutoff | | 结果过多 | 阈值过松 | 使用qvalue或更严格的p.adjust | | 通路名称不显示 | 网络问题 | 检查KEGG数据库连接 | | 基因数过少 | 差异基因少 | 放宽差异基因筛选阈值 | --- ## 10. 本讲小结 | 要点 | 内容 | |------|------| | 功能数据库 | **GO**（CC/MF/BP）、**KEGG**是主要资源 | | ID转换 | 使用`bitr()`进行SYMBOL↔ENTREZID↔UNIPROT转换 | | ORA原理 | **超几何检验**评估基因集是否过度代表 | | GSEA优势 | 不依赖阈值，可检测协调的微小变化 | | 可视化 | 条形图、点图、**cnetplot**、GSEA mountain plot | | 结果解读 | 关注显著性、GeneRatio、生物学意义 | --- ## 思考题 1. **为什么富集分析需要设置背景基因集？** <details><summary>点击查看答案</summary>背景基因集决定了"预期"的随机富集水平。使用所有检测到的基因作为背景，而非全基因组，因为未检测到的基因不应该参与计算。</details> 2. **ORA和GSEA各自适用于什么场景？** <details><summary>点击查看答案</summary>ORA适用于快速筛选已有明确差异基因列表的场景；GSEA适用于不预设阈值、希望检测协调微小变化的深入分析场景。</details> 3. **富集分析p值<0.05就一定有意义吗？** <details><summary>点击查看答案</summary>不一定。需要：1) 使用校正后的p值(p.adjust)；2) 考虑GeneRatio（富集比）；3) 结合生物学知识判断；4) 需要文献支持。</details> 4. **如何判断一个富集结果是否可靠？** <details><summary>点击查看答案</summary>1) p.adjust < 0.05且qvalue < 0.2；2) GeneRatio > 2；3) 富集基因数Count >= 5；4) 与文献报道一致；5) 多个数据库(GO/KEGG)一致。</details> --- ## 参考文献 1. Yu G, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284-287. 2. Subramanian A, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. PNAS. 2005;102(43):15545-15550. 3. Ashburner M, et al. Gene ontology: tool for the unification of biology. Nat Genet. 2000;25(1):25-29. 4. Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27-30. 5. Yu G, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation. 2021;2(3):100141.