做单细胞分析的人手里基本都攥着一长串差异基因可能是某个细胞亚群的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我的理解是这俩不是一个替代关系而是互补关系。GOGene 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富集大致有两条路线。第一是经典的ORAOver-Representation Analysis就是你拿一组差异基因或marker基因去富集看哪些通路在这个基因集合里显著富集。这是绝大多数文章的做法门槛低、结果好解释我下面讲的实操也以这个为主。第二是打分法比如用AUCell、AddModuleScore这类方法先算出每个细胞在某条通路的活性得分再做细胞亚群间的差异比较。这种方法不需要预先筛差异基因能保留单细胞的连续信息但解释起来更绕审稿人有时候会问你这个通路活性得分是怎么算的阈值怎么定的。我的建议是如果是为了快速出图、清晰表达先用ORA如果要把通路活性跟拟时序、细胞状态转化结合再考虑打分法。不要上来就整花活先把ORA跑明白。2. 富集计算的原理到底在算什么2.1 超几何分布和Fisher精确检验的通俗版本很多人第一次跑enrichKEGG看到底层用的是超几何检验瞬间头大。我换个说法你就懂了。假设全校有10000个学生其中200个是参加过数学竞赛集训的你随手抓了50个人发现里面有20个参加过集训。你会觉得这不是巧合因为你随手抓的50人里按照全校比例只能期待抓到1个200/10000×501结果实际出现了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) # 筛选显著的markeravg_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时注意设置bgwhite别留下透明背景。圈图可以先用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。场景二结果出来全是NApvalue全空。八成是基因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之后再单独加一列手动备注记录每个富集条目的信息来源或初步解读。后续写文章讨论部分这列备注就是你现成的素材。你用六个月后回头写论文时看着当时的备注比对着几百行通路名硬回忆靠谱得多。这条流程你最好按自己的数据过一遍别只看代码。富集分析是单细胞流程里最出故事的一环但也是最能暴露你是否真正理解自己数据的一环——一张图漂不漂亮倒在其次能不能跟你的生物学问题对得上才决定这篇分析有没有价值。