企业官网建设流程全解析

单细胞测序跑到第十篇，手头终于有了那一批差异基因，但这种时候最容易被卡住：几百个基因摆在表格里，怎么告诉别人它们到底在干什么？我的答案很直接——做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 ID

clusterProfiler的enrichKEGG识别基因ID比较挑剔，默认keyType是kegg，实际上需要的是Entrez Gene ID或者KEGG ID，不是人类容易读的Symbol。所以拿到单细胞流程输出的差异基因表后，第一件事是转换ID：

library(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.05，q值控制在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富集结果表长这样：

GeneRatio和BgRatio是很多人忽略的重点。GeneRatio=15/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 ID，logFC表里也要对应使用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存PDF，base绘图用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.1，qvalueCutoff从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报错排查速查表

单细胞流程跑到KEGG这一步，真正的分水岭不是会不会跑代码，而是能不能把富集表和可视化图结合着讲故事。圈图之所以值得花时间调，就是因为它能在一张图里同时承载通路活性和基因方向两个关键信息。

我个人做完富集之后，一定会做一件事：回到差异结果里，把最显著两三条通路的基因再人工核对一遍，看看它们在其他cluster里的表达趋势是不是和富集结论一致。这个复核很费时间，但能避免被单次富集的假阳性带偏方向。另外，把这次的clusterProfiler、GOplot版本和sessionInfo存成文本留档，比圈图画好看更重要，审稿人一旦问起来，你当场就能给出可复现环境。下一期可以做GSEA和GSVA的基因集联合分析，也可以聊聊通路互作网络图，到时候接着写。