企业官网建设流程全解析

做转录组项目的人都知道，拿到差异基因列表的那一刻，既兴奋又心虚。兴奋的是，前面经历了质控、比对、定量、差异分析层层筛选，终于在几百上千个基因里看到了可能跟表型变化相关的候选；心虚的是，如果接下来这一步跑不好，之前积累的所有可信度都会打折扣——这一步就是GO/KEGG富集分析。说它是转录组学RNA-Seq测序数据生信分析里最被高频使用、同时也最容易被误用的环节，一点都不夸张。网上铺天盖地的教程，基本就是"把差异基因粘贴进去、点一下提交、下载个表格"，但真正的难点在于：你的输入数据准备好了吗？背景基因选对了吗？为什么KEGG富集经常只有寥寥几条通路？为什么同一个项目用不同参数跑出来的结论完全相反？

这篇文章是转录组学RNA-Seq测序数据生信分析系列的第5篇，我会围绕GO和KEGG富集分析，把从输入数据规范化、背景基因集选择、R包实操、结果解读到可视化出图的完整链路拆开讲，重点放在那些常规教程里不会写、但你实际跑项目时百分之百会踩的问题上。适合已经把差异基因列表拿到手、准备开始做注释和富集分析的同学，也适合已经跑过一次但被结果弄得一头雾水的人。

1. 富集分析的第一道门槛：输入数据格式与背景基因的选择

很多人以为富集分析就是把一个Excel里的基因名列表导进去就行，这个想法会带来两个很隐蔽的问题：第一，基因ID不统一，有的写的是Symbol（如TP53），有的是ENTREZID（如7157），有的甚至是Ensembl ID，工具根本认不全；第二，没有指定背景基因集，部分在线工具默认用全数据库或某物种的所有注释基因当背景，最后算出来的富集条目，经不起推敲。

先解决ID的问题。富集分析本质上是一个"在给定注释库中做映射"的过程，无论GO还是KEGG，后台都是把基因ID关联到功能条目或通路上。clusterProfiler等R包通常支持多种ID类型，但为了保证效率和减少歧义，我习惯统一转成ENTREZID再往下走。转换方法很简单：

library(clusterProfiler) library(org.Hs.eg.db) # 假设deg_symbol是差异基因的Symbol列表 entrez_ids <- bitr(deg_symbol, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)

这里有个容易忽视的细节：bitr()转出来如果丢失了一部分基因，不要急着往下跑，先看丢失比例。丢失的原因通常有两个。一个是Symbol本身就是历史别名，比如某些基因在最新注释版本里已经换了官方名；另一个是格式问题，比如大小写、空格、行尾回车符。处理方式是先unique()去重，再用toupper()统一大小写核对一遍，实在查不到的就放弃，或者等做KEGG的时候用旧的别名桥接。丢失超过20%时，我会觉得差异分析那一步的注释版本可能就有问题，会回去检查基因注释文件。

再说背景基因集，这个是重灾区。背景基因（universe）指的是你"本实验里有机会被检测到的所有基因"，不是当前物种所有注释基因，更不是全数据库基因。为什么？举个例子，如果你的测序方案是表达谱芯片或某个靶向RNA捕获panel，那很多基因压根没有探针或捕获探针，它们在筛选阶段就不会出现在表达矩阵里，也就不可能成为差异基因。用全基因组注释当背景，等价于把大量"不可能被检出"的基因塞进分母，这会虚高富集显著性。正确做法是：用你定量矩阵中所有检测到的基因id作为背景基因集。实操中，我会把DESeq2结果表里所有基因（无论是否差异）的ENTREZID列提取出来当背景，而差异基因列表则用筛选阈值过滤后的子集。

# 结果表res_filtered按|log2FC|>1且padj<0.05筛选 deg_entrez <- bitr(rownames(res_filtered[res_filtered$padj < 0.05 & abs(res_filtered$log2FoldChange) > 1, ]), fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db) # 背景：所有检测到的基因 background_entrez <- bitr(rownames(res_all), fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)

提示：如果用的是edgeR或limma，思路完全一样，核心原则就一条——差异基因必须是你背景基因的子集。否则富集结果没有任何生物学意义，审稿人一眼就能看穿。

