ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

R语言一键生成KEGG气泡图与桑基图:生信结果可视化完整工作流

R语言一键生成KEGG气泡图与桑基图:生信结果可视化完整工作流 如果你在生物信息学分析中已经完成了差异表达基因的筛选和KEGG通路富集分析面对一长串富集结果下一步该怎么做直接贴一个满是P值和基因数的表格到论文里审稿人大概率会皱眉头。用Excel手动做个柱状图不仅费时费力而且信息呈现单一难以同时展示通路的统计学显著性和基因的流向关系。这正是很多生信初学者和科研工作者遇到的真实瓶颈分析做完了但不知道如何高效、专业地将结果可视化做出既能上得了台面又能讲清楚故事的图表。今天要解决的问题就是如何用同一套R语言代码和数据一键生成两种科研级图表——直观展示富集显著性的KEGG气泡图和清晰揭示基因-通路关联的桑基图Sankey Diagram。本文不会停留在“调用某个函数出图”的层面。我们将深入一套完整的工作流从富集结果表格的整理开始到使用ggplot2绘制可高度定制化的气泡图再到利用networkD3包构建交互式桑基图最后实现两图联动用数据讲述一个完整的生物学故事。无论你是刚接触R语言的生信新手还是想优化分析流程的老手这套方法都能直接提升你的成果展示效率与专业度。1. 核心价值为什么你需要掌握“一码两图”在生信数据分析中可视化不是分析的终点而是沟通的起点。KEGG富集分析的结果至少包含两个维度的信息通路的显著性如P值、Q值、富集因子和基因与通路的归属关系。单一图表往往只能强调其中一个维度。KEGG气泡图擅长展示“哪些通路最显著”。通过点的大小基因数和颜色P值一眼就能锁定关键通路适用于结果摘要和初步结论展示。桑基图擅长展示“基因具体是如何分配到各个通路中的”。它能清晰呈现基因与通路之间的多对多关系直观显示哪些基因是“多面手”参与多个通路哪些通路共享关键基因非常适合深入机制探讨。传统做法是分别用不同工具如R的ggplot2画气泡图再用Origin或在线工具画桑基图处理流程割裂且当数据更新时需要重复劳动。而“一码两图”的核心优势在于效率提升基于同一份清洁数据运行一次脚本同时获得两种视角的图表避免重复数据整理。一致性保证两图来源于同一数据源确保在汇报或论文中数据口径绝对一致。深度挖掘气泡图帮你筛选出TOP通路桑基图则帮你深入解读这些通路背后的基因联系形成分析闭环。灵活性高全部流程在R中完成从数据预处理到图形美化每一步都可控、可复现。接下来我们将从零开始拆解整个流程。2. 环境准备构建你的R生信绘图工作流工欲善其事必先利其器。本节将确保你拥有一个可稳定运行的环境。2.1 R与RStudio的安装与配置首先你需要安装R语言和RStudio集成开发环境IDE。R是引擎RStudio是方向盘和仪表盘能极大提升编码效率。安装R访问R官网https://www.r-project.org/选择与操作系统对应的CRAN镜像下载安装。建议选择较新的稳定版本如4.3.x。安装RStudio访问RStudio官网https://posit.co/download/rstudio-desktop/下载免费的Desktop版本并安装。安装后打开RStudio界面通常分为四个窗格脚本编辑器、控制台、环境/历史、文件/图/帮助。2.2 必需R包的安装与加载我们将主要依赖以下三个包ggplot2绘图语法之王用于绘制静态气泡图。dplyr/tidyr数据整理神器用于清洗和转换富集结果。networkD3用于创建交互式桑基图D3.js的R接口。在RStudio的控制台Console中一次性安装并加载它们# 安装包如果尚未安装 install.packages(c(ggplot2, dplyr, tidyr, networkD3)) # 加载包到当前会话 library(ggplot2) library(dplyr) library(tidyr) library(networkD3)注意networkD3包生成的桑基图是HTML交互式图表可在浏览器中查看支持鼠标悬停查看详细信息这比静态图更具探索性。2.3 准备示例数据为了演示我们模拟一份典型的KEGG富集分析结果。在实际项目中你只需将你的结果文件通常是.csv或.txt格式读入即可。# 模拟创建一份富集结果数据框 set.seed(123) # 确保结果可重复 kegg_result - data.frame( Pathway c(Metabolic pathways, Biosynthesis of secondary metabolites, Microbial metabolism in diverse environments, Carbon metabolism, Biosynthesis of amino acids, 2-Oxocarboxylic acid metabolism, Fatty acid metabolism, Degradation of aromatic compounds, ABC transporters, Quorum sensing), GeneRatio c(120/2500, 85/2500, 78/2500, 65/2500, 54/2500, 48/2500, 42/2500, 38/2500, 35/2500, 30/2500), BgRatio c(1500/8000, 1200/8000, 1100/8000, 900/8000, 800/8000, 700/8000, 600/8000, 500/8000, 450/8000, 400/8000), pvalue c(1.2e-12, 3.5e-09, 2.1e-07, 8.7e-06, 4.5e-05, 0.00012, 0.00045, 0.0012, 0.0033, 0.0088), p.adjust c(1.2e-10, 1.8e-07, 7.0e-06, 2.2e-04, 9.0e-04, 0.0020, 0.0060, 0.0150, 0.0330, 0.0660), qvalue c(1.0e-10, 1.5e-07, 5.8e-06, 1.8e-04, 7.5e-04, 0.0017, 0.0050, 0.0125, 0.0275, 0.0550), geneID c(geneA,geneB,geneC,geneD,geneE,geneF,geneG, geneH,geneI,geneJ,geneK,geneL, geneM,geneN,geneO,geneP,geneQ,geneR, geneS,geneT,geneU,geneV, geneW,geneX,geneY,geneZ,geneAA, geneB,geneF,geneL,geneP,geneV, geneC,geneI,geneO,geneU,geneAA,geneBB, geneD,geneJ,geneP,geneV,geneCC, geneE,geneK,geneQ,geneW,geneDD, geneG,geneR,geneX,geneBB,geneCC,geneDD) ) # 查看数据结构 head(kegg_result)这份数据包含了通路名称、富集基因比例、背景基因比例、P值、校正P值、Q值以及富集到的基因ID列表以逗号分隔。geneID列是连接气泡图和桑基图的关键。3. 数据预处理从原始结果到绘图就绪数据原始数据通常不能直接用于绘图尤其是geneID列。我们需要将其转换为适合ggplot2和networkD3的格式。3.1 清洗与筛选富集结果首先我们通常只关注最显著的一些通路比如选择校正P值p.adjust小于0.05的并按P值排序。# 筛选显著通路并按p.adjust升序排序 kegg_sig - kegg_result %% filter(p.adjust 0.05) %% arrange(p.adjust) # 为了演示我们选取前8条最显著的通路 top_pathways - head(kegg_sig, 8) print(top_pathways)3.2 关键步骤拆分基因列表构建“基因-通路”关联表这是生成桑基图的核心准备步骤。我们需要将geneID列一个包含多个基因的字符串拆分成多行每行代表一个基因与其所属通路的关系。# 拆分geneID列构建长格式的基因-通路关联表 gene_pathway_df - top_pathways %% select(Pathway, geneID) %% # 选择需要的列 separate_rows(geneID, sep ,) %% # 按逗号拆分一行变多行 rename(Gene geneID) %% # 重命名列 mutate(Gene trimws(Gene)) # 去除基因名两端的空格 # 查看关联表的前几行 head(gene_pathway_df)现在gene_pathway_df数据框的每一行都是一个明确的Gene属于某个Pathway的记录。这个格式完美契合桑基图对“源-目标”链接数据的要求。3.3 为绘图准备衍生数据对于气泡图我们可能希望用-log10(p.adjust)来表示显著性使得值越大点颜色越深越显著同时计算富集因子Enrichment Factor。# 为气泡图准备数据计算 -log10(p.adjust) 和富集因子 bubble_data - top_pathways %% mutate( -log10(p.adjust) -log10(p.adjust), EnrichmentFactor (GeneRatio) / (BgRatio) # 简化计算实际需注意格式转换 ) %% # 重新调整GeneRatio的格式便于理解 mutate(GeneCount as.numeric(sub(/.*, , GeneRatio))) # 提取基因数 head(bubble_data)至此我们得到了两个核心数据对象用于气泡图的bubble_data和用于桑基图的gene_pathway_df。它们同源但形态各异。4. 核心图表绘制ggplot2气泡图实战我们将使用ggplot2的图层语法逐步构建一个出版级的气泡图。4.1 基础气泡图绘制最基本的映射关系是X轴富集因子或GeneRatioY轴通路名称点大小基因数点颜色显著性。# 基础气泡图 p_bubble - ggplot(bubble_data, aes(x EnrichmentFactor, y reorder(Pathway, -log10(p.adjust)), # 按显著性排序通路 size GeneCount, color -log10(p.adjust))) geom_point(alpha 0.8) # 添加点图层设置透明度 scale_size_area(name Gene Count, max_size 12) # 控制点大小范围 scale_color_gradient(low blue, high red, name -log10(p.adjust)) # 设置颜色渐变 labs(x Enrichment Factor, y Pathway, title KEGG Pathway Enrichment Analysis) theme_bw(base_size 14) # 使用黑白主题设置基础字体大小 theme(axis.text.y element_text(size 10, color black), plot.title element_text(hjust 0.5, face bold)) # 标题居中加粗 print(p_bubble)这段代码会生成一个可用的气泡图。但我们可以做得更好。4.2 高级美化与定制科研图表讲究清晰、准确、美观。以下是一些常见的优化技巧# 高级美化版气泡图 p_bubble_enhanced - p_bubble # 1. 优化图例 guides(size guide_legend(order 1), color guide_colorbar(order 2)) # 2. 调整坐标轴和网格线 theme( panel.grid.major.y element_line(linetype dashed, color grey90), # 横向虚线网格 panel.grid.major.x element_blank(), # 去除纵向主网格线 panel.grid.minor element_blank(), # 去除次要网格线 axis.line.x element_line(color black), # X轴线 axis.ticks.y element_blank() # 去除Y轴刻度线 ) # 3. 扩展颜色标度让极值更突出 scale_color_gradientn( colours c(#4393C3, #FFD700, #D73027), name -log10(p.adjust) ) # 4. 添加数值标签可选在点上显示基因数 geom_text(aes(label GeneCount), color white, size 3, show.legend FALSE) print(p_bubble_enhanced)通过theme()函数你可以精细控制图表的每一个元素。现在你的气泡图已经具备了投稿期刊的潜力。5. 核心图表绘制networkD3桑基图实战桑基图需要一种特定的数据格式一个包含“链接”Links和“节点”Nodes的数据框。5.1 构建桑基图数据格式“链接”数据框需要三列source源节点索引、target目标节点索引、value链接权重通常为1。“节点”数据框需要一列name按顺序列出所有唯一的节点名。# 1. 准备节点列表包含所有唯一的基因和通路 # 注意节点顺序至关重要它将决定索引号。 all_genes - unique(gene_pathway_df$Gene) all_pathways - unique(gene_pathway_df$Pathway) node_names - c(all_genes, all_pathways) # 通常将“源”基因放在前面“目标”通路放在后面 # 创建节点数据框 nodes - data.frame(name node_names, stringsAsFactors FALSE) # 2. 准备链接数据框 # 为每个链接找到对应的源节点和目标节点索引 links - gene_pathway_df %% mutate( source match(Gene, node_names) - 1, # networkD3索引从0开始 target match(Pathway, node_names) - 1, value 1 # 每个链接的权重设为1 ) %% select(source, target, value) # 查看链接数据前几行 head(links)关键点networkD3要求索引从0开始。match()函数返回的是在node_names向量中的位置R索引从1开始所以需要减1。5.2 绘制交互式桑基图使用sankeyNetwork()函数传入链接和节点数据。# 绘制基础桑基图 sankey_plot - sankeyNetwork(Links links, Nodes nodes, Source source, Target target, Value value, NodeID name, units genes, # 链接的单位 fontSize 12, nodeWidth 20, height 600, width 900) # 在RStudio的Viewer窗格中显示交互式 sankey_plot # 保存为独立的HTML文件可在浏览器中打开并交互 saveNetwork(sankey_plot, file KEGG_Sankey_Diagram.html)运行后你会在RStudio的Viewer窗口看到一个可交互的流程图。鼠标悬停在节点基因或通路上会高亮显示所有与之相连的流悬停在链接上会显示详细信息。这极大地便利了数据探索。5.3 桑基图的美化与问题处理默认的桑基图可能颜色单一、节点拥挤。我们可以进行优化# 美化桑基图为基因和通路节点设置不同颜色 # 假设我们想用蓝色系表示基因橙色系表示通路 node_colors - c(rep(steelblue, length(all_genes)), # 基因节点颜色 rep(darkorange, length(all_pathways))) # 通路节点颜色 sankey_plot_enhanced - sankeyNetwork(Links links, Nodes nodes, Source source, Target target, Value value, NodeID name, units genes, fontSize 14, nodeWidth 25, nodePadding 15, # 增加节点间距 height 700, width 1000, colourScale JS(d3.scaleOrdinal().range([#999])), # 链接颜色 NodeGroup name, # 用于分组的列这里我们用名字但配合自定义颜色 LinkGroup source, # 链接按源节点分组着色 sinksRight TRUE) # 右侧节点对齐 # 直接修改HTML对象的样式是复杂的更简单的办法是在保存后编辑HTML或使用更高级的包。 # 一个实用技巧通过调整nodeWidth, nodePadding, height, width来改善布局。如果基因和通路数量过多桑基图会变得非常复杂难以阅读。最佳实践是先用气泡图筛选出最显著的少数几个通路如5-10个再用这些通路对应的基因子集来绘制桑基图。这正是“一码两图”工作流的精髓气泡图用于筛选桑基图用于深挖。6. 流程整合与自动化脚本将以上步骤整合到一个R脚本或函数中即可实现“一键出两图”。# 文件名kegg_dual_plot.R # 功能输入KEGG富集结果数据框输出气泡图和桑基图 generate_kegg_plots - function(enrichment_df, p_adjust_cutoff 0.05, top_n 8) { # 加载必要库如果在函数外未加载 library(ggplot2) library(dplyr) library(tidyr) library(networkD3) # 1. 数据筛选与排序 sig_data - enrichment_df %% filter(p.adjust p_adjust_cutoff) %% arrange(p.adjust) %% head(top_n) if (nrow(sig_data) 0) { stop(No significant pathways found with the given cutoff.) } # 2. 准备气泡图数据 bubble_data - sig_data %% mutate(-log10(p.adjust) -log10(p.adjust), GeneCount as.numeric(sub(/.*, , GeneRatio)), EnrichmentFactor GeneCount / as.numeric(sub(.*/, , BgRatio))) # 3. 绘制气泡图 p_bubble - ggplot(bubble_data, aes(x EnrichmentFactor, y reorder(Pathway, -log10(p.adjust)), size GeneCount, color -log10(p.adjust))) geom_point(alpha 0.7) scale_size_area(max_size 10, name Gene Count) scale_color_gradient(low lightblue, high red, name -log10(p.adjust)) labs(x Enrichment Factor, y NULL, title Top KEGG Enriched Pathways) theme_minimal(base_size 12) theme(axis.text.y element_text(size 10), plot.title element_text(hjust 0.5, face bold), legend.position right) # 4. 准备桑基图数据 gene_pathway_long - sig_data %% select(Pathway, geneID) %% separate_rows(geneID, sep ,) %% mutate(geneID trimws(geneID)) %% rename(Gene geneID) all_genes - unique(gene_pathway_long$Gene) all_pathways - unique(gene_pathway_long$Pathway) node_names - c(all_genes, all_pathways) nodes - data.frame(name node_names) links - gene_pathway_long %% mutate(source match(Gene, node_names) - 1, target match(Pathway, node_names) - 1, value 1) %% select(source, target, value) # 5. 绘制桑基图 sankey_plot - sankeyNetwork(Links links, Nodes nodes, Source source, Target target, Value value, NodeID name, fontSize 10, nodeWidth 20, height 500, width 800, sinksRight TRUE) # 6. 返回结果 return(list(bubble_plot p_bubble, sankey_plot sankey_plot, bubble_data bubble_data, sankey_data list(links links, nodes nodes))) } # 使用函数 # 假设你的富集结果在 my_kegg_results 数据框中 plots - generate_kegg_plots(my_kegg_results, p_adjust_cutoff 0.05, top_n 6) # 查看气泡图 print(plots$bubble_plot) # 查看并保存桑基图 plots$sankey_plot saveNetwork(plots$sankey_plot, file My_Analysis_Sankey.html)这个函数封装了核心流程你只需要提供自己的enrichment_df调整p_adjust_cutoff和top_n参数即可快速生成图表。7. 常见问题与排查指南在实际操作中你可能会遇到以下问题问题现象可能原因排查方式解决方案气泡图点的大小或颜色映射错误用于映射的列如GeneCount不是数值型。使用str(bubble_data)检查数据类型。用as.numeric()转换列或检查字符串提取逻辑。桑基图节点重叠布局混乱节点基因通路数量过多。检查length(node_names)。务必先筛选。只保留最显著的少数通路如Top 5-8及其基因。桑基图链接不显示或显示错误source/target索引错误或node_names顺序与链接不匹配。检查links数据框的source/target值是否在合理范围0到nrow(nodes)-1。确保node_names顺序是先所有基因后所有通路。仔细检查match()函数。separate_rows报错geneID列分隔符不一致或存在NA。使用unique(gene_pathway_df$geneID[1:5])查看分隔符。检查是否有NA。统一分隔符如全部改为逗号。用drop_na()删除包含NA的行。桑基图保存为HTML后无法交互用saveNetwork保存但用文本编辑器打开可能丢失依赖。在浏览器中打开保存的HTML文件。确保在浏览器中打开。如果网络受限可尝试保存为包含所有依赖的单个HTMLselfcontained TRUE。图形主题或字体不生效theme_*()设置被后续代码覆盖。检查代码顺序确保主题设置在最后。将theme()调整放在绘图语句的最后部分。最重要的建议始终从一个小规模的、可控的测试数据集开始比如只选3个通路确保每一步的代码都按预期工作然后再应用到全数据集上。8. 最佳实践与进阶技巧掌握了基础流程后这些技巧能让你的分析更上一层楼数据溯源与可复现性在脚本开头使用set.seed()保证随机过程可重复。使用sessionInfo()记录R和包的版本。将原始数据、处理脚本和最终图表放在同一个项目目录下。桑基图性能优化当基因数过多时桑基图会变得极其复杂。考虑在桑基图中只展示差异最显著的基因如logFC绝对值最大的Top 50而不是所有富集基因。使用networkD3的nodePadding和margin参数调整布局避免节点挤在一起。气泡图的美学定制使用scale_color_gradient2()可以设置中间色如白色让颜色对比更柔和。通过theme(legend.position “bottom”)将图例放在底部节省纵向空间。使用ggsave(“bubble_plot.png”, width8, height6, dpi300)导出高清图片用于投稿。结果的生物学解读气泡图帮你找到“什么”通路重要。桑基图帮你回答“为什么”重要——通过展示哪些核心基因同时参与了多个关键通路提示潜在的调控枢纽。将两图并列放在报告或论文中并配文说明左图展示了显著性排名靠前的通路右图揭示了这些通路之间通过共享基因形成的功能网络。扩展到其他富集分析本流程不仅适用于KEGG稍作修改即可用于GO、Reactome、MSigDB等任何提供基因列表的富集分析结果。关键在于结果表中需包含geneID这类列。从混乱的富集结果表格到直观的气泡图和揭示内在联系的桑基图你不仅完成了一次可视化升级更构建了一个从宏观显著性筛选到微观基因网络探查的完整分析叙事。这套基于R的“一码两图”工作流其价值在于将固定的分析模式转化为可复用的自动化脚本。下次当你拿到新的测序数据并完成富集分析后只需将结果文件路径指向这个脚本几分钟内就能获得可用于组会、报告或论文插图的专业图表。
返回列表