基因本体GO富集分析¶

一句话概述：GO富集分析用clusterProfiler包检验差异基因是否在特定功能（BP/CC/MF）或KEGG通路上"扎堆"，是转录组分析的必备后续步骤。

核心知识点速览¶

概念	白话解释
Gene Ontology（GO）	用统一术语描述基因功能的"字典"，分三大类
BP（Biological Process）	生物过程，基因参与了什么事件（如"细胞凋亡"）
CC（Cellular Component）	细胞组分，基因产物在细胞哪里（如"线粒体"）
MF（Molecular Function）	分子功能，基因产物能干什么（如"激酶活性"）
KEGG	代谢通路数据库，画出基因参与的信号通路全景图
ORA（过表示分析）	超几何检验，看差异基因列表里某功能的基因是否"超标"
背景基因集	"全班同学"，用来和"考满分的同学"对比
p.adjust	多重检验校正后的p值，常用BH方法，<0.05才算显著
Rich Factor	差异基因中属于某通路的比例，越大说明富集越强
clusterProfiler	余光创开发的R包，GO/KEGG分析的"瑞士军刀"，2026年版本4.19

一、clusterProfiler安装¶

# 安装BiocManager（如果没有的话）
if (!require("BiocManager")) install.packages("BiocManager")  # 检查并安装

# 安装clusterProfiler及相关包
BiocManager::install("clusterProfiler")  # 核心包（v4.19，2026最新）
BiocManager::install("org.Hs.eg.db")     # 人类基因注释数据库
BiocManager::install("enrichplot")        # 富集分析可视化
BiocManager::install("DOSE")              # 疾病本体富集

# 加载包
library(clusterProfiler)  # 富集分析核心
library(org.Hs.eg.db)     # 人类基因注释
library(enrichplot)        # 可视化
library(ggplot2)           # 画图基础

二、基因ID转换¶

# 差异基因列表（通常来自DESeq2或edgeR的结果）
gene_symbols <- c("TP53", "BRCA1", "EGFR", "MYC", "CDK2",
                   "RB1", "PTEN", "AKT1", "MTOR", "PIK3CA")  # 基因Symbol列表

# Symbol转Entrez ID（clusterProfiler很多函数需要Entrez ID）
gene_entrez <- bitr(
  gene_symbols,                # 输入基因列表
  fromType = "SYMBOL",         # 输入ID类型
  toType = "ENTREZID",         # 输出ID类型
  OrgDb = org.Hs.eg.db         # 物种注释数据库
)
head(gene_entrez)  # 查看转换结果

# 也可以一次转多种ID
gene_multi <- bitr(
  gene_symbols,
  fromType = "SYMBOL",
  toType = c("ENTREZID", "ENSEMBL", "UNIPROT"),  # 同时转多种ID
  OrgDb = org.Hs.eg.db
)

三、GO富集分析（ORA方法）¶

3.1 标准GO富集¶

# 运行GO富集分析
go_enrich <- enrichGO(
  gene = gene_entrez$ENTREZID,   # 差异基因的Entrez ID
  OrgDb = org.Hs.eg.db,          # 物种注释数据库
  keyType = "ENTREZID",          # 输入基因ID类型
  ont = "ALL",                    # GO类别：BP/CC/MF/ALL
  pAdjustMethod = "BH",          # 多重检验校正方法（Benjamini-Hochberg）
  pvalueCutoff = 0.05,           # p值阈值
  qvalueCutoff = 0.05,           # q值（FDR）阈值
  readable = TRUE                 # 结果中显示基因Symbol而非ID
)

# 查看结果
head(go_enrich, 10)  # 显示前10个富集条目

# 分别做三个子类
go_bp <- enrichGO(gene = gene_entrez$ENTREZID, OrgDb = org.Hs.eg.db,
                  ont = "BP", pvalueCutoff = 0.05, readable = TRUE)  # 生物过程
go_cc <- enrichGO(gene = gene_entrez$ENTREZID, OrgDb = org.Hs.eg.db,
                  ont = "CC", pvalueCutoff = 0.05, readable = TRUE)  # 细胞组分
