单细胞测序跑到第十篇手头终于有了那一批差异基因但这种时候最容易被卡住几百个基因摆在表格里怎么告诉别人它们到底在干什么我的答案很直接——做KEGG通路富集分析然后把结果画成圈图。KEGG是单细胞测序分析里绕不开的一步。无论是细胞类型注释后找marker基因还是拟时序分析之后找关键调控模块最终都要回到“这些基因落在哪些通路上”这个问题。这篇就把KEGG通路富集分析和可视化圈图这条线完整走一遍从概念讲到实操从clusterProfiler跑富集到GOplot画圈图最后附上我踩过的坑和排查思路。适合已经会跑单细胞基础分析、想补上通路解释这一环的初学者也适合想把手里的富集表画得更有说服力的老手。1. KEGG富集分析开始前先把这几个概念捋清楚1.1 KEGG数据库不是“基因列表”是功能模块地图很多人第一次接触KEGG时会把它理解成一种类似GO注释的基因功能标签这个理解不准确。KEGG全称是Kyoto Encyclopedia of Genes and Genomes核心组织方式不是给每个基因打标签而是把基因、酶、化合物、反应组织成一张又一张通路图。整库分几个层级通路层级Pathway比如hsa04110是细胞周期、KO层级KEGG Orthology比如K06630是CDK1、模块和反应层级。kegg注释跑得顺不顺利取决于你是否理解这个层级关系。我们自己测序得到的是一个个基因名但KEGG做富集时先把基因映射到KO条目上再通过KO对应到具体通路里去数个数。所以你会看到富集结果表里第一列永远是hsa04110这种带物种前缀的通路ID而不是一排排基因名堆在前面。单细胞流程里这一步尤其容易出问题因为上游降维聚类之后拿到的基因名经常是Symbol格式如果没有正确转换KEGG注释率会低得让人怀疑人生。1.2 KEGG通路里的一个成员是一个基因吗很多人问过这个这个热搜词对应的疑问我在带学生的时候被问过好多次“通路图里那个方框是不是代表一个基因”严格说不是。KEGG通路图上方框标的是KO条目或EC编号不是具体某个基因名。比如细胞周期通路的图里标着K06630的那个方框对应的是CDK1的蛋白产物不同物种中执行同样功能的同源基因都会归到这一个KO下。理解这一点后面的富集逻辑就通了。所谓“差异基因在通路中”本质是“差异基因的蛋白产物归属的KO条目出现在这条通路的节点集合里”。这就像坐地铁KO是站台各物种的基因是乘客乘客虽然来自四面八方但进了同一个站台就能坐同一条线路。所以画圈图时一个基因可能出现在多条通路里一条通路也可以收容多个差异基因这种多对多关系恰恰是KEGG富集分析和圈图可视化最想表达的信息。1.3 ORA和GSEA两种思路单细胞场景下怎么选KEGG富集分析有两条技术路线实践中经常被混为一谈。一条是过表达分析ORA核心是对差异基因列表做超几何检验看目标通路是否比随机更富集另一条是基因集富集分析GSEA不要求先筛出差异基因而是把全部基因按表达变化排序再看每条通路在排序顶端的富集程度。我的建议是如果单细胞数据里差异基因数量适中比如一个cluster的marker基因有100到1000个用ORA效率高、结果直观如果差异基因只有二三十个或者你担心硬阈值把弱信号滤掉了用GSEA补充一次会稳妥很多。单细胞项目里我通常两条线都跑一下ORA结果用来出主图GSEA结果用来验证那些看起来显著但基因数很少的通路。2. 实操第一步用clusterProfiler把富集结果跑出来2.1 环境准备和包安装R里做KEGG富集的主流选择是clusterProfiler配套需要OrgDb注释包和后续画图的GOplot。安装代码很简单但要注意版本兼容if (!require(BiocManager, quietly TRUE)) { install.packages(BiocManager) } BiocManager::install(c(clusterProfiler, org.Hs.eg.db, GOplot)) install.packages(circlize)这里的org.Hs.eg.db是人类基因组注释包其他物种要换对应的比如小鼠用org.Mm.eg.db大鼠用org.Rn.eg.db。Windows环境下建议R版本不低于4.2否则部分依赖包的编译会遇到麻烦。2.2 差异基因ID转换Symbol到Entrez IDclusterProfiler的enrichKEGG识别基因ID比较挑剔默认keyType是kegg实际上需要的是Entrez Gene ID或者KEGG ID不是人类容易读的Symbol。所以拿到单细胞流程输出的差异基因表后第一件事是转换IDlibrary(clusterProfiler) library(org.Hs.eg.db) # 假设deg_symbols是你从单细胞差异分析中拿到的基因Symbol列表 deg_symbols - c(TP53, CDK1, CCNB1, EGFR, MYC) # Symbol转Entrez ID entrez - bitr(deg_symbols, fromType SYMBOL, toType ENTREZID, OrgDb org.Hs.eg.db) # 去重很重要同一个Symbol可能映射到多个Entrez ID entrez - unique(entrez$ENTREZID) length(entrez) # 看一眼成功转换了多少个bitr是clusterProfiler里非常好用的ID转换函数。转换之后一定检查转换率如果低于70%先回头查一查差异基因表里的基因名格式是不是有问题比如带上了版本号或者Ensembl ID后缀。2.3 enrichKEGG核心参数逐个说ID转换完成后跑ORA版本的KEGG富集就是一行函数的事ekegg - enrichKEGG(gene entrez, organism hsa, keyType kegg, pvalueCutoff 0.05, qvalueCutoff 0.2)organism填的是NCBI的物种三字母缩写人类是hsa小鼠是mmu大鼠是rno斑马鱼是dre猪是ssc鸡是gga。很多初学者在这里填成“human”直接报错。pvalueCutoff控制的是富集检验的原始p值阈值qvalueCutoff控制的是多重假设检验校正后的FDR阈值。一般p值放宽到0.05q值控制在0.2单细胞数据里如果差异基因不多可以再把p值放到0.1看看结果轮廓。需要额外说明的是背景基因集。enrichKEGG默认用的是该物种全部KEGG注释基因做背景不是你自己传入的那个差异基因列表。这意味着你的差异基因再怎么全面富集结果也只能反映“这群基因在整个物种注释背景下的相对富集程度”。有些人会误把传入的gene列表当成背景导致结果看起来奇怪。2.4 GSEA版本怎么写GSEA版KEGG富集不需要先定义差异基因阈值直接把所有基因的表达变化排序传进去# res是差异分析结果表需要包含log2FoldChange和entrez ID geneList - res$log2FoldChange names(geneList) - res$entrez # 去掉缺失值按表达变化从大到小排序 geneList - geneList[!is.na(names(geneList))] geneList - sort(geneList, decreasing TRUE) # 跑GSEA gkegg - gseKEGG(geneList geneList, organism hsa, pvalueCutoff 0.05)这里有个高频报错点gseKEGG要求geneList必须是已经排序的命名数值向量如果名字里有NA或者没有按降序排好函数会直接报错。单细胞流程里建议把上游差异分析输出的所有基因都拿来做排序不要只选显著性基因否则GSEA会失去随机排序的背景意义。2.5 结果表字段看不明白怎么办富集结果跑完很多人拿到结果表只盯着p值看忽视了其他字段。一个标准的KEGG富集结果表长这样字段含义使用建议IDKEGG通路ID如hsa04110后续画图和查询全凭它Description通路名称如Cell cycle写论文时用这个GeneRatio差异基因中命中该通路的比例如15/200用来判断富集强度BgRatio背景基因中命中该通路的比例如120/8000与GeneRatio对比才有意义pvalue超几何检验原始p值越小越显著p.adjustBH校正后p值推荐看这个qvalue基于FDR的q值某些审稿人会问geneID命中的Entrez ID斜杠分隔画圈图时需要拆开用Count命中该通路的差异基因个数注意不是通路基因总数GeneRatio和BgRatio是很多人忽略的重点。GeneRatio15/200表示200个差异基因里有15个落在细胞周期通路里如果BgRatio是120/8000说明背景里只有1.5%的注释基因属于该通路而你的差异基因里有7.5%落在这里这个富集就是有说服力的。看表的时候把两个比例放一起读比单独盯p值有用得多。3. 把结果变圈图GOplot从数据整理到成品图3.1 圈图凭什么比气泡图更能打KEGG富集分析最常见的可视化成套动作是气泡图横轴是GeneRatio纵轴是通路名点大小是基因数颜色是p值。气泡图信息量没问题但它在单细胞场景里有一个天然短板——看不出基因表达方向也没法呈现基因和通路之间的多对多关系。圈图正是冲这个需求来的。它把通路级别和基因级别的信息叠在同一张图里外圈告诉你哪些通路显著内圈告诉你每个通路里的基因在上调还是下调基因点的大小还能反映富集到的基因数量。审稿人看到圈图的第一反应往往是“这张图信息量很足但又不乱”这适合作为文章主图。我用它展示cluster特异的marker基因通路时曾在一张图里同时讲清楚了“细胞周期通路显著富集”和“该通路以CDC20、CCNB1等上调基因为主”两个结论。3.2 数据重排从enrichKEGG结果到GOplot能认的格式GOplot包画圈图时需要一个特定结构的数据框通常是从DAVID导出的格式每一行是一个基因在某条通路中的记录必须包含通路ID、通路名称、基因名和该基因的logFC值。用enrichKEGG结果转换的代码如下library(GOplot) # 把富集结果转成数据框 kegg_df - as.data.frame(ekegg) # 逐条通路拆开geneID列生成“通路-基因”长表 gene_pairs - lapply(seq_len(nrow(kegg_df)), function(i) { genes - strsplit(kegg_df$geneID[i], /)[[1]] data.frame( ID kegg_df$ID[i], Term kegg_df$Description[i], Genes genes, logFC deg_logFC[genes] # 从差异分析里按Entrez ID匹配logFC ) }) david - do.call(rbind, gene_pairs) # 转成GOplot内部格式 circ - circle_dat(david, term Term)这里最容易翻车的坑是logFC匹配。基因名用Entrez IDlogFC表里也要对应使用Entrez ID如果你差异分析表里的行名是Symbol需要先转换再合并。另一个坑是重复行同一个基因在同一条通路里只保留一行不要在拆geneID时把重复的拼进去。3.3 GOCircle参数调优数据准备妥当画图本体反而简单GOCircle(circ, nsub 12, lfc.col c(cornflowerblue, firebrick), label.size 4, rad1 0.5, rad2 1.5, rad3 2.0)nsub控制展示多少条通路我建议控制在8到15之间超过15条会变成一圈密密麻麻的色带内圈基因点也会挤到分不清。lfc.col是内圈基因点的颜色向量默认是红配绿但我的实际体验是红配蓝在彩色打印和色盲友好度上都更好所以通常手动改为蓝红组合。rad1/rad2/rad3控制三个圆环的半径位置文字重叠时把rad3调大就能把标签往外推。还有几个参数值得调。table.legend默认会附带一个图例表格投稿时建议关掉改在正文里描述图例避免挤占主图空间。内圈基因点的纵轴范围可以通过lfc.min和lfc.max控制如果某一通路的基因logFC跨度特别大不设置上下限会把同一条通路里的点拉变形。3.4 配色与排序的审美细节圈图的观感很大程度上取决于细节。通路排序上我通常按调整后p值从小到大排列让最显著的通路排在最上方这样第一眼就能抓到重点。配色方面外圈通路色带可以考虑用RColorBrewer的Set1或Dark2内圈基因点颜色要和通路色带区分开避免撞色。基因点自身的大小默认代表该通路富集到的基因数也就是Count列。这个设计很实用但要注意一个陷阱如果两条通路的Count相差十几倍小点的变化会被压缩到几乎不可见这时候可以考虑对Count取对数再映射到点大小信息损失不大但可视差异明显得多。4. 进阶玩法GOChord和弦图和基于circlize的自定义圈图4.1 GOChord展示基因与通路的多对多关系如果你觉得GOCircle的信息密度还不够或者想换个角度强调“基因同时参与多个通路”这个事实GOChord是更好的选择。它呈现的是和弦图左侧排列基因右侧排列通路连线的有无表示归属关系基因块的颜色按logFC方向渐变通路块颜色区分通路。GOChord的输入是0/1矩阵加一列logFC# 构建基因与通路的0/1矩阵 # 行名是基因列名是通路值为1表示该基因属于该通路 # 最后一列必须叫logFC存基因的表达变化 chord - matrix(0, nrow length(genes), ncol length(pathways)) rownames(chord) - genes colnames(chord) - pathways for (i in seq_len(nrow(david))) { chord[david$Genes[i], david$Term[i]] - 1 } chord - as.data.frame(chord) chord$logFC - deg_logFC[rownames(chord)] GOChord(chord, space 0.02, gene.order logFC, lfc.col c(blue, white, red))GOChord对输入基因数量很敏感一次塞100个基因大概率会画成一团乱麻。实际操作时应该手动筛选只保留最核心的20到40个基因和5到10条通路。筛选逻辑可以考虑出现频次高于阈值、logFC绝对值大、在重点通路中反复出现的基因优先保留。4.2 用circlize自由定制圈图GOplot的圈图是封装好的想大改布局反而不容易。当你需要画一张完全按自己逻辑组织的圈图比如外层是通路分类、中层是富集显著性、内层是基因表达、最里层是平均表达趋势时我建议直接用circlize从零搭。大致思路是用chordDiagram画基因到通路的连接带再用circos.track叠加logFC柱状条和显著性色带。library(circlize) # 彩色富集结果 # 这一步通常需要把通路按显著性排好序 # 然后用circos.par设定起始角度 circos.clear() circos.par(start.degree 90, gap.degree 2) chordDiagram(gene_pathway_df, transparency 0.5)circlize胜在灵活代价是学习成本高。对多数期刊级图表来说GOplot已经够用了circlize更多用于汇报展示或需要批量定制配色方案的场景。我一般建议先把GOplot吃透再考虑上circlize。4.3 输出格式与尺寸圈图画完导出是个不能省的环节。发文章建议输出成PDF或SVG矢量格式这样文字和点都不会糊。如果期刊非要位图用tiff输出300dpi以上宽度按期刊要求设置一般是单栏8.5cm左右、双栏17.5cm左右。我用过的稳妥方式是ggplot对象用ggsave存PDFbase绘图用pdf()函数保存前先用dev.size确认画布尺寸。单细胞项目里如果需要同时展示多个cluster或分组的富集圈图我通常把每张图分别导出再用AI或PPT排版时统一字体字号而不是在一张R图里硬塞多个panel那样字体会变小变挤反而不好看。5. 常见问题与排查实录5.1 enrichKEGG总提示下载失败怎么办这个问题出现频率极高。clusterProfiler的enrichKEGG默认会去KEGG官方API拉取最新通路数据官方接口不稳定或者本机网络受限时就会在运行中途卡住或者提示download failed。这种情况我会放弃实时联网改用本地KEGG注释缓存来跑。做法是先把KEGG数据下载到本地library(clusterProfiler) # 将KEGG数据保存到本地后续设置use_internal_data TRUE即可 downloadKEGG(species hsa)之后在enrichKEGG里加一行use_internal_data TRUE就能直接读取本地注释不依赖远程接口。我自己的经验是跑之前先试一次在线版失败就立刻切本地缓存不纠结省时间要紧。5.2 富集结果全是空集或显著通路少得可怜单细胞数据分析中最常见的挫败感来源就是跑了半天得到一个空结果表。排查顺序记住一条链差异基因数量太少、ID转换率太低、物种代码错了、阈值太严。差异基因少于30个时ORA基本没有统计功效建议改上GSEA。如果是转换率问题检查Symbol大小写和是否有基因名里带有“-AS1”这类容易转换失败的lncRNA命名。阈值问题最好解决把pvalueCutoff从0.05放到0.1qvalueCutoff从0.2放到0.3先看看轮廓再收紧。5.3 GOCircle报错和文字重叠问题GOCircle最常见的报错是“Error in circle_dat”九成原因是输入数据框列名不规范。GOplot对列名有严格约定必须包含ID、Term、Genes、logFC四列其中logFC列名不能改。另一个高频问题是图上文字重叠通路标签和基因点挤作一团。解决路径按优先级排序调大rad3把标签外推调小label.size最后减少nsub。5.4 多个cluster的通路结果怎么比较单细胞项目里很少有人只做一个cluster的KEGG。当你有四五个cluster各自跑出富集结果时不建议把四五张圈图简单拼在一起。我更推荐的做法是每个cluster分别跑富集然后取几个核心通路的显著性做热图横轴是通路、纵轴是cluster颜色是-log10(p.adjust)这样一眼就能看出来哪个cluster处于增殖状态、哪个cluster在走炎症通路。圈图留给最核心的cluster做深描不要平均用力。5.5 gseKEGG报错排查速查表报错关键词原因处理办法geneList is not sorted传入向量没有按数值降序排列用sort(decreasing TRUE)NA in geneList基因名或数值含NA过滤后再传入Error in download.KEGG.Path远程接口不通安装KEGG.db或使用本地缓存No gene can be mappedID格式不正确确认使用Entrez ID重新走bitr转换organism not found物种缩写错误用NCBI三字母代码单细胞流程跑到KEGG这一步真正的分水岭不是会不会跑代码而是能不能把富集表和可视化图结合着讲故事。圈图之所以值得花时间调就是因为它能在一张图里同时承载通路活性和基因方向两个关键信息。我个人做完富集之后一定会做一件事回到差异结果里把最显著两三条通路的基因再人工核对一遍看看它们在其他cluster里的表达趋势是不是和富集结论一致。这个复核很费时间但能避免被单次富集的假阳性带偏方向。另外把这次的clusterProfiler、GOplot版本和sessionInfo存成文本留档比圈图画好看更重要审稿人一旦问起来你当场就能给出可复现环境。下一期可以做GSEA和GSVA的基因集联合分析也可以聊聊通路互作网络图到时候接着写。