你是不是也遇到过这样的困境做KEGG通路富集分析时生成了富集结果但总感觉用气泡图展示不够“带感”想换个更直观、能体现通路间关联的桑基图却发现自己得从头写另一套代码费时费力今天这篇文章就是来解决这个痛点的。我们将打破常规教你如何用一套核心R语言代码同时生成两种风格迥异但都极具价值的可视化图表经典的KEGG气泡图和能揭示通路间基因流动关系的桑基流向图。这不仅仅是“画两个图”而是让你掌握一种高效、灵活的数据可视化工作流一次分析双重视角最大化挖掘富集结果的价值。读完本文你将能理解KEGG富集分析结果的核心数据结构。掌握使用ggplot2绘制标准KEGG气泡图的核心技巧与美化要点。学会将富集结果数据转换为桑基图所需的“源-目标-流量”格式。使用ggalluvial或networkD3等包一键生成揭示基因在通路间共享关系的桑基图。获得一套可复用的R脚本轻松适配你自己的数据。1. 为什么需要“一图两出”气泡图与桑基图的互补价值在生信分析中可视化从来不只是为了“好看”而是为了更高效地“传达信息”。KEGG气泡图和桑基图服务于不同的洞察目的。KEGG气泡图是富集分析的“标准答案”。它通过点的大小如基因数量、颜色如p值或FDR和位置富集因子在一个二维平面上直观地展示了哪些通路最显著、影响最大。它的优势在于排序和筛选能让你快速锁定Top N的重要通路。但它有一个天生的局限每个通路是孤立的点无法展示基因在不同通路间的重叠与交叉关系。而桑基图Sankey Diagram正是为此而生。它通过流动的线条流连接不同的节点通路线条的粗细代表共享基因的数量。这能直观回答“我的差异基因除了富集在最显著的A通路是否也大量参与了B、C通路”、“这些显著通路之间是否存在共同的调控基因群”。这对于理解复杂的生物学过程、发现潜在的核心调控网络至关重要。因此“一套代码出两图”的本质是让你从“找重点通路”气泡图深入到“理通路关系”桑基图用最低的边际成本获得双倍的数据洞察。下面我们就从数据准备开始一步步实现。2. 环境准备与核心R包安装在开始之前请确保你的R环境已经就绪。本文假设你已经完成了差异表达分析并得到了KEGG富集分析的结果。我们主要需要以下可视化相关的包。2.1 必需R包安装与加载打开你的RStudio或R控制台执行以下命令安装所需的包。如果已安装library()命令会直接加载。# 安装必要的包如果尚未安装 install.packages(c(ggplot2, dplyr, tidyr, stringr)) install.packages(ggalluvial) # 用于绘制静态桑基图 install.packages(networkD3) # 用于绘制交互式桑基图依赖htmlwidgets # 注意如果进行KEGG富集分析你可能还需要clusterProfiler、org.Hs.eg.db等此处聚焦可视化。 # 加载所有需要的包 library(ggplot2) library(dplyr) library(tidyr) library(stringr) library(ggalluvial) library(networkD3)2.2 示例数据准备为了演示我们模拟一个典型的KEGG富集分析结果数据框kegg_result。你的真实数据格式应与此类似。# 模拟KEGG富集结果数据 set.seed(123) # 确保结果可重复 kegg_result - data.frame( ID c(hsa04110, hsa03010, hsa04114, hsa05200, hsa04010, hsa04630, hsa04150, hsa05212), Description c(Cell cycle, Ribosome, Oocyte meiosis, Pathways in cancer, MAPK signaling pathway, Jak-STAT signaling pathway, mTOR signaling pathway, Pancreatic cancer), GeneRatio c(25/200, 18/200, 15/200, 22/200, 20/200, 12/200, 10/200, 8/200), BgRatio c(120/8000, 90/8000, 80/8000, 150/8000, 130/8000, 70/8000, 60/8000, 50/8000), pvalue c(1.2e-08, 3.5e-06, 7.2e-05, 9.8e-05, 1.5e-04, 3.1e-04, 5.0e-04, 8.0e-04), p.adjust c(1.5e-06, 2.1e-04, 3.0e-03, 3.8e-03, 4.5e-03, 8.0e-03, 1.1e-02, 1.6e-02), qvalue c(1.0e-06, 1.5e-04, 2.2e-03, 2.8e-03, 3.3e-03, 6.0e-03, 8.0e-03, 1.2e-02), geneID c(CDC20/CDK1/CCNB1/..., RPSA/RPS2/RPLP1/..., AURKA/PLK1/CDC25C/..., TP53/MYC/EGFR/..., MAPK1/MAPK3/DUSP1/..., STAT1/STAT3/JAK2/..., MTOR/RPTOR/MLST8/..., KRAS/TP53/CDKN2A/...), Count c(25, 18, 15, 22, 20, 12, 10, 8) ) # 查看数据结构 head(kegg_result)这段代码创建了一个包含通路ID、描述、富集基因数、p值、校正p值以及基因列表的模拟数据框。geneID列是以“/”分隔的基因符号字符串这是clusterProfiler等包的常见输出格式。3. 核心数据预处理一套代码的基石“一套代码”的关键在于前期构建一个干净、规整的数据集。预处理做得好后续画图就是调用不同函数的事情。3.1 数据清洗与格式转换我们需要从原始结果中提取画图所需的关键信息通路名、富集基因数、显著性指标并为桑基图准备基因-通路的对应关系。# 1. 为气泡图准备数据选择Top N通路并计算富集因子(Enrichment Factor) # 富集因子 ≈ (GeneRatio) / (BgRatio)这里用Count和背景基因总数近似计算 top_n - 10 # 选择最显著的前10条通路或根据你的数据调整 kegg_for_bubble - kegg_result %% arrange(p.adjust) %% # 按校正p值排序 slice_head(n top_n) %% # 取前top_n行 mutate( GeneRatio_num sapply(strsplit(GeneRatio, /), function(x) as.numeric(x[1]) / as.numeric(x[2])), BgRatio_num sapply(strsplit(BgRatio, /), function(x) as.numeric(x[1]) / as.numeric(x[2])), EnrichmentFactor GeneRatio_num / BgRatio_num, # 将p.adjust转换为-log10用于颜色映射值越大越显著 log10_padj -log10(p.adjust) ) %% # 对通路描述进行排序因子化保证图中顺序 mutate(Description factor(Description, levels rev(unique(Description)))) # 查看处理后的气泡图数据 print(kegg_for_bubble[, c(Description, Count, p.adjust, EnrichmentFactor, log10_padj)])3.2 为桑基图构建“边-节点”数据这是从气泡图数据延伸到桑基图的关键一步。我们需要把geneID列拆开形成“每个基因属于哪些通路”的长格式数据。# 2. 为桑基图准备数据构建基因-通路的关联矩阵长数据格式 # 假设我们只使用前6条最显著的通路来制作桑基图避免过于复杂 top_pathways_for_sankey - kegg_for_bubble$Description[1:6] sankey_edges - kegg_result %% filter(Description %in% top_pathways_for_sankey) %% # 拆分基因列表 separate_rows(geneID, sep /) %% # 清理可能的空字符串或NA filter(!is.na(geneID), geneID ! ) %% # 选择需要的列基因和通路 select(geneID, Description) %% # 去重每个基因在一条通路中只出现一次即使原始数据重复 distinct() %% # 统计每个基因出现的通路数用于后续筛选例如只保留属于至少2条通路的基因以简化图形 group_by(geneID) %% mutate(pathway_count n()) %% ungroup() # 查看桑基图边数据 head(sankey_edges) table(sankey_edges$pathway_count) # 查看基因的通路分布现在我们有了两个核心数据集kegg_for_bubble用于气泡图sankey_edges用于桑基图。预处理完成接下来就是可视化。4. 绘制标准KEGG气泡图ggplot2进阶版我们将绘制一个信息丰富且美观的气泡图包含颜色、大小、坐标轴和标签的优化。# 使用ggplot2绘制气泡图 bubble_plot - ggplot(kegg_for_bubble, aes(x EnrichmentFactor, y Description, size Count, color log10_padj)) # 1. 绘制点 geom_point(alpha 0.8) # 设置透明度 # 2. 颜色标尺使用连续色系越红越显著 scale_color_gradient(low blue, high red, name -log10(Adj.P), guide guide_colorbar(reverse FALSE)) # 3. 大小标尺设置点的大小范围 scale_size_continuous(range c(3, 10), name Gene Count, breaks seq(min(kegg_for_bubble$Count), max(kegg_for_bubble$Count), length.out 3) %% round()) ) # 4. 坐标轴与标签 labs(x Enrichment Factor, y NULL, title KEGG Pathway Enrichment Analysis, subtitle paste(Top, top_n, most significant pathways)) # 5. 主题美化 theme_bw(base_size 12) theme( axis.title.x element_text(face bold, size 11), axis.text.y element_text(face bold, color black, size 10), axis.text.x element_text(color black, size 10), plot.title element_text(hjust 0.5, face bold, size 14), plot.subtitle element_text(hjust 0.5, size 10), legend.position right, legend.box vertical, panel.grid.major.y element_line(linetype dashed, color grey90), panel.grid.minor element_blank() ) # 6. 可选为点添加文本标签基因数 geom_text(aes(label Count), color white, size 3, fontface bold) # 显示图形 print(bubble_plot) # 保存图形高分辨率适合发表 ggsave(filename KEGG_Bubble_Plot.png, plot bubble_plot, width 10, height 7, dpi 300)代码关键点解读aes映射将富集因子映射给X轴体现富集程度通路名称给Y轴基因数给点大小显著性给颜色。这是气泡图的灵魂。scale_*函数精细控制颜色和大小标尺的图例让图例清晰易懂。theme调整theme_bw提供干净背景调整网格线、字体、标题位置等细节提升专业度。geom_text在点上叠加白色粗体数字直观显示基因数避免读者在大小和颜色图例间来回比对。至此一张出版级的气泡图已经生成。接下来我们利用已处理好的sankey_edges数据绘制桑基图。5. 绘制桑基流向图揭示通路间的基因共享网络桑基图有两种主流画法静态的ggalluvial和交互式的networkD3。我们将分别展示。5.1 方法一使用ggalluvial绘制静态桑基图ggalluvial基于ggplot2风格统一易于保存为PDF/PNG。# 首先我们需要将边数据转换为ggalluvial需要的格式每个基因作为一条观测通路作为分层变量。 # 但ggalluvial通常用于展示分类变量间的流动我们需要先创建一个从“基因”到“通路”的“流动”。 # 一个简单的转换是假设所有基因从一个虚拟的“Gene Pool”流向它们所属的各个通路。 # 创建ggalluvial数据格式 sankey_data_alluvial - sankey_edges %% # 筛选至少在2条通路中出现的基因使图形更清晰 filter(pathway_count 2) %% mutate(source Gene Set, # 虚拟的源节点 target Description) %% # 目标节点为通路 group_by(source, target) %% summarise(value n(), .groups drop) # 流量为共享基因数 # 查看转换后的数据 print(sankey_data_alluvial) # 使用ggalluvial绘制桑基图 library(ggalluvial) alluvial_plot - ggplot(sankey_data_alluvial, aes(axis1 source, axis2 target, y value)) geom_alluvium(aes(fill target), width 1/12, alpha 0.7, knot.pos 0) geom_stratum(width 1/12, fill grey90, color grey50) geom_text(stat stratum, aes(label after_stat(stratum)), size 3) scale_x_discrete(limits c(Source, Pathway), expand c(0.05, 0.05)) scale_fill_brewer(palette Set3, name KEGG Pathway) labs(title Gene Flow to KEGG Pathways (Sankey Diagram), subtitle Showing genes shared among top significant pathways, y Number of Shared Genes) theme_minimal() theme(legend.position none, axis.text.y element_blank(), axis.title.y element_text(angle 0, vjust 0.5), panel.grid element_blank()) print(alluvial_plot) ggsave(KEGG_Sankey_Alluvial.png, alluvial_plot, width 9, height 6, dpi300)这种方法直观展示了基因集到不同通路的“分配”但对于展示通路-通路间的基因共享不够直接。更经典的方法是构建通路-通路的共基因网络。5.2 方法二使用networkD3绘制交互式通路-通路桑基图推荐交互式桑基图能让你用鼠标悬停查看细节更适合探索复杂关系。我们需要构建节点列表和边列表。# 构建节点列表每个唯一的通路名称就是一个节点 nodes - data.frame(name unique(sankey_edges$Description)) # 为每个节点添加ID从0开始networkD3的要求 nodes$id - 0:(nrow(nodes) - 1) # 构建边列表计算每对通路之间共享的基因数量 # 这是一个组合问题找出所有通路对 pathway_pairs - t(combn(unique(sankey_edges$Description), 2)) links - data.frame(source pathway_pairs[,1], target pathway_pairs[,2], value 0) # 计算共享基因数 for(i in 1:nrow(links)) { genes_in_source - sankey_edges$geneID[sankey_edges$Description links$source[i]] genes_in_target - sankey_edges$geneID[sankey_edges$Description links$target[i]] shared_genes - intersect(genes_in_source, genes_in_target) links$value[i] - length(shared_genes) } # 过滤掉共享基因数为0的边可选使图形更简洁 links - links[links$value 0, ] # 将通路名称转换为节点ID links - links %% left_join(nodes, by c(source name)) %% rename(source_id id) %% left_join(nodes, by c(target name)) %% rename(target_id id) %% select(source_id, target_id, value) # 使用networkD3绘制交互式桑基图 sankeyNetwork(Links links, # 边数据框 Nodes nodes, # 节点数据框 Source source_id, Target target_id, Value value, NodeID name, fontSize 14, nodeWidth 20, nodePadding 10, sinksRight FALSE) # 节点不强制右对齐运行这段代码会在RStudio的Viewer面板或浏览器中生成一个交互式桑基图。你可以鼠标悬停在节点上高亮显示所有与该通路相连的流。鼠标悬停在流上显示该连接通路对共享的基因数量。拖动节点重新布局。保存为HTML使用htmlwidgets::saveWidget()保存为独立网页分享。这才是真正的“通路关系桑基图”。它清晰展示了哪些通路之间共享大量基因流粗暗示这些通路可能在功能上紧密相关受共同的核心基因调控。6. 将两图整合进一个分析流程现在我们将所有步骤整合到一个连贯的、可复用的R脚本中。你只需要将kegg_result替换成自己的富集结果数据框。# 完整流程从KEGG结果到双图 # 作者你的名字 # 功能输入KEGG富集结果输出气泡图和交互式桑基图 # # 0. 加载包 library(ggplot2) library(dplyr) library(tidyr) library(networkD3) # 1. 数据预处理函数 prepare_kegg_data - function(kegg_df, top_n 15, min_shared_genes 2) { # 气泡图数据 bubble_df - kegg_df %% arrange(p.adjust) %% slice_head(n top_n) %% mutate( GeneRatio_num sapply(strsplit(GeneRatio, /), function(x) as.numeric(x[1]) / as.numeric(x[2])), BgRatio_num sapply(strsplit(BgRatio, /), function(x) as.numeric(x[1]) / as.numeric(x[2])), EnrichmentFactor GeneRatio_num / BgRatio_num, log10_padj -log10(p.adjust), Description factor(Description, levels rev(unique(Description))) ) # 桑基图边数据基因-通路关联 edges - kegg_df %% filter(Description %in% bubble_df$Description) %% separate_rows(geneID, sep /) %% filter(!is.na(geneID), geneID ! ) %% select(geneID, Description) %% distinct() %% group_by(geneID) %% mutate(pathway_count n()) %% ungroup() %% filter(pathway_count min_shared_genes) # 过滤只出现在一条通路的基因 return(list(bubble bubble_df, edges edges)) } # 2. 绘制气泡图函数 plot_bubble - function(bubble_df, title_suffix ) { p - ggplot(bubble_df, aes(x EnrichmentFactor, y Description, size Count, color log10_padj)) geom_point(alpha 0.8) scale_color_gradient(low blue, high red, name -log10(Adj.P)) scale_size_continuous(range c(3, 10), name Gene Count) labs(x Enrichment Factor, y NULL, title paste(KEGG Pathway Enrichment, title_suffix), subtitle paste(Top, nrow(bubble_df), significant pathways)) theme_bw() theme(axis.text.y element_text(face bold), plot.title element_text(hjust 0.5, face bold)) return(p) } # 3. 绘制桑基图函数 plot_sankey - function(edge_df) { if(nrow(edge_df) 0) { message(No genes shared across multiple pathways. Sankey diagram skipped.) return(NULL) } nodes - data.frame(name unique(edge_df$Description), stringsAsFactors FALSE) nodes$id - 0:(nrow(nodes) - 1) pathway_pairs - t(combn(nodes$name, 2)) links - data.frame(source pathway_pairs[,1], target pathway_pairs[,2], value 0, stringsAsFactors FALSE) for(i in 1:nrow(links)) { g1 - edge_df$geneID[edge_df$Description links$source[i]] g2 - edge_df$geneID[edge_df$Description links$target[i]] links$value[i] - length(intersect(g1, g2)) } links - links[links$value 0, ] if(nrow(links) 0) { message(No pathway pairs share genes. Sankey diagram skipped.) return(NULL) } links - links %% left_join(nodes, by c(source name)) %% rename(source_id id) %% left_join(nodes, by c(target name)) %% rename(target_id id) %% select(source_id, target_id, value) sankey - sankeyNetwork(Links links, Nodes nodes, Source source_id, Target target_id, Value value, NodeID name, fontSize 14, nodeWidth 20, nodePadding 10, sinksRight FALSE) return(sankey) } # 4. 主执行流程 # 假设你的KEGG结果已经在 my_kegg_result 数据框中 # my_kegg_result - read.csv(your_kegg_enrichment.csv) # 或从clusterProfiler结果转换 # 使用模拟数据演示 data_prepared - prepare_kegg_data(kegg_result, top_n 8, min_shared_genes 1) # 生成气泡图 bubble_p - plot_bubble(data_prepared$bubble, title_suffix (Demo Data)) print(bubble_p) ggsave(My_KEGG_Bubble.png, bubble_p, width10, height6, dpi300) # 生成交互式桑基图 sankey_p - plot_sankey(data_prepared$edges) if(!is.null(sankey_p)) { print(sankey_p) # 在RStudio中查看 # 保存为HTML文件 library(htmlwidgets) saveWidget(sankey_p, file My_KEGG_Sankey.html) }这个脚本提供了完整的函数封装你可以轻松地将其整合到自己的分析流程中。7. 常见问题与排查思路在实际操作中你可能会遇到一些问题。下表列出了常见问题及解决方法。问题现象可能原因排查方式解决方案气泡图点的大小或颜色异常数据列Count或p.adjust包含非数值如字符、NA。使用str(kegg_for_bubble)检查数据结构用summary()查看数值列范围。确保用于映射的列是数值型。用as.numeric()转换或用na.omit()处理缺失值。桑基图没有显示或节点缺失links数据框为空或value全为0。打印links数据框检查nrow(links)和head(links)。检查min_shared_genes参数是否设置过高。降低min_shared_genes阈值例如设为1或检查原始geneID列分隔符是否正确不一定是“/”。separate_rows报错geneID列不是字符型或分隔符与代码中sep参数不匹配。使用class(kegg_result$geneID)查看类型。查看kegg_result$geneID[1]的实际格式。用mutate(geneID as.character(geneID))转换类型。根据实际情况修改sep参数如sep /或sep 、。交互式桑基图在浏览器中不显示未安装或加载htmlwidgets包或浏览器安全策略阻止。确保已安装并加载library(htmlwidgets)。尝试在RStudio的Viewer面板中查看。使用saveWidget()保存为HTML后用浏览器直接打开该文件。确保网络环境允许加载D3.js库通常离线也可用。图形过于拥挤选择的通路数 (top_n) 过多或共享基因的过滤条件 (min_shared_genes) 太宽松。观察生成的节点和边数量。超过15个节点或50条边就会显得杂乱。减少top_n如8-12增加min_shared_genes如2或3。在桑基图函数中可添加links - links[links$value 1, ]过滤弱连接。富集因子计算为NA或InfGeneRatio或BgRatio格式非“a/b”或分母为0。打印kegg_result$GeneRatio和kegg_result$BgRatio的前几行。确保比率字符串格式正确。使用更稳健的解析函数例如sapply(strsplit(GeneRatio, /), function(x) ifelse(length(x)2, as.numeric(x[1])/as.numeric(x[2]), NA))。8. 最佳实践与进阶技巧掌握了基础流程后以下几点能让你的分析更上一层楼数据源适配本文假设geneID列以“/”分隔。如果你的数据来自clusterProfiler的enrichResult对象可以使用as.data.frame()直接转换geneID列默认即为该格式。如果是其他工具输出请先确认分隔符。结果筛选策略绘制桑基图前对通路进行智能筛选。基于显著性只选择p.adjust 0.05的通路。基于基因数选择Count大于某个阈值如5的通路。基于功能手动挑选你感兴趣的通路子集如所有癌症相关通路。桑基图的美化与定制节点颜色在nodes数据框中增加一列group可以在sankeyNetwork()中通过NodeGroup group参数对节点按类着色。链接颜色设置LinkGroup source_name可以让流颜色与源节点一致。布局优化sankeyNetwork()的iterations参数默认32控制布局算法迭代次数增加此值可能使布局更优。从桑基图到网络分析桑基图揭示了通路间的关联你可以进一步将links数据导入igraph或Cytoscape进行更复杂的网络分析如计算中心性、识别模块等。自动化与报告生成将整个流程封装进一个R函数或R Markdown文档中。每次更新富集分析结果只需运行一个命令或点击“Knit”即可自动生成包含双图的HTML或PDF报告。9. 总结通过本文我们完成了一次从数据预处理到双图可视化的完整旅程。核心收获在于理解气泡图是纵向深度的展示告诉我们每条通路自身的显著性强度和规模。桑基图是横向关联的展示揭示不同通路之间通过共享基因构成的功能网络。“一套代码出两图”的精髓不在于代码本身有多短而在于构建了一个可复用、可扩展的数据处理管道。你投入一次时间进行数据清洗和转换就能轻松获得两种不同维度的洞察极大地提升了分析效率和结果的解读深度。下次当你完成KEGG富集分析不妨在生成标准气泡图之后多花几分钟运行一下本文的桑基图代码。你可能会发现那些共享基因最多的通路往往指向了你课题中最核心、最交织的生物学故事。
R语言实现KEGG通路富集分析：气泡图与桑基图双重视角可视化
你是不是也遇到过这样的困境做KEGG通路富集分析时生成了富集结果但总感觉用气泡图展示不够“带感”想换个更直观、能体现通路间关联的桑基图却发现自己得从头写另一套代码费时费力今天这篇文章就是来解决这个痛点的。我们将打破常规教你如何用一套核心R语言代码同时生成两种风格迥异但都极具价值的可视化图表经典的KEGG气泡图和能揭示通路间基因流动关系的桑基流向图。这不仅仅是“画两个图”而是让你掌握一种高效、灵活的数据可视化工作流一次分析双重视角最大化挖掘富集结果的价值。读完本文你将能理解KEGG富集分析结果的核心数据结构。掌握使用ggplot2绘制标准KEGG气泡图的核心技巧与美化要点。学会将富集结果数据转换为桑基图所需的“源-目标-流量”格式。使用ggalluvial或networkD3等包一键生成揭示基因在通路间共享关系的桑基图。获得一套可复用的R脚本轻松适配你自己的数据。1. 为什么需要“一图两出”气泡图与桑基图的互补价值在生信分析中可视化从来不只是为了“好看”而是为了更高效地“传达信息”。KEGG气泡图和桑基图服务于不同的洞察目的。KEGG气泡图是富集分析的“标准答案”。它通过点的大小如基因数量、颜色如p值或FDR和位置富集因子在一个二维平面上直观地展示了哪些通路最显著、影响最大。它的优势在于排序和筛选能让你快速锁定Top N的重要通路。但它有一个天生的局限每个通路是孤立的点无法展示基因在不同通路间的重叠与交叉关系。而桑基图Sankey Diagram正是为此而生。它通过流动的线条流连接不同的节点通路线条的粗细代表共享基因的数量。这能直观回答“我的差异基因除了富集在最显著的A通路是否也大量参与了B、C通路”、“这些显著通路之间是否存在共同的调控基因群”。这对于理解复杂的生物学过程、发现潜在的核心调控网络至关重要。因此“一套代码出两图”的本质是让你从“找重点通路”气泡图深入到“理通路关系”桑基图用最低的边际成本获得双倍的数据洞察。下面我们就从数据准备开始一步步实现。2. 环境准备与核心R包安装在开始之前请确保你的R环境已经就绪。本文假设你已经完成了差异表达分析并得到了KEGG富集分析的结果。我们主要需要以下可视化相关的包。2.1 必需R包安装与加载打开你的RStudio或R控制台执行以下命令安装所需的包。如果已安装library()命令会直接加载。# 安装必要的包如果尚未安装 install.packages(c(ggplot2, dplyr, tidyr, stringr)) install.packages(ggalluvial) # 用于绘制静态桑基图 install.packages(networkD3) # 用于绘制交互式桑基图依赖htmlwidgets # 注意如果进行KEGG富集分析你可能还需要clusterProfiler、org.Hs.eg.db等此处聚焦可视化。 # 加载所有需要的包 library(ggplot2) library(dplyr) library(tidyr) library(stringr) library(ggalluvial) library(networkD3)2.2 示例数据准备为了演示我们模拟一个典型的KEGG富集分析结果数据框kegg_result。你的真实数据格式应与此类似。# 模拟KEGG富集结果数据 set.seed(123) # 确保结果可重复 kegg_result - data.frame( ID c(hsa04110, hsa03010, hsa04114, hsa05200, hsa04010, hsa04630, hsa04150, hsa05212), Description c(Cell cycle, Ribosome, Oocyte meiosis, Pathways in cancer, MAPK signaling pathway, Jak-STAT signaling pathway, mTOR signaling pathway, Pancreatic cancer), GeneRatio c(25/200, 18/200, 15/200, 22/200, 20/200, 12/200, 10/200, 8/200), BgRatio c(120/8000, 90/8000, 80/8000, 150/8000, 130/8000, 70/8000, 60/8000, 50/8000), pvalue c(1.2e-08, 3.5e-06, 7.2e-05, 9.8e-05, 1.5e-04, 3.1e-04, 5.0e-04, 8.0e-04), p.adjust c(1.5e-06, 2.1e-04, 3.0e-03, 3.8e-03, 4.5e-03, 8.0e-03, 1.1e-02, 1.6e-02), qvalue c(1.0e-06, 1.5e-04, 2.2e-03, 2.8e-03, 3.3e-03, 6.0e-03, 8.0e-03, 1.2e-02), geneID c(CDC20/CDK1/CCNB1/..., RPSA/RPS2/RPLP1/..., AURKA/PLK1/CDC25C/..., TP53/MYC/EGFR/..., MAPK1/MAPK3/DUSP1/..., STAT1/STAT3/JAK2/..., MTOR/RPTOR/MLST8/..., KRAS/TP53/CDKN2A/...), Count c(25, 18, 15, 22, 20, 12, 10, 8) ) # 查看数据结构 head(kegg_result)这段代码创建了一个包含通路ID、描述、富集基因数、p值、校正p值以及基因列表的模拟数据框。geneID列是以“/”分隔的基因符号字符串这是clusterProfiler等包的常见输出格式。3. 核心数据预处理一套代码的基石“一套代码”的关键在于前期构建一个干净、规整的数据集。预处理做得好后续画图就是调用不同函数的事情。3.1 数据清洗与格式转换我们需要从原始结果中提取画图所需的关键信息通路名、富集基因数、显著性指标并为桑基图准备基因-通路的对应关系。# 1. 为气泡图准备数据选择Top N通路并计算富集因子(Enrichment Factor) # 富集因子 ≈ (GeneRatio) / (BgRatio)这里用Count和背景基因总数近似计算 top_n - 10 # 选择最显著的前10条通路或根据你的数据调整 kegg_for_bubble - kegg_result %% arrange(p.adjust) %% # 按校正p值排序 slice_head(n top_n) %% # 取前top_n行 mutate( GeneRatio_num sapply(strsplit(GeneRatio, /), function(x) as.numeric(x[1]) / as.numeric(x[2])), BgRatio_num sapply(strsplit(BgRatio, /), function(x) as.numeric(x[1]) / as.numeric(x[2])), EnrichmentFactor GeneRatio_num / BgRatio_num, # 将p.adjust转换为-log10用于颜色映射值越大越显著 log10_padj -log10(p.adjust) ) %% # 对通路描述进行排序因子化保证图中顺序 mutate(Description factor(Description, levels rev(unique(Description)))) # 查看处理后的气泡图数据 print(kegg_for_bubble[, c(Description, Count, p.adjust, EnrichmentFactor, log10_padj)])3.2 为桑基图构建“边-节点”数据这是从气泡图数据延伸到桑基图的关键一步。我们需要把geneID列拆开形成“每个基因属于哪些通路”的长格式数据。# 2. 为桑基图准备数据构建基因-通路的关联矩阵长数据格式 # 假设我们只使用前6条最显著的通路来制作桑基图避免过于复杂 top_pathways_for_sankey - kegg_for_bubble$Description[1:6] sankey_edges - kegg_result %% filter(Description %in% top_pathways_for_sankey) %% # 拆分基因列表 separate_rows(geneID, sep /) %% # 清理可能的空字符串或NA filter(!is.na(geneID), geneID ! ) %% # 选择需要的列基因和通路 select(geneID, Description) %% # 去重每个基因在一条通路中只出现一次即使原始数据重复 distinct() %% # 统计每个基因出现的通路数用于后续筛选例如只保留属于至少2条通路的基因以简化图形 group_by(geneID) %% mutate(pathway_count n()) %% ungroup() # 查看桑基图边数据 head(sankey_edges) table(sankey_edges$pathway_count) # 查看基因的通路分布现在我们有了两个核心数据集kegg_for_bubble用于气泡图sankey_edges用于桑基图。预处理完成接下来就是可视化。4. 绘制标准KEGG气泡图ggplot2进阶版我们将绘制一个信息丰富且美观的气泡图包含颜色、大小、坐标轴和标签的优化。# 使用ggplot2绘制气泡图 bubble_plot - ggplot(kegg_for_bubble, aes(x EnrichmentFactor, y Description, size Count, color log10_padj)) # 1. 绘制点 geom_point(alpha 0.8) # 设置透明度 # 2. 颜色标尺使用连续色系越红越显著 scale_color_gradient(low blue, high red, name -log10(Adj.P), guide guide_colorbar(reverse FALSE)) # 3. 大小标尺设置点的大小范围 scale_size_continuous(range c(3, 10), name Gene Count, breaks seq(min(kegg_for_bubble$Count), max(kegg_for_bubble$Count), length.out 3) %% round()) ) # 4. 坐标轴与标签 labs(x Enrichment Factor, y NULL, title KEGG Pathway Enrichment Analysis, subtitle paste(Top, top_n, most significant pathways)) # 5. 主题美化 theme_bw(base_size 12) theme( axis.title.x element_text(face bold, size 11), axis.text.y element_text(face bold, color black, size 10), axis.text.x element_text(color black, size 10), plot.title element_text(hjust 0.5, face bold, size 14), plot.subtitle element_text(hjust 0.5, size 10), legend.position right, legend.box vertical, panel.grid.major.y element_line(linetype dashed, color grey90), panel.grid.minor element_blank() ) # 6. 可选为点添加文本标签基因数 geom_text(aes(label Count), color white, size 3, fontface bold) # 显示图形 print(bubble_plot) # 保存图形高分辨率适合发表 ggsave(filename KEGG_Bubble_Plot.png, plot bubble_plot, width 10, height 7, dpi 300)代码关键点解读aes映射将富集因子映射给X轴体现富集程度通路名称给Y轴基因数给点大小显著性给颜色。这是气泡图的灵魂。scale_*函数精细控制颜色和大小标尺的图例让图例清晰易懂。theme调整theme_bw提供干净背景调整网格线、字体、标题位置等细节提升专业度。geom_text在点上叠加白色粗体数字直观显示基因数避免读者在大小和颜色图例间来回比对。至此一张出版级的气泡图已经生成。接下来我们利用已处理好的sankey_edges数据绘制桑基图。5. 绘制桑基流向图揭示通路间的基因共享网络桑基图有两种主流画法静态的ggalluvial和交互式的networkD3。我们将分别展示。5.1 方法一使用ggalluvial绘制静态桑基图ggalluvial基于ggplot2风格统一易于保存为PDF/PNG。# 首先我们需要将边数据转换为ggalluvial需要的格式每个基因作为一条观测通路作为分层变量。 # 但ggalluvial通常用于展示分类变量间的流动我们需要先创建一个从“基因”到“通路”的“流动”。 # 一个简单的转换是假设所有基因从一个虚拟的“Gene Pool”流向它们所属的各个通路。 # 创建ggalluvial数据格式 sankey_data_alluvial - sankey_edges %% # 筛选至少在2条通路中出现的基因使图形更清晰 filter(pathway_count 2) %% mutate(source Gene Set, # 虚拟的源节点 target Description) %% # 目标节点为通路 group_by(source, target) %% summarise(value n(), .groups drop) # 流量为共享基因数 # 查看转换后的数据 print(sankey_data_alluvial) # 使用ggalluvial绘制桑基图 library(ggalluvial) alluvial_plot - ggplot(sankey_data_alluvial, aes(axis1 source, axis2 target, y value)) geom_alluvium(aes(fill target), width 1/12, alpha 0.7, knot.pos 0) geom_stratum(width 1/12, fill grey90, color grey50) geom_text(stat stratum, aes(label after_stat(stratum)), size 3) scale_x_discrete(limits c(Source, Pathway), expand c(0.05, 0.05)) scale_fill_brewer(palette Set3, name KEGG Pathway) labs(title Gene Flow to KEGG Pathways (Sankey Diagram), subtitle Showing genes shared among top significant pathways, y Number of Shared Genes) theme_minimal() theme(legend.position none, axis.text.y element_blank(), axis.title.y element_text(angle 0, vjust 0.5), panel.grid element_blank()) print(alluvial_plot) ggsave(KEGG_Sankey_Alluvial.png, alluvial_plot, width 9, height 6, dpi300)这种方法直观展示了基因集到不同通路的“分配”但对于展示通路-通路间的基因共享不够直接。更经典的方法是构建通路-通路的共基因网络。5.2 方法二使用networkD3绘制交互式通路-通路桑基图推荐交互式桑基图能让你用鼠标悬停查看细节更适合探索复杂关系。我们需要构建节点列表和边列表。# 构建节点列表每个唯一的通路名称就是一个节点 nodes - data.frame(name unique(sankey_edges$Description)) # 为每个节点添加ID从0开始networkD3的要求 nodes$id - 0:(nrow(nodes) - 1) # 构建边列表计算每对通路之间共享的基因数量 # 这是一个组合问题找出所有通路对 pathway_pairs - t(combn(unique(sankey_edges$Description), 2)) links - data.frame(source pathway_pairs[,1], target pathway_pairs[,2], value 0) # 计算共享基因数 for(i in 1:nrow(links)) { genes_in_source - sankey_edges$geneID[sankey_edges$Description links$source[i]] genes_in_target - sankey_edges$geneID[sankey_edges$Description links$target[i]] shared_genes - intersect(genes_in_source, genes_in_target) links$value[i] - length(shared_genes) } # 过滤掉共享基因数为0的边可选使图形更简洁 links - links[links$value 0, ] # 将通路名称转换为节点ID links - links %% left_join(nodes, by c(source name)) %% rename(source_id id) %% left_join(nodes, by c(target name)) %% rename(target_id id) %% select(source_id, target_id, value) # 使用networkD3绘制交互式桑基图 sankeyNetwork(Links links, # 边数据框 Nodes nodes, # 节点数据框 Source source_id, Target target_id, Value value, NodeID name, fontSize 14, nodeWidth 20, nodePadding 10, sinksRight FALSE) # 节点不强制右对齐运行这段代码会在RStudio的Viewer面板或浏览器中生成一个交互式桑基图。你可以鼠标悬停在节点上高亮显示所有与该通路相连的流。鼠标悬停在流上显示该连接通路对共享的基因数量。拖动节点重新布局。保存为HTML使用htmlwidgets::saveWidget()保存为独立网页分享。这才是真正的“通路关系桑基图”。它清晰展示了哪些通路之间共享大量基因流粗暗示这些通路可能在功能上紧密相关受共同的核心基因调控。6. 将两图整合进一个分析流程现在我们将所有步骤整合到一个连贯的、可复用的R脚本中。你只需要将kegg_result替换成自己的富集结果数据框。# 完整流程从KEGG结果到双图 # 作者你的名字 # 功能输入KEGG富集结果输出气泡图和交互式桑基图 # # 0. 加载包 library(ggplot2) library(dplyr) library(tidyr) library(networkD3) # 1. 数据预处理函数 prepare_kegg_data - function(kegg_df, top_n 15, min_shared_genes 2) { # 气泡图数据 bubble_df - kegg_df %% arrange(p.adjust) %% slice_head(n top_n) %% mutate( GeneRatio_num sapply(strsplit(GeneRatio, /), function(x) as.numeric(x[1]) / as.numeric(x[2])), BgRatio_num sapply(strsplit(BgRatio, /), function(x) as.numeric(x[1]) / as.numeric(x[2])), EnrichmentFactor GeneRatio_num / BgRatio_num, log10_padj -log10(p.adjust), Description factor(Description, levels rev(unique(Description))) ) # 桑基图边数据基因-通路关联 edges - kegg_df %% filter(Description %in% bubble_df$Description) %% separate_rows(geneID, sep /) %% filter(!is.na(geneID), geneID ! ) %% select(geneID, Description) %% distinct() %% group_by(geneID) %% mutate(pathway_count n()) %% ungroup() %% filter(pathway_count min_shared_genes) # 过滤只出现在一条通路的基因 return(list(bubble bubble_df, edges edges)) } # 2. 绘制气泡图函数 plot_bubble - function(bubble_df, title_suffix ) { p - ggplot(bubble_df, aes(x EnrichmentFactor, y Description, size Count, color log10_padj)) geom_point(alpha 0.8) scale_color_gradient(low blue, high red, name -log10(Adj.P)) scale_size_continuous(range c(3, 10), name Gene Count) labs(x Enrichment Factor, y NULL, title paste(KEGG Pathway Enrichment, title_suffix), subtitle paste(Top, nrow(bubble_df), significant pathways)) theme_bw() theme(axis.text.y element_text(face bold), plot.title element_text(hjust 0.5, face bold)) return(p) } # 3. 绘制桑基图函数 plot_sankey - function(edge_df) { if(nrow(edge_df) 0) { message(No genes shared across multiple pathways. Sankey diagram skipped.) return(NULL) } nodes - data.frame(name unique(edge_df$Description), stringsAsFactors FALSE) nodes$id - 0:(nrow(nodes) - 1) pathway_pairs - t(combn(nodes$name, 2)) links - data.frame(source pathway_pairs[,1], target pathway_pairs[,2], value 0, stringsAsFactors FALSE) for(i in 1:nrow(links)) { g1 - edge_df$geneID[edge_df$Description links$source[i]] g2 - edge_df$geneID[edge_df$Description links$target[i]] links$value[i] - length(intersect(g1, g2)) } links - links[links$value 0, ] if(nrow(links) 0) { message(No pathway pairs share genes. Sankey diagram skipped.) return(NULL) } links - links %% left_join(nodes, by c(source name)) %% rename(source_id id) %% left_join(nodes, by c(target name)) %% rename(target_id id) %% select(source_id, target_id, value) sankey - sankeyNetwork(Links links, Nodes nodes, Source source_id, Target target_id, Value value, NodeID name, fontSize 14, nodeWidth 20, nodePadding 10, sinksRight FALSE) return(sankey) } # 4. 主执行流程 # 假设你的KEGG结果已经在 my_kegg_result 数据框中 # my_kegg_result - read.csv(your_kegg_enrichment.csv) # 或从clusterProfiler结果转换 # 使用模拟数据演示 data_prepared - prepare_kegg_data(kegg_result, top_n 8, min_shared_genes 1) # 生成气泡图 bubble_p - plot_bubble(data_prepared$bubble, title_suffix (Demo Data)) print(bubble_p) ggsave(My_KEGG_Bubble.png, bubble_p, width10, height6, dpi300) # 生成交互式桑基图 sankey_p - plot_sankey(data_prepared$edges) if(!is.null(sankey_p)) { print(sankey_p) # 在RStudio中查看 # 保存为HTML文件 library(htmlwidgets) saveWidget(sankey_p, file My_KEGG_Sankey.html) }这个脚本提供了完整的函数封装你可以轻松地将其整合到自己的分析流程中。7. 常见问题与排查思路在实际操作中你可能会遇到一些问题。下表列出了常见问题及解决方法。问题现象可能原因排查方式解决方案气泡图点的大小或颜色异常数据列Count或p.adjust包含非数值如字符、NA。使用str(kegg_for_bubble)检查数据结构用summary()查看数值列范围。确保用于映射的列是数值型。用as.numeric()转换或用na.omit()处理缺失值。桑基图没有显示或节点缺失links数据框为空或value全为0。打印links数据框检查nrow(links)和head(links)。检查min_shared_genes参数是否设置过高。降低min_shared_genes阈值例如设为1或检查原始geneID列分隔符是否正确不一定是“/”。separate_rows报错geneID列不是字符型或分隔符与代码中sep参数不匹配。使用class(kegg_result$geneID)查看类型。查看kegg_result$geneID[1]的实际格式。用mutate(geneID as.character(geneID))转换类型。根据实际情况修改sep参数如sep /或sep 、。交互式桑基图在浏览器中不显示未安装或加载htmlwidgets包或浏览器安全策略阻止。确保已安装并加载library(htmlwidgets)。尝试在RStudio的Viewer面板中查看。使用saveWidget()保存为HTML后用浏览器直接打开该文件。确保网络环境允许加载D3.js库通常离线也可用。图形过于拥挤选择的通路数 (top_n) 过多或共享基因的过滤条件 (min_shared_genes) 太宽松。观察生成的节点和边数量。超过15个节点或50条边就会显得杂乱。减少top_n如8-12增加min_shared_genes如2或3。在桑基图函数中可添加links - links[links$value 1, ]过滤弱连接。富集因子计算为NA或InfGeneRatio或BgRatio格式非“a/b”或分母为0。打印kegg_result$GeneRatio和kegg_result$BgRatio的前几行。确保比率字符串格式正确。使用更稳健的解析函数例如sapply(strsplit(GeneRatio, /), function(x) ifelse(length(x)2, as.numeric(x[1])/as.numeric(x[2]), NA))。8. 最佳实践与进阶技巧掌握了基础流程后以下几点能让你的分析更上一层楼数据源适配本文假设geneID列以“/”分隔。如果你的数据来自clusterProfiler的enrichResult对象可以使用as.data.frame()直接转换geneID列默认即为该格式。如果是其他工具输出请先确认分隔符。结果筛选策略绘制桑基图前对通路进行智能筛选。基于显著性只选择p.adjust 0.05的通路。基于基因数选择Count大于某个阈值如5的通路。基于功能手动挑选你感兴趣的通路子集如所有癌症相关通路。桑基图的美化与定制节点颜色在nodes数据框中增加一列group可以在sankeyNetwork()中通过NodeGroup group参数对节点按类着色。链接颜色设置LinkGroup source_name可以让流颜色与源节点一致。布局优化sankeyNetwork()的iterations参数默认32控制布局算法迭代次数增加此值可能使布局更优。从桑基图到网络分析桑基图揭示了通路间的关联你可以进一步将links数据导入igraph或Cytoscape进行更复杂的网络分析如计算中心性、识别模块等。自动化与报告生成将整个流程封装进一个R函数或R Markdown文档中。每次更新富集分析结果只需运行一个命令或点击“Knit”即可自动生成包含双图的HTML或PDF报告。9. 总结通过本文我们完成了一次从数据预处理到双图可视化的完整旅程。核心收获在于理解气泡图是纵向深度的展示告诉我们每条通路自身的显著性强度和规模。桑基图是横向关联的展示揭示不同通路之间通过共享基因构成的功能网络。“一套代码出两图”的精髓不在于代码本身有多短而在于构建了一个可复用、可扩展的数据处理管道。你投入一次时间进行数据清洗和转换就能轻松获得两种不同维度的洞察极大地提升了分析效率和结果的解读深度。下次当你完成KEGG富集分析不妨在生成标准气泡图之后多花几分钟运行一下本文的桑基图代码。你可能会发现那些共享基因最多的通路往往指向了你课题中最核心、最交织的生物学故事。
