KEGG通路富集分析可视化：气泡图与桑基图组合方案详解

在生信分析中KEGG通路富集分析是解读基因功能与生物过程的关键步骤。然而如何将富集结果以更直观、更具信息量的方式呈现常常让分析者感到困扰。传统的条形图或表格虽然能展示富集程度但难以同时体现通路间的层级关系、基因流向以及显著性水平。本文将介绍一种高效的可视化组合方案使用同一套富集分析数据同时生成KEGG气泡图和桑基流向图。这套流程不仅能让你的分析报告更具专业性和美观度还能从不同维度揭示数据背后的生物学故事无论是用于论文图表还是项目汇报都能显著提升表现力。本文将从R语言环境搭建开始逐步讲解如何获取KEGG注释数据、进行富集分析、使用ggplot2绘制标准气泡图并最终利用ggalluvial或networkD3包将富集结果与基因映射关系转化为精美的桑基图。整个过程代码完整、可复现适合有一定R语言基础的生物信息学初学者和希望提升可视化技能的分析人员。1. 背景与核心概念为何需要组合图表在深入代码之前我们有必要理解这两种图表各自的价值以及组合使用的意义。KEGG气泡图 (Bubble Plot/ Dot Plot)是富集分析结果最常用的可视化方法之一。它的每个气泡代表一个富集的通路或GO条目。X轴通常表示富集因子Enrichment Factor或基因比率Gene RatioY轴是通路名称。气泡的大小代表映射到该通路的基因数量或差异基因数量颜色则用于表示富集的显著性水平如P值或校正后的Q值。气泡图能一目了然地展示哪些通路最显著、影响最大。桑基图 (Sankey Diagram)是一种流图它通过“流”的宽度来显示数据在多个维度或节点之间的转移或分配情况。在生信语境下我们可以将“差异基因”作为源节点将“富集的KEGG通路”作为目标节点连接线的粗细代表有多少个基因共同映射到某个通路上。桑基图能清晰揭示基因的多功能性一个基因可能参与多个通路桑基图可以展示这种复杂的多对多关系。通路的核心基因哪些基因是多个关键通路共有的枢纽。数据的整体流向从基因集合到功能模块的宏观分布。组合使用的优势气泡图擅长展示通路的“重要性排序”而桑基图擅长展示基因与通路之间的“归属网络”。将两者结合一份分析既能回答“哪些通路最重要”气泡图也能回答“是哪些基因驱动了这些重要通路它们之间有何关联”桑基图。这实现了从宏观统计到微观关联的全方位解读。2. 环境准备与R包安装本教程基于R语言进行请确保你已安装R建议版本4.0以上和RStudio。我们将使用一系列强大的R包请按顺序安装。# 设置CRAN镜像加速安装国内用户建议使用 options(repos c(CRAN https://mirrors.tuna.tsinghua.edu.cn/CRAN/)) # 安装生物信息学核心分析包 if (!requireNamespace(BiocManager, quietly TRUE)) install.packages(BiocManager) # 通过BiocManager安装所需生物信息学包 BiocManager::install(c(clusterProfiler, org.Hs.eg.db, DOSE, enrichplot)) # 说明 # - clusterProfiler: 富集分析核心包支持GO/KEGG等。 # - org.Hs.eg.db: 人类基因注释数据库。如果你是其他物种请替换如 org.Mm.eg.db小鼠。 # - DOSE: 用于语义相似性计算和可视化。 # - enrichplot: 提供丰富的富集结果可视化函数包括气泡图。 # 安装数据处理与可视化包 install.packages(c(tidyverse, ggalluvial, networkD3, viridis)) # 说明 # - tidyverse: 包含dplyr, tidyr, ggplot2等数据处理和绘图的基础。 # - ggalluvial: 基于ggplot2的桑基图绘制包语法与ggplot2一致易于上手。 # - networkD3: 生成交互式桑基图可输出为HTML。 # - viridis: 提供美观且色盲友好的配色方案。安装完成后加载所有必要的包。library(clusterProfiler) library(org.Hs.eg.db) library(enrichplot) library(tidyverse) library(ggalluvial) library(networkD3) library(viridis)3. 数据准备模拟差异基因列表与KEGG富集分析为了演示我们首先生成一个模拟的差异表达基因列表。在实际项目中这通常来自RNA-seq等分析得到的DESeq2或limma结果。3.1 生成模拟基因列表我们使用人类基因的Entrez ID进行模拟。# 从数据库中获取所有Entrez ID并随机抽取300个作为我们的“差异基因” all_genes - keys(org.Hs.eg.db, keytype ENTREZID) set.seed(123) # 设置随机种子保证结果可重复 diff_genes - sample(all_genes, 300) head(diff_genes) # 查看前几个基因ID3.2 进行KEGG通路富集分析使用clusterProfiler的enrichKEGG函数进行分析。注意在线分析需要网络连接。# 执行KEGG富集分析 kegg_enrich - enrichKEGG( gene diff_genes, # 差异基因列表 organism hsa, # 人类其他物种如小鼠‘mmu’ keyType kegg, # 输入的ID类型我们用的是Entrez ID与KEGG内部ID一致 pvalueCutoff 0.05, # P值阈值 pAdjustMethod BH, # P值校正方法常用BHBenjamini-Hochberg qvalueCutoff 0.2, # Q值阈值 minGSSize 10, # 通路最小基因集大小 maxGSSize 500 # 通路最大基因集大小 ) # 查看富集结果摘要 head(kegg_enrich, n10) # 结果是一个enrichResult对象包含ID、Description、GeneRatio、BgRatio、pvalue、p.adjust、qvalue、geneID等列。重要提示如果你的基因列表是Symbol基因符号需要先转换为Entrez ID。可以使用bitr函数。# 假设你的基因列表是Symbol gene_symbols - c(TP53, “BRCA1”, “MYC”, “EGFR”, ...) # 转换 gene_df - bitr(gene_symbols, fromType “SYMBOL, toType c(“ENTREZID), OrgDb org.Hs.eg.db) diff_genes - gene_df$ENTREZID4. 绘制标准KEGG气泡图富集分析完成后我们首先用enrichplot包快速绘制气泡图。# 方法1使用enrichplot包的dotplot函数最简单 p_bubble - dotplot(kegg_enrich, showCategory 20, # 显示最显著的20个通路 font.size 10, title “KEGG Pathway Enrichment Analysis”, color “p.adjust”, # 按校正后P值着色 size “Count”) # 按基因计数决定点大小 print(p_bubble)dotplot函数非常便捷但自定义程度有限。为了获得更精美的图表并与后续桑基图数据衔接我们使用ggplot2手动绘制。4.1 提取并整理富集结果数据# 将富集结果转换为数据框 kegg_result_df - as.data.frame(kegg_enrich) # 查看数据结构 str(kegg_result_df) # 为了绘图美观我们通常需要整理数据 # 1. 计算富集因子 (GeneRatio) kegg_result_df - kegg_result_df %% separate(GeneRatio, into c(“GeneInPathway”, “TotalGenes”), sep “/”) %% separate(BgRatio, into c(“PathwaySize”, “BackgroundGenes”), sep “/”) %% mutate(across(c(GeneInPathway, TotalGenes, PathwaySize, BackgroundGenes), as.numeric)) %% mutate(EnrichmentFactor (GeneInPathway / TotalGenes) / (PathwaySize / BackgroundGenes)) %% arrange(p.adjust) # 按校正P值排序 # 2. 选择Top N个通路用于绘图例如Top 15 top_n - 15 plot_data - kegg_result_df[1:top_n, ] # 3. 对通路描述进行排序因子化保证绘图时顺序正确 plot_data$Description - factor(plot_data$Description, levels rev(plot_data$Description)) # rev()使最重要的在顶部4.2 使用ggplot2绘制自定义气泡图p_custom_bubble - ggplot(plot_data, aes(x EnrichmentFactor, y Description)) geom_point(aes(size Count, color -log10(p.adjust))) # 颜色用-log10(p.adjust)值越大越显著 scale_size_continuous(range c(3, 8), name “Gene Count”) # 控制气泡大小范围 scale_color_viridis(option “C”, begin 0.3, end 0.9, name “-log10(Adj.P)”) # 使用viridis配色 labs(x “Enrichment Factor”, y NULL, title “Top Enriched KEGG Pathways”, subtitle “Bubble size: number of genes; Color: significance level”) theme_minimal(base_size 12) theme(axis.text.y element_text(color “black”, size 10), panel.grid.major.y element_line(linetype “dashed”, color “grey90”), legend.position “right”) print(p_custom_bubble) # 保存图片 ggsave(“KEGG_bubble_plot.png”, p_custom_bubble, width 10, height 7, dpi 300)至此我们得到了一张高度定制化的KEGG气泡图。接下来我们将利用同一份kegg_enrich对象中的数据来构建桑基图。5. 从富集结果到桑基图数据重构桑基图需要一种特定的“长格式”数据包含“源节点”基因、“目标节点”通路以及连接它们的“流”通常用基因计数表示。我们需要从enrichResult对象中提取基因-通路的映射关系。5.1 提取基因-通路关联矩阵# 从kegg_enrich对象中提取基因与通路的对应关系 # geneID列包含了映射到每个通路的基因列表以‘/’分隔 gene_pathway_list - strsplit(kegg_enrich$geneID, “/”) # 为每个通路创建数据框记录基因与通路的对应 edges - data.frame() for (i in 1:length(gene_pathway_list)) { pathway - kegg_enrich$ID[i] genes - gene_pathway_list[[i]] temp_df - data.frame(Gene genes, Pathway pathway, stringsAsFactors FALSE) edges - rbind(edges, temp_df) } # 查看前几行 head(edges) # 输出类似 # Gene Pathway # 1 1030 hsa04110 # 2 1962 hsa04110 # 3 2072 hsa041105.2 为桑基图准备节点和连接数据桑基图需要两个数据框一个描述所有节点nodes一个描述所有连接links。# 创建节点数据框 # 节点包括所有唯一的基因和通路 gene_nodes - unique(edges$Gene) pathway_nodes - unique(edges$Pathway) # 我们需要将通路ID转换为描述以便在图中显示 pathway_names - setNames(kegg_enrich$Description, kegg_enrich$ID) pathway_nodes_named - pathway_names[pathway_nodes] # 构建节点数据框包含节点ID和名称 nodes_df - data.frame( node_id 0:(length(gene_nodes) length(pathway_nodes) - 1), # 从0开始编号 node_name c(gene_nodes, pathway_nodes_named), group c(rep(“Gene”, length(gene_nodes)), rep(“Pathway”, length(pathway_nodes))) ) # 创建连接数据框 # 需要将基因名和通路名转换为对应的节点ID links_df - edges %% mutate(source match(Gene, nodes_df$node_name) - 1, # networkD3要求索引从0开始 target match(Pathway, nodes_df$node_name) - 1) %% group_by(source, target) %% summarise(value n(), .groups ‘drop’) # value代表连接权重这里每个基因-通路对记为1汇总后即为基因数 head(links_df)6. 绘制桑基流向图我们介绍两种方法静态的ggalluvial和交互式的networkD3。6.1 方法一使用ggalluvial绘制静态桑基图ggalluvial语法与ggplot2一致易于集成到图形组合中。# 首先我们需要将edges数据转换为ggalluvial需要的格式每个基因-通路对为一行的数据框 # 同时我们需要通路描述而不是ID sankey_data - edges %% left_join(kegg_result_df[, c(“ID”, “Description”)], by c(“Pathway” “ID”)) %% select(Gene, Pathway Description) %% # 为了图形可读性通常只展示与Top通路相关的基因 filter(Pathway %in% plot_data$Description) # 计算每个基因连接到Top通路的总数用于排序可选 gene_freq - sankey_data %% count(Gene, name “Freq”) %% arrange(desc(Freq)) # 将基因按频率排序使图形更有序 sankey_data$Gene - factor(sankey_data$Gene, levels gene_freq$Gene) # 绘制桑基图 p_sankey_static - ggplot(sankey_data, aes(axis1 Gene, axis2 Pathway)) geom_alluvium(aes(fill Pathway), width 1/12, alpha 0.7) geom_stratum(width 1/12, fill “grey80”, color “grey”) geom_text(stat “stratum”, aes(label after_stat(stratum)), size 3) scale_x_discrete(limits c(“Genes”, “Pathways”), expand c(0.05, 0.05)) scale_fill_viridis_d(option “plasma”) labs(title “Gene-Pathway Sankey Diagram (Top Pathways)”, subtitle “Shows the flow of genes from the differential list into enriched KEGG pathways”) theme_minimal() theme(legend.position “none”, # 图例可能太复杂隐藏 axis.text.y element_blank(), axis.ticks element_blank(), panel.grid element_blank()) print(p_sankey_static) ggsave(“Sankey_static_plot.png”, p_sankey_static, width 14, height 10, dpi 300)6.2 方法二使用networkD3绘制交互式桑基图交互式桑基图允许鼠标悬停查看详细信息非常适合在网页报告中展示。# 使用之前准备好的nodes_df和links_df sankeyNetwork(Links links_df, # 连接数据框 Nodes nodes_df, # 节点数据框 Source “source”, Target “target”, Value “value”, NodeID “node_name”, NodeGroup “group”, # 按组基因/通路着色 units “Genes”, fontSize 12, nodeWidth 30, height 600, width 900, sinksRight FALSE) # 让最右侧的节点通路左对齐布局更合理运行这行代码会在RStudio的Viewer窗口生成一个交互式图表。你可以使用saveNetwork函数将其保存为独立的HTML文件。# 将图表保存为HTML文件 sn - sankeyNetwork(Links links_df, Nodes nodes_df, Source “source”, Target “target”, Value “value”, NodeID “node_name”, NodeGroup “group”, fontSize 12, nodeWidth 30) saveNetwork(sn, “Interactive_Sankey.html”)7. 组合与优化将两图整合到分析报告中在实际报告中我们通常将气泡图和桑基图并列展示。可以使用patchwork或cowplot包轻松实现。# 安装并加载patchwork包 # install.packages(“patchwork”) library(patchwork) # 组合图表气泡图在上桑基图在下 combined_plot - p_custom_bubble / p_sankey_static plot_layout(heights c(1, 1.5)) # 调整上下两部分的高度比例 combined_plot # 保存组合图 ggsave(“Combined_KEGG_Visualization.png”, combined_plot, width 16, height 14, dpi 300)8. 常见问题与排查思路在实践过程中你可能会遇到以下问题问题现象可能原因解决思路enrichKEGG报错“InternetOpenUrl failed:...”或长时间无响应网络连接问题或KEGG官网访问不稳定。1. 检查网络。2. 使用use_internal_data TRUE参数但数据可能非最新。3. 考虑使用clusterProfiler的enricher函数配合自定义的KEGG背景基因集需提前下载。气泡图中通路描述文字过长重叠显示通路名称太长。1. 使用stringr::str_wrap函数在绘图数据中截断或换行。2. 调整图形尺寸width,height。3. 调整theme中的axis.text.y的size和hjust。桑基图节点过多图形杂乱无法辨认展示的基因和通路太多。1.严格过滤只展示p.adjust最显著的前10-15个通路及其相关基因。2. 在ggalluvial中可以过滤掉只映射到1个通路的基因简化图形。3. 使用交互式networkD3图通过缩放和拖拽查看细节。ggalluvial绘图时提示“No alluvia to plot”数据格式不符合ggalluvial要求。确保数据框的每一行代表一个独立的“流动单元”如一个基因-通路对。检查aes中的axis1和axis2列名是否正确。networkD3桑基图节点名称显示不全节点名称太长或图形区域太小。1. 在准备nodes_df时对过长的通路描述进行缩写。2. 调整fontSize和图形width/height参数。3. 鼠标悬停时可以看到完整名称。基因ID转换失败 (bitr函数返回空)输入的基因标识符类型错误或不在数据库中。1. 确认fromType参数是否正确如“SYMBOL”,“ENSEMBL”。2. 检查基因标识符的版本和格式是否与数据库匹配。3. 使用library(org.XX.eg.db)后用columns(org.XX.eg.db)查看支持的ID类型。9. 最佳实践与工程建议数据溯源与记录始终在R脚本开头使用set.seed()保证随机过程可重复。使用sessionInfo()记录所有包版本这对于生信分析的可复现性至关重要。结果过滤策略不要盲目展示所有富集结果。结合生物学意义、显著性p.adjust和富集因子进行综合筛选。通常关注p.adjust 0.05且Count 2的通路。图形定制与美化配色坚持使用色盲友好配色如viridis、RColorBrewer的Set2、Set3。标签桑基图中若基因过多可考虑只标注高频基因或枢纽基因或用“Other Genes”聚合低频基因。布局静态桑基图ggalluvial中可以通过对基因和通路节点进行排序如按连接度、按字母顺序来改善可读性。性能考虑当差异基因数量很大1000时桑基图的边和节点会急剧增加导致渲染缓慢或图形混乱。务必在数据准备阶段进行有效聚合和过滤。输出格式对于出版物输出PDF或高分辨率PNGdpi600或更高。对于网页报告或演示交互式HTMLnetworkD3是更好的选择。代码模块化将数据准备、富集分析、气泡图绘制、桑基图数据重构、图形绘制分别写成函数。这样当你分析新的数据集时只需更换输入基因列表即可快速生成全套图表。掌握这套“一套数据双图呈现”的流程不仅能提升你的生信数据分析效率更能让你的研究成果以更专业、更深刻的视觉形式展现出来。从理解富集分析的统计结果到洞察基因与通路间的复杂网络R语言提供了强大而灵活的工具链。建议读者在理解本文代码的基础上尝试将自己的差异表达分析结果代入并进一步探索ggplot2和ggalluvial的主题theme系统定制出具有个人或实验室风格的图表模板。