本节概览：
1.获取DEG结果的上下调差异基因
2.bitr()函数转化基因名为entrez ID
3.利用clusterProfiler进行KEGG与GO富集
4.用enrichplot进行富集结果可视化：pathview goplot barplot dotplot cnetplot emapplot treeplot heatplot upsetplot
ps:本次分析示例文件可在关注后私信领取~
承接上期文章：2026版RNA-seq入门实战（六）：差异分析——DESeq2 edgeR limma的使用与比较 - 简书
1. 获取DEG结果的上下调差异基因
载入上节中保存的三种差异分析结果数据，这里示范选取DESeq2的结果数据，进行筛选条件设置后获取上下调基因名
# 环境准备
rm(list=ls()) # 清空当前工作环境中的所有对象
options(stringsAsFactors = FALSE) # 禁止字符串自动转换为因子
# 加载富集分析相关R包
library(clusterProfiler) # 核心富集分析工具
library(enrichplot) # 富集结果可视化
library(tidyverse) # 数据处理与可视化套件
library(org.Mm.eg.db) # 小鼠基因注释数据库（本示例为小鼠）
library(DOSE) # 疾病本体语义富集分析
library(pathview) # KEGG通路可视化
#数据载入与工作目录设置
# setwd("目的文件夹") # 根据实际项目路径设置工作目录（示例被注释掉）
# 加载基因表达计数矩阵
load(file = '1.counts.Rdata')
# 自动查找并加载指定文件夹下的差异分析结果文件（文件名包含"DEG_results.Rdata"）
load(list.files(path = "./3.DEG", pattern = 'DEG_results.Rdata', full.names = TRUE))
# 创建富集分析输出文件夹，若已存在则跳过
dir.create("4.Enrichment_KEGG_GO")
setwd("4.Enrichment_KEGG_GO") # 进入该文件夹，后续结果将保存在此
# 差异基因筛选阈值设置
log2FC_cutoff <- log2(10) # log2FoldChange阈值（对应原始差异倍数10倍）
pvalue_cutoff <- 0.05 # 原始p值阈值
padj_cutoff <- 0.05 # 多重检验校正后的p值（padj）阈值
# 提取上调和下调的差异基因
# 选择DESeq2差异分析结果中的关键列：log2FoldChange、pvalue、padj
# 根据实际情况可替换为其他工具结果（如DEG_limma_voom[,c(1,4,5)] 或 DEG_edgeR[,c(1,4,5)]）
need_DEG <- DEG_DEseq2[, c(2,5,6)]
head(need_DEG) # 预览前几行，检查数据格式
# 将列名统一为易于理解的名称
colnames(need_DEG) <- c('log2FoldChange', 'pvalue', 'padj')
# 筛选上调基因：log2FC > 阈值 且 p值 < 阈值 且 矫正p值 < 阈值
gene_up <- rownames(need_DEG[
with(need_DEG, log2FoldChange > log2FC_cutoff &
pvalue < pvalue_cutoff &
padj < padj_cutoff),
])
# 筛选下调基因：log2FC < -阈值 且 p值 < 阈值 且 矫正p值 < 阈值
gene_down <- rownames(need_DEG[
with(need_DEG, log2FoldChange < -log2FC_cutoff &
pvalue < pvalue_cutoff &
padj < padj_cutoff),
])
# 注：
# gene_up / gene_down 为基因ID向量，可用于后续GO/KEGG富集分析
# 此处阈值设置较为宽松（FC>10, padj<0.05），可根据生物学问题调整
2. 转化基因名为entrez ID
在执行差异基因富集分析之前，需先将基因标识符统一为 Entrez ID。
进行转换前，要加载对应的物种注释数据库（如 human 用 org.Hs.eg.db，mouse 用 org.Mm.eg.db），这些库整合了多个主流数据库的标识信息，包含 Entrez ID、Ensembl 等类型。可通过 keytypes(org.Mm.eg.db) 查看该库支持的全部可转换 ID 类型。
通常使用 clusterProfiler 包中的 bitr() 或 mapIds() 函数来完成标识转换，典型用法如下：
# 基因ID转换：SYMBOL → ENTREZID
# 使用clusterProfiler的bitr函数将基因名从SYMBOL转换为ENTREZID
# 根据物种选择合适的OrgDb包：
# 人类 -> org.Hs.eg.db
# 小鼠 -> org.Mm.eg.db
# keytypes(org.Hs.eg.db) # 可查看该数据库支持的所有ID类型，如 ENTREZID, ENSEMBL, SYMBOL 等
# 上调基因ID转换
gene_up_entrez <- as.character(na.omit(
bitr(gene_up, # 待转换的基因向量（SYMBOL格式）
fromType = "SYMBOL", # 输入ID类型
toType = "ENTREZID", # 目标ID类型
OrgDb = "org.Mm.eg.db" # 小鼠基因组注释数据库
)[, 2] # 提取转换后的ENTREZID列
))
# 下调基因ID转换（流程同上）
gene_down_entrez <- as.character(na.omit(
bitr(gene_down,
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = "org.Mm.eg.db"
)[, 2]
))
# 合并全部差异基因（上、下调）的ENTREZID，并去重
gene_diff_entrez <- unique(c(gene_up_entrez, gene_down_entrez))
# 注：gene_diff_entrez 将用于后续富集分析，包含所有显著差异基因（不论方向）
3. 利用clusterProfiler进行KEGG与GO富集
在拿到差异基因（上调和下调）的 Entrez ID 后，就可以借助 clusterProfiler 包中的 enrichKEGG() 和 enrichGO() 来进行功能富集。这里只以上调基因为例演示富集分析的具体步骤；在实际研究里，建议对上调、下调以及全体差异基因分别进行富集分析，以便全方位考察结果。在分析过程中，需要留意以下几点：
- pvalueCutoff 和 qvalueCutoff 的默认值分别是 0.05 和 0.2，实际分析中可根据需求灵活调整，以控制显著性和多重检验矫正的严格程度。
- 在使用 enrichGO() 时，必须通过 OrgDb 指定物种注释库（如小鼠用 "org.Mm.eg.db"）；设置 readable = TRUE 可自动将 Entrez ID 转换为基因 Symbol；ont 参数可选择 "BP"（生物过程）、"MF"（分子功能）、"CC"（细胞组分）或 "ALL"（全部三类），分别对感兴趣的本体类别进行富集。
- enrichKEGG() 中需要用 organism 指定物种（小鼠为 "mmu"），它没有 readable 选项，因此后续需借助 DOSE 包的 setReadable() 将结果中的 Entrez ID 转换成更可读的基因 Symbol。
#### KEGG、GO富集 ####
# KEGG 通路富集分析（上调基因）
kegg_enrich_results <- enrichKEGG(
gene = gene_up_entrez, # 输入ENTREZID格式的上调基因
organism = "mmu", # 物种缩写：人类 "hsa"，小鼠 "mmu"
pvalueCutoff = 0.05, # p值阈值
qvalueCutoff = 0.05 # 多重检验校正后的q值阈值
)
# 将ENTREZID替换为基因Symbol，方便结果解读
kk_read <- DOSE::setReadable(
kegg_enrich_results, # 富集结果对象
OrgDb = "org.Mm.eg.db", # 小鼠基因组注释数据库
keyType = "ENTREZID" # 指明当前基因ID类型为ENTREZID
) # 结果中会在ENTREZID列基础上新增SYMBOL等可读信息
# 输出KEGG富集结果表格
write.csv(kk_read@result, 'KEGG_gene_up_enrichresults.csv')
# 保存原始结果对象，便于后续可视化
save(kegg_enrich_results, file = 'KEGG_gene_up_enrichresults.Rdata')
#GO 富集分析（上调基因）
go_enrich_results <- enrichGO(
gene = gene_up_entrez, # 使用相同的上调基因
OrgDb = "org.Mm.eg.db", # 小鼠数据库
ont = "ALL", # GO本体类型："BP"（生物过程），"MF"（分子功能），
# "CC"（细胞组分），"ALL"（包含全部三类）
pvalueCutoff = 0.05, # p值阈值
qvalueCutoff = 0.05, # q值阈值
readable = TRUE # 自动将ENTREZID转换为基因Symbol
)
# 输出GO富集结果表格（文件名中“BP”为示例，实际包含全部本体）
write.csv(go_enrich_results@result, 'GO_gene_up_BP_enrichresults.csv')
# 保存GO富集结果对象
save(go_enrich_results, file = 'GO_gene_up_enrichresults.Rdata')
4. 富集结果可视化：pathview goplot barplot dotplot cnetplot emapplot treeplot heatplot upsetplot
富集通路后就要进行可视化展示了，参见clusterprofiler说明书 Introduction | Biomedical Knowledge Mining using GOSemSim and clusterProfiler (yulab-smu.top)，其中enrichplot包可以对富集结果进行超级丰富的可视化。
下面就来总结介绍一下clusterprofiler说明书介绍的各种富集结果可视化方法，其中4.1pathview 只针对KEGG， 4.2goplot只针对GO结果，其他可视化方法则对两者都通用
4.1 pathview ——KEGG通路可视化
将KEGG通路进行可视化一般有以下三种方法：
*使用函数 browseKEGG(enrich_results, select_pathway)进行网页查看相关通路，如
browseKEGG(kegg_enrich_results, 'mmu04310') #网页查看通路
-
网页版的pathview https://pathview.uncc.edu/guest-home，可以上传数据进行在线可视化
pathview包：
在使用 pathview() 时，要给出含 log2 倍变化信息的 gene.data、目标通路的 pathway.id 以及 species，函数据此生成一幅在通路上标注了基因上、下调方向的通路图。
需要特别留意参数 kegg.native：当设为 TRUE 时，会生成一幅包含完整 KEGG 通路的 PNG 图片；若设为 FALSE，则仅输出一张仅展示输入基因信息的 PDF 文件。
参数limit可以对图例color bar范围（即log2FC范围）进行调整。
其他参数使用详见官方说明：pathview.pdf (bioconductor.org)，再推荐一篇简要中文教程：可视化kegg通路-pathview包 | KeepNotes blog (bioinfo-scrounger.com)
下面我们选取富集排名第3（富集结果是按pvalue由小到大排序的）进行示范。
# 构建含 log2FoldChange 信息的基因列表
# 从差异分析结果中提取 log2FoldChange 值
# DEG_DEseq2 为之前获得的 DESeq2 差异分析完整结果，第2列为 log2FoldChange
genelist <- as.numeric(DEG_DEseq2[, 2]) # 提取 log2FC 值作为数值向量
names(genelist) <- row.names(DEG_DEseq2) # 用基因SYMBOL命名向量元素
# 基因列表ID转换：SYMBOL → ENTREZID
genelist_entrez <- genelist # 保留原始 log2FC 值
# 使用 bitr 将基因名从 SYMBOL 转换为 ENTREZID，作为新向量的名称
names(genelist_entrez) <- bitr(
names(genelist), # 输入基因SYMBOL
fromType = "SYMBOL", # 输入ID类型
toType = "ENTREZID", # 目标ID类型
OrgDb = "org.Mm.eg.db" # 小鼠数据库
)[, 2] # 提取 ENTREZID 列
# 移除转换失败（NA）的基因
genelist_entrez <- genelist_entrez[!is.na(names(genelist_entrez))]
# 此时 genelist_entrez 为一个命名数值向量，names 为 ENTREZID，值为 log2FC
# 查看并选择目标通路
# 需要先运行上一步的 KEGG 富集分析，得到 kk_read 对象
kk_read@result$Description[1:10] # 查看前10条富集通路的描述
i <- 3 # 选择第3条通路作为演示
select_pathway <- kk_read@result$ID[i] # 获取该通路的 KEGG ID（如 mmu04512）
# 使用 pathview 绘制 KEGG 通路图
# 方式一：输出完整 KEGG 通路图（png格式），基因表达变化以颜色标注在通路图上
pathview(
gene.data = genelist_entrez, # 含 log2FC 的命名向量
pathway.id = select_pathway, # 目标通路ID
species = 'mmu', # 物种：小鼠 mmu，人类 hsa
kegg.native = TRUE, # TRUE：输出原生KEGG通路PNG图
new.signature = FALSE, # PDF相关参数，此处无效
limit = list(gene = 2.5, cpd = 1) # 设置颜色标尺范围，gene上限2.5，化合物上限1
)
# 方式二：输出基因列表的 PDF 文件（仅含差异基因在通路上位置，非完整通路）
pathview(
gene.data = genelist_entrez,
pathway.id = select_pathway,
species = 'mmu',
kegg.native = FALSE, # FALSE：输出基因/化合物的简化PDF展示
new.signature = FALSE, # PDF中是否添加自定义签名/标注
limit = list(gene = 2.5, cpd = 1)
)
参数kegg.native = T，所得如下：
参数kegg.native = F，所得如下：
4.2 goplot—— GO富集结果的有向无环图 directed acyclic graph
注意当enrichGO()的ont为'ALL'时不能运行该图，会报错。以下 goplot展示的是enrichGO()的ont='BP'时的go_enrich_results
### goplot : Gene Ontology is organized as a directed acyclic graph.有向无环图
gop <- goplot(go_enrich_results, showCategory = 10)
ggsave(gop, filename = "goplot.pdf", width=10, height=10)
4.3 barplot与dotplot——最常用的可视化图形
barplot 与 dotplot 主要用于展示富集通路的校正 p 值（p.adj）、基因比率（GeneRatio）以及计数（count）等信息。
如果将 enrichGO 的 ont 设为 "ALL"，还可以在这两类图中按不同本体（ONTOLOGY）将通路分面，使 BP、CC、MF 的结果分别罗列，便于对比。
#绘制 GO 富集分析条形图（barplot）
### barplot
# 基础条形图：展示前10个最显著的GO条目，字体大小14
barp <- barplot(go_enrich_results, font.size=14, showCategory=10)+
theme(plot.margin=unit(c(1,1,1,1),'lines')) # 调整图形四周边距，防止文字被截断
# 如果enrichGO的ont为'ALL'（即同时包含BP、MF、CC三类本体），则使用分面展示
if (T) { # 此处始终为TRUE，因为上游设定 ont="ALL"，强制采用分面
barp <- barplot(go_enrich_results, split= "ONTOLOGY", font.size =14)+
facet_grid(ONTOLOGY~., scale="free")+ # 按ONTOLOGY纵向分面，每类本体x轴范围自适应
theme(plot.margin=unit(c(1,1,1,1),'lines')) #调整图形四周留白大小
}
# 保存条形图为PDF，文件名“barplot.pdf”，尺寸宽10高15英寸
ggsave(barp,filename = paste0("barplot",'.pdf'), width=10, height=15)
# 绘制 GO 富集分析气泡图（dotplot）
### dotplot
# 基础气泡图：展示所有显著富集的GO条目，字体14，图例键尺寸10pt
dotp <- enrichplot::dotplot(go_enrich_results,font.size =14)+
theme(legend.key.size = unit(10, "pt"),#调整图例大小，使其更精致
plot.margin=unit(c(1,1,1,1),'lines'))#调整四周留白大小，避免标签拥挤
# 同样，若GO结果包含多个ONTOLOGY，则按ONTOLOGY分面展示
if (T) {
dotp <- enrichplot::dotplot(go_enrich_results,font.size =8,split = 'ONTOLOGY')+ # 分面后字体可适当调小
facet_grid(ONTOLOGY~., scale="free")+ # 纵向分面，x轴（GeneRatio等）范围独立
theme(legend.key.size = unit(10, "pt"), #调整图例大小
plot.margin=unit(c(1,1,1,1),'lines'))#调整四周留白大小
}
# 保存气泡图为PDF，文件名“dotplot.pdf”，尺寸宽10高12英寸
ggsave(dotp,filename = paste0("dotplot",'.pdf'),width =10,height =12)
4.4 cnetplot——Gene-Concept Network
cnetplot展示了各通路下的基因网络。绘制cnetplot有两种展现方式， 更改参数circular 为 F（默认）或T可以分别得到散布状和圈状分布的cnetplot；cnetplot还可以输入含log2FC信息的genelist ，会将log2FC信息展现在图中
# 基因-概念网络图（cnetplot）
# cnetplot 用于展示基因与富集GO条目（概念）之间的关联网络，结合 log2FoldChange 着色
#--- 构建含 log2FoldChange 信息的基因列表 ---
# 从 DESeq2 差异分析结果中提取 log2FoldChange，用于对网络图中基因节点着色
genelist <- as.numeric(DEG_DEseq2[, 2]) # 提取 log2FoldChange 数值
names(genelist) <- row.names(DEG_DEseq2) # 用基因 SYMBOL 命名
#--- 绘制网络图：所有节点均显示标签 ---
cnetp1 <- cnetp1 <- cnetplot( # 重复赋值（通常多余，保留原样）
go_enrich_results, # GO 富集分析结果对象
showCategory = 6, # 展示前6个最显著的GO条目
foldChange = genelist, # 提供基因的差异倍数，用于基因节点着色
node_label = "all" # 所有节点（基因与GO条目）均显示标签
)
#--- 绘制环形布局网络图：仅基因节点显示标签 ---
cnetp2 <- cnetp2 <- cnetplot(
go_enrich_results,
showCategory = 6,
foldChange = genelist,
node_label = "gene", # 仅对基因节点显示标签，概念节点不标
layout = "circle" # 使用圆形布局
)
# 保存网络图到 PDF 文件
ggsave(cnetp1, filename = 'cnetplot.pdf', # 默认布局的网络图
width = 12, height = 10)
ggsave(cnetp2, filename = 'cnetplot_cir.pdf', # 圆形布局的网络图
width = 15, height = 10)
4.5 emapplot ——Enrichment Map
emapplot富集图将被富集的术语组织成一个边缘连接重叠基因集的网络，相互重叠的基因集往往会聚集在一起，因此有助于我们识别功能模块。
注意使用emapplot前还需要用pairwise_termsim()处理富集结果
#富集网络图（Enrichment Map）
# 计算GO条目的两两语义相似度（默认方法为“jaccard”），构建用于网络图的相似矩阵
pt <- pairwise_termsim(go_enrich_results)
# 绘制富集网络图，节点代表GO条目，边表示相似性
emapp <- emapplot(pt,
showCategory = 30, # 显示前30个最显著的GO条目
node_label = 'category') # 节点标签显示类别名称（可选：'category', 'group', 'all', 'none'）
# 保存富集网络图为PDF，尺寸10x10英寸
ggsave(emapp, filename = 'emapplot.pdf', width = 10, height = 10)
4.6 heatplot: Heatmap-like functional classification
与cnetplot展示内容类似，但是将不同富集通路的相同基因聚在了一起，有助于我们更好识别基因模块
# 基因-功能热图（Heatplot）
# heatplot 通过热图形式展示富集条目与相关基因的关联，并用颜色反映基因的差异表达水平
#--- 构建含 log2FoldChange 信息的基因列表（用于着色）---
genelist <- as.numeric(DEG_DEseq2[, 2]) # 提取 DESeq2 结果的 log2FoldChange 值
names(genelist) <- row.names(DEG_DEseq2) # 用基因 SYMBOL 命名，确保与富集结果匹配
#--- 绘制热图 ---
heatp <- heatplot(go_enrich_results, # GO 富集分析结果对象
foldChange = genelist, # 提供基因变化倍数，用于热图颜色映射
showCategory = 5) # 显示前5个最显著的GO条目
# 保存热图为 PDF，宽度加倍以适应基因名称等细节
ggsave(heatp, filename = 'heatplot.pdf', width = 20, height = 10)
4.7 upsetplot——不同富集通路间的重叠基因数量
upsetplot展示不同富集通路间的重叠基因数量。
# UpSet 图：展示不同基因集间的基因重叠
# UpSet 图可直观显示多个基因集（此处为不同 GO 条目）所包含基因的交叉情况，
# 比传统韦恩图更适合展示大量集合间的重叠关系。
upsetp <- upsetplot(go_enrich_results, n = 10) + # 对前10个最显著富集的GO条目绘制UpSet图
theme(plot.margin = unit(c(1, 1, 1, 1), 'lines')) # 调整图形四周边距，避免文字溢出
ggsave(upsetp, filename = 'upsetplot.pdf', width = 15, height = 10) # 将UpSet图保存为PDF
4.7 treeplot——富集结果术语的层次聚类与高频词标记
treeplot 会对富集术语进行层次聚类，并用高频词加以标注，便于从繁杂的富集结果中迅速捕捉关键信息。
绘制 treeplot 之前，也要先用 pairwise_termsim() 对富集结果进行处理。
# 树图（Tree plot）：展示GO条目间的层次聚类关系
# 再次计算GO条目间的成对语义相似度（与前文富集网络图类似）
pt <- pairwise_termsim(go_enrich_results)
# 基于语义相似度矩阵绘制树状图，直观呈现富集条目的聚类结构
treep <- treeplot(pt,
showCategory = 30) # 显示前30个GO条目
# 保存树图为PDF，宽度18英寸以适应标签
ggsave(treep, filename = 'treeplot.pdf', width = 18, height = 10)
GO、KEGG富集与可视化到此就结束啦
如果所得显著差异基因很少，或无法富集到有生物学意义的通路时又该怎么办呢，下节将介绍一种不依赖于差异基因筛选的富集分析方法——GSEA
ps:本次分析示例文件可在关注后私信领取~
这里是两栖生物手册，中科院生物医学在读博士，正在努力成为贴心负责的RNA-seq分析者和实验技术分享者~
