在生信分析中，KEGG通路富集分析是解读基因功能与生物过程的关键步骤。然而，如何将富集结果以更直观、更具信息量的方式呈现，常常让分析者感到困扰。传统的条形图或表格虽然能展示富集程度，但难以同时体现通路间的层级关系、基因流向以及显著性水平。本文将介绍一种高效的可视化组合方案：使用同一套富集分析数据，同时生成KEGG气泡图和桑基流向图。这套流程不仅能让你的分析报告更具专业性和美观度，还能从不同维度揭示数据背后的生物学故事，无论是用于论文图表还是项目汇报，都能显著提升表现力。

本文将从R语言环境搭建开始，逐步讲解如何获取KEGG注释数据、进行富集分析、使用ggplot2绘制标准气泡图，并最终利用ggalluvial或networkD3包将富集结果与基因映射关系转化为精美的桑基图。整个过程代码完整、可复现，适合有一定R语言基础的生物信息学初学者和希望提升可视化技能的分析人员。

1. 背景与核心概念：为何需要组合图表？

在深入代码之前，我们有必要理解这两种图表各自的价值以及组合使用的意义。

KEGG气泡图 (Bubble Plot/ Dot Plot)是富集分析结果最常用的可视化方法之一。它的每个气泡代表一个富集的通路（或GO条目）。X轴通常表示富集因子（Enrichment Factor）或基因比率（Gene Ratio），Y轴是通路名称。气泡的大小代表映射到该通路的基因数量（或差异基因数量），颜色则用于表示富集的显著性水平（如P值或校正后的Q值）。气泡图能一目了然地展示哪些通路最显著、影响最大。

桑基图 (Sankey Diagram)是一种流图，它通过“流”的宽度来显示数据在多个维度（或节点）之间的转移或分配情况。在生信语境下，我们可以将“差异基因”作为源节点，将“富集的KEGG通路”作为目标节点，连接线的粗细代表有多少个基因共同映射到某个通路上。桑基图能清晰揭示：

1. 基因的多功能性：一个基因可能参与多个通路，桑基图可以展示这种复杂的多对多关系。
2. 通路的核心基因：哪些基因是多个关键通路共有的枢纽。
3. 数据的整体流向：从基因集合到功能模块的宏观分布。

组合使用的优势：气泡图擅长展示通路的“重要性排序”，而桑基图擅长展示基因与通路之间的“归属网络”。将两者结合，一份分析既能回答“哪些通路最重要？”（气泡图），也能回答“是哪些基因驱动了这些重要通路，它们之间有何关联？”（桑基图）。这实现了从宏观统计到微观关联的全方位解读。

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) # 查看前几个基因ID

3.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值校正方法，常用BH（Benjamini-Hochberg） qvalueCutoff = 0.2, # Q值阈值 minGSSize = 10, # 通路最小基因集大小 maxGSSize = 500 # 通路最大基因集大小 ) # 查看富集结果摘要 head(kegg_enrich, n=10) # 结果是一个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$ENTREZID

4. 绘制标准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 hsa04110

5.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. 常见问题与排查思路

在实践过程中，你可能会遇到以下问题：

9. 最佳实践与工程建议

1. 数据溯源与记录：始终在R脚本开头使用set.seed()保证随机过程可重复。使用sessionInfo()记录所有包版本，这对于生信分析的可复现性至关重要。
2. 结果过滤策略：不要盲目展示所有富集结果。结合生物学意义、显著性（p.adjust）和富集因子进行综合筛选。通常关注p.adjust < 0.05且Count > 2的通路。
3. 图形定制与美化：
   • 配色：坚持使用色盲友好配色（如viridis、RColorBrewer的Set2、Set3）。
   • 标签：桑基图中若基因过多，可考虑只标注高频基因或枢纽基因，或用“Other Genes”聚合低频基因。
   • 布局：静态桑基图（ggalluvial）中，可以通过对基因和通路节点进行排序（如按连接度、按字母顺序）来改善可读性。
4. 性能考虑：当差异基因数量很大（>1000）时，桑基图的边和节点会急剧增加，导致渲染缓慢或图形混乱。务必在数据准备阶段进行有效聚合和过滤。
5. 输出格式：对于出版物，输出PDF或高分辨率PNG（dpi=600或更高）。对于网页报告或演示，交互式HTML（networkD3）是更好的选择。
6. 代码模块化：将数据准备、富集分析、气泡图绘制、桑基图数据重构、图形绘制分别写成函数。这样，当你分析新的数据集时，只需更换输入基因列表即可快速生成全套图表。

掌握这套“一套数据，双图呈现”的流程，不仅能提升你的生信数据分析效率，更能让你的研究成果以更专业、更深刻的视觉形式展现出来。从理解富集分析的统计结果，到洞察基因与通路间的复杂网络，R语言提供了强大而灵活的工具链。建议读者在理解本文代码的基础上，尝试将自己的差异表达分析结果代入，并进一步探索ggplot2和ggalluvial的主题（theme）系统，定制出具有个人或实验室风格的图表模板。