企业官网建设流程全解析

做单细胞分析的人，手里基本都攥着一长串差异基因，可能是某个细胞亚群的marker，也可能是疾病组对比对照组筛出来的上调下调基因。但这串基因列在表格里、火山图上，撑死了只能证明"它们确实有差异"。审稿人问一句"这些基因变化到底影响了哪些生物学过程"——回答不上来，后面就没法聊了。这时候就得靠KEGG通路富集分析把基因翻译成机制。这篇是单细胞测序流程系列的第十篇，重点讲清楚两件事：KEGG富集分析从原理到实操怎么跑通，以及富集结果怎么用圈图做出能直接放进文章的可视化。

先说清楚这篇适合谁看：正在做单细胞转录组课题的研究生、已经拿到Seurat对象但卡在"拿到marker基因之后不知道下一步干什么"的人，以及需要给文章补一张漂亮的通路富集图的从业者。默认你有R语言基础、跑过Seurat的标准流程，但KEGG富集没系统做过。读完这篇文章，你能用clusterProfiler跑出KEGG富集结果，再基于circlize或GOplot画出两种不同风格的圈图，顺便知道哪些坑是前人已经踩烂的。

1. 为什么单细胞分析绕不开KEGG富集这一步

1.1 从一长串marker基因到生物学通路的思维转变

说实话，我第一次跑完FindAllMarkers，看到输出的几百上千个基因，第一反应是兴奋，第二反应是懵。兴奋是因为终于有了自己这个细胞亚群的"身份证"，懵是因为如果要从几百个基因里逐个编故事解释细胞功能，写出来的东西发出去是要被同行笑话的。

KEGG富集分析解决的就是这个问题。它把你手上的基因列表，放到一个"通路数据库"里去比对，看这些基因是不是特别集中在某些已知的通路上。比如你筛出一批T细胞相关的marker基因，富集结果大概率会出现T细胞受体信号通路；如果你分析的是肿瘤相关巨噬细胞，结果里很可能有趋化因子信号通路、抗原加工提呈这些条目。这相当于给基因列表装了一个"语义翻译器"，把"这些基因都变了"翻译成"这些基因变化指向了哪几条核心的生物学通路"。

在单细胞场景下，这一步还有特殊价值。单细胞数据本身稀疏性高、噪声大，单个基因的表达变化往往很不稳定，但通路层面的变化是相对稳定的。一个基因在某个细胞里counts掉到0可能只是dropout，但整条通路上多个基因协同变化的信号就可靠得多。所以KEGG富集不只是为了文章好看，它在一定程度上是帮你在低信噪比的数据里捞出稳健的生物学信号。

1.2 KEGG和GO到底该看哪个

很多新手上来就问：我要做富集分析，是选GO还是选KEGG？我的理解是，这俩不是一个替代关系，而是互补关系。

GO（Gene Ontology）分三大类：生物过程（BP）、分子功能（MF）、细胞组分（CC）。它的覆盖面最广，任何基因基本都能注释到，问题是条目太多，容易富集出一大堆"泛泛"的结果，比如"信号传导""蛋白结合"这种，看起来很全但没有重点。

KEGG不一样，它收录的是经过人工整理的代谢通路、信号转导通路、疾病通路，是带有"方向感"的。KEGG的条目明显更少，但每一条都是一个可以被讲故事的整体机制图。比如KEGG里有一条通路叫"NF-kappa B signaling pathway"，你看名字就知道它讲的是什么，再去关联自己的基因，讲起机制来顺理成章。所以行业里普遍的做法是两个都做，GO看大方向、KEGG看具体机制。

单细胞文章里，最常见的图是两个：一个GO富集气泡图，一个KEGG富集圈图，或者KEGG集中用气泡图展示top通路。这篇咱们重点把KEGG讲透，因为KEGG的数据库结构、ID体系、可视化逻辑都比GO更"挑人"，踩坑概率大得多。

1.3 单细胞场景下富集的两种玩法

