☰
基因富集分析全解析:从超几何分布到GO/KEGG实操
2026/10/6 3:49:15 网站建设 项目流程

1. 从“拿到一堆基因”到“看懂生物学故事”,富集分析到底在解决什么问题

做组学数据分析的人,心里基本都有过这么一段经历:好不容易把转录组、蛋白组甚至ChIP-seq的数据跑完,差异表达基因列表也筛出来了,几千个基因上调下调摆在面前。然后呢?盯着这个列表,除了能看出几个眼熟的明星基因,剩下的一大半完全不知道它们在干什么。

这不是能力问题,而是单基因逐个解释这件事,在规模上根本不现实。人类的基因组有两万多个蛋白编码基因,一次差异分析筛出几百上千个基因再正常不过。你不可能一个个去查文献,也没必要——因为生物学过程本身就不是靠单个基因独立完成的,而是靠一整套基因协同运作。某个通路激活,一组基因一起上调;某个生物学过程被抑制,另一组基因整体沉默。如果你只看单个基因,看到的是一片树叶;富集分析做的事情,是让你退后一步,看到整片森林的轮廓。

我当时第一次跑完RNA-seq拿到2500多个差异基因时,整个人是懵的。老板问“这些基因主要涉及什么通路”,我只能支支吾吾说不出个所以然。后来把富集分析跑完,输出结果的第一屏就是几十条显著富集的通路,那种“豁然开朗”的感觉,做生信的人一定都能共鸣——原来这一堆基因背后,是炎症反应、T细胞活化、细胞因子受体通路这几个核心故事。之前的一堆散沙,瞬间变成了几张有逻辑的图。

所以富集分析这个词,翻译成大白话就是:把你手里的基因列表,放到一个预先定义好的“基因功能注释数据库”里,去统计哪些功能分类、哪些信号通路、哪些疾病关联在你的基因列表里被显著富集。这里的核心思想是一个超几何分布的统计学判断——如果某个通路里的基因在你的列表里出现了远超随机期望的数目,那这个通路就和你的实验条件密切相关。

这篇文章我打算从原理讲到实操,再讲到图表解读和排坑。适合刚拿到测序结果却不知道怎么往下走的新手,也适合已经跑过富集分析但总觉得结果解释得不够透的人。我会以最常见的GO和KEGG分析为主线,穿插我个人这几年在不同项目里积累下来的经验,不绕弯子,直接讲干货。

2. 原理必须理解到位:Q值、背景基因和超几何分布,三者决定结果可信度

很多人跑富集分析的时候,其实就是点两下鼠标或者敲几行代码,出来一个表格就算完事。但如果你不理解结果里那些数字到底是怎么算出来的,你根本没法判断这个结果靠不靠谱,更没法在审稿人问“你这个富集分析参数怎么设的”的时候给出一个得体的回答。

2.1 超几何分布:富集分析背后的统计引擎

富集分析的经典统计学基础是超几何分布。这玩意儿听起来吓人,但可以用一个特别生活化的例子讲清楚。

假设你面前有一个大缸,缸里有10000个球,其中500个是红球(代表某个通路里的全部基因,我们称为通路基因集),9500个是白球(剩余的所有基因)。现在你闭上眼睛从缸里随机抓取200个球(代表你做实验筛选出来的差异基因),理论上你应该抓到 (200 \times 500 / 10000 = 10) 个红球。如果你实际抓到了35个红球,你就得思考一个问题:这纯属手气好,还是缸里的球分布本来就有猫腻?

富集分析的逻辑完全一样。它的零假设是:你的差异基因列表是从所有基因中随机抽样得到的,通路基因在差异基因里出现的比例,应当与通路基因在整个基因组中的比例一致。如果实际比例显著高于期望比例,我们就拒绝零假设,认为这个通路在你的条件下被富集了。超几何分布模型算出来的p值,就是“手气这么好”的概率有多大。p值越小,说明这种富集越不可能是随机碰巧。

这里有个常见的理解误区:很多人以为富集分析是“用我的基因去做分类”,其实不是。它本质上是一个抽样检验,检验的是你的基因列表和数据库里预先定义好的通路基因集之间是否存在显著重叠。

2.2 基因数最少的通路往往最容易“显著”——先搞懂分母逻辑

