企业官网建设流程全解析

如果你在生物信息学分析中，已经完成了差异表达基因的筛选，下一步通常会做什么？很多人的第一反应是去做GO和KEGG富集分析，然后生成一张气泡图或柱状图，把富集到的通路展示出来。这没错，但问题在于，这张静态图往往成了分析的终点——我们看到了哪些通路显著，然后呢？

真正的分析价值，往往藏在“然后”之后。一张气泡图告诉你“What”（是什么通路富集了），但它很难清晰地展示“How”（这些基因是如何在不同通路间流动和协作的），以及“So What”（这对理解生物学过程意味着什么）。这时，你需要的不只是一张图，而是一套能够串联起“发现-解释-呈现”的分析叙事。

今天要讨论的，就是如何用R语言，从同一套KEGG富集分析结果出发，生成两种互补的、极具信息量的可视化图形：经典的KEGG气泡图和揭示基因流向的桑基图。这不仅仅是“一套代码出两图”的技术操作，其核心在于通过两种视角的交叉验证与叙事补充，将你的数据分析从“罗列结果”提升到“讲述故事”的层次。气泡图负责呈现显著性，是分析的“锚点”；桑基图负责揭示关联与流动，是叙事的“桥梁”。下面，我们就从为什么需要这种组合开始，一步步拆解全流程。

1. 为什么是气泡图+桑基图？超越单图呈现的分析叙事

在深入代码之前，我们必须先理解这两种图形各自承担的角色，以及它们组合起来产生的“1+1>2”效应。这决定了我们整个分析流程的设计思路，而不仅仅是机械地执行绘图命令。

1.1 气泡图：展示富集分析的“基本面”

KEGG气泡图是富集分析结果最直观的呈现方式。它的X轴通常是基因比例（GeneRatio）或富集因子（Fold Enrichment），Y轴是通路名称，点的大小代表基因数目，颜色代表显著性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, showCategory=15) + 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")

代码解读与常见问题：

1. GeneRatio的计算：GeneRatio是Count（该通路中富集到的基因数）除以BgRatio（该通路在背景基因组中的总基因数）。它比单纯的Count更能反映富集强度。
2. 排序：使用reorder(Description, p.adjust)让Y轴通路按校正P值从小到大（即显著性从高到低）排列，这是最合理的排序方式。
3. 颜色映射：用-log10(p.adjust)作为颜色映射值，使得P值越小（越显著）的点颜色越红（假设梯度是蓝到红）。
4. 图形尺寸：scale_size_continuous中的range参数控制点的大小范围，可根据通路数量调整。

至此，你已经得到了一张信息丰富、可用于发表的KEGG气泡图。它标识出了本次分析中最值得关注的Top N条通路。接下来，我们将以这些通路为焦点，构建桑基图。

3. 构建桑基图：可视化基因在通路间的共享网络

桑基图需要三类数据：节点（Nodes）和连接（Links）。在我们的场景中，节点有两层：第一层是“差异基因”，第二层是“KEGG通路”。连接表示某个基因被富集到了某个通路。

3.1 数据准备：从富集结果中提取连接关系

我们需要从kegg_enrich对象中，提取出“基因-通路”的对应关系。

# 提取富集结果中的基因-通路对应关系 # kegg_enrich@geneSets 存储了每个通路对应的基因列表（Entrez ID） # kegg_enrich@result 存储了结果表格 # 方法：遍历每条通路，获取其富集到的基因 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_enrich@result$ID == path_id) if (length(idx) > 0) { # 获取该通路富集到的基因Entrez ID genes_in_path <- kegg_enrich@geneSets[[path_id]] # 获取通路名称 path_name <- kegg_enrich@result$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 整合解读框架

不要孤立地看两张图。建议按以下流程进行整合解读：

1. 定位核心通路：在气泡图中，锁定P值最小、GeneRatio最高的几个通路（如Top 3-5）。这些是你的“一级焦点”。
2. 识别枢纽基因：在桑基图中，找到那些从“DEGs”流出、同时连接到多个“一级焦点”通路的基因。这些基因是潜在的关键调控因子。将鼠标悬停在交互式桑基图的连接上（或从静态图数据中筛选），记录下这些基因符号。
3. 验证与深挖：
   • 回到你的原始差异表达数据，查看这些枢纽基因的log2FoldChange。它们是上调还是下调？
   • 使用STRING数据库或类似工具，检查这些枢纽基因编码的蛋白质之间是否存在已知的相互作用（PPI网络）。一个紧密连接的枢纽基因簇说服力更强。
   • 查阅文献，确认这些枢纽基因在你所研究的生物学背景（如特定癌症、发育阶段）中是否已被报道具有核心作用。
4. 形成叙事：将发现串联起来。例如：“我们的分析发现，X信号通路和Y代谢通路在YY条件下被显著激活（气泡图）。进一步分析揭示，基因A和基因B是连接这两条通路的关键节点（桑基图）。已知基因A编码的蛋白能调控Y通路中的关键酶，而本研究中基因A显著上调，这可能是导致X和Y通路协同变化的核心机制。”

4.2 进阶应用与定制

基础流程跑通后，可以根据需求进行深度定制：

• 筛选连接：桑基图可能因为基因太多而显得杂乱。可以只保留那些连接到至少2个通路的基因，或者只保留表达变化最显著（如|log2FC| > 2）的基因来绘图，使图形更清晰，重点更突出。
• 分层桑基图：如果你的分析涉及多个比较组（如处理vs对照，时间点T1 vs T2），可以构建更复杂的分层桑基图，展示基因在不同条件、不同通路间的动态变化。
• 与其它数据整合：将桑基图中识别出的枢纽基因，与其拷贝数变异（CNV）、甲基化状态或生存分析数据关联，进行多组学层面的整合分析。
• 自动化脚本：将上述流程封装成一个函数，输入差异基因列表，自动输出气泡图、桑基图的数据和图形，以及一个包含枢纽基因的表格，大大提高重复分析效率。

4.3 常见问题排查

1. 桑基图节点过多，图形混乱：

   • 解决：严格限制气泡图中通路的数量（如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)

2. clusterProfiler富集分析结果为空：

   • 检查：基因ID类型是否正确？organism参数是否正确？输入的基因列表是否有效（成功转换为Entrez ID）？P值阈值是否太严格？
   • 尝试：放宽pvalueCutoff和qvalueCutoff；确保使用的注释数据库（OrgDb）与物种匹配。

3. ggalluvial绘图时出现警告或错误：

   • 检查：数据格式是否正确？aes中的axis1,axis2是否对应数据框中的列名？所有用于分组的变量（如fill）是否已转换为因子（factor）并设置了合理的水平（levels）？

4. 交互式桑基图无法显示或节点错位：

   • 检查：links数据框中的source_id和target_id是否与nodes数据框中的id正确对应（从0开始计数）。确保Value列的值是数值型。

通过这套组合拳，你的KEGG富集分析将不再是一张孤立的、陈述性的图片，而是一个包含发现、关联和假设的完整分析故事。气泡图提供了故事的目录和摘要，而桑基图则揭示了章节之间隐藏的人物关系和情节线索。掌握这种方法，意味着你掌握了将生信数据转化为生物学洞察的更强大工具。