单细胞数据做KEGG富集，大致有两条路线。第一是经典的ORA（Over-Representation Analysis），就是你拿一组差异基因或marker基因去富集，看哪些通路在这个基因集合里显著富集。这是绝大多数文章的做法，门槛低、结果好解释，我下面讲的实操也以这个为主。

第二是打分法，比如用AUCell、AddModuleScore这类方法，先算出每个细胞在某条通路的活性得分，再做细胞亚群间的差异比较。这种方法不需要预先筛差异基因，能保留单细胞的连续信息，但解释起来更绕，审稿人有时候会问"你这个通路活性得分是怎么算的，阈值怎么定的"。

我的建议是：如果是为了快速出图、清晰表达，先用ORA；如果要把通路活性跟拟时序、细胞状态转化结合，再考虑打分法。不要上来就整花活，先把ORA跑明白。

2. 富集计算的原理，到底在算什么

2.1 超几何分布和Fisher精确检验的通俗版本

很多人第一次跑enrichKEGG，看到底层用的是超几何检验，瞬间头大。我换个说法你就懂了。假设全校有10000个学生，其中200个是"参加过数学竞赛集训"的，你随手抓了50个人，发现里面有20个参加过集训。你会觉得这不是巧合，因为你随手抓的50人里，按照全校比例只能期待抓到1个（200/10000×50=1），结果实际出现了20个，这显然"富集"了。

KEGG富集是一个道理。全校学生就是你的背景基因库（整个物种的注释基因），参加过集训的学生就是某一条KEGG通路上的基因，你抓的50个人就是你筛出来的差异基因。如果差异基因里落在某条通路上的数目，明显高于随机抓取能出现的数目，就说这条通路在你的差异基因里"显著富集"。

这个"明显高于"，统计学上就是用超几何分布或者等价的Fisher精确检验算出来的p值。p值越小，说明这个通路跟你的基因列表的关系越不可能是碰巧。

2.2 富集结果表里的每个数字代表什么

等你跑完enrichKEGG，得到一个结果表，里面有几列你必须看懂。GeneRatio，是差异基因中落在该通路的基因数占你提交的总基因数的比例，等于前面例子里"20/50"这个数。BgRatio，是注释到该通路的基因总数占背景基因库总数的比例，对应"200/10000"。pvalue上面说了，是这个富集程度的显著性。p.adjust是经BH方法校正后的p值，qvalue是Storey方法做的校正后值。Count是落在通路里的基因个数。

还有一个更直观的指标叫RichFactor，很多在线工具喜欢用，公式是差异基因在通路中的个数/该通路背景基因的总数，比如有10个差异基因落在某通路，这条通路背景里有100个基因，RichFactor就是0.1。这个值越高，说明该通路在你的差异基因里比例越大，富集强度越高。

判断显著性的红线通常取p.adjust < 0.05，但如果你做的是单细胞这种基因列表很长、通路又多的分析，我建议可以放宽到p.adjust < 0.1甚至p.adjust < 0.2，配合qvalue一起看。毕竟通路富集是一个探索性分析，过严的阈值会漏掉有线索的方向。

3. 实操：R语言做KEGG富集全流程

3.1 准备环境、数据，以及一个重要的ID转换

先看代码环境。我用的R版本是4.3.x，核心包是Bioconductor的clusterProfiler，版本4.10以上。建议用TRUE默认参数。这个包依赖org.Hs.eg.db等注释包做ID转换，还有pathview等做通路图，但KEGG富集本身主要用enrichKEGG函数。

你的输入数据来源有两种：要么从Seurat对象里拿FindAllMarkers的结果，要么自己准备一个两列的差异基因表。我这里用一个Seurat对象的实际例子。