我在带新人做分析的时候,经常被问到这样一个问题:为什么我富集出来的结果里,排在最前面的通路都只有几十个基因,而一些上千个基因的大通路(比如MAPK信号通路、PI3K-Akt通路)反而排名不高?

原因是:一个通路里的基因总数越少,越容易在随机情况下产生看起来“显著”的富集。举个极端例子,如果一个通路只有3个基因,你的差异列表里恰好包含这3个基因,那么p值会小得惊人,因为随机情况下同时抽中这3个基因的概率极低。但一个包含1200个基因的通路,你从中抽出30个基因,这在随机抽样下可能并不算罕见,p值反而可能不显著。

这就是为什么富集结果排在前面的通路,不见得是生物学上最重要的通路,而只是统计学上最容易“撞中”的通路。我个人的习惯是:结果出来之后不看排名,先看富集基因的数目。如果一个通路的p值很显著但富集的基因数只有2到3个,一般不会优先去讨论它——除非这2个基因本身就是文献里反复验证过的重要角色。相反,那些富集了三四十个基因的大通路,即使p值排名不是第一,往往才是你实验表型里真正的主线通路。

2.3 p值必须换算成Q值,做多重检验校正

另一个绕不开的概念是多重检验校正。你做一次GO分析,富集到的术语条目可能是几百上千条,每个条目都有一个p值。其中每条p值小于0.05的条目,都有约5%的概率是假阳性。几百条检验同时做,假阳性的数量会相当可观——甚至可以说,如果不做校正,你的结果表里可能有十几条甚至几十条是“撞大运”撞出来的。

所以论文里通常要求报告Q值(校正后的p值)。常见的校正方法有Benjamini-Hochberg(BH)、Bonferroni、Holm等。其中BH法控制的是错误发现率(FDR),简单说就是:在所有被判断为显著的结果里,期望的假阳性比例是多少。做GO/KEGG富集分析,BH是默认选择,它比Bonferroni温和得多,能在控制假阳性的同时尽量保留真实信号,也是大多数主流工具自带的方法。

我看到不少人在结果里只写p值不写Q值,或者画图的时候只用p<0.05做筛选阈值。如果你投稿的是专业生信期刊,审稿人大概率会质疑这一点。我的建议是:看图或看表时,除非你真的只需要粗略看个趋势,否则永远以q值为准。比较严格的筛选可以用q<0.01,常规的用q<0.05,尽量不要用原始p值直接下结论。

2.4 背景基因集设不对,结果废一半

这是富集分析里最容易被忽视、但影响最大的一个参数。所谓背景基因(background genes),就是超几何分布公式里的“10000个球”这一项——你用来计算“随机期望比例”的基因全集。

如果背景设得不对,结果会出现系统性偏差。我举一个很典型的场景:你的实验用的是只包含免疫相关基因的芯片,但做富集分析时背景基因选成了“人类全部蛋白编码基因”。由于你的芯片本身只测了免疫基因,免疫相关通路在你的差异列表里自然占据了极高的比例,计算结果会告诉你“全部通路都和免疫相关”——这就是假象,不能写进论文。

反过来,如果你做的是全转录组测序,背景用只包含几千个基因的芯片注释集,那也会严重低估很多通路的显著性。正确做法是:背景基因 = 你这项实验实际检测到的所有基因(或所有表达量达到检测阈值的基因)。用R跑clusterProfiler的时候,universe参数就是干这个用的,很多人没有设置,工具就会默认用数据库的全部基因。对于全转录组数据来说这影响不是太大,但对于靶向测序、芯片数据或者某些表达量过滤严格的转录组数据,不设背景就是给自己挖坑。

3. GO和KEGG的分工关系,以及主流富集工具的真实水平对比

知道了原理之后,下一步就是选工具。但选工具之前先把GO和KEGG这两个概念理清楚,因为这是所有人刚开始接触富集分析时最容易混淆的地方。

3.1 GO的三个维度,以及为什么“分子功能”结果最不好解释

GO(Gene Ontology,基因本体论)是一个结构化的功能注释体系,它把所有基因的功能描述组织成一个有向无环图。GO注释分三大类:分子功能(Molecular Function,MF)、生物学过程(Biological Process,BP)和细胞组分(Cellular Component,CC)。

  • BP:描述基因参与的生物学过程,比如“炎症反应”“DNA修复”“细胞增殖的负调控”。BP层面的富集结果最直观、最好讲故事,也是我在实际分析中最先看的一类。
  • MF:描述基因的分子活性,比如“ATP结合”“蛋白激酶活性”。MF层面的结果通常比较枯燥,而且很多MF条目在不同通路里都会出现,解释余地不大。
  • CC:描述基因产物所在的细胞位置,比如“线粒体”“细胞膜”“核小体”。CC富集和实验的处理方式(比如是否做了亚细胞组分分离)关系很大,有一定的参考价值,但通常不是论文讨论的重点。

