ARTICLE DETAIL

资讯详情

深耕网站建设、视觉设计与SEO优化的一线实战洞察。

R语言生信分析:KEGG富集结果可视化与桑基图绘制实战

R语言生信分析:KEGG富集结果可视化与桑基图绘制实战 1. 先搞清楚这套流程到底能解决什么问题如果你正在做转录组、蛋白组或者代谢组学分析做完差异表达或富集分析后下一步通常就是可视化。KEGG富集分析的结果最常见的就是画个气泡图或者柱状图展示一下富集到的通路和显著性。但有时候你可能会想这些通路之间有什么关系我的差异基因是怎么在不同通路间“流动”的这时候桑基图Sankey Diagram就能很直观地展示这种“从基因到通路”的归属关系。所以这个主题的核心价值是用同一套数据通常是KEGG富集分析结果和基因-通路映射关系一次性生成两种互补的图表——展示整体富集概况的气泡图和展示具体归属关系的桑基流向图。这比单独画两张图更高效也更能从宏观哪些通路重要和微观基因如何分布两个层面讲好数据故事。适合谁看主要是刚入门R语言生信分析已经能跑通差异分析和富集分析但在结果可视化和深度解读上想更进一步的研究者。你不用是R语言高手但需要能理解data.frame、ggplot2基本操作和富集分析结果的基本结构。最关键的一点是我建议你不要把这两个图看成完全独立的步骤。它们的核心是共享同一份“基因-通路”关联数据。理解了这一点后面的代码组织就会清晰很多。2. 环境准备与核心数据理解别急着写代码在动手敲代码之前先把环境和数据搞清楚。这一步做扎实了后面能省掉80%的报错。2.1 R环境与包管理别让包安装卡住你首先确保你的R版本不要太老R 4.0以上比较稳妥。新手最容易卡住的地方就是包安装失败。根据搜索热词里提到的“causalweight包为何装不上 r语言”这提醒我们有些包的安装依赖特定系统库或者需要从特定源安装。对于我们要用到的可视化包主要依赖如下ggplot2: 画气泡图的核心基本都会装。ggsankey/networkD3/plotly: 用于绘制桑基图。ggsankey生成的是静态图集成在ggplot2体系里风格统一适合放入论文。networkD3和plotly生成的是交互式HTML图表可以在浏览器里拖动节点适合探索性分析和汇报。这里我以ggsankey为例因为它和ggplot2语法一致学习成本低。dplyr/tidyr: 用于数据整理和转换几乎是现代R数据分析的标配。clusterProfiler: 如果你是用这个神包做的KEGG富集分析那你的结果对象直接可以用。它也是数据来源的关键。stringr: 处理通路名称、基因ID等文本信息非常有用。安装命令很简单但要注意网络问题# 设置CRAN镜像国内用户必备能极大提升安装速度和成功率 options(repos c(CRAN https://mirrors.tuna.tsinghua.edu.cn/CRAN/)) # 安装必要包 install.packages(c(ggplot2, dplyr, tidyr, stringr, ggsankey)) # 如果是Bioconductor的包如clusterProfiler if (!require(BiocManager, quietly TRUE)) install.packages(BiocManager) BiocManager::install(clusterProfiler)如果某个包比如热词里提到的causalweight安装失败先别慌。错误信息通常会提示缺少什么系统依赖比如在Linux上或者尝试从GitHub安装(remotes::install_github)。但对于ggsankey这类流行包直接从CRAN安装通常没问题。2.2 理解你的数据从富集结果到绘图数据框这是最核心的一步。你的输入数据长什么样决定了代码怎么写。通常你有两种起点clusterProfiler的富集结果对象这是最理想的情况。假设你的结果对象叫kegg_result。一个包含富集分析结果的data.frame通常是CSV或Excel文件读入很多在线工具或其它软件的分析结果会导出成表格。无论哪种来源你最终都需要整理出两个核心数据框用于气泡图的数据框 (bubble_df)需要至少包含以下几列Description: KEGG通路描述。GeneRatio: 富集到该通路的基因数 / 背景基因数例如10/100。注意clusterProfiler的结果里GeneRatio列是字符型需要转换。BgRatio: 通路中总基因数 / 背景基因数。pvalue/p.adjust/qvalue: 显著性P值或校正后的P值。Count: 富集到该通路的基因数目。这个通常由GeneRatio计算得来但有时也会直接有一列。用于桑基图的数据框 (sankey_df)需要整理成一个“长格式”数据框至少包含两列或三列source: 源节点通常是基因ID比如gene1,gene2。target: 目标节点通常是通路ID或通路描述比如hsa04110,Cell cycle。value: (可选) 边的权重如果每个基因-通路连接权重为1这列可以省略或统一设为1。但有些桑基图函数需要。关键经验桑基图的数据准备比气泡图麻烦。你需要从富集结果中把“哪个基因属于哪个通路”的映射关系提取出来。clusterProfiler的结果对象中有一个geneID列里面是用/分隔的基因列表这就是你的原始材料。3. 实战第一步从数据整理到气泡图生成我们先从相对简单的气泡图开始这个过程也会帮我们整理出后续桑基图需要的数据。3.1 数据整理与清洗假设我们有一个clusterProfiler生成的kegg_result对象类型是enrichResult。我们首先把它转换成数据框并整理出气泡图所需列。# 加载包 library(clusterProfiler) library(dplyr) library(tidyr) library(stringr) library(ggplot2) # 1. 将富集结果转为数据框并筛选显著通路例如 padj 0.05 bubble_df - as.data.frame(kegg_result) %% filter(p.adjust 0.05) %% # 根据你的阈值调整 # 2. 计算富集基因数量(Count)和基因比例(GeneRatio数值) mutate( Count as.numeric(str_split(GeneRatio, /, simplify TRUE)[,1]), GeneRatio_num Count / as.numeric(str_split(GeneRatio, /, simplify TRUE)[,2]), # 3. 通常我们按富集因子(Enrichment Factor)或p值排序这里按p.adjust升序越小越显著 Description factor(Description, levels rev(Description[order(p.adjust)])) ) %% # 4. 选择需要的列并可能限制显示通路的数量如前20条 select(Description, Count, GeneRatio_num, p.adjust) %% head(20) # 只展示最显著的20条通路避免图太拥挤 # 查看整理好的数据框 head(bubble_df)为什么这么操作str_split用于拆分GeneRatio如”10/100″得到分子和分母。将通路描述(Description)转换为因子并排序是为了让气泡图里的Y轴通路顺序是固定的最显著的在顶部或底部。限制显示数量是因为通路太多会导致图上的点重叠严重可读性变差。前20条是一个常用选择。3.2 绘制基础气泡图有了bubble_df用ggplot2画气泡图就非常直接了。p_bubble - ggplot(bubble_df, aes(x GeneRatio_num, y Description)) geom_point(aes(size Count, color -log10(p.adjust))) # 点的大小代表基因数颜色代表显著性 scale_color_gradient(low blue, high red, name -log10(adj.P)) # 颜色梯度 scale_size_continuous(range c(3, 8), name Gene Count) # 点大小范围 labs( x Gene Ratio, y KEGG Pathway, title KEGG Pathway Enrichment Analysis ) theme_bw() theme( axis.text.y element_text(size 10), axis.title element_text(size 12, face bold), plot.title element_text(hjust 0.5, size 14, face bold) ) # 显示图形 print(p_bubble) # 保存图形 ggsave(KEGG_bubble_plot.png, p_bubble, width 10, height 8, dpi 300)参数解释与避坑点aes(size Count, color -log10(p.adjust)): 这是气泡图的精髓。size映射到基因数量直观显示通路规模color映射到-log10(p.adjust)使得P值越小越显著的颜色越“热”如红色。scale_size_continuous(range c(3, 8)): 调整点的大小范围。如果Count值跨度很大可以调整这个范围让图更美观。theme_bw(): 经典的白底黑线主题适合出版。常见问题如果通路名称太长Y轴的标签会重叠。可以用stringr::str_wrap来截断或换行或者调整theme(axis.text.y element_text(...))中的size和margin。4. 实战第二步准备桑基图数据并绘制气泡图告诉我们“哪些通路重要”桑基图则要展示“重要的基因具体流向了哪些通路”。所以我们需要从原始数据中提取基因-通路的对应关系。4.1 从富集结果中提取基因-通路映射这是最关键且稍显繁琐的一步。我们需要把geneID列包含多个基因拆分成多行。# 继续使用kegg_result对象 # 1. 同样先转为数据框并筛选显著通路 sankey_raw_df - as.data.frame(kegg_result) %% filter(p.adjust 0.05) %% select(ID, Description, geneID, Count) %% head(10) # 桑基图节点不宜过多先取前10条显著通路演示 # 2. 拆分geneID列将用/分隔的基因字符串拆分成多行 sankey_links - sankey_raw_df %% # 使用separate_rows将一行的多个基因拆成多行 separate_rows(geneID, sep /) %% # 重命名列形成 source (基因) - target (通路) 的连接 rename(gene geneID, pathway Description) %% select(gene, pathway) %% # 每个连接基因-通路对的value设为1 mutate(value 1) # 查看连接数据 head(sankey_links)现在sankey_links数据框里每一行代表一个“基因-通路”连接。例如genepathwayvaluegeneACell cycle1geneAp53 signaling pathway1geneBCell cycle1注意一个基因可能富集到多个通路如上例的geneA这在桑基图中表现为一个源节点连接到多个目标节点这正是桑基图要展示的“分流”效果。4.2 使用ggsankey绘制静态桑基图ggsankey包提供了geom_sankey和geom_sankey_label等函数可以无缝融入ggplot2。library(ggsankey) # 为了绘图需要将数据框转换为ggsankey需要的格式 # 使用make_long函数指定节点列这里gene是x pathway是next_x sankey_data_for_plot - sankey_links %% make_long(gene, pathway) # 这个函数会将数据转换成x, next_x, node, next_node格式 # 绘制桑基图 p_sankey - ggplot(sankey_data_for_plot, aes(x x, next_x next_x, node node, next_node next_node, fill node, # 按节点填充颜色 label node)) # 节点标签 geom_sankey(flow.alpha 0.5, # 流线的透明度 node.color black, # 节点边框颜色 show.legend FALSE) # 通常节点太多图例没意义 geom_sankey_label(size 3, color black, fill white) # 添加节点标签 theme_void() # 清空背景和坐标轴桑基图通常不需要 theme(plot.margin unit(c(1, 1, 1, 1), cm)) # 增加边距防止标签被切 # 显示图形 print(p_sankey) ggsave(KEGG_sankey_plot.png, p_sankey, width 14, height 10, dpi 300)参数解释与避坑点make_long(): 是ggsankey的关键函数它把“源-目标”格式的数据转换成绘图需要的长格式。flow.alpha: 调整流线连接线的透明度当流线很多时适当调低透明度如0.5可以避免画面过脏。theme_void(): 桑基图通常不需要坐标轴用这个主题清空。最大的坑节点过多。如果你把成百上千个基因和几十条通路都放进去图形会变成一团乱麻根本无法阅读。务必筛选只保留最显著的前N条通路比如前10或前15并且只保留富集到这些通路里的基因。这是绘图前最重要的数据过滤步骤。标签重叠节点标签(geom_sankey_label)可能会重叠。ggsankey对标签位置的处理有时不够智能。如果重叠严重可以考虑使用geom_sankey_text并手动调整或者转向交互式绘图库如plotly它允许手动拖动节点。4.3 进阶使用plotly绘制交互式桑基图如果你需要探索性分析交互式图表更合适。plotly库功能强大。library(plotly) # 为plotly准备数据需要定义节点列表和连接列表 # 1. 获取所有唯一的节点基因名通路名 nodes - unique(c(sankey_links$gene, sankey_links$pathway)) # 创建节点数据框plotly需要索引 node_df - data.frame(name nodes, id 0:(length(nodes)-1)) # 2. 创建连接数据框将基因和通路名称映射到索引 links_df - sankey_links %% left_join(node_df, by c(“gene” “name”)) %% rename(source id) %% left_join(node_df, by c(“pathway” “name”)) %% rename(target id) %% select(source, target, value) # 3. 绘制交互式桑基图 fig - plot_ly( type “sankey”, orientation “h”, # 水平流向 node list( label node_df$name, pad 15, thickness 20, line list(color “black”, width 0.5) ), link list( source links_df$source, target links_df$target, value links_df$value ) ) # 显示图形在RStudio的Viewer或浏览器中 fig # 保存为独立的HTML文件 htmlwidgets::saveWidget(as_widget(fig), “KEGG_sankey_interactive.html”)交互式图的优势你可以用鼠标悬停查看每个节点或连接的详细信息拖动节点来重新布局这对于理解复杂的关系网络非常有帮助。生成的HTML文件可以单独打开方便在报告或网页中展示。5. 将两图生成流程封装与通用化上面是分步演示。在实际项目中我们肯定希望有一个函数或一套脚本输入富集结果就能输出两张图。这里提供一个简化的流程框架。5.1 创建一个整合绘图函数你可以创建一个R脚本例如plot_kegg_dual.R里面包含一个主函数。# 生成KEGG气泡图和桑基图 # # param enrich_obj clusterProfiler的富集结果对象 # param padj_cutoff 显著性阈值默认0.05 # param top_n_pathway 用于绘图的前N条通路气泡图默认20桑基图默认10 # param bubble_outfile 气泡图输出文件名 # param sankey_outfile 桑基图输出文件名静态 # param sankey_interactive_outfile 交互式桑基图输出文件名可选 # # return 一个列表包含气泡图对象和桑基图连接数据 generate_kegg_dual_plots - function(enrich_obj, padj_cutoff 0.05, top_n_bubble 20, top_n_sankey 10, bubble_outfile “kegg_bubble.png”, sankey_outfile “kegg_sankey.png”, sankey_interactive_outfile NULL) { library(dplyr); library(tidyr); library(stringr); library(ggplot2); library(ggsankey) # 1. 数据准备 df_raw - as.data.frame(enrich_obj) # 2. 生成气泡图数据并绘图 bubble_df - df_raw %% filter(p.adjust padj_cutoff) %% mutate( Count as.numeric(str_split(GeneRatio, “/“, simplify TRUE)[,1]), GeneRatio_num Count / as.numeric(str_split(GeneRatio, “/“, simplify TRUE)[,2]), Description factor(Description, levels rev(Description[order(p.adjust)])) ) %% arrange(p.adjust) %% head(top_n_bubble) p_bubble - ggplot(bubble_df, aes(x GeneRatio_num, y Description)) geom_point(aes(size Count, color -log10(p.adjust))) scale_color_gradient(low “blue”, high “red”, name “-log10(adj.P)”) scale_size_continuous(range c(3, 8), name “Gene Count”) labs(x “Gene Ratio”, y “”, title “KEGG Pathway Enrichment”) theme_bw() theme(axis.text.y element_text(size 9), plot.title element_text(hjust 0.5)) ggsave(bubble_outfile, p_bubble, width 9, height 6, dpi 300) message(“Bubble plot saved to: “, bubble_outfile) # 3. 生成桑基图数据并绘图 sankey_raw_df - df_raw %% filter(p.adjust padj_cutoff) %% arrange(p.adjust) %% head(top_n_sankey) %% select(ID, Description, geneID) sankey_links - sankey_raw_df %% separate_rows(geneID, sep “/“) %% rename(gene geneID, pathway Description) %% select(gene, pathway) %% mutate(value 1) # 如果连接数据为空则跳过桑基图绘制 if (nrow(sankey_links) 0) { warning(“No significant pathways found for Sankey diagram after filtering.”) return(list(bubble_plot p_bubble, sankey_links NULL)) } sankey_data - sankey_links %% make_long(gene, pathway) p_sankey - ggplot(sankey_data, aes(x x, next_x next_x, node node, next_node next_node, fill node, label node)) geom_sankey(flow.alpha 0.5, node.color “black”, show.legend FALSE) geom_sankey_label(size 2.5, color “black”, fill “white”) theme_void() theme(plot.margin unit(c(1, 1, 1, 1), “cm”)) ggsave(sankey_outfile, p_sankey, width 12, height 8, dpi 300) message(“Static Sankey plot saved to: “, sankey_outfile) # 4. 可选生成交互式桑基图 if (!is.null(sankey_interactive_outfile)) { library(plotly) # … (此处插入前面plotly的代码逻辑生成fig并保存为HTML) … message(“Interactive Sankey plot saved to: “, sankey_interactive_outfile) } # 返回结果 return(list(bubble_plot p_bubble, sankey_plot p_sankey, links_data sankey_links)) } # 使用示例 # result - generate_kegg_dual_plots(kegg_result, # top_n_sankey 8, # bubble_outfile “my_bubble.pdf”, # 支持pdf, png等格式 # sankey_outfile “my_sankey.pdf”)5.2 通用化建议处理非clusterProfiler的输入数据如果你的数据不是来自clusterProfiler而是一个普通的data.frame比如列名是Pathway,PValue,Genes基因列表用逗号分隔。你需要调整数据整理的代码。核心修改点在数据提取部分# 假设你的数据框叫 my_df bubble_df - my_df %% filter(PValue 0.05) %% mutate( # 计算Count可能需要从Genes列数逗号 Count str_count(Genes, “,”) 1, GeneRatio_num Count / background_gene_count, # 你需要知道背景基因总数 Pathway factor(Pathway, levels rev(Pathway[order(PValue)])) ) # 桑基图连接数据 sankey_links - my_df %% filter(PValue 0.05) %% head(10) %% separate_rows(Genes, sep “,”) %% # 按逗号拆分 rename(gene Genes, pathway Pathway) %% select(gene, pathway) %% mutate(value 1)关键确保你能从数据中准确提取出“基因列表”和“通路”的对应关系以及计算或获取GeneRatio和Count。6. 常见问题排查与图形优化在实际运行中你可能会遇到以下问题。按照这个顺序排查大部分都能解决。6.1 图形渲染或保存问题报错Error in geom_sankey()或could not find function “make_long”原因ggsankey包没有正确安装或加载。解决确认已安装(install.packages(“ggsankey”))并加载(library(ggsankey))。注意包名是ggsankey不是sankey。桑基图节点/标签重叠严重看不清原因数据太多超过了静态图的表达能力。解决严格筛选这是最有效的办法。将top_n_sankey参数调小比如只画最显著的5-8条通路。桑基图不适合展示太多节点。调整图形尺寸保存时增加width和height如width16, height12给标签更多空间。改用交互式图这是解决重叠问题的终极方案。用plotly生成HTML然后手动拖动和缩放查看。气泡图通路名称太长Y轴标签重叠解决在theme(axis.text.y …)中使用element_text(angle 0, hjust 1)调整或者用stringr::str_wrap(Description, width 40)在绘图前将长通路名自动换行。6.2 数据相关问题桑基图数据sankey_links为空画不出图原因top_n_sankey设置过大但显著通路没那么多或者padj_cutoff太严格没有显著通路。解决检查df_raw %% filter(p.adjust padj_cutoff) %% nrow()有多少行。确保筛选后有数据。可以先画气泡图确认有多少条显著通路。geneID列是空的或是NA原因有些富集分析结果可能不包含具体的基因列表或者你在运行enrichKEGG时没有设置universe或gene参数导致。解决确保你的富集分析步骤包含了基因ID信息。对于clusterProfiler检查enrichKEGG()函数的参数是否正确结果对象的geneID列是否有内容。一个基因出现在过多通路中导致桑基图像“扫把”现象某个基因节点伸出大量流线连接几乎所有通路图形失去重点。原因该基因可能是一个广泛表达的“管家基因”或者富集分析本身不够特异。解决这更多是生物学问题。可以从分析角度在富集前过滤掉低表达或变化不显著的基因。从绘图角度可以尝试在桑基图中过滤掉连接数超过某个阈值比如5的基因但需谨慎因为这可能掩盖真实生物学信息。6.3 性能与美化数据量很大时绘图速度慢解决对于气泡图限制top_n_bubble。对于桑基图必须限制top_n_sankey。这是保证可读性和性能的关键。交互式plotly图在处理数百个节点时也可能变慢需要权衡。想自定义颜色气泡图使用scale_color_gradient2或scale_color_viridis_c等函数可以更换颜色渐变方案。桑基图(ggsankey)aes(fill node)会根据节点自动分配颜色。你可以通过scale_fill_manual(values my_colors)来手动指定颜色向量my_colors但需要颜色数量与节点数匹配操作较复杂。一个更简单的方法是aes(fill x)这样只根据层级基因层或通路层来分配两种颜色。这套“一码出两图”的流程其核心思想在于数据驱动的可视化。一旦你整理好了标准的富集结果数据框生成气泡图是水到渠成。而桑基图所需的核心映射关系基因-通路其实已经蕴含在同一个数据源里只需要一次separate_rows的转换就能提取出来。我建议你在自己的项目上先确保富集分析结果本身是可靠的然后用前5-10条最显著的通路来尝试桑基图。先跑通再调整最后考虑封装成可复用的函数。这样你的下一次KEGG可视化就不仅仅是发一张图而是讲一个更有层次的数据故事了。
返回列表