library(clusterProfiler) library(org.Hs.eg.db) library(Seurat) library(dplyr) # 读取Seurat对象，提取marker基因 pbmc <- readRDS("pbmc_final.rds") markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.5) # 筛选显著的marker，avg_log2FC和p_val_adj是硬条件 sig_markers <- markers %>% filter(p_val_adj < 0.05, avg_log2FC > 0.5) # 提取基因名，转成ENTREZID genes <- unique(sig_markers$gene) gene_entrez <- bitr(genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)

这里有一个新手必踩的坑：KEGG的经典富集需要的是ENTREZID，不是基因symbol。原因在于KEGG数据库是以ENTREZID为索引来标记每个基因在该通路中的位置的，你用symbol去跑enrichKEGG，会丢失大量映射关系。bitr这个函数就是专门做ID转换的，转换后尽量检查一下匹配率，正常人类基因的转换覆盖率在90%以上，如果匹配率极低，先检查你的基因名格式对不对，是不是带着 Ensembl ID 或者其他奇怪后缀。

3.2 enrichKEGG的参数逐个解释

来看核心代码：

# KEGG富集 # organism参数查下面这张表，hsa=人类，mmu=小鼠，rno=大鼠 ekegg <- enrichKEGG( gene = gene_entrez$ENTREZID, organism = "hsa", keyType = "kegg", pvalueCutoff = 0.05, qvalueCutoff = 0.2, use_internal_data = FALSE )

organism是物种的三字母缩写，不是随便填的。实验室最常见的就是人和小鼠，分别对应hsa和mmu，其他物种建议去KEGG官网Genome条目里查一下。你的实验物种如果没有KEGG数据库，比如某些植物或模式生物没有收录，那KEGG富集就跑不了，只能用GO或者其他数据库替代。

keyType默认是"kegg"，这个参数告诉函数你的基因ID是什么类型。如果你的基因ID是ENTREZID，保留默认即可。如果你非要用symbol跑，keyType改为"ncbi-geneid"是不行的，你需要先把symbol转成ENTREZID再用。这是一条死路，别走。

use_internal_data这个参数值得提一下。clusterProfiler做KEGG注释时默认联网到KEGG官网下载最新数据，如果你的服务器不能上网，或者网速不稳导致每次运行都卡在下载环节，将use_internal_data = TRUE改为使用包内置数据。但注意，内置数据的注释版本可能比官网滞后一到两年，对结果影响不大，但你在方法里最好写明版本。

3.3 按细胞类型批量富集compareCluster

单细胞数据天然就是多个细胞类型并列的，如果每个亚群都手动跑一遍enrichKEGG，代码会非常臃肿。clusterProfiler提供了compareCluster，可直接按分组批量做富集分析，再配合后续可视化，一整套流程非常流畅。

# 把marker基因列和细胞类型列整理成数据框 marker_df <- sig_markers %>% select(gene, cluster) %>% # cluster就是Seurat的细胞亚群列 filter(gene %in% gene_entrez$SYMBOL) %>% left_join(gene_entrez, by = c("gene" = "SYMBOL")) # 每个cluster取一个唯一ENTREZID，避免重复 marker_df <- marker_df %>% distinct(ENTREZID, .keep_all = TRUE) # 批量富集 ckegg <- compareCluster( ENTREZID ~ cluster, data = marker_df, fun = "enrichKEGG", organism = "hsa", pvalueCutoff = 0.05, qvalueCutoff = 0.2 ) # 转成数据框查看结果 ekegg_df <- as.data.frame(ckegg) head(ekegg_df)

compareCluster输出的对象可以直接用dotplot画多组别气泡图，这也是单细胞文章里最常见的KEGG展示方式。你拿到的结果里会有cluster列，对应每个细胞亚群的富集结果。这里有个细节：如果你的某个cluster标记基因特别少，可能跑出来一条显著通路都没有，这是正常现象，不代表你数据有问题，只是这个亚群的转录特征不集中，或者你筛选marker的阈值太严了。

3.4 富集结果的保存和整理

富集结果一定记得导出CSV存底，后面做圈图、做PPT汇报、给审稿人回复都需要反复查。