需要提醒的是,GO条目的层级结构问题会被很多人忽略。GO里有层级关系,上层术语很宽泛(比如“生物学过程”这种级别),下层术语很具体(比如“T细胞受体信号通路的正调控”)。富集分析得到的条目如果层级太粗,解释起来等于没有信息量;如果太细,又可能因为覆盖基因太少而不稳定。我的习惯是优先关注富集基因数目在10到数百之间、且术语描述中有明确生物学行为动词的BP条目,这种层级的信息熵刚刚好。

3.2 KEGG通路图和简单的“通路是基因的协作网络”认知

如果说GO是一本“功能词典”,那KEGG就是一张“通路地图”。KEGG(Kyoto Encyclopedia of Genes and Genomes)最核心的资源是KEGG PATHWAY数据库,它把基因之间的相互作用关系整理成了可视化的通路图——比如Toll样受体信号通路、糖酵解通路、p53信号通路等。

KEGG通路图里的每个节点是基因(或基因产物),节点之间的连线代表激活、抑制、磷酸化、转录调控等关系。富集分析告诉你哪些通路在你的条件下被激活了,通路图则告诉你这条通路里到底发生了什么事——比如你的差异基因是富集在通路的起始受体端,还是在末端的转录因子,位置不同生物学含义天差地别。

我在做炎症相关项目的时候,就遇到过KEGG富集出“IL-17信号通路”、但读图发现差异基因全部集中在负调控因子那一侧的情况。如果不看通路图,直接说“该通路被激活”,会被同行一眼看出问题。所以我的原则是:GO富集用来看方向,KEGG富集用来讲故事,而KEGG通路图用来核对故事讲得对不对。

3.3 工具实测:clusterProfiler、DAVID、Enrichr、g:Profiler的优缺点对比

市面上富集分析工具多得离谱,我挑四个最常用的对比一下。

工具平台优点缺点适合场景
clusterProfilerR支持数十种物种和注释包;支持GO、KEGG、Reactome、MSigDB等;结果对象可自由绘图;支持GSEA需要R基础;依赖Bioconductor注释包;本地跑较大的GSEA时占内存常规转录组/蛋白组/甲基化,推荐首选
DAVID网页操作简单;自带通路图跳转;输出格式友好数据库更新滞后;基因ID格式要求死板;物种覆盖有限;无法批量自定义背景想快速看结果、不想写代码的入门用户
Enrichr网页内置海量基因集库(含ChEA、TRRUST等);支持多种基因ID格式;界面清爽注释库质量参差;对多样化背景基因支持弱;不适合正式发表级分析快速验证思路、看转录因子富集、单细胞数据初步探索
g:Profiler网页/R包支持多物种;内置g:SCS多重校正;操作简单可定制性较低;通路资源以GO、KEGG及少量其他为主快速批量分析多组基因列表时效率极高

用这么久下来,我最认可的方案是:主力用R的clusterProfiler做正式分析,结果需要看一眼趋势或者跟审稿人快速解释时,拿Enrichr或g:Profiler补一个交叉验证。交叉验证的价值在于,不同工具背后的数据库版本和注释策略不同,如果一个富集结果能在两个独立工具中得到一致的结论,那这个结论的底气就足很多。如果两边的结果不一致,也别急着说谁错了——先检查基因ID转换是否正确,再看数据库版本,绝大多数不一致都是这两个原因造成的。

顺便提一句,DAVID虽然对新用户友好,但数据库版本更新缓慢是硬伤,尤其在做非模式生物分析时,注释覆盖率会明显不够。能用R解决的,尽量别把时间耗在网页工具上,因为一旦基因数量多了、比较组多了,网页操作的低效会让你怀疑人生。

4. 实战演示:用R完成从基因列表到富集图表的一整套流程

这一节进入实操环节。我用R的clusterProfiler包,给你完整跑一遍基因富集分析的流程。这套流程我在多个项目里反复验证过,直接复制改文件路径就能用。