go_mf <- enrichGO(gene = gene_entrez$ENTREZID, OrgDb = org.Hs.eg.db,
                  ont = "MF", pvalueCutoff = 0.05, readable = TRUE)  # 分子功能

3.2 去冗余（simplify）¶

# GO术语有层级关系，结果中会有很多重复/冗余条目
# 用simplify去除冗余
go_simplified <- simplify(
  go_enrich,
  cutoff = 0.7,          # 语义相似性阈值（越小去除越多）
  by = "p.adjust",       # 保留p值更小的条目
  select_fun = min       # 选择函数
)

四、KEGG通路富集分析¶

# KEGG富集分析（在线查询，需要联网）
kegg_enrich <- enrichKEGG(
  gene = gene_entrez$ENTREZID,   # 基因Entrez ID
  organism = "hsa",               # 物种缩写（hsa=人类，mmu=小鼠）
  keyType = "kegg",               # ID类型
  pAdjustMethod = "BH",           # 校正方法
  pvalueCutoff = 0.05,            # p值阈值
  qvalueCutoff = 0.05             # q值阈值
)

# 查看结果
head(kegg_enrich)  # 显示前几条

# KEGG结果基因ID转换为Symbol（方便阅读）
kegg_readable <- setReadable(kegg_enrich, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")

常用物种代码¶

人类: hsa | 小鼠: mmu | 大鼠: rno | 果蝇: dme
斑马鱼: dre | 拟南芥: ath | 酿酒酵母: sce

五、可视化¶

5.1 条形图（Barplot）¶

# 条形图 — 最基础的展示方式
barplot(go_bp,
        showCategory = 15,     # 显示前15个条目
        title = "GO-BP Enrichment") +  # 标题
  theme_minimal()              # 简洁主题

5.2 气泡图（Dotplot）¶

# 气泡图 — 同时展示基因比例、p值和基因数
dotplot(go_bp,
        showCategory = 15,     # 显示前15个条目
        title = "GO-BP Dotplot") +
  theme(axis.text.y = element_text(size = 10))  # 调整Y轴文字大小

5.3 基因-概念网络图（cnetplot）¶

# 网络图 — 展示基因和GO条目之间的关联
cnetplot(go_bp,
         showCategory = 5,         # 显示前5个条目
         categorySize = "pvalue",   # 节点大小按p值
         foldChange = NULL,         # 可选：传入log2FC上色
         colorEdge = TRUE)          # 边着色

5.4 富集地图（emapplot）¶

# 先计算GO条目之间的相似性
go_bp_sim <- pairwise_termsim(go_bp)  # 计算术语间相似性

# 富集地图 — 展示GO条目之间的关系
emapplot(go_bp_sim,
         showCategory = 20,    # 显示前20个
         cex.params = list(category_label = 0.8))  # 标签大小

5.5 KEGG通路图¶

# 安装pathview包
BiocManager::install("pathview")  # 安装pathview
library(pathview)                  # 加载

# 准备基因表达数据（需要Entrez ID为名字的向量）
gene_fc <- gene_entrez$ENTREZID  # 实际使用时用log2FC值
names(gene_fc) <- gene_entrez$ENTREZID

# 绘制特定通路图（如 hsa04110 = 细胞周期）
pathview(
  gene.data = gene_fc,     # 基因表达数据
  pathway.id = "hsa04110", # KEGG通路ID
  species = "hsa",         # 物种
  out.suffix = "cellcycle"  # 输出文件后缀
)
# 结果保存为 hsa04110.cellcycle.png

六、高级用法¶

6.1 比较不同组的富集结果¶

# 准备多组基因列表
gene_list <- list(
  upregulated = up_genes_entrez,    # 上调基因
  downregulated = down_genes_entrez  # 下调基因
)

# 比较GO富集
compare_go <- compareCluster(
  geneCluster = gene_list,  # 多组基因
  fun = "enrichGO",         # 使用enrichGO函数
  OrgDb = org.Hs.eg.db,
  ont = "BP",
  pvalueCutoff = 0.05
)

# 可视化比较结果
dotplot(compare_go, showCategory = 10)  # 分面气泡图

6.2 自定义基因集富集¶

# 使用自定义基因集（GMT格式）
# 适用于非模式生物或自定义通路
library(gson)  # 加载gson包

# 读取GMT文件
custom_gs <- read.gmt("my_pathways.gmt")  # 读取自定义基因集

# 运行富集分析
custom_enrich <- enricher(
  gene = gene_entrez$ENTREZID,  # 差异基因
  TERM2GENE = custom_gs         # 自定义基因集
)

常见报错与解决¶

报错信息	原因	解决方案
`No gene can be mapped`	基因ID类型不对	确认keyType参数，用bitr()先转换
`No enrichment found`	基因太少或阈值太严	放宽pvalueCutoff，增加基因数量
`KEGG website not available`	网络问题	检查网络，或使用本地KEGG数据
`org.Hs.eg.db not found`	未安装注释包	BiocManager::install("org.Hs.eg.db")
`Error in bitr: duplicated`	输入基因有重复	unique()去重后再运行
`Column 'ONTOLOGY' not found`	ont参数写错	用 "BP"/"CC"/"MF"/"ALL"

速查表¶

# 安装
BiocManager::install("clusterProfiler")

# GO富集三步走
bitr() → enrichGO() → dotplot()

# KEGG富集三步走
bitr() → enrichKEGG() → dotplot()

# 常用可视化函数
barplot()   — 条形图
dotplot()   — 气泡图（最常用）
cnetplot()  — 基因-概念网络
emapplot()  — 富集地图
heatplot()  — 热图
pathview()  — KEGG通路图

# 结果表关键列
ID          — GO/KEGG条目ID
Description — 功能描述
GeneRatio   — 差异基因中命中比例
BgRatio     — 背景基因中命中比例
pvalue      — 原始p值
p.adjust    — BH校正后p值
count       — 命中基因数

面试高频问题¶

Q1：GO富集分析的统计原理是什么？¶

答：ORA（Over-Representation Analysis）基于超几何检验（Fisher精确检验）。假设背景有N个基因，其中M个属于某GO条目；你的差异基因有n个，其中k个属于该条目。超几何检验计算"从N中随机抽n个，恰好≥k个属于该条目"的概率。p值越小说明富集越显著（不是随机碰巧的）。

Q2：为什么需要多重检验校正？¶

答：GO数据库有上万个条目，同时做上万次假设检验，按5%的假阳性率也会产生几百个假阳性结果。BH（Benjamini-Hochberg）方法控制FDR（False Discovery Rate），确保报告的显著结果中假阳性比例不超过设定阈值（通常5%）。所以要看p.adjust而非pvalue。

Q3：背景基因集怎么选？¶

答：背景基因集（universe）应该是"所有被检测到的基因"，而不是"全基因组所有基因"。比如RNA-seq分析时，应该用所有表达的基因（如FPKM>1的）作为背景。如果用全基因组做背景，会引入偏差——很多未表达的基因永远不可能出现在差异基因列表中，相当于"不考试的学生也算入总人数"。

Q4：enrichGO和enrichKEGG的主要区别？¶

答：①数据来源不同——GO用本地注释数据库（org.Hs.eg.db），KEGG实时在线查询；②分类体系不同——GO有BP/CC/MF三个子类的层级结构，KEGG按代谢/信号/疾病等通路分类；③GO更偏功能描述，KEGG更偏通路机制；④KEGG支持pathview画彩色通路图。

Q5：clusterProfiler 4.0有哪些重要更新？¶

答：①统一接口——支持上千物种的GO/KEGG/Reactome/WikiPathway等数据库；②2026新增interpret()函数——利用大语言模型（LLM）辅助解读富集结果；③优化了KEGG缓存机制——不再打包KEGG数据，每次在线获取最新版；④新增gson_GO_local()函数支持本地GO注释。Nature Protocols 2024年发表了详细使用教程。