# 保存完整结果 write.csv(ekegg_df, file = "KEGG_enrichment_results.csv", row.names = FALSE) # 查看显著结果有多少条 sum(ekegg_df$p.adjust < 0.05) # 单独提取某个cluster的显著通路，后面画圈图要用 cluster1_kegg <- ekegg_df %>% filter(cluster == "0", p.adjust < 0.05) %>% arrange(desc(GeneRatio))

注意看GeneRatio这一列，它是计算过的字符串比例，排序的时候需要转成数值再排，或者直接用order函数按Count排也凑合。真要精确排序，我通常把"10/250"这种拆一下。另外，富集结果的Description列是通路名，务必检查一下有没有重复条目，有的版本会把同一个通路名拆成pathway和module分开显示，实际是同一件事，画画之前顺手去重。

4. 可视化：从气泡图到圈图实战

4.1 先看最经典的气泡图dotplot

如果你时间紧，只想快速出一张能放进文章的图，直接用clusterProfiler的dotplot。

# 对enrichKEGG的结果对象直接出图 dotplot(ekegg, showCategory = 15, title = "KEGG Pathway Enrichment")

气泡图里，横轴是GeneRatio，纵轴是通路名，点的大小表示Count，颜色表示p.adjust。这张图信息密度高、读者理解成本低，是同行审稿时最买账的图之一。你在实际使用中，showCategory不要放太多，15条以内最佳，题目尽量写清楚物种、分组、比较条件。不然图一多，审稿人根本分不清哪张是哪组的。

4.2 圈图到底画什么，很多人理解偏了

说回到题目里的重头戏：圈图。首先我提醒一个概念问题。很多人以为"KEGG圈图"是这个样子：一个圆形结构，最外圈标通路，内圈连线到基因。实际上，你如果去查KEGG官网，它的通路图是长方形的，有各种框和箭头。而"圈图"这种环形展示，常见的来源是两个：一个是GOplot包画的"富集圈图"，一个是circlize包画的弦图。

GOplot的圈图通常适合展示GO富集结果，它能同时呈现Z-score和基因表达变化；而circlize的弦图适合展示通路和基因之间的隶属关系，你在文章里说"KEGG通路与基因关联圈图"，指的往往就是这种弦图。另外，GOplot画的圈图也可以用在KEGG上，但因为KEGG结果字段和GOplot默认输入格式不完全兼容，需要做一点数据变换。所以这里我把两种方案都讲给你，按自己需求选。

4.3 方案一：circlize画“通路-基因”弦图

弦图是我更推荐的一种圈图形式，特别是你要展示"某几条通路包含了哪些基因"的时候，它非常直观。外圈是通路条目，内圈是基因，弧线连接说明这个基因属于那条通路。这里的可视化逻辑类似于"这个基因被富集到了哪条通路"。实现方式如下：

library(circlize) library(stringr) # 取富集结果里显著的top通路，示例取前8条 terms <- ekegg_df %>% filter(p.adjust < 0.05) %>% head(8) # 把geneID这一列拆开，geneID一般是ENTREZID用/分隔的字符串 # 可以先转成symbol，方便看图 idx_list <- str_split(terms$geneID, "/", simplify = FALSE) gene_syms <- lapply(idx_list, function(x) { mapIds(org.Hs.eg.db, keys = x, column = "SYMBOL", keytype = "ENTREZID") }) # 生成两列数据：第一列通路名，第二列基因list里的符号 chord_df <- lapply(seq_along(terms$Description), function(i) { data.frame( pathway = rep(terms$Description[i], length(gene_syms[[i]])), gene = as.character(gene_syms[[i]]), stringsAsFactors = FALSE ) }) %>% bind_rows() # 去掉NA基因 chord_df <- chord_df %>% filter(!is.na(gene)) # 画弦图 pdf("KEGG_chord.pdf", width = 10, height = 10) circos.clear() circos.par(start.degree = 90, gap.degree = 3) chordDiagram( chord_df, transparency = 0.3, directional = -1, direction.type = c("diffHeight", "arrows"), link.arr.type = "big.arrow", annotationTrack = c("grid", "name"), big.gap = 8 ) dev.off()