4.1 事前准备:安装R包与确认基因ID格式

首先你需要安装几个R包。如果已经装过,可以跳过这一步。

if (!requireNamespace("BiocManager", quietly = TRUE)) { install.packages("BiocManager") } BiocManager::install(c( "clusterProfiler", "org.Hs.eg.db", # 人类注释包,其他物种替换为对应的org包 "enrichplot", # 绘图扩展包 "AnnotationDbi", "pathview" # KEGG通路图可视化 ))

安装之后,别急着跑,先想清楚你的基因ID是什么格式。clusterProfiler对ID格式是有明确要求的——bitr()函数主要支持Entrez ID、Ensembl ID、SYMBOL(基因符号)等格式。如果你手里的文件是“TP53”“EGFR”这种SYMBOL格式,后续必须先用bitr()转换成Entrez ID,因为KEGG数据库的注释逻辑是以Entrez ID为核心的。

这里我把转换代码直接放出来:

library(clusterProfiler) library(org.Hs.eg.db) # 假设你有一个差异基因列表,里面是SYMBOL gene_symbols <- c("TP53", "EGFR", "AKT1", "IL6", "TNF", "MYC") # SYMBOL转Entrez ID gene_entrez <- bitr(gene_symbols, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db) head(gene_entrez)

需要注意:bitr()转换之后,可能有部分基因因为注释不完整而丢掉了,这是正常的。关键是要确认转换率没有低到离谱。如果1000个基因只转出200个,那八成是你的输入ID格式写错了(比如把Ensembl ID当成了SYMBOL),先回头检查数据。

4.2 GO富集分析代码:三个ont参数一次跑全

以下是GO富集分析的核心代码,我顺序标注了参数的意图。

# GO富集分析:BP维度 go_bp <- enrichGO( gene = gene_entrez$ENTREZID, # 差异基因Entrez ID universe = rownames(expr_matrix), # 背景基因,建议换成你自己数据的全部基因ID OrgDb = org.Hs.eg.db, keyType = "ENTREZID", ont = "BP", # 可选 BP / MF / CC / ALL pAdjustMethod = "BH", # 多重检验校正方法 pvalueCutoff = 0.05, qvalueCutoff = 0.05, readable = TRUE # 结果表显示SYMBOL而非Entrez ID,阅读友好 ) # 查看结果概览 head(as.data.frame(go_bp))

这里有几个参数值得多说几句:

  • universe是背景基因。我刚才强调了背景的重要性,如果你手头是全转录组数据,可以把所有检测到的基因的Entrez ID放进去。如果你没有提供,enrichGO会用OrgDb里的全部基因做背景,对于全转录组来说问题不大,但如果你过滤掉了大量低表达基因,建议还是显式指定。
  • readable = TRUE会把结果表里的geneID列从Entrez ID转换成基因符号,这对后续解释结果来说太重要了。不然你看到一串数字ID,还得手动去查。
  • ont = "BP"是可以换的。需要三维度全覆盖时,可以设置ont = "ALL",结果里会有一列ONTOLOGY标记每条记录的层次。

4.3 KEGG富集分析:记得处理物种缩写

KEGG富集分析相对GO更依赖物种注释。

# KEGG富集分析,注意hsa是人(Homo sapiens)的KEGG缩写 kegg <- enrichKEGG( gene = gene_entrez$ENTREZID, organism = "hsa", # 小鼠是mmu,大鼠是rno,斑马鱼是dre,果蝇是dme keyType = "kegg", pvalueCutoff = 0.05, qvalueCutoff = 0.05, use_internal_data = FALSE ) head(as.data.frame(kegg))

KEGG数据库的物种覆盖远不如GO宽,如果你在跑一个冷门物种,很可能跑出来结果寥寥甚至直接报错。此时可以考虑换用KEGG的在线下载方案,或者干脆只用GO加Reactome分析。这里我不展开讲Reactome了,但对非模式生物,ReactomePA包是个不错的替补。

跑完后你会得到一个enrichResult类的对象。平时我们最关心的列包括:ID、Description、GeneRatio、BgRatio、pvalue、p.adjust、qvalue、geneID、Count。其中GeneRatio是差异基因中属于该通路的比例,BgRatio是全部背景基因中属于该通路的比例,这两个值结合起来看,能帮你快速判断富集的强度。

4.4 结果可视化的标准套餐:条形图、气泡图、通路图、富集网络图

有了富集结果,下一步就是画图。clusterProfiler配合enrichplot,可以很轻松地完成几种标准图。

条形图(Bar Plot)和气泡图(Dot Plot)

# 气泡图,推荐优先使用 library(enrichplot) dotplot(kegg, showCategory = 20, title = "KEGG Pathway Enrichment")

气泡图是富集分析里信息密度最高的一种图。纵轴是通路名称,横轴是GeneRatio(富集比例),点的大小代表富集到的基因数目,颜色代表p.adjust或q值的大小。我建议默认就画气泡图,因为条形图只能表达p值和富集基因数目两个维度,而气泡图可以同时表达三个维度。

通路图(Pathway View)

KEGG富集出通路名之后,可以用pathview把差异基因标到通路图上:

# 选出你想查看的通路ID,比如hsa04657(IL-17 signaling pathway) pathview( gene.data = logFC_vector, # 需要名称是Entrez ID且有logFC值的向量 pathway.id = "04657", species = "hsa", out.suffix = "IL17" )

跑完以后工作目录下会生成一张带基因颜色标记的通路图。红色代表上调,绿色代表下调。这张图在写论文时价值极高,也是我看数据时必做的一步。

富集网络图(Enrichment Map)

如果富集出的通路太多,可以直接做一张网络图来看通路之间的重叠关系:

ego <- enrichplot::pairwise_termsim(kegg) emapplot(ego, showCategory = 30, cex_category = 0.5)

网络图里每个节点是一条通路,节点之间的连线表示两条通路共享了大量基因。看到连成一大片的通路网络,通常意味着你的差异基因集中在某几个核心生物程序里——这种“抱团”现象本身就是重要的生物学信号。

5. 富集结果到底怎么读?图表解读中常见的误判和陷阱

跑出图和表之后,最后一公里的解读才是真正区分新手和老手的地方。这一节把最常见的失败姿势逐个拆解。

5.1 富集到的通路越“多”越好吗?警惕碎片化的显著条目

我第一次给学生看结果时,学生兴奋地告诉我:“老师,我富集了三百多条显著通路!”我当时的回复是:“恭喜你,但这并不一定是好事。”

富集到的通路数量多,可能说明你的处理效应强烈,但也可能说明你的筛选阈值太宽松,或者背景基因设置有问题。更麻烦的是条目碎片化——比如你看到“细胞对干扰素的反应”“干扰素介导的信号通路”“对I型干扰素的应答”同时出现在列表里,它们其实是高度相关的同一件事,只是GO术语把同一生物过程拆到了不同层级。这种情况下,如果只看表头按p值排序,你以为发现了三个不同的发现,实际上是同一个信号重复计数了三次。

我的建议是:先看top通路,再对富集条目做冗余归并。clusterProfiler里可以用simplify()函数基于语义相似度对GO结果去冗余,保留代表性的条目。KEGG分析里也可以用enrichplot::cnetplot()把基因-通路关系画成网络图,帮助你看清楚哪些通路本质上是同一批基因在支撑。

5.2 “显著”和“重要”是两码事:基因数、表达趋势和网络位置综合判断

这是我在无数场合强调过的点:统计显著的富集不等于生物学重要。一个p值只有1e-10但只富集了2个基因的通路,和一个p值0.001但富集了50个基因的通路,我通常更重视后者。原因很简单:50个基因共同参与一个通路,说明这个通路的整体活性发生了系统性的变化,这种变化才有更大概率驱动表型。

另外要看表达趋势。如果一个通路富集了30个基因,其中25个是一致的上调或者一致的下调,说明这个通路确实是整体激活或抑制;如果上下调基因在通路里参半分布,那这个“富集”可能反映的是通路内部的复杂性,而非简洁的激活/抑制关系。此时最好结合pathview的通路图,看关键节点到底是哪一侧的基因发生了变化。

还想提醒一点:不要把富集分析的结果当文献综述来用。有些人拿到富集结果后,喜欢把所有显著通路挨个在论文里介绍一遍,看起来内容丰富,实际上稀释了主线。好的文章通常只讲2到3个最核心的生物学故事,富集分析的真正作用就是帮你找出这2到3个故事,其余的通路可以放在补充材料里。

5.3 没富集到显著通路,是先看数据还是先改阈值?

这是个让很多人崩溃的场景:差异基因明明很多,但富集分析结果全部不显著。这时候怎么办?

我的排查顺序是:

  1. 检查背景基因(最优先)。你是不是把全基因组当背景,但实际检测到的基因只有几千个?这会让几乎所有通路的期望比例严重偏高,导致p值普遍不显著。换成实际检出基因做背景后,问题可能瞬间解决。
  2. 检查差异基因数目。如果差异基因只有二十几个,不显著是正常的——统计功效不够,神仙来了也救不了。这时候可以适当放宽差异筛选阈值(比如logFC从1降到0.5),先多纳入一些基因看趋势。
  3. 检查ID转换。基因注释丢失比例过高,会直接削弱富集的信号。
  4. 换用GSEA试试。ORA方法本质上只用了差异基因的名单,把表达量的强弱信息全部丢掉了。GSEA使用全基因组的表达排序做分析,往往比ORA更灵敏。这一步经常会带来意外之喜。

综合我的经验,ORA不显著时,先别急着调p值阈值。先修背景、再修ID、最后还是不行就换GSEA。调阈值是万不得已的下策,因为审稿人看到p<0.1这类阈值时,观感会很差。

5.4 GSEA、ORA和GSA的区别,以及什么时候该用哪种

最后简单讲一下富集分析的方法论家族,因为很多人被缩写搞晕了。

  • ORA(Over-Representation Analysis,过表达分析):就是前文讲的超几何分布方法,输入是差异基因列表,输出是哪些通路在这个列表里显著富集。优点是简单直观,缺点是丢失定量信息,且对阈值敏感。
  • GSEA(Gene Set Enrichment Analysis,基因集富集分析):输入是所有基因的表达量或统计量(如logFC、t值),不需要预先设定显著基因阈值。它回答的问题是:某个通路里的全部基因,整体上是倾向于在实验组高表达还是在对照组高表达。GSEA不需要武断地设置logFC阈值,统计功效更高。
  • GSA(Gene Set Analysis):是比GSEA更广的一个家族统称,涵盖了GSEA、ssGSEA、GSVA等一大堆方法。GSVA是一种不需要分组标签、可以按单个样本计算通路活性的方法,单细胞分析里用的尤其多。

如何选?如果你的差异基因不多,且你关心的是“这群基因在哪些通路集中”,用ORA;如果你的差异基因很多,或者你的处理效应是温和但广泛的那种(比如每个基因只变化20%,但一整条通路都在变),用GSEA;如果你做单细胞,想比较不同细胞亚群之间的通路状态,优先考虑GSVA或AUCell这一类单样本通路活性打分方法。

GSEA在clusterProfiler里跑也很方便,我这里给一个最小示例:

# geneList需要是从高到低排序的统计量向量,名称是Entrez ID geneList <- sort(logFC_vector, decreasing = TRUE) gsea_res <- gseKEGG( geneList = geneList, organism = "hsa", pvalueCutoff = 0.05, seed = 42 ) head(as.data.frame(gsea_res))

GSEA的结果里有个关键概念叫normalized enrichment score(NES),正负号代表通路在实验组中整体上调还是下调,绝对值大小反映富集强度。看GSEA结果时,我一般按NES绝对值排序,再结合p.adjust筛选,效果比只看p值排序合理得多。

6. 非模式生物、多组学数据和批量分析时,必须变通的那些细节

走到这一步,常规流程你已经能跑通了。但我再往下挖几层,聊聊那些真正影响你项目效率和结果质量的进阶场景。

6.1 非模式生物怎么破:Ortholog映射和通用注释包方案

模式生物的富集分析是很无脑的——人类的org.Hs.eg.db、小鼠的org.Mm.eg.db装了就能用。但真到了做斑马鱼、猪、牛、甚至一些冷门物种(比如我做过一段时间的蜜蜂)时,注释包就可能不存在了。

我踩过最深的坑是:某物种在NCBI里有基因注释,但Bioconductor上并没有对应的org包。这个时候有几个备选方案:

  1. 用Ortholog映射。先把你的物种基因序列通过序列比对或已有的同源映射工具(比如Ensembl BioMart的getLDS()函数、biomaRt包)映射到人类或小鼠的同源基因上,然后拿人类的注释包做富集分析,最后在结果解释时把基因ID翻译回去。这是目前处理非模式生物最通用的方案。
  2. 用go_enrichment类的在线工具。有些网页工具内置了较多物种的GO注释,但数据库更新和维护情况参差不齐,用之前要确认最近有没有更新。
  3. 自建注释。如果你所在的实验室长期做某个非模式物种,最彻底的办法是基于该物种的参考基因组注释文件(GTF/GFF)自己构建GO/KEGG注释映射表,然后用clusterProfiler的enricher()函数配合自定义的TERM2GENE格式数据去做分析。这个方法能让你彻底摆脱第三方注释包的束缚,下次再分析直接复用。代价是前期需要投入一两天时间做注释清洗。

不管用哪种方案,我都建议把“基因ID转换率”作为质控指标写进你的分析流程里。转换率低于60%时,不要继续往下分析,回头查数据才是正事。

6.2 多组学数据(转录组+蛋白组+磷酸化组)联合富集的思路

如果你手头同时有转录组和蛋白组的数据,想做一个联动分析,我推荐一个简单的思路:先分别做各自的富集分析,再做交集比较。重点不是看两个组学富集出多少相同的通路,而是看它们不同在哪里——

转录组富集的通路代表的是转录调控层面的变化,蛋白组富集代表的更接近功能蛋白层面的最终结果。如果某条通路在转录组显著富集但蛋白组没有,可能是转录后调控(比如microRNA、蛋白降解)在其中起作用;如果蛋白组显著富集但转录组没有,那就要考虑翻译效率调控或者蛋白稳定性变化。这种差异本身就是生物学故事的一个切入点。

代码层面的实现,主要是把两个组学各自的差异基因列表整理好,分别跑enrichKEGG,然后用compareCluster()函数合并:

# gene_lists是一个命名列表,比如 list(transcriptome = genes_rna, proteome = genes_prot) comp <- compareCluster( gene_lists, fun = "enrichKEGG", organism = "hsa", pvalueCutoff = 0.05 ) dotplot(comp, showCategory = 15)

这样一张图就能看出两个组学在通路上的一致性和差异性,做汇报或者写文章都很好用。

6.3 批量跑多个比较组时的效率管理:函数封装和结果汇总

最后分享一个实操层面的小技巧。现实项目里你很少只分析一组比较——比如一个包含3个时间点、2个处理组的设计,可能会有五六个比较组合。如果每次清洗数据、调参、画图都手动来一遍,既容易出错又浪费时间。我的做法是写一个包装函数,把差异基因列表加对比组名传进去,函数内部自动跑GO、KEGG、GSEA,并把结果导出成统一的Excel文件。

用clusterProfiler的对象结构可以很方便地做这一步:因为每个富集结果本质都是enrichResult对象,可以直接用as.data.frame()转成数据框,再聚合到一个总表里。给每行加上比较组标签之后,可以用dplyr::filter()快速筛选出多个比较组里都显著的共同通路——这种“核心通路”通常就是项目最需要强调的结论。

导出到Excel我习惯用writexl包,它不需要Java环境,比xlsx省心很多。输出文件可以按比较组分别建sheet,每个sheet里放GO和KEGG的结果,再加一个简短的参数记录sheet(比如筛选阈值、背景基因数),这样后续整理数据时不会留下说不清的参数漏洞。

我在实际项目里,光是把富集分析流程函数化封装这件事,就让组里新人每次跑新数据的耗时从一天压缩到了半小时。自动化不是为了让分析变得“玄学”,而是让你把精力花在真正的生物学解释上,而不是反复复制粘贴代码。

7. 写在最后:富集分析只是起点,不是终点

做富集分析这几年,我觉得它教会我最重要的一个思维方式是:实验数据不会主动开口说话,你得用对的问题去问它。富集分析就是那个提问的工具,工具本身不产生结论,结论来自你对生物学问题的理解和判断。

如果非得给刚入门的读者一条最值得记住的经验,我想说是这个:富集分析最大的坑从来都不是代码不会写,而是对原理理解不透彻导致的误用。背景基因设错了、ID格式转错了、读图时分不清“显著”和“重要”——这些才是让一堆分析结果变成废纸的真正原因。不论用哪个工具、跑多少条通路,先确保你手里的基因列表是干净的,背景是合理可辩护的,你才有底气去解释任何一条通路背后的故事。把上面的环节都做到位之后,剩下的就是尽情享受那个“豁然开朗”的时刻了。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询