做转录组分析的同行应该都有过这个经历:差异表达跑完,挑出几百上千个基因,表格已经生成,但下一步怎么走很容易懵。拿着一张基因列表,老板问“这批基因和哪些生物学过程有关,富集在哪条通路里”,这时候就需要给这批基因做GO和KEGG注释。我见过不少人在这个环节卡住,有人卡在ID转换,有人选错工具,有人富集结果跑出来一大堆条目却不知道怎么筛。这篇东西就把我从拿到基因列表到出富集结果再到解读的全套流程讲清楚,工具选型、数据准备、代码、判读标准都覆盖到,适合刚接触组学分析的科研入门者,也适合做惯了RNA-seq但没系统整理过注释流程的从业者。
1. 从基因列表到生物学结论:GO与KEGG注释的核心逻辑
1.1 基因注释的本质是什么
先说清楚一个概念。GO和KEGG注释,本质上是把你手里的基因列表放到一个“标准词典”里去查含义。GO全称是Gene Ontology,它把基因功能拆成三个维度:生物学过程(Biological Process)、分子功能(Molecular Function)、细胞组分(Cellular Component)。打个比方,你拿到一份公司员工名单,想知道这家公司主要干什么,GO就是按“每个人在哪个部门、负责什么业务、有什么技能”给每个人贴标签,然后统计哪个部门人手最多、哪个业务最热门。
KEGG则是从“通路”的角度看问题。它把已知的代谢与信号通路整理成一张张地图,比如细胞凋亡通路、糖酵解通路、TNF信号通路,每个基因在地图上都有一个位置。你要做的事就是把基因列表里的名字丢到地图上,看它们在哪些通路上密集出现。KEGG的核心价值在于它自带通路图,能直观看到基因在通路中的上下游关系,这是GO做不到的。
我接触过不少刚入行的同学,最容易混淆的一点是把“注释”和“富集”当成两件事。其实它们是同一件事的两个侧面。单纯注释是每个基因逐个查它的功能标签,而富集分析是统计整批基因在某个功能类别或通路上的集中程度。通常我们说的“做GO和KEGG注释”,实际操作里绝大多数是指GO/KEGG富集分析,也就是检验你手里的基因列表相对基因组背景,是否显著富集在某类功能或通路上。
1.2 为什么要做富集而不是挨个查基因
有人问过我很实际的问题:几百个基因我挨个去数据库查不就行了?遇到一千个以上基因时,这个方法基本不可行。更关键的是,单个基因的功能注释给不了你系统层面的结论。比如你拿到200个差异基因,挨个看可能看到TP53、EGFR、AKT1这些熟人,但老板问“这批基因整体上跟什么通路相关”,你回答“有细胞增殖相关的成员”,这个结论太弱。
富集分析的原理是超几何分布检验,简单理解就是:从整个基因组里随机抓200个基因,某个通路有50个节点,理论上可能抓到3个;你现在实际抓到15个,这远远超出随机水平,于是得到一个很小的p值,说明这批基因在这个通路上发生了明显的“扎堆”。这个“扎堆”往往对应着样本处理前后真正被扰动的核心生物学过程。比如某实验室用某药物处理细胞后,差异基因显著富集在“细胞周期”和“p53信号通路”,那你就有了一个可以往下挖的生物学猜想,而不是面对一堆基因名无从下手。
1.3 在线工具与本地工具的大分流
明确了什么是要做的,接下来就要选工具。市面上的工具大致分两个阵营:一类是网页在线提交的,像DAVID、KOBAS、Enrichr、Metascape;另一类是本地运行的程序包,最常用的是R语言环境下的clusterProfiler。两个阵营用起来完全是两个思路。
在线工具的好处是没有编程门槛,上传基因列表、选物种、点提交就能拿到结果。它的代价是你能控制的东西很少——背景基因往往是工具内置的,数据库版本是固定的,p值校正方法不一定能自由选。本地工具的好处是流程可控、可复现,参数可以按需调整,缺点是至少得会基本R操作。从我多年的使用经验看,如果你只是临时分析一批基因,想快速看结果,在线工具够了;但如果这批基因是课题的核心数据,之后要反复调整参数、做对比、出发表级图表,建议直接上clusterProfiler,一步到位。
2. 在线注释还是本地R包:四类主流工具选型对比
2.1 四个我实际用过的工具横向对比
先把我实际用过的四类工具放在一张表里,再逐个讲我的使用感受和选择标准。
| 工具 | 运行方式 | 支持物种 | 默认背景 | 可定制性 | 典型用途 |
|---|---|---|---|---|---|
| DAVID | 网页提交 | 数十个模式物种 | 内置全基因组 | 中 | 快速出注释报告 |
| KOBAS | 网页/命令行 | 数百个物种 | 内置全基因组 | 中 | 无参物种支持好 |
| Enrichr | 网页/API | 多个模式物种 | 内置文库 | 低 | 快速多文库富集 |
| clusterProfiler | 本地R包 | 依赖OrgDb/在线KEGG | 自定或默认真实基因 | 高 | 发表级图表与批量流程 |
2.2 DAVID的使用心得与细节
DAVID是公认的老牌工具,也是很多人第一个接触的在线注释平台。它的操作流程很直接:打开页面把基因列表粘贴进去,选择基因ID类型,比如官方基因符号(Official Gene Symbol)或探针ID,再选物种和背景,提交后等一会儿就能看结果。
但DAVID有几个让我不怎么舒服的地方。第一,它的KEGG注释依赖KEGG数据库的API更新,而KEGG官方后来调整了规则,导致DAVID里的KEGG通路更新频率并不固定。第二,它一次提交的基因数量有限制,如果你做的是大规模转录组分析,有三千个以上差异基因,拆批次提交就很麻烦,因为拆开之后背景变了,结果可能对不上。第三,DAVID的GO/KEGG结果导出格式比较老旧,重新整理成发表需要的表格,得花不少功夫。
2.3 KOBAS与Enrichr的差异化优势
KOBAS在非模式物种和植物研究里口碑很好,原因很简单:KOBAS收录的物种数远超DAVID,它不仅能注释基因ID,还能对无参考基因组物种做在线的序列比对注释,这在很多动辄几十个基因但参考注释不全的项目里救过我的命。它的在线版操作也不复杂,粘贴基因列表,选物种,选注释范围(GO、KEGG、KEGG Orthology),提交后能得到一个比较规范的富集表。
Enrichr则是另一个极端——它更加“轻”,集成了大量基因集文库,除了GO和KEGG之外,还有各种转录因子靶标、细胞类型标记基因等。Enrichr的交互式可视化做得比较漂亮,适合快速探查基因列表的潜在功能方向。它的缺点也很明显:GO和KEGG的数据库更新相对慢,而且背景基因集不透明,不建议直接用它出最终发表结果。
2.4 为什么最终推荐本地R包
我的原则是:初筛用在线工具,出结果用本地工具。本地工具里首推clusterProfiler。逻辑很简单:一是可复现,脚本存下来,下次换一个基因列表,改一个变量就能重跑,数据管理上有据可查;二是背景可控,你可以自己指定表达矩阵中所有检测到表达的基因作为背景,避免因为背景差异导致结果失真;三是可视化组件完整,barplot、dotplot、cnetplot这些函数一条命令就能画出符合期刊投稿要求的图。
clusterProfiler入门也有一些前置成本:需要会装R包、处理依赖、理解几段代码。针对这个痛点,下文第三、四章专门把从数据准备到完整跑通的流程写出来,照着操作即可。
3. 基因ID转换与文件整理:注释前最容易被卡住的环节
3.1 三种最常见基因ID类型
我见过太多人卡在ID转换这个环节。RNA-seq的定量结果一般提供多种基因ID:基因符号(SYMBOL,比如TP53、BRCA1),Entrez ID(纯数字,比如7157是TP53),Ensembl ID(以ENSG开头的长字符串)。在线工具和R包要求的输入类型各不相同,有的直接认SYMBOL,有的非要Entrez ID,还有的得用Ensembl ID加版本号后缀。不搞清楚这个,你提交1000个基因,工具只匹配上300个,剩下的全部丢弃,做出来的富集结果自然失真。
一个比较稳妥的做法是:拿到基因列表后先判断自己手头是哪类ID,再根据目标工具的要求去做转换。如果目标工具是DAVID,它能自动识别多种ID类型,但识别率不是百分之百;如果目标工具是clusterProfiler的enrichGO,它最好用Entrez ID;enrichKEGG在最新版本里也推荐统一到Entrez ID再做。
3.2 用R包bitr做ID转换实操
R环境下做ID转换最常用的是AnnotationHub体系下的OrgDb包。人类用org.Hs.eg.db,小鼠用org.Mm.eg.db,其他物种就要去Bioconductor找对应的OrgDb。转换函数是clusterProfiler包的bitr。
我实际跑通的一段转换代码如下:
library(clusterProfiler) library(org.Hs.eg.db) # 这是示例基因列表,实际分析时从差异表达表格读取 gene_list <- read.csv("deg_symbols.csv", stringsAsFactors = FALSE) symbols <- gene_list$symbol # SYMBOL 转 ENTREZID conversion <- bitr(symbols, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db) # 检查转换成功率 cat("转换成功比例:", nrow(conversion) / length(symbols) * 100, "%\n")这里有个细节值得留意:转换后结果里可能出现重复行,比如一个符号对应两个Entrez ID。实际操作中我会先按“是否有重复”检查一遍,原则上保留第一次出现的行,再用unique去重。比例太低时,我会回头检查原始基因列表里是否有旧符号或过时编号,必要时用通用数据库的别名表批量修正。还有一个小经验:如果一部分基因是来自商业芯片平台,会附带Illumina探针ID或Affymetrix探针ID,这类ID需要额外的注释包来转换,bitr本身不直接处理。
3.3 在线工具的输入格式注意事项
走在线工具路线时,格式问题更琐碎。DAVID要求删掉表头,一列一个基因ID;KOBAS接受带表头的CSV,但要求指定哪一列是基因ID;Enrichr对纯文本格式最宽容,直接粘贴一列符号即可。我习惯把所有输入文件统一整理成如下形式:单列、无表头、无重复、无空行、无空格,保存为纯文本或CSV。这个“五无”标准能避免绝大多数在线工具的解析报错。
还有一个容易被忽略的坑:有些基因符号在Excel打开和复制粘贴时会被自动转为日期格式,比如“MARCH1”“SEPT1”这种符号,进Excel再出来就变成了“3月1日”“9月1日”。预防办法是导入Excel前把那一列设为文本格式,或者干脆全程用文本编辑器处理基因列表,不轻易经过Excel。
4. 本地实操:用clusterProfiler完成GO和KEGG注释的完整流程
4.1 环境配置与包安装
这一步是很多人第一个掉坑的地方。安装clusterProfiler需要先有R环境,建议直接用R 4.2以上的版本。默认的CRAN源安装不了这个包,必须用Bioconductor的安装命令。
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("clusterProfiler") BiocManager::install("org.Hs.eg.db") # 人的OrgDb包,小鼠用 org.Mm.eg.db安装过程中遇到缺什么依赖包就跟着装什么,比较常见的是igraph、tidyverse等。如果你所在网络环境下载较慢,可以把Bioconductor源换到镜像站,比如使用上海交大或清华的镜像,能明显提速。注意新版clusterProfiler依赖KEGG在线API,跑enrichKEGG时网络不稳定会导致连接失败,但这属于运行时问题,下一章细讲。
4.2 GO富集分析的完整代码与参数解读
数据准备完成、环境就绪后,GO富集的核心代码非常简洁。下面这段是我常用的标准流程。
library(clusterProfiler) library(org.Hs.eg.db) # 1. 读取基因列表并转换ID gene_symbols <- readLines("gene_list.txt") # 一列SYMBOL entrez_result <- bitr(gene_symbols, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db) gene_entrez <- unique(entrez_result$ENTREZID) # 2. 做GO富集,ont可选BP/MF/CC,ALL表示三个都跑 ego <- enrichGO(gene = gene_entrez, OrgDb = org.Hs.eg.db, keyType = "ENTREZID", ont = "BP", # BP/MF/CC/ALL pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.2, readable = TRUE) # 3. 结果保存 write.csv(as.data.frame(ego), "GO_BP_results.csv", row.names = FALSE)这里每个参数背后都有讲究。pAdjustMethod选择“BH”是Benjamini-Hochberg校正,这是组学分析的事实标准,控制的是错误发现率,比Bonferroni温和得多。pvalueCutoff默认是0.05,qvalueCutoff默认是0.2,后者更严格,因为qvalue是经过校正和控制假发现率之后的值。我在多数课题里会先把qvalueCutoff定在0.2,先看整体富集面貌,再根据后续验证需要收紧到0.05。
ont参数如果你不设为“ALL”则会只跑一个维度。我经常遇到的问题是:用户只跑BP,结果看起来通路杂乱;后来把MF和CC也跑出来,发现真正的信号在MF里。我自己的习惯是三个维度都先跑一遍,保存后合并到同一份Excel工作簿里,分别命名为BP、MF、CC,这样后续翻看非常方便。
readable = TRUE这个参数很实用,它会把结果里的Entrez ID自动替换回基因符号,生成的表格直接能贴进文章补充材料。
4.3 KEGG富集分析的完整代码与常见bug修复
KEGG富集的代码在clusterProfiler新版本里发生过参数变动,这也是很多旧教程失效的原因。新版用法如下:
# KEGG 富集需要物种缩写,人类是 hsa,小鼠是 mmu ekegg <- enrichKEGG(gene = gene_entrez, organism = "hsa", keyType = "kegg", pvalueCutoff = 0.05, qvalueCutoff = 0.2) # 保存结果 write.csv(as.data.frame(ekegg), "KEGG_results.csv", row.names = FALSE)新版有个变化需要注意:旧版本里有个use_internal_data参数,现在在线查询KEGG REST API获取注释数据;如果调用时提示连接超时,解决思路分三步:首先检查网络,确认域名可访问;其次在代码里加上Sys.sleep(3)之类的间隔请求逻辑,把大量基因拆成几批去查;最后才是考虑换用旧版本离线注释数据包。
4.4 让GO结果更整洁:去冗余与多列表比较
GO富集结果有个常见毛病:条目大量重叠。比如“细胞增殖的正调控”“细胞增殖的调控”“细胞增殖”同时出现,因为它们共享大量基因,导致p值都很小。这时可以用simplify函数按语义相似度去冗余,保留代表性条目。
ego_simplified <- simplify(ego, cutoff = 0.7, by = "p.adjust", select_fun = min)如果你同时比较多个基因组的富集情况,比如对照组和给药组分别两批基因,clusterProfiler还提供compareCluster函数。先将两组基因打好分组标签,然后一条语句实现分组富集,输出结果直接用于后续分组柱状图。此部分设计我一般会在项目开始前就规划好,因为后期再补标签,工作量远大于跑一次富集。
4.5 可视化出图参数推荐
我最后出发表图一般只用三种图。barplot展示显著条目数量与p值,适合看出大致全貌;dotplot以气泡大小表示富集基因数,颜色表示p值,信息量最大;cnetplot把基因与通路的网络关系画出来,一眼能看出核心基因。
barplot(ego, showCategory = 20) dotplot(ego, showCategory = 20) cnetplot(ego, showCategory = 5) # 保存为专业出版格式 ggsave("GO_dotplot.pdf", width = 10, height = 8) ggsave("GO_dotplot.png", width = 10, height = 8, dpi = 300)出图时要控制showCategory数量。我很少超过20个,太多会压扁看不清;太少会遗漏重要条目。另外dotplot默认按p值排序,如果你的结果里有一个通路基因数特别多,它在图中可能不突出,可以手动调整排序参数。
5. 富集结果判读与可视化:别让p值成为唯一的答案
5.1 富集表里每个字段的现实含义
跑完的结果是一个DataFrame或在线表格,通常包含几列:ID、Description、GeneRatio、BgRatio、pvalue、p.adjust、qvalue、geneID、Count。我看到很多人只盯p.adjust,其实前几列信息量非常大。
GeneRatio的含义是“富集到该条目下的输入基因数 / 输入总基因数”。举个例子,输入200个差异基因,富集到“DNA修复”的基因有30个,GeneRatio就是30/200。BgRatio则是“该条目在背景基因中的基因数 / 背景总基因数”。两者相除能得到Fold Enrichment,这个值大于1说明该条目在输入列表里比随机背景更集中。我在实际解读中会特别关注那些p值很小但Count只有2~3的条目,这类结果虽然显著,生物学意义往往单薄;相反Count达到几十的大条目,即使p值稍大,更值得深挖。
5.2 忽视背景基因会导致什么恶果
这是我最想强调的一点。很多人跑出结果后抱怨“怎么富集出来的全是些基础代谢通路,跟我的研究方向毫无关系”,十有八九是没管背景基因。clusterProfiler默认以模式物种全部基因作为背景,但更合理的背景是你样本中实际检测到表达的基因。比如你在某个细胞系里做了敲除实验,本身这个细胞系不上调“神经发育”通路,但你用全基因组背景做检验时,可能因为输入基因里有几个边缘相关的成员,就“显著”富集到无关通路。
clusterProfiler里有两种方式指定背景:buildAnno函数(已废弃的旧方式)或使用compareCluster的universe参数。更常见的做法是先从表达数据里筛出有表达的基因集合,然后作为gene参数传入fect并设置universe。这里放一个对比示例:
ego_with_bg <- enrichGO(gene = gene_entrez, universe = background_entrez, # 表达矩阵中过滤后的背景基因 OrgDb = org.Hs.eg.db, keyType = "ENTREZID", ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.2, readable = TRUE)用不用合理背景,结果差别非常大。某课题组曾经因为默认背景,一个药物处理实验富集出一堆心肌收缩相关条目,换用实际检测表达的背景后,这批信号完全消失,换成了跟药物靶点代谢直接相关的通路。所以你在论文方法部分写“背景基因为所有检测到的蛋白质编码基因”是有实际意义的,不是套话。
5.3 多组学结果放一起时怎么确定优先级
当你同时拥有mRNA、miRNA或蛋白组的多个基因列表时,富集结果需要横向比较。最常见做法是把富集分析结果按p值排序,取每组Top10条目,做韦恩图看交集。这里要小心一点:KEGG通路ID在不同组件里都叫hsa04150之类,直接ID对齐没问题;GO条目则可能因为BP/MF/CC维度不同,ID前缀也不同。
优先级判断上,我给一个比较实用的原则:优先关注那些在多个分子层面同时富集到的通路,其次是关键差异基因直接参与的通路,最后才是单纯p值小的通路。很多人做多组学只追求p值交集,忽略了基因层面的证据,结果筛选出的通路跟实际表达变化方向对不上。我一般会把富集到的通路再拿去做基因层面的趋势一致性检查,比如某通路被富集,通路内的关键基因在组间的表达变化方向是否一致。
5.4 生成报告时的常见消费习惯
跑完富集后把结果拷进Word就完事,这样容易让导师抓细节漏洞。我通常会在方法里写清楚:使用clusterProfiler进行了GO功能富集分析,以BH方法校正p值,显著富集阈值设定为p.adjust小于0.05;KEGG通路富集使用KEGG数据库,以qvalue小于0.2为显著。这比只写“做了富集分析”要严谨得多。这一行方法描述放在文章里,审稿人通常不会追问过多。
6. 高频问题排查与避坑实录
6.1 ID匹配率低:绝大多数问题的源头
症状:提交几百个基因,富集结果里只出现二三十条ID,提示大量基因无法映射。原因无非三类:ID类型给错、基因列表里混入了版本号或无法识别的特殊字符、物种的注释包没装对。我排错的顺序是:先检查前10个基因是不是正常格式,再确认OrgDb包和物种匹配,最后检查基因列表是否有多余空格或中文标点。有一个很隐蔽的情况:很多下载的公开数据里基因符号全是小写,比如tp53而非TP53,bitr在这种情况下的匹配率会瞬间降到一半以下。解决方法是先统一大小写,遵守官方基因命名规则。
6.2 enrichKEGG连不上服务器
新版clusterProfiler跑enrichKEGG会去访问KEGG官方API,网络不稳定时报错很频繁,报错信息往往让人以为是自己代码写错了。解决办法我前面提过,这里补充详细操作:先把基因列表拆成每200个一组,每组查询之间加Sys.sleep(2);如果依然频繁失败,可以考虑从KEGG FTP下载物种全基因组的KO注释文件离线完成注释,再用clusterProfiler配合离线文件跑富集。这个操作需要配置KEGG文件路径,相对繁琐,非必要不推荐普通用户处理。日常分析遇到这类问题,我最实用的做法是在早晨网络情况较稳定的时段重跑,成功率会高不少。
6.3 富集结果大量重复冗长条目
这种情况最常见于GO的BP维度。用simplify去冗余是首选方案,但要注意cutoff的度。我试过cutoff从0.5一路试到0.9,发现0.7是一个平衡点:既能去掉明显重复的大类,又不至于把相关但不同的过程合并掉。去冗余之后再看结果,通常能从三十个纠缠不清的条目里清理出七八个明晰的主题,报告写起来也舒服很多。
6.4 设定阈值的黄金法则
p值调得很严,结果几乎为空;调得很松,出来几百条全是无效信息。我见过新手在pvalueCutoff=0.01时跑GO,跑完一个条目都没有,直接怀疑代码出问题。但基因数差异也很大的情况下,阈值本来就应该是浮动的。我的经验法则:输入基因在500个以上时,pvalueCutoff用0.05比较合理,qvalueCutoff先用0.2探路;输入基因不足100个时,可以适当放宽到pvalueCutoff=0.1,重点看富集趋势而不是严格显著性。当然,这些阈值在最终方法描述里必须如实写明,不能做了0.2的报告转头对外宣称是0.05。
6.5 在线工具结果格式输出与再整理
在线工具给的结果一般是一个表格,包含Term、PValue、Benjamini等列。直接拿去发表的话格式不够规范。我通常把在线结果的Term列拆成ID和描述两部分,再手动补上GeneRatio或观察到的基因名称,统一成与clusterProfiler输出兼容的表格结构。这里偷懒不得,因为期刊编辑越来越喜欢看统计数据齐全的富集表。另一个细节是引用问题,在线工具对应的数据库和算法均有对应引用文献,交稿前务必核对引用列表是否有遗漏。
6.6 运行速度慢的优化策略
大基因列表跑GO富集其实很快,慢多半是KEGG在线注释在拖后腿。本地优化策略有三个:一是在bitr转换时就去掉所有无法映射的基因,二是把universe背景基因尽量精简到实际检测的蛋白编码基因,三是多利用并行运算。clusterProfiler本身支持并行,在enrichGO中设置BPPARAM参数为多核配置即可。一个一万个背景基因的样本,开八个线程,整体耗时能从十几分钟降到三分钟左右,效果明显。
7. 几个真实场景的实战复盘
7.1 八百个DEGs的标准处理流程复盘
某次我为一份接受某干扰处理后的转录组数据做分析,差异基因有八百多个。我的处理路线是:先从差异表格提取SYMBOL,用bitr转出792个Entrez ID,再取表达矩阵所有检测到且表达量过滤后的约一万两千个基因作为背景。GO三个维度各跑一遍,KEGG跑一遍。结果里BP方向最明显的条目是“对外界刺激响应的调控”,KEGG里最显著的是“TNF信号通路”,这两个结论和该条件已知的生物学效应完全吻合,最终报告只花了半小时就整理完成。
如果我当时图省事直接用默认背景去跑,部分与刺激响应相关的基因背景比例会膨胀,最终富集出来的通路可能截然不同,后续验证实验的方向也会被带偏。这个案例再次印证合理指定universe的重要性。
7.2 非模式物种基因列表的注释自救方案
非模式物种最恶心的地方在于没有对应的OrgDb包,clusterProfiler常规流程直接跑不起来。这种情况下我有两条路:如果该物种基因组有GFF注释,先自己用AnnotationHub或自建TxDb从GFF生成OrgDb,之后流程完全一样;如果连GFF都不太完整,就退回到KOBAS在线工具,它能接受蛋白序列提交,用序列同源比对的方式注释到KEGG通路,虽然麻烦,但胜在多一个解决方案。
7.3 两批基因的联合分析实操
一个常见需求:比较敲除前后两批差异基因的功能变化。用compareCluster写起来很清晰:
gene_list_by_group <- list( treatment_up = c(...), treatment_down = c(...) ) ck <- compareCluster(geneClusters = gene_list_by_group, fun = "enrichKEGG", organism = "hsa", pvalueCutoff = 0.05) dotplot(ck)结果会自动按分组列出各自的显著通路。这种分组的优势是能看出药物处理中哪些通路上调、哪些下调,而不是所有差异基因混在一起平均化。我习惯配合组间表达方向再看一遍原有数据,避免因为上下调抵消隐藏了真实信号。
8. 我的最终选择与给你的建议
这些年下来,我做基因注释的固定动作已经形成肌肉记忆:先检查ID类型,再做ID转换,随后以表达检测背景跑clusterProfiler的GO和KEGG,最后用dotplot和cnetplot出图。遇到在线工具好像更快的诱惑,我心里会提醒自己一句:在线工具节约的是几分钟的点击时间,但损失的是背景可控性、可复现性和数据的标准化程度,对正式课题来说得不偿失。
如果你只做一次探索性分析,直接上DAVID或KOBAS没有问题;但如果这个基因列表之后还要跟其他数据比较,或要反复出现在论文里,请一定从一开始就走上本地R包的路线。这个技术路线入门多花的半天时间,会在后面无数个日夜里帮你省下成倍的沟通成本和时间成本。
还有一个小建议,把这个流程中所有用于ID转换和富集的参数,包括注释包版本、R版本,一并写进项目的记录文件里。我吃过版本更新后结果无法复现的亏,当时项目已经过了大半年,重新安装旧版本包费了很大劲。如果你把这些信息从一开始就记好,未来回看会轻松得多。