这段代码核心就是把富集结果里的geneID按通路拆开，整理成一个两列的长表，再输入给chordDiagram。这里注意几个细节：第一，通路名不要太长，太长的话外圈标签会重叠，我一般精简通路名，比如"NF-kappa B signaling pathway"改成"NF-κB"；第二，如果基因太多，弦图会变成一团乱麻，建议控制每条通路只展示显著性最高的前10-15个基因；第三，如果只想突出一个细胞亚群的KEGG结果，建议单独取该亚群的数据做，不要所有亚群混在一张图里，颜色根本没法区分。

4.4 方案二：GOplot风格富集圈图

另一种圈图是GOplot包里的GOCircle，它长这样：中心是一个圆形，最外圈一圈关键词，中间扇形区域用颜色表示Z-score，内圈散点表示基因的表达倍数变化。这个图的好处是能同时展示富集显著性、基因上下调方向和富集分数，文章里用它来做主图非常有冲击力。

GOplot原本是给GO富集设计的，但如果你的KEGG结果也包含geneID、logFC这些字段，完全可以改造。GOCircle的输入要求有一个"zscore"列，需要根据通路里上调下调基因数目算一个分数。基本思路是：针对每条通路，把落在通路内的基因，按照avg_log2FC是否大于0分成上调组和下调组，然后算Z-score = (up - down) / sqrt(count)。这个分数在GOplot的文档里有明确说明。

# 构建GOplot需要的输入格式 # 准备基因的logFC信息 gene_logfc <- sig_markers %>% select(gene, avg_log2FC) gene_logfc$gene <- mapIds(org.Hs.eg.db, keys = gene_logfc$gene, column = "SYMBOL", keytype = "ENTREZID") # 对每个通路计算up/down和zscore terms_z <- terms %>% rowwise() %>% mutate( genes = list(str_split(geneID, "/")[[1]]), up = sum(gene_logfc$avg_log2FC[gene_logfc$gene %in% unlist(genes)] > 0), down = sum(gene_logfc$avg_log2FC[gene_logfc$gene %in% unlist(genes)] < 0), zscore = (up - down) / sqrt(Count) ) %>% ungroup() # 用GOplot时，还需要一个包含基因ID和logFC的独立对象，之后调用 # GOCircle(circ) 即可

这个方案稍微绕一点，但画出来的图确实比单纯的弦图信息量大，因为它是按"通路富集方向"来组织的，既展示了通路之间的差异，也展示了基因在通路中的角色。不过我要提醒你，GOplot的图对数据量很敏感，通路数量最好控制在10条以内，基因数量50-100个之间，太多了中心区域会完全糊掉。

4.5 配色、导出，以及投稿前的最后检查

圈图的配色其实很影响审稿印象。circlize默认的高饱和度用色比较"原始"，我建议自己调色。简单做法：指定chordDiagram的col参数，或者设置circos.par的track.height。如果嫌麻烦，用RColorBrewer的Set2或者Paired配色，比默认颜色耐看得多。

导出格式上，文章排版大图优先用PDF矢量图，方便后期编辑。分辨率要求：如果期刊要求TIFF 300dpi，用ggsave导出PNG或者TIFF时注意设置bg="white"，别留下透明背景。圈图可以先用PDF导出，然后用Adobe Illustrator或Inkscape转存成TIFF，这样最保险。

还有一个细节，做富集图时，基因名尽量保证能对应到KEGG结果里的ENTREZID，否则你圈图里画了一大堆基因，回头审稿人问"你这个基因在通路里对应的是哪个节点"，你答不上来又得重新查。

5. 常见问题和避坑实录

5.1 enrichKEGG跑不出来，先排查这几种情况