2. GO富集：BP/MF/CC三分法、OrgDb注释体系与自建注释库的玩法

2.1 GO三个ontology到底在讲什么

Gene Ontology（GO）由三个独立的本体组成：生物学过程（Biological Process, BP）、分子功能（Molecular Function, MF）、细胞组分（Cellular Component, CC）。BP描述的是"基因参与的一系列事件"，例如"炎症反应""DNA修复"；MF描述的是"基因产物在分子层面执行的任务"，例如"ATP结合""蛋白激酶活性"；CC描述的是"基因产物在细胞里待的位置"，例如"线粒体膜""核内体"。

做富集时，我倾向于把BP、MF、CC分开跑。原因很简单，三者粒度完全不同。BP条目往往层级深、词条多，一次能出来上百条；MF通常比较稳定但数量少；CC结果容易跟"细胞类型标志物"混淆，需要结合实验背景判断。如果三个一起去跑，结果表会非常臃肿，而且不同ontology之间的条目完全无法横向比较。分开跑，后续绘图时也更容易聚焦生物学问题。

2.2 enrichGO的完整参数解读

用clusterProfiler跑GO富集，核心代码是这样的：

ego_bp <- enrichGO(gene = deg_entrez$ENTREZID, universe = background_entrez$ENTREZID, OrgDb = org.Hs.eg.db, keyType = "ENTREZID", ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.2, readable = TRUE)

这里面的参数没有一个可以随便填，逐一说明：

• gene：差异基因ID向量，必须是字符型，且来自universe子集。
• universe：背景基因。不写，R包会默认使用OrgDb里的全部基因，这就是我前面强调的坑。
• OrgDb：物种对应的注释包，人类是org.Hs.eg.db，小鼠是org.Mm.eg.db，大鼠是org.Rn.eg.db。没有现成物种包时需要用AnnotationHub或自建，后面展开讲。
• keyType：输入基因的ID类型，bitr()转成ENTREZID后，这里就是ENTREZID。
• pAdjustMethod：多重假设检验校正方法，一般选BH（Benjamini-Hochberg）。为什么要校正？因为你要同时检验成百上千个GO条目，不做多重检验校正，假阳性会爆炸。
• pvalueCutoff和qvalueCutoff：这两个不是一回事。pvalue是校正前的原始P值过滤，qvalue是校正后的FDR过滤。实际操作中我会保留pvalueCutoff=0.05，qvalueCutoff=0.2，主要看qvalue，它更严格。
• readable=TRUE：把ENTREZID映射回Symbol，方便后续人工查阅基因名称。但如果后面还要做cnetplot或emapplot，我一般先设FALSE，画完图再映射，避免某些工具对Symbol的兼容性问题。

跑完之后，ego_bp实际上是一个enrichResult对象，转换成data.frame就能看了：

ego_df <- as.data.frame(ego_bp) head(ego_df[, c("ID", "Description", "GeneRatio", "BgRatio", "pvalue", "p.adjust", "qvalue")])

看见GeneRatio和BgRatio没？GeneRatio是"落在该GO条目的差异基因数/总差异基因数"，BgRatio是"落在该GO条目的背景基因数/总背景基因数"。富集因子就是这两个比值的比，具体表现为条目是否"过度代表"。解释结果时，GeneRatio高不代表一定有生物学意义，关键要看p.adjust和qvalue。

2.3 OrgDb缺失时自建GO注释库

非模式物种经常会遇到一个尴尬情况：org.xx.eg.db不存在。解决办法有几个，按优先级排序：

第一，去AnnotationHub搜：

library(AnnotationHub) ah <- AnnotationHub() # 按物种名搜索 query(ah, c("GFF", "Aedes aegypti"))

AnnotationHub上的注释资源一般会转录成OrgDb对象，但质量参差不齐，需要自己检查下版本日期和源数据库。

