单细胞KEGG富集分析全流程：从基因ID转换到专业圈图可视化
在单细胞转录组数据分析中，功能富集分析是揭示差异基因生物学意义的关键步骤。许多研究者习惯性地只进行GO分析，却忽略了KEGG通路富集的重要性。实际上，KEGG分析能够提供基因在代谢通路和信号转导网络中的具体作用，与GO分析的泛功能注释形成互补。本文将详细介绍如何利用R语言和Perl脚本构建完整的KEGG富集分析流程，解决实际分析中的常见痛点，并生成高质量的圈图可视化结果。
1. KEGG与GO分析的互补价值
功能富集分析是单细胞研究从数据挖掘到生物学解释的桥梁。GO（Gene Ontology）和KEGG（Kyoto Encyclopedia of Genes and Genomes）作为两大主流数据库，分别从不同角度阐释基因功能：
- GO分析：提供基因在细胞组分（cellular component）、分子功能（molecular function）和生物过程（biological process）三个层面的系统描述
- KEGG分析：聚焦基因在代谢通路、信号转导和疾病相关通路中的具体角色
两者对比关键差异如下表所示：
| 特征 | GO分析 | KEGG分析 |
|---|---|---|
| 注释维度 | 功能分类体系 | 通路网络图 |
| 结果呈现 | 层级结构 | 通路图谱 |
| 适用场景 | 广泛功能注释 | 具体通路机制研究 |
| 可视化方式 | 柱状图/气泡图 | 通路图/圈图 |
| 基因关系体现 | 功能相似性 | 通路上下游关系 |
在实际项目中，我们推荐同时进行GO和KEGG分析以获得更全面的生物学洞见。特别是在研究疾病机制或药物靶点时，KEGG通路分析能够直观展示基因在特定病理生理过程中的相互作用网络。
2. 分析前的数据准备与环境配置
完整的KEGG富集分析流程需要准备以下输入文件和环境：
必需输入文件
id.txt：包含基因ID和对应logFC值的差异基因列表（通常来自Seurat等单细胞分析流程）- 文件格式要求：
gene avg_logFC entrezID CD4 1.25 920 IL2 0.87 3558 ...
R环境配置
需要安装以下关键R包：
install.packages(c("clusterProfiler", "org.Hs.eg.db", "DOSE", "GOplot"))Perl环境准备
Perl脚本主要用于基因ID转换，确保系统已安装Perl解释器（通常Linux/macOS已预装，Windows需单独安装Strawberry Perl等发行版）
注意：如果
id.txt中缺少entrezID列，需要先使用bitr函数进行基因ID转换：
library(clusterProfiler) gene_df <- bitr(rt$gene, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")3. KEGG富集分析核心步骤
3.1 执行KEGG富集分析
使用clusterProfiler包的enrichKEGG函数进行富集分析：
library(clusterProfiler) # 读取差异基因列表 rt <- read.table("id.txt", sep="\t", header=TRUE, check.names=FALSE) rt <- rt[!is.na(rt$entrezID), ] # 去除无entrezID的基因 # 执行KEGG富集 kk <- enrichKEGG( gene = rt$entrezID, organism = "hsa", # 人类使用hsa，小鼠用mmu pvalueCutoff = 0.05, qvalueCutoff = 0.05, use_internal_data = FALSE ) # 保存富集结果 write.table(kk, file="KEGGId.txt", sep="\t", quote=FALSE, row.names=FALSE)3.2 基础可视化
生成标准的柱状图和气泡图：
# 柱状图 pdf("barplot.pdf", width=10, height=7) barplot(kk, drop=TRUE, showCategory=20, title="KEGG Pathway Enrichment", font.size=8) dev.off() # 气泡图 pdf("bubble.pdf", width=10, height=7) dotplot(kk, showCategory=20, title="KEGG Pathway Enrichment", font.size=8) dev.off()4. 基因ID转换的Perl脚本解决方案
KEGG富集结果中的基因列表通常只包含Entrez ID，缺乏直观的基因符号。以下Perl脚本可将ID转换为基因名：
#!/usr/bin/perl use strict; use warnings; # 建立基因ID到符号的映射哈希 my %id2symbol = (); open(my $id_fh, "<", "id.txt") or die "无法打开id.txt: $!"; while(my $line = <$id_fh>) { chomp $line; my @fields = split(/\t/, $line); $id2symbol{$fields[2]} = $fields[0]; # entrezID => gene symbol } close($id_fh); # 处理KEGG结果文件 open(my $kegg_fh, "<", "KEGGId.txt") or die "无法打开KEGGId.txt: $!"; open(my $out_fh, ">", "kegg.txt") or die "无法创建kegg.txt: $!"; while(my $line = <$kegg_fh>) { chomp $line; if($. == 1) { # 保留标题行 print $out_fh "$line\n"; next; } my @cols = split(/\t/, $line); my @entrez_ids = split(/\//, $cols[-2]); # 分割基因ID列表 # 转换每个Entrez ID为基因符号 my @symbols = (); foreach my $id (@entrez_ids) { if(exists $id2symbol{$id}) { push @symbols, $id2symbol{$id}; } } $cols[-2] = join("/", @symbols); # 用基因符号替换原始ID print $out_fh join("\t", @cols) . "\n"; } close($kegg_fh); close($out_fh);将上述代码保存为convert_ids.pl并运行：
perl convert_ids.pl5. 高级可视化：GOplot圈图绘制
GOplot包提供了丰富的圈图可视化功能，能直观展示基因-通路关系：
library(GOplot) # 准备输入数据 ego <- read.table("kegg.txt", header=TRUE, sep="\t", check.names=FALSE) go_data <- data.frame( Category = "KEGG", ID = ego$ID, Term = ego$Description, Genes = gsub("/", ", ", ego$geneID), adj_pval = ego$p.adjust ) id_fc <- read.table("id.txt", header=TRUE, sep="\t", check.names=FALSE) genelist <- data.frame(ID = id_fc$gene, logFC = id_fc$avg_logFC) rownames(genelist) <- genelist[,1] # 生成圈图数据 circ <- circle_dat(go_data, genelist) # 设置可视化参数 term_num <- 5 # 展示top5通路 gene_num <- 50 # 展示top50基因 # 生成弦图 chord <- chord_dat(circ, genelist[1:gene_num,], go_data$Term[1:term_num]) pdf("circ.pdf", width=11, height=10) GOChord(chord, space=0.001, gene.order="logFC", gene.space=0.25, gene.size=5, border.size=0.1, process.label=8) dev.off() # 生成聚类图 term_colors <- c("#1F78B4", "#33A02C", "#E31A1C", "#FF7F00", "#6A3D9A") pdf("cluster.pdf", width=11, height=10) GOCluster(circ, go_data$Term[1:term_num], lfc.space=0.2, lfc.width=1, term.col=term_colors, term.space=0.2, term.width=1) dev.off()6. 结果解读与优化技巧
6.1 圈图元素解读
弦图（circ.pdf）：展示基因与通路的关联关系
- 内圈：基因，按logFC值排序
- 外圈：KEGG通路
- 连接线：表示基因在对应通路中富集
- 颜色深浅：反映logFC值大小
聚类图（cluster.pdf）：展示基因表达模式与通路的关系
- 左侧：基因表达变化（logFC）
- 右侧：通路富集情况
- 颜色：代表不同通路
6.2 常见问题解决
基因ID转换失败：
- 检查
id.txt中entrezID列是否完整 - 确保Perl脚本中的列索引正确
- 检查
富集结果不显著：
- 调整pvalueCutoff和qvalueCutoff参数
- 检查差异基因分析阈值是否合适
可视化效果不佳：
# 调整弦图参数示例 GOChord(chord, space=0.005, # 增加基因间距 gene.size=6, # 增大基因标签 process.label=10) # 增大通路标签
6.3 高级定制技巧
通路筛选：在可视化前筛选特定通路
# 只展示免疫相关通路 immune_pathways <- grep("immune|infection|virus", ego$Description, ignore.case=TRUE) go_data <- go_data[immune_pathways, ]颜色自定义：使用RColorBrewer创建专业配色
library(RColorBrewer) term_colors <- brewer.pal(term_num, "Set1")
在实际项目中，KEGG圈图能够直观展示关键基因在重要通路中的分布情况，特别是在研究疾病机制或药物靶点时，这种可视化方式能清晰呈现基因-通路网络关系。一个典型的应用场景是比较不同细胞亚群间的通路活性差异，通过并排展示多个圈图，可以直观发现各亚群的特异性通路特征。