场景一：控制台报错"Failed to download KEGG data"。这是最常见的联网失败，尤其是服务器在防火墙后面。解决办法把use_internal_data改成TRUE先用旧版数据跑通，或者手动下载KEGG的本地注释文件，再通过clusterProfiler的read.gmt思路读入，但这套操作比较折腾，不建议新手折腾，我一般直接建议用use_internal_data。

场景二：结果出来全是NA，pvalue全空。八成是基因ID类型搞错了。你传进去的gene参数需要是ENTREZID的字符向量，如果你直接塞了一堆symbol进去，很多基因在KEGG里匹配不上，就没法做统计。检查办法是跑完enrichKEGG后打印一下对象，看里面geneID列是否为NULL，如果是，说明输入基因无法映射到KEGG上。

场景三：p.adjust全大于0.05，一条显著通路都没有。这种情况在单细胞数据里不少见，原因可能是你的差异基因个数太少（少于50个就很难富集出显著通路），或者你用了all markers做富集而不是差异最显著的基因。建议检查一下过滤条件，avg_log2FC阈值放宽到0.25试试，或者把多个cluster的基因合并成大集合再做。

5.2 解读KEGG结果时的三个提醒

第一个提醒是"富集到某条通路，不等于你的细胞'激活'了这条通路"。ORA只是一种统计学上的过表达分析，说明这些基因在通路里的比例异常高，不代表通路活性真的上升了。真要证明通路激活，你需要回到表达量层面看这些基因到底是上调还是下调。这也是我建议你用GOplot之类的工具展示logFC的原因，它至少让人能看到方向。

第二个提醒是通路名自带"光环"问题。比如富集到"Pathways in cancer"，这个条目在KEGG里是一个整合型通路，包含了很多跟肿瘤相关的子通路基因，如果你做的不是肿瘤课题，富集到它往往只是因为这组基因覆盖面太广，生物学特异性弱。汇报的时候如果硬拿这条去讲机制，容易被懂行的专家挑刺。这种条目可以保留显示，但解释时点到为止。

第三个提醒是不要只挑P值最小的几个通路讲。富集分析本质是筛选假设，最显著的通路不一定跟你的生物学问题最相关。我见过有人富集出一堆"核糖体"通路，这是因为样本处理过程中细胞应激，核糖体蛋白基因大规模变化导致的，跟研究目标八竿子打不着。这时候要认真审视数据质量，而不是硬编故事。

5.3 送给新手的实操节奏建议

如果从零开始做单细胞的KEGG富集，我的节奏建议是：第一步先把常规的差异基因筛选跑通，明白哪些基因留哪些基因去，这一步比富集本身更影响结果。第二步是选几个细胞亚群练手，跑一次性价比最高的ORA，看看结果是否符合预期。第三步再追求可视化，从dotplot做起，等dotplot已经能讲清楚故事了，再挑战圈图。不要一上来就直奔圈图，因为圈图对输入数据格式的要求更高，你还没摸清富集结果里各个字段的关系就去做图，大概率会卡在莫名其妙的地方。

我在实际项目中，KEGG富集部分真正花时间最多的不是跑程序，而是"反复拷问结果"：这条通路是不是合理、这个基因簇是不是有批次效应、这个亚群的marker基因是不是真的亚群特异的。跑代码十分钟，讨论结果一整天，这是单细胞数据分析的常态。

最后再分享一个小技巧：富集分析的结果表，我习惯在导出CSV之后，再单独加一列手动备注，记录每个富集条目的"信息来源"或"初步解读"。后续写文章讨论部分，这列备注就是你现成的素材。你用六个月后回头写论文时，看着当时的备注，比对着几百行通路名硬回忆靠谱得多。

这条流程你最好按自己的数据过一遍，别只看代码。富集分析是单细胞流程里最"出故事"的一环，但也是最能暴露你是否真正理解自己数据的一环——一张图漂不漂亮倒在其次，能不能跟你的生物学问题对得上，才决定这篇分析有没有价值。