第二，如果你有该物种GFF注释文件，并且上面有GO注释信息，可以选择自建OrgDb。这里说的自建，在clusterProfiler里对应的函数是makeOrgPackage，它接受三个表：基因ID和GO term的对应关系表、基因ID和基因名的对应表、以及基因本身的描述信息表。最麻烦的是第一张表，本质上就是把GFF第三列上的Ontology_term属性拆出来，转成长表格式。我之前做过一个无参考物种的转录组项目，用的就是NCBI上导出的基因-GO关联表，构建流程大约长这样：

suppressMessages({ library(AnnotationForge) }) # gene2go：两列，gene_id, go_id（可以有多行，相当于一个基因多个GO） # gene_info：gene_id, gene_symbol, description # gene2pubmed不是必须的，可以给空表 makeOrgPackage(gene_info = gene_info, go = gene2go, version = "1.0.0", maintainer = "your name <your@email.com>", author = "your name", outputDir = "orgdb_pkg", tax_id = "XXXXX", genus = "YourSpecies", species = "YourSpecies_name")

生成的是一个标准R包目录，用install.packages("orgdb_pkg", repos = NULL, type = "source")装进去后，就能像org.Hs.eg.db一样被enrichGO()读取了。这条路第一次走比较折腾，但一劳永逸，后续所有基因注释、ID转换、富集分析都会顺畅很多。

2.4 GO结果里的冗余条目处理

GO条目本身是DAG结构（有向无环图），父子条目之间有包含关系。直接跑出来的结果里经常出现"父条目显著、子条目也显著"的冗余现象，比如"免疫反应"和"适应性免疫反应"同时都高亮。用enrichGO()之后接simplify()做一步去冗余：

ego_simplified <- simplify(ego_bp, cutoff = 0.7, by = "p.adjust")

cutoff是相似度阈值，默认0.7，越高越激进，去得越狠。处理完后，条目数量通常会缩水20%到40%，但它让结果更干净、更利于报告撰写。

3. KEGG富集：通路数据库的注释逻辑、物种代码坑与无参考物种兜底方案

3.1 KEGG是怎么把基因连到通路上的

KEGG跟GO最大的差别，在于它研究的不是"功能分类"，而是"基因在代谢和信号通路网络中的位置关系"。KEGG底层有一条隐式的映射链：基因 → KO（KEGG Orthology） → pathways。一个基因先通过序列/功能相似性归到某个KO编号（例如K04534对应TNF），然后这个KO编号被画到某个通路图谱上（例如hsa04668）。因此，做KEGG富集时，工具本质上是帮你判断"给定的一组差异基因，是不是在某些KO节点上扎堆"。

这个机制直接衍生出了几个实操层面的坑。坑一：KO编号和通路图谱一直在更新，同一个基因在不同版本的KEGG数据库里可能归属不同，所以分析时最好固定数据库版本，写论文时把KEGG的发布日期标上。坑二：KEGG里不是所有基因都有KO号，有些物种注释不全会导致大量差异基因映射不到通路，结果里通路数量少、且富集到的基因比率低，这很正常。坑三：enrichKEGG默认连在线KEGG API，网络不稳定时容易报错，后面专门讲。

3.2 enrichKEGG参数与物种代码表

kk <- enrichKEGG(gene = deg_entrez$ENTREZID, organism = "hsa", keyType = "kegg", pAdjustMethod = "BH", pvalueCutoff = 0.05)

organism参数是KEGG自己要的三字母/四字母物种代码，比如人类hsa、小鼠mmu、大鼠rno、斑马鱼dre、拟南芥ath、水稻osa。这个代码跟NCBI的taxid不一样，也别指望用常见拉丁名，写错直接报错或者返回空结果。完整代码表可以到KEGG官网的Organism页面查，但更快的办法是直接按我这份常用表对照：

keyType这里要特别提醒：虽然我们的基因是ENTREZID，但enrichKEGG里默认的keyType="kegg"对应的其实是"KEGG内部基因ID"（通常是ncbi-geneid形式，跟ENTREZID基本一致）。只要你用的是人类、小鼠这类注释完善的物种，这个参数不用动；换到别的物种，比如某些真菌、植物，基因ID可能不是ENTREZID体系，就要先通过KEGG的keggConv做映射，否则输入进去也查不到。

3.3 KEGG在线下载失败：换个思路接入本地KEGG数据

