如果你在生物信息学分析中已经完成了差异表达基因的筛选下一步通常会做什么很多人的第一反应是去做GO和KEGG富集分析然后生成一张气泡图或柱状图把富集到的通路展示出来。这没错但问题在于这张静态图往往成了分析的终点——我们看到了哪些通路显著然后呢真正的分析价值往往藏在“然后”之后。一张气泡图告诉你“What”是什么通路富集了但它很难清晰地展示“How”这些基因是如何在不同通路间流动和协作的以及“So What”这对理解生物学过程意味着什么。这时你需要的不只是一张图而是一套能够串联起“发现-解释-呈现”的分析叙事。今天要讨论的就是如何用R语言从同一套KEGG富集分析结果出发生成两种互补的、极具信息量的可视化图形经典的KEGG气泡图和揭示基因流向的桑基图。这不仅仅是“一套代码出两图”的技术操作其核心在于通过两种视角的交叉验证与叙事补充将你的数据分析从“罗列结果”提升到“讲述故事”的层次。气泡图负责呈现显著性是分析的“锚点”桑基图负责揭示关联与流动是叙事的“桥梁”。下面我们就从为什么需要这种组合开始一步步拆解全流程。1. 为什么是气泡图桑基图超越单图呈现的分析叙事在深入代码之前我们必须先理解这两种图形各自承担的角色以及它们组合起来产生的“112”效应。这决定了我们整个分析流程的设计思路而不仅仅是机械地执行绘图命令。1.1 气泡图展示富集分析的“基本面”KEGG气泡图是富集分析结果最直观的呈现方式。它的X轴通常是基因比例GeneRatio或富集因子Fold EnrichmentY轴是通路名称点的大小代表基因数目颜色代表显著性P值或校正后的Q值。它的核心价值在于快速定位哪些通路最显著颜色最深哪些通路涉及的基因最多点最大富集的效果有多强点在X轴上的位置气泡图是一个优秀的“总结者”和“筛选者”。在报告或文章中它通常是展示富集分析结果的第一个也是必有的图形。它回答了“发生了什么”的基础问题。1.2 桑基图揭示基因与通路的“网络关系”桑基图是一种流图通过流的宽度来显示数据在不同节点间转移的数量。在生信分析中我们将它用于展示同一个基因集合在不同通路间的共享关系。它的核心价值在于揭示隐藏模式核心基因节点是哪些哪些基因参与了多个显著通路这些基因可能是调控多个生物学过程的关键枢纽。通路之间如何关联哪些通路共享了大量基因这暗示了这些通路在功能上可能存在协同或上下游关系。数据的流向是什么从“差异基因”这个源节点流向各个“通路”目标节点流的宽度直观展示了每个通路贡献的基因数量同时保留了基因多归属的信息。桑基图回答的是“如何发生”以及“内在联系是什么”的深层问题。它能帮你从一堆显著通路中识别出那个连接多个通路的、功能可能更核心的基因子集。1.3 组合叙事从“是什么”到“为什么”的进阶单独看气泡图你只知道通路A和通路B都重要。但结合桑基图你可能会发现通路A和通路B共享了30%的关键基因而这些基因恰好都受某个核心转录因子调控。这个洞察是任何单张图都无法提供的。因此我们的流程设计原则是以气泡图确定分析范围和重点通路以桑基图深入挖掘这些重点通路之间的功能联系和关键基因。代码实现上我们将使用clusterProfiler进行富集分析并绘制气泡图然后提取结果中的数据用ggalluvial或networkD3来构建桑基图。两者共享同一数据源确保叙事的一致性。2. 从差异基因到KEGG气泡图奠定分析基石我们的旅程从一组差异表达基因DEGs的列表开始。假设你已经通过DESeq2、edgeR等工具得到了deg_df这个数据框其中包含基因ID如ENTREZID或SYMBOL和显著性指标。2.1 环境准备与数据转换首先加载必要的R包。clusterProfiler是富集分析的核心org.Hs.eg.db是人类基因注释数据库请根据你的物种替换如org.Mm.eg.db对应小鼠。# 安装如果尚未安装 # BiocManager::install(c(clusterProfiler, org.Hs.eg.db, DOSE, enrichplot)) # 加载 library(clusterProfiler) library(org.Hs.eg.db) library(DOSE) library(ggplot2) library(enrichplot) # 用于增强可视化 # 假设你的差异基因数据框 deg_df 有一列叫 gene_symbol 存放基因符号 # 另有一列 log2FoldChange 和 pvalue head(deg_df)进行KEGG富集分析需要基因的Entrez ID。我们需要将基因符号SYMBOL转换为Entrez ID。# 基因ID转换 gene_ids - bitr(deg_df$gene_symbol, fromType SYMBOL, toType ENTREZID, OrgDb org.Hs.eg.db) # 查看转换结果可能会有部分基因无法映射而被剔除 head(gene_ids)2.2 执行KEGG富集分析使用转换后的Entrez ID列表进行富集分析。enrichKEGG函数会计算每个KEGG通路中富集到的基因情况。# 提取Entrez ID向量 gene_list - gene_ids$ENTREZID # 执行KEGG富集分析 kegg_enrich - enrichKEGG(gene gene_list, organism hsa, # ‘hsa’为人 ‘mmu’为鼠 keyType kegg, pvalueCutoff 0.05, pAdjustMethod BH, # Benjamini-Hochberg校正 qvalueCutoff 0.2) # 查看富集分析结果摘要 head(kegg_enrich)关键参数理解pvalueCutoff: 原始P值的阈值。通常设为0.05。pAdjustMethod: 多重检验校正方法。“BH”是最常用的FDR校正方法之一。qvalueCutoff: 校正后Q值的阈值。通常比P值阈值宽松一些如0.1或0.2以捕获更多有意义的通路。organism: 必须指定正确否则会报错或得到空结果。2.3 绘制与美化气泡图clusterProfiler提供了dotplot函数来快速绘制气泡图但我们可以通过ggplot2进行深度定制使其更符合出版要求。基础气泡图dotplot(kegg_enrich, showCategory15) ggtitle(KEGG Pathway Enrichment Analysis)进阶美化与定制通常我们需要对结果进行排序、筛选并调整颜色、尺寸等美学特征。# 1. 将富集结果转换为数据框便于操作 kegg_result - as.data.frame(kegg_enrich) # 2. 按显著性p.adjust排序并选择Top N条通路 top_n - 15 kegg_result_sorted - kegg_result[order(kegg_result$p.adjust), ][1:top_n, ] # 3. 自定义绘图 library(ggplot2) ggplot(kegg_result_sorted, aes(x GeneRatio, y reorder(Description, p.adjust))) geom_point(aes(size Count, color -log10(p.adjust))) scale_color_gradient(low blue, high red, name -log10(Adj.P), guide guide_colorbar(reverse FALSE)) scale_size_continuous(range c(3, 8), name Gene Count) labs(x Gene Ratio, y NULL, title Top Enriched KEGG Pathways, subtitle paste(Based on, length(gene_list), input genes)) theme_bw(base_size 12) theme(axis.text.y element_text(size 10, color black), plot.title element_text(face bold, hjust 0.5), legend.position right)代码解读与常见问题GeneRatio的计算GeneRatio是Count该通路中富集到的基因数除以BgRatio该通路在背景基因组中的总基因数。它比单纯的Count更能反映富集强度。排序使用reorder(Description, p.adjust)让Y轴通路按校正P值从小到大即显著性从高到低排列这是最合理的排序方式。颜色映射用-log10(p.adjust)作为颜色映射值使得P值越小越显著的点颜色越红假设梯度是蓝到红。图形尺寸scale_size_continuous中的range参数控制点的大小范围可根据通路数量调整。至此你已经得到了一张信息丰富、可用于发表的KEGG气泡图。它标识出了本次分析中最值得关注的Top N条通路。接下来我们将以这些通路为焦点构建桑基图。3. 构建桑基图可视化基因在通路间的共享网络桑基图需要三类数据节点Nodes和连接Links。在我们的场景中节点有两层第一层是“差异基因”第二层是“KEGG通路”。连接表示某个基因被富集到了某个通路。3.1 数据准备从富集结果中提取连接关系我们需要从kegg_enrich对象中提取出“基因-通路”的对应关系。# 提取富集结果中的基因-通路对应关系 # kegg_enrichgeneSets 存储了每个通路对应的基因列表Entrez ID # kegg_enrichresult 存储了结果表格 # 方法遍历每条通路获取其富集到的基因 library(dplyr) library(tidyr) # 获取我们之前选定的Top N条通路的ID top_pathway_ids - kegg_result_sorted$ID # 初始化一个空列表来收集数据 edge_list - list() for (path_id in top_pathway_ids) { # 获取该通路在富集结果对象中的索引 idx - which(kegg_enrichresult$ID path_id) if (length(idx) 0) { # 获取该通路富集到的基因Entrez ID genes_in_path - kegg_enrichgeneSets[[path_id]] # 获取通路名称 path_name - kegg_enrichresult$Description[idx] # 构建数据框基因ID - 通路名称 df - data.frame(source genes_in_path, target path_name, stringsAsFactors FALSE) edge_list[[path_id]] - df } } # 合并所有数据框 gene_pathway_edges - bind_rows(edge_list) # 查看连接关系的前几行 head(gene_pathway_edges)现在gene_pathway_edges数据框包含了source基因Entrez ID和target通路名称两列。但桑基图通常需要基因是可读的符号所以我们再转换一次。# 将Entrez ID转换回基因符号以便于识别 gene_pathway_edges$source - mapIds(org.Hs.eg.db, keys gene_pathway_edges$source, column SYMBOL, keytype ENTREZID) # 移除可能转换失败的NA值 gene_pathway_edges - na.omit(gene_pathway_edges) # 为了桑基图结构清晰我们引入一个虚拟的“DEGs”源节点 # 这意味着所有连接都将是 “DEGs” - “Gene” - “Pathway” sankey_links - data.frame( source c(rep(DEGs, nrow(gene_pathway_edges)), gene_pathway_edges$source), target c(gene_pathway_edges$source, gene_pathway_edges$target), value 1 # 每个连接的权重这里设为1 )3.2 使用ggalluvial绘制静态桑基图ggalluvial是基于ggplot2的扩展包可以绘制漂亮的静态桑基图或称冲击图。# 安装并加载 # install.packages(ggalluvial) library(ggalluvial) # 准备节点因子确保顺序可选但有助于图形美观 node_levels - unique(c(sankey_links$source, sankey_links$target)) sankey_links$source - factor(sankey_links$source, levels node_levels) sankey_links$target - factor(sankey_links$target, levels node_levels) # 绘制桑基图 ggplot(sankey_links, aes(axis1 source, axis2 target, y value)) geom_alluvium(aes(fill target), # 按目标节点填充颜色 alpha 0.7, width 1/8) geom_stratum(width 1/8, fill grey80, color grey) geom_text(stat stratum, aes(label after_stat(stratum)), size 3) scale_x_discrete(limits c(Source, Target), expand c(0.05, 0.05)) scale_fill_viridis_d(option C, guide none) # 使用viridis配色 labs(title Gene-Pathway Sankey Diagram, subtitle Flow from DEGs to Top Enriched KEGG Pathways) theme_void() theme(plot.title element_text(hjust 0.5, face bold), legend.position none)静态图优缺点优点输出为PDF/SVG等矢量格式印刷质量高与ggplot2生态无缝集成样式控制灵活。缺点当节点和连接非常多时图形会变得拥挤难以交互查看。3.3 使用networkD3绘制交互式桑基图对于更复杂的网络交互式图表体验更好。networkD3包可以生成基于D3.js的HTML交互图。# 安装并加载 # install.packages(networkD3) library(networkD3) # 为networkD3准备数据需要节点列表和连接列表 # 节点列表所有唯一节点的索引和名称 nodes - data.frame(name unique(c(sankey_links$source, sankey_links$target)), stringsAsFactors FALSE) nodes$id - 0:(nrow(nodes) - 1) # 连接列表用节点索引代替名称 links - sankey_links links$source_id - match(links$source, nodes$name) - 1 # 索引从0开始 links$target_id - match(links$target, nodes$name) - 1 # 绘制交互式桑基图 sankeyNetwork(Links links, Nodes nodes, Source source_id, Target target_id, Value value, NodeID name, fontSize 12, nodeWidth 20, sinksRight FALSE) # 让最右侧的节点也左对齐布局更平衡运行这段代码会在RStudio的Viewer面板或浏览器中打开一个交互式图表。你可以用鼠标拖动节点悬停查看连接细节。这对于探索哪些基因是连接多个通路的核心枢纽非常直观。4. 整合、解读与进阶应用从图形到生物学故事两张图都生成了但工作只完成了一半。更重要的是如何整合解读并基于此进行深入分析。4.1 整合解读框架不要孤立地看两张图。建议按以下流程进行整合解读定位核心通路在气泡图中锁定P值最小、GeneRatio最高的几个通路如Top 3-5。这些是你的“一级焦点”。识别枢纽基因在桑基图中找到那些从“DEGs”流出、同时连接到多个“一级焦点”通路的基因。这些基因是潜在的关键调控因子。将鼠标悬停在交互式桑基图的连接上或从静态图数据中筛选记录下这些基因符号。验证与深挖回到你的原始差异表达数据查看这些枢纽基因的log2FoldChange。它们是上调还是下调使用STRING数据库或类似工具检查这些枢纽基因编码的蛋白质之间是否存在已知的相互作用PPI网络。一个紧密连接的枢纽基因簇说服力更强。查阅文献确认这些枢纽基因在你所研究的生物学背景如特定癌症、发育阶段中是否已被报道具有核心作用。形成叙事将发现串联起来。例如“我们的分析发现X信号通路和Y代谢通路在YY条件下被显著激活气泡图。进一步分析揭示基因A和基因B是连接这两条通路的关键节点桑基图。已知基因A编码的蛋白能调控Y通路中的关键酶而本研究中基因A显著上调这可能是导致X和Y通路协同变化的核心机制。”4.2 进阶应用与定制基础流程跑通后可以根据需求进行深度定制筛选连接桑基图可能因为基因太多而显得杂乱。可以只保留那些连接到至少2个通路的基因或者只保留表达变化最显著如|log2FC| 2的基因来绘图使图形更清晰重点更突出。分层桑基图如果你的分析涉及多个比较组如处理vs对照时间点T1 vs T2可以构建更复杂的分层桑基图展示基因在不同条件、不同通路间的动态变化。与其它数据整合将桑基图中识别出的枢纽基因与其拷贝数变异CNV、甲基化状态或生存分析数据关联进行多组学层面的整合分析。自动化脚本将上述流程封装成一个函数输入差异基因列表自动输出气泡图、桑基图的数据和图形以及一个包含枢纽基因的表格大大提高重复分析效率。4.3 常见问题排查桑基图节点过多图形混乱解决严格限制气泡图中通路的数量如Top 10。在构建桑基图连接数据前过滤掉只出现在一个通路中的基因即非共享基因。代码示例gene_count - table(gene_pathway_edges$source); shared_genes - names(gene_count[gene_count 1]); gene_pathway_edges_filtered - subset(gene_pathway_edges, source %in% shared_genes)clusterProfiler富集分析结果为空检查基因ID类型是否正确organism参数是否正确输入的基因列表是否有效成功转换为Entrez IDP值阈值是否太严格尝试放宽pvalueCutoff和qvalueCutoff确保使用的注释数据库OrgDb与物种匹配。ggalluvial绘图时出现警告或错误检查数据格式是否正确aes中的axis1,axis2是否对应数据框中的列名所有用于分组的变量如fill是否已转换为因子factor并设置了合理的水平levels交互式桑基图无法显示或节点错位检查links数据框中的source_id和target_id是否与nodes数据框中的id正确对应从0开始计数。确保Value列的值是数值型。通过这套组合拳你的KEGG富集分析将不再是一张孤立的、陈述性的图片而是一个包含发现、关联和假设的完整分析故事。气泡图提供了故事的目录和摘要而桑基图则揭示了章节之间隐藏的人物关系和情节线索。掌握这种方法意味着你掌握了将生信数据转化为生物学洞察的更强大工具。