你的clusterProfiler结果只用了4维?试试这个桑吉气泡图R包/代码复现教程
突破四维限制:用R打造桑吉气泡图实现KEGG富集五维可视化
在基因功能富集分析领域,气泡图一直是展示KEGG通路或GO富集结果的经典选择。传统的气泡图通过四个维度传递信息:Y轴的通路名称、X轴的富集程度、点的大小表示基因数量以及颜色深浅对应p值。但当我们面对clusterProfiler等工具生成的富集结果时,往往会发现一个关键信息被遗漏——具体的基因列表。这些基因名称通常被压缩在表格的geneID列中,以斜杠分隔的字符串形式存在,无法在标准可视化中直接展现。
1. 为什么需要五维桑吉气泡图
生物信息分析的核心价值在于建立基因与功能之间的关联网络。传统的四维气泡图虽然能够展示通路的重要性排序和统计显著性,但却切断了基因与通路之间的具体联系。这种信息断裂在实际研究中可能导致几个关键问题:
- 验证困难:当发现某个通路显著富集时,研究者需要额外步骤查看具体包含哪些基因
- 模式识别障碍:无法直观发现多个通路共享的关键基因
- 结果解释局限:难以评估富集结果是否由少数高频基因驱动
桑吉气泡图通过引入第五维度——基因名称及其与通路的隶属关系,完美解决了这些问题。这种可视化技术结合了桑吉图的网络属性和气泡图的统计表达能力,呈现出更完整的生物学故事。
技术对比:与微生信等在线工具相比,R语言实现具有三大优势:
- 完全可定制:从字体大小到颜色映射均可精细调整
- 可重复流程:代码化分析确保结果可复现
- 扩展性强:可轻松整合到自动化分析流程中
2. 数据准备与结构重塑
要实现桑吉气泡图,首先需要正确处理clusterProfiler的输出结果。标准的富集分析结果通常包含以下关键列:
library(clusterProfiler) # 示例数据加载 data(geneList, package = "DOSE") data(geneSets, package = "DOSE") ego <- enrichGO(gene = names(geneList)[1:100], OrgDb = org.Hs.eg.db, keyType = "ENTREZID", ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.01, qvalueCutoff = 0.05) enrich_df <- as.data.frame(ego)原始数据需要转换为适合桑吉图的长格式。以下代码演示如何拆分geneID列并重构数据:
library(tidyr) library(dplyr) # 转换数据为桑吉图所需格式 sankey_data <- enrich_df %>% separate_rows(geneID, sep = "/") %>% select(Description, geneID, GeneRatio, pvalue, Count) %>% mutate(link_strength = 1) # 桑吉图连接强度关键步骤说明:
separate_rows将geneID列按分隔符拆分为多行- 选择需要的列并添加连接强度变量
- 确保每个基因-通路对都成为独立观测
提示:对于大型数据集,建议在拆分前过滤掉不显著的通路以减少计算负担
3. 构建基础气泡图
在添加桑吉元素前,我们需要先创建标准的四维气泡图作为基础。ggplot2提供了高度灵活的图形构建能力:
library(ggplot2) library(scales) # 计算GeneRatio的数值形式 enrich_df <- enrich_df %>% mutate(GeneRatio_num = sapply(strsplit(GeneRatio, "/"), function(x) as.numeric(x[1])/as.numeric(x[2]))) # 创建基础气泡图 base_plot <- ggplot(enrich_df, aes(x = GeneRatio_num, y = reorder(Description, GeneRatio_num))) + geom_point(aes(size = Count, color = p.adjust)) + scale_color_gradient(low = "red", high = "blue", name = "Adjusted p-value", limits = c(0, 0.05)) + scale_size(range = c(3, 10), name = "Gene count") + labs(x = "Gene Ratio", y = "", title = "KEGG Pathway Enrichment") + theme_minimal(base_size = 12) + theme(axis.text.y = element_text(size = 10), panel.grid.major.y = element_line(linetype = "dotted"))这段代码实现了:
- GeneRatio从"a/b"格式转换为数值
- 通路按富集程度排序
- 点的大小和颜色分别映射到基因数和p值
- 精心调整的视觉元素提升可读性
4. 集成桑吉图元素
桑吉图的核心是展示基因与通路之间的流动关系。我们将使用ggalluvial包在气泡图左侧添加这一维度:
library(ggalluvial) # 创建桑吉图层 sankey_layer <- ggplot(sankey_data, aes(y = link_strength, axis1 = geneID, axis2 = Description)) + geom_alluvium(aes(fill = Description), alpha = 0.6, width = 1/8) + geom_stratum(width = 1/8, fill = "grey90", color = "grey50") + geom_text(stat = "stratum", aes(label = after_stat(stratum)), size = 3, angle = 45, hjust = 0) + scale_x_discrete(limits = c("Genes", "Pathways"), expand = c(0.05, 0.05)) + theme_void() + theme(legend.position = "none")关键参数说明:
geom_alluvium绘制流动曲线,展示基因-通路关联geom_stratum和geom_text添加分层标签width参数控制桑吉元素的宽度比例theme_void创建干净背景便于组合
5. 图形组合与美化
使用patchwork包将两个图形无缝组合:
library(patchwork) final_plot <- sankey_layer + base_plot + plot_layout(widths = c(1, 2)) + plot_annotation(title = "Five-dimensional Enrichment Visualization", subtitle = "Combining Sankey relationships with bubble plot statistics") # 输出图形 ggsave("sankey_bubble.png", final_plot, width = 16, height = 10, dpi = 300)组合技巧:
- 调整widths参数平衡桑吉图与气泡图的比例
- 添加整体标题说明图形含义
- 控制输出尺寸确保所有元素清晰可读
注意:当通路数量较多时,可能需要调整图形高度或缩小字体以避免标签重叠
6. 高级定制技巧
基础图形完成后,我们可以通过多种方式提升其信息量和美观度:
颜色方案优化:
# 使用viridis色盲友好调色板 base_plot <- base_plot + scale_color_viridis_c(option = "plasma", direction = -1, name = "-log10(p.adjust)", trans = "log10")交互式探索:
library(plotly) ggplotly(final_plot) %>% layout(hoverlabel = list(bgcolor = "white"), margin = list(l = 150)) # 增加左侧边距基因频率标记:
# 在桑吉图中高亮高频基因 gene_counts <- sankey_data %>% count(geneID) %>% arrange(desc(n)) sankey_data <- sankey_data %>% mutate(gene_importance = case_when( geneID %in% gene_counts$geneID[1:5] ~ "High", geneID %in% gene_counts$geneID[6:15] ~ "Medium", TRUE ~ "Low" ))实际项目中,我经常发现这种可视化能够揭示传统方法忽略的模式。例如,在一次癌症差异表达分析中,桑吉气泡图清晰显示了几个关键通路实际上共享相同的调控基因,这一发现在标准气泡图中完全被掩盖。