enrichKEGG()默认联网访问KEGG REST API。如果你用的是公司服务器、校园网、或国外数据库连接不畅的环境，批量跑多个物种时经常弹出403 Forbidden或无法连接到rest.kegg.jp。这个时候我有两个备选方案。

方案一：加参数use_internal_data=FALSE不行，那是对enrichWP的。对enrichKEGG来说，可以尝试从KEGG官网下载注释文件到本地，构建自定义的KEGG term2gene表，用enricher()函数跑：

# 假设已经下载好对应物种的kegg基因-通路对应表 # 格式是两列：pathway_id, gene_id kegg_term2gene <- read.table("kegg_hsa_term2gene.txt", header = TRUE, sep = "\t") kk <- enricher(gene = deg_entrez$ENTREZID, TERM2GENE = kegg_term2gene, pAdjustMethod = "BH", pvalueCutoff = 0.05)

这条路的好处是彻底离线可跑，重复性极高；缺点是通路表需要定期自己更新，而且没有KO层级信息，下游画通路图会麻烦。

方案二：用clusterProfiler配套的download_KEGG类函数先把KEGG数据缓存到本地，再跑。比较老的方式是：

library(KEGG.db) # 已不再更新，仅做参考

说实话，KEGG.db这个包已经不太行了。我更推荐直接用enricher()配合自己维护的term2gene文件，数据来源可以选KEGG官网的link和conv接口，也可以用gage包里的kegg.gsets函数生成物种通路的GSET列表。脚本化之后，公司服务器上脱离外网也能稳定复现，这在很多生产环境里是硬需求。

3.4 无参考物种或无GO/KEGG注释信息时的策略

这一点必须单独拎出来讲，因为它太常见了。很多人做的是无参考基因组的转录组，拿到的是Trinity或StringTie拼接出来的transcript序列，没有现成的GO/KEGG注释。此时直接跑enrichGO是不可行的。

我个人的建议流程是：先做"功能注释"再跑"富集分析"。功能注释的首选工具是eggNOG-mapper，它会基于eggNOG数据库的直系同源关系，把转录本序列一次性注释出GO term、KEGG KO、COG分类等。输出结果里通常有三列非常有用：query序列名、GO term列表、KEGG KO列表。拿到这个表之后，自己把"转录本ID到GO/KEGG"的两列关系拆出来，然后照样用enricher()做富集。整个流程前期工作量集中在序列注释这一步，后面反而很顺。

这里有个经验分享：eggNOG-mapper的--target_orthologs参数建议选one2one或best，能有效减少多拷贝同源基因造成的假注释；跑完之后，建议对注释结果做一个覆盖度统计（注释到的基因/总基因的比例），低于50%要长个心眼，后续富集结果可能严重偏向注释完好的方向。

4. 从差异表到富集结果：批量处理多比较组、参数选择与灵敏度分析

4.1 一个标准化的批量处理流程

真实项目很少只有一个比较组。疾病对照组、药物治疗组、不同时间点……我要面对的是十几个比较组，一个一个手工跑会累死，而且容易参数不一致。我的做法是写一个函数封装富集全流程，输入是一张DE结果表和比较组名单，输出是每个组的GO/KEGG富集结果对象。

run_enrichment <- function(res_table, orgdb, kegg_code, p_filter = 0.05, lfc_filter = 1, ont = "BP") { # 差异基因 sig_genes <- rownames(res_table[res_table$padj < p_filter & abs(res_table$log2FoldChange) > lfc_filter, ]) if (length(sig_genes) < 10) return(NULL) sig_entrez <- bitr(sig_genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = orgdb) # 背景 bg_entrez <- bitr(rownames(res_table), fromType = "SYMBOL", toType = "ENTREZID", OrgDb = orgdb) # GO ego <- enrichGO(gene = sig_entrez$ENTREZID, universe = bg_entrez$ENTREZID, OrgDb = orgdb, ont = ont, pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.2) # KEGG kk <- enrichKEGG(gene = sig_entrez$ENTREZID, organism = kegg_code, pvalueCutoff = 0.05) list(GO = ego, KEGG = kk) }

批量调用时用lapply按组遍历，输出存成RDS文件。这一步的经验是：每跑完一组就保存一次RDS，别全跑完一次性保存，原因很简单，中途服务器断连或内存崩了不至于全军覆没。我还会顺手记录每个组输入了多少差异基因、富集到多少条目，存储在一个Excel里，做QC用。

4.2 富集统计的底层逻辑和阈值到底怎么选

富集分析用的统计模型是超几何分布（或Fisher精确检验）。我把它的核心逻辑用最朴素的话说：在所有背景基因里，某个通路相关的基因占比是X；在你的差异基因列表里，这个通路相关基因的占比是Y；如果Y显著大于X，就说这个通路在差异基因中富集了。这里的"显著"，就是pvalue。做上千次这样的检验后，用BH方法校正FDR，得到qvalue。

参数选择方面，比较稳妥的是：pvalueCutoff在0.01到0.05之间，qvalueCutoff在0.1到0.2之间。需要说明一点：不要把pvalue和qvalue的阈值都设得很松（比如都是0.05），两者都设为0.05时实际上是在套两层过滤，结果几乎全是极显著条目，看似干净，实则过于严格，可能会丢掉一些真正的低频生物学信号。我一般是qvalueCutoff=0.2，给后续的人工筛选留空间。

4.3 sensitive analysis：结论稳不稳，换个阈值看看

"灵敏度分析"这个词在生信里用得不多，但几乎每个成熟的生信分析流程里都会做：把关键阈值上下浮动，看看富集结论是否稳定。尤其是差异基因筛选这一步的阈值，对富集结果影响极大。同一个项目，padj<0.05、|log2FC|>1筛选出来的基因可能有三四百个，换到padj<0.01、|log2FC|>1.5可能只剩几十个。富集到的通路可能完全变样。

操作上，我会至少跑三组：严格阈值、默认阈值、宽松阈值，对比每组排名前十的GO条目和KEGG通路。如果某个通路在三组中都稳定出现在前列，那这个结论才敢写进报告里。如果只有严格阈值下显著、宽松阈值下完全不见，那它很可能是被少数几个极端基因驱动的伪富集，需要警惕。

5. 结果可视化：barplot、dotplot、cnetplot、emapplot的选用逻辑与出图规范

5.1 常用图形的信息表达差异

结果可视化不是随便挑一张顺眼的图，每一种图回答的问题不同：

• barplot()：查看富集条目的显著性和数量对比，x轴是pvalue或qvalue，y轴是GO/KEGG条目名，适合展示Top富集列表。
• dotplot()：在barplot基础上多了一个维度——GeneRatio（差异基因中落在该通路的比例），用点大小表示，信息量更丰富，是目前论文里最主流的展示方式。
• cnetplot()：展示"差异基因-富集条目"的关系网络，能直观看到哪些基因同时参与多个通路，适合把笼统的"通路富集"落到具体基因层面。
• emapplot()：展示"富集条目之间"的相似性网络，通过共有的基因连接两个通路，适合找功能模块，在多个条目高度雷同时尤其好使。

作图时我会先看数据本身决定用图，而不是千篇一律dotplot。如果差异基因主要集中在少数几条通路上，且和实验处理强相关，我会优先用cnetplot，这张图在组会上讲故事最有效；如果通路数量多、想整体展示功能布局，dotplot更清晰。

barplot(ego_bp, showCategory = 15) dotplot(ego_bp, showCategory = 15) cnetplot(ego_bp, showCategory = 5, categorySize = "pvalue", colorEdge = TRUE) emapplot(pairwise_termsim(ego_bp), showCategory = 15)

注意emapplot()在R新版本里可能会要求先运行pairwise_termsim()计算条目相似性矩阵，否则报错。小细节，但容易卡壳。

5.2 出图与排版：论文级图形的保存方式

富集图的文字经常很长，GO条目描述动辄十几个单词，默认尺寸下必然拥挤。我的保存习惯是：

ggsave("dotplot_BP.png", plot = last_plot(), width = 10, height = 8, dpi = 300)

GO的BP条目如果显示前15个，8英寸高度勉强够；如果显示20个，高度调到10英寸以上。文字重叠时优先考虑拉宽画布而不是缩小字号，因为投稿时图片会被缩放到单栏宽度，字号太小容易糊。

多组结果放一起时，我用cowplot或patchwork拼图，右侧加统一的qvalue图例。还有个小技巧：GO三个ontology（BP/MF/CC）分别画一张dotplot，横过来拼成一个三合一大图，整张图能非常紧凑地展示功能倾向，这也是很多论文里的常见排版方式。

6. 实战避坑录：实际项目中反复踩过的五个问题

6.1 ID转换丢失大量基因，但误以为是正常情况

bitr()转换时，会存在一部分基因转不过去。尤其是从Ensembl ID转ENTREZID时，因为Ensembl ID的版本后缀可能和注释库对不上，或者某些基因确实没有对应的ENTREZID记录。我遇到过一次，一个四百基因的列表，转完只剩两百多个，当时差点就直接跑了，后来检查发现是Ensembl版本和org.Hs.eg.db版本跨了三个Release，部分新基因的ID在旧版本里还没收录。解决思路是定期更新Bioconductor注释包，或者用biomaRt::getBM()去取最新版本映射，不要用bitr()一把梭。

6.2 enrichKEGG返回空结果的系统性排查

这里有个没写在doc里的经验：当enrichKEGG()返回的结果里pvalue全为NA或者行数只有一两条时，多半不是你的基因没富集到通路，而是输入的基因ID和KEGG数据库的映射出错了。按这个顺序排查：

1. 确认organism代码正确。
2. 用keggConv()检查目标物种的ID类型是否和输入一致。
3. 对输入的每个基因，用keggFind()单独查一次，确认是否有映射。

在人类转录组里，排查效率最高的是第2步。因为keyType="kegg"默认情况下，enrichKEGG其实是参考了KEGGREST包的映射关系，有时候DESeq2结果表里的行名带着Ensembl版本号（例如ENSG000001... .14），直接提取去跑就会丢一大批。

6.3 富集结果大量重复，同一个通路出现多次不同ID

这个在我跑植物转录组时遇到过。KEGG里不同物种的同源通路，ID前缀不同但内容几乎一样，比如ko04110（KEGG直系同源通路）和hsa04110（人类特异通路）。当使用多物种整合注释时，同一个基因可能同时映射到多个不同前缀的通路ID上。解决方法是：在下游分析时统一用KO层的通路ID（ko前缀），不要区分物种。这一点在比较不同物种间的富集结果时特别重要，避免在比对时把同一生物学通路当成不同通路。

6.4 为了追求"全"把pvalueCutoff调成1

我理解有些同学想先看看全部富集条目再筛选，于是把pvalueCutoff=1跑，结果得到一张几千行的表，大部分是pvalue=0.9的垃圾条目。我不推荐这种做法。富集分析本身依赖多重检验校正，如果你在filter阶段就放进大量无差异信号，反而会稀释校正强度，让真正显著的条目排位下降。更好的做法是先用pvalueCutoff=0.05跑，如果觉得结果太少，可以稍微放宽到0.1，但不要一下放到1。

6.5 多组批量跑富集时的内存与缓存失控

用lapply批量跑几十个比较组时，每个enrichGO对象里都保存着庞大的背景基因和映射数据，如果不清理，R全局环境的内存占用会迅速攀升到十几GB。我常用两个手段：一是每跑完一组gc()一次；二是只保留必要的列和结果对象，把中间过程对象及时清除。遇到需要反复比较的组，直接存成RDS文件，要用的时候readRDS()单独载入，不要在一个会话里堆几十个结果对象。

回看整个GO/KEGG富集分析，最核心的心智模型就是"完备的注释+正确的背景+统一的ID+合理的阈值"。这四个条件全满足时，富集分析能给你提供非常扎实的生物学线索；任何一个环节偷懒，结果就只是统计学噪音的另一种表现形式。我见过不少项目最后被复核时，仅仅因为背景基因选择不当就被要求重跑全部富集，非常可惜。希望这篇文章能让你少走这些弯路，拿到差异基因列表后稳稳当当地把富集分析跑完，然后有底气地在组会上说："这批差异基因确实富集在XX通路里，而且换了阈值也稳。"