做生信分析这几年,我见过太多人拿到差异基因列表后,直接丢进某个在线工具里点一下“富集”,然后截个气泡图就完事了。结果文章被审稿人一问“你这个背景基因集是什么”“为什么用这个数据库版本”,就答不上来。基因的富集分析表面上看就是一次超几何分布检验或者一个GSEA跑分,但实际用起来,里面全是细节。这篇东西我想把富集分析从原理到实操完整捋一遍,重点讲清楚ORA和GSEA这两条主流路径各自解决什么问题、参数怎么选、结果怎么解读,以及我踩过的那些坑。无论是刚接触生信的湿实验同学,还是想系统梳理富集分析流程的分析人员,这篇文章都能给你一套可以直接上手的思路。
1. 富集分析解决的到底是什么问题
1.1 从“差异基因列表”到“生物学解释”的鸿沟
测序做完,差异表达分析跑完,你会拿到一张表格,里面有几百上千个上调或者下调的基因。这些基因名字你大部分都认识,但你很难回答一个最简单的问题:这些基因整体上在干什么?
单个基因看,A基因跟免疫相关,B基因跟增殖相关,C基因似乎参与代谢。你可不能因为看到几个眼熟的基因就下结论说这项研究的机制是免疫加增殖加代谢,这是典型的“挑着看数据”,审稿人一眼就能识破。
富集分析要解决的,就是从“散装基因列表”到“系统性生物学结论”这一步。它做的事情本质上是:拿你的差异基因列表去跟一个个已经注释好的基因集合(比如某条代谢通路、某个生物学过程)做比对,看看你列表里的基因是否在某个已知功能集合中显著富集。如果某个通路里的基因在你列表里扎堆出现,那么即使单个基因的fold change并不惊人,你也有理由相信这条通路确实在生物学过程中扮演了重要角色。
这也是为什么我强调富集分析不是“跑一下就有结果”那么简单。你的基因列表质量、你选的注释数据库、你设定的统计检验参数,每一步都在影响最终结论。后面你会看到,同样的差异基因,用不同的背景基因集跑GO分析,结果可能完全不同。
1.2 富集分析的核心逻辑与统计思想
先不急着上代码,我用一个极其朴素的例子把富集分析背后的统计思想讲明白。
想象你面前有10000个球(代表全部基因,这里理解为你的背景基因集),其中100个是红球(代表属于“细胞增殖”这个功能类别的基因)。你随机抓了200个球(代表你的差异基因列表),发现里面有50个红球。按照背景比例,200个球里你期望只看到2个红球(100/10000×200),结果你看到了50个。这明显高于随机预期,于是统计学告诉你:“红球在这个抓取结果里是显著富集的”。
上述思想落到实际算法上,常见的有两类。一类叫过表达分析(ORA,Over-Representation Analysis),问的是“我列表里的基因在某个功能集合里是不是过多”,用的统计模型是超几何分布或Fisher精确检验;另一类叫功能类别打分(FCS,Functional Class Scoring),其中最典型的是GSEA,它不需要你事先圈定基因列表,而是看全部基因的表达排序后,某个功能集合的基因是否整体偏向排序的一端。
这两类思路差异巨大。ORA的输入是“阈值切出来的基因列表”,GSEA的输入是“全部基因的排序信息”。很多初学者混淆这两者,结果该用GSEA的场景用了ORA,白白丢掉信息。后面我分别展开讲。
2. ORA过表达分析:最经典的富集策略
2.1 ORA的基本流程与统计公式
ORA是最早出现、也是目前在线工具用得最多的富集分析方法。它的流程很清晰,四步:
- 确定目标基因集(比如差异表达基因列表,up或down分开)。
- 确定背景基因集(这个极其关键,稍后单独说)。
- 获取基因与功能注释的对应关系(GO、KEGG、Reactome等)。
- 对每一个功能条目做统计检验,判断目标基因集中属于该条目的基因数量是否显著高于随机期望。
统计检验的核心是超几何分布。假设背景基因集规模为N,其中属于某功能条目K的基因数为M;目标基因集规模为n,其中恰好有k个基因落在该功能条目里。那么“随机情况下看到k个或更多基因落入该条目”的概率由超几何分布给出:
P(X ≥ k) = Σ (C(M, i) × C(N-M, n-i)) / C(N, n),i从k到min(M,n)。
实际操作中,很多工具会用Fisher精确检验来等效计算这个概率。两者数学上等价或近似,你不需要手动算,但你要理解它在做什么——它在算的是“富集到这个程度纯靠运气发生的可能性”。
算完P值还不够。你同时检验了几千个GO条目,每个条目都有一个P值,这就带来了多重检验问题。我一向主张用BH(Benjamini-Hochberg)方法校正得到FDR(False Discovery Rate),选显著阈值时看padj或q-value小于0.05。只报原始P值不校正,在几千次检验里必然有假阳性,审稿人必定会问。
2.2 实操中三个最常见的坑
ORA看起来简单,但我审阅过大量结果文件,最常见的错误集中在三处。
第一个坑是背景基因集设错。很多人做差异基因富集分析时,背景基因集随便选了“全部注释基因”。但你的差异基因是从“检测到的基因”里筛出来的,RNA-seq实际检测到的基因可能只有两万个中的一万八千个,被表达量过低过滤掉的那两千个基因根本没有机会进入差异列表。假如你用全部两万个基因做背景,那相当于把一批不可能出现的球也算进了池子,会低估或高估富集显著性。规范做法是:背景基因集应为“表达矩阵里实际参与差异检验的基因全集”。这一点几乎每个工具里都有参数可以改,但很多人从没注意过。
第二个坑是基因ID体系不统一。Ensembl ID、Entrez ID、Gene Symbol在工具里混着用,或者转换的时候掉了一大批基因。比如你在差异表里用的是Symbol,注释库用的是Entrez,直接跑GSEA时就会提示匹配率过低。正确做法是分析前统一ID类型,并在结果里报告匹配率。我给自己定的及格线是匹配率不低于80%,低于这个数就要回头查转换过程。
第三个坑是数据库版本。GO注释每个月都在更新,KEGG通路也会调整,不同版本跑同一份数据结果会有出入。文章里必须写明用的是哪个数据库的哪个版本、什么时间下载的。这不仅是规范问题,也直接关系到结果能否被复现。我一般在方法部分会写类似“GO enrichment was performed using clusterProfiler (v4.6.0) with org.Hs.eg.db (v3.16)”这样的字样,读者跟着跑一遍就能还原。
3. GSEA富集分析:如何看功能趋势而不是单个基因
3.1 GSEA做了哪些ORA做不到的事情
ORA有一个天生的局限:它要求你先把基因列表用阈值切成“显著”和“不显著”两组。这个切法的问题在于,生物体内很多通路的改变是细微而协调的——单个基因的表达变化可能都达不到差异显著的P值门槛,但整条通路的基因都朝同一方向发生了小幅变化。这种“趋势性”的信号,ORA完全看不见。
GSEA全称Gene Set Enrichment Analysis,核心思想是抛弃阈值,保留排序。你把所有基因按照某种指标(比如log2 fold change,或signal-to-noise ratio)从高到低排一列,再去检查事先定义好的每个基因集中所有基因在这个排序里是否均匀分布。如果某个基因集的成员整体偏向排序的顶端(上调端)或底端(下调端),就说这个基因集在这个条件下被富集了。
打个比方,ORA是看“你选的50个人里有没有10个都来自某个公司”,GSEA是看“整个会场入场时,某公司的人是不是特别早到或者特别晚到”。前者在乎命中比例,后者在乎整体分布趋势。
GSEA输出里有两个核心统计量需要理解。一个是富集分数(ES,Enrichment Score),它衡量基因集成员在排序中的聚集程度,实际上是一个加权Kolmogorov-Smirnov-like统计量;另一个是归一化富集分数(NES,Normalized Enrichment Score),它把ES按照基因集大小做了归一化,方便不同大小的基因集之间比较。显著性的判据是用排列检验得到的FDR q-value,一般取小于0.25作为阈值,这比ORA的0.05宽松,原因是GSEA的检验本身更保守,而且通常用于发现趋势而不是确定单点结论。
3.2 GSEA实操:排序文件怎么构建,参数怎么设
跑GSEA前,最关键的准备工作是构建基因排序文件(.rnk文件)。每一行是两个字段:基因ID和排序得分。排序得分用什么值,直接影响分析结果。
我见过有人直接拿P值取负对数当排序分,这不是不可以,但会扭曲生物学含义——P值只反映统计显著性,不反映变化方向和幅度。我更推荐的做法是:有生物学重复时用signal-to-noise ratio (S2N)或limma的t值,没有重复或者只想看表达量变化方向时,直接用log2 fold change。注意,用log2FC排序时一定要保证上调基因在前面(正值),下调在后面(负值),别把方向搞反。
接下来是GSEA运行时的关键参数,我逐个说。
- 排列次数(permutations):建议至少1000次,少于这个数,最小可达P值都会被限制在0.001左右,这对多重检验校正非常不利。样本量很小的时候可以适当增加,但计算时间会变长。
- 加权指数(weight exponent, p):默认值是1,意思是基因在排序中位置越靠前,其对ES的贡献用该基因排序得分的p次方加权。想突出高排位基因的作用可以设p=2,想做无加权版本设p=0。默认1在绝大多数场景下表现稳定,没有特殊理由不要改。
- 基因集大小过滤:一般过滤掉少于15个和多于500个成员的基因集。太小的基因集统计不稳定,太大的基因集太宽泛,结论没有针对性。
- 基因集数据库选择:做GSEA常用MSigDB,里面有H(hallmark gene sets)以及C2(curated gene sets,含KEGG、Reactome等)等几个大类。我的经验是,第一轮看hallmark结果快筛方向,锁定整体生物学主题后再去医院化的C2子集深入。
运行完你会得到一个富集结果表格和一堆可视化图。最值得仔细看的是running enrichment score图(就是那个爬坡状的折线图)、热图以及leading edge子集。leading edge指的是基因集中真正推动富集分数的那部分基因,也就是出现在最大ES峰值之前的成员,它们才是你后续做机制研究要优先关注的候选基因。
4. 实操细节:注释库、工具链与可视化
4.1 常用注释数据库盘点
富集分析的结果质量上限由注释数据库决定。程序再花哨,数据库烂就得不出好结论。我常用的数据库有以下几类,按使用频率排个序。
GO(Gene Ontology)是覆盖面最广的功能注释体系,分三个子本体:生物学过程(BP)、细胞组分(CC)、分子功能(MF)。做富集分析时三个子本体要分开跑,否则混合在一起结果非常难解读。实践中BP条目通常最受关注,但冗余度高,聚类后看会清晰很多;CC和MF的信号相对集中,经常能提示你关注亚细胞定位或具体分子活性。不同物种要用对应的org包,人用org.Hs.eg.db,小鼠org.Mm.eg.db,大鼠org.Rn.eg.db,其他物种可以试试AnnotationHub或biomaRt在线注释。
KEGG通路是最常被引用的通路数据库,胜在“通路图”直观、审稿人熟悉。但KEGG数据有版权限制,有些工具已经不再更新或采用通过API获取的策略。最新的KEGG富集分析最好用clusterProfiler配合KEGG REST API做在线查询,用之前检查网络连通性。
Reactome是一个人工注释的通路数据库,结构层次比KEGG更细,覆盖也广,尤其适合信号通路相关研究。MSigDB则是GSEA官方推荐的基因集来源,不只包含通路,还包括各种带生物学主题的基因集(如特定细胞类型标志基因、癌症相关signature等),做GSEA时能打开很多新视角。
4.2 从差异基因列表到富集结果的完整命令链
具体到代码实现,我日常主力是R语言的clusterProfiler包。下面这段是我做ORA的标准流程,可以直接复制改路径用。
library(clusterProfiler) library(org.Hs.eg.db) # 假设你的差异基因列表是数据框deg,其中包含基因列和log2FC列 deg <- read.csv("deg_results.csv", stringsAsFactors = FALSE) # 统一ID:用Symbol转Entrez,注意去掉版本号 gene_symbols <- bitr(deg$gene, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db) merged <- merge(deg, gene_symbols, by.x = "gene", by.y = "SYMBOL") # 上面bitr没匹配到的基因会被丢弃,建议统计一下匹配率 cat("matched:", nrow(merged), "total:", nrow(deg), "\n") # 分离上调和下调基因(阈值自己定,通常log2FC绝对值>1, padj<0.05) up_genes <- merged$ENTREZID[merged$log2FoldChange > 1 & merged$padj < 0.05] down_genes <- merged$ENTREZID[merged$log2FoldChange < -1 & merged$padj < 0.05] # 关键:background必须是所有参与差异检验的基因的Entrez ID background <- bitr(deg$gene, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)$ENTREZID # GO富集分析(BP子本体) ego_up <- enrichGO(gene = up_genes, universe = background, OrgDb = org.Hs.eg.db, ont = "BP", keyType = "ENTREZID", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.2, readable = TRUE) head(ego_up@result) # KEGG富集分析 ekegg_up <- enrichKEGG(gene = up_genes, universe = background, organism = "hsa", keyType = "kegg", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.2) head(ekegg_up@result)这段代码里有一个我特别想强调的点:universe参数。很多人enrichGO时漏掉它,默认会用OrgDb里所有基因当背景,结果常常高估显著性。我在2.2里埋的坑就在这里暴露——把universe设成实际参与检验的基因集,富集结果才真正可信。
输出结果的解读,我一般看三列:GeneRatio(目标基因中命中该条目的比例)、BgRatio(背景基因中属于该条目的比例)、p.adjust(校正后的显著性)。如果GeneRatio很高但BgRatio也很高,说明这个条目本来就大而全,富集到它不值得大写特写;反之,BgRatio很低但命中集中在你的列表里,这个条目才是真正有故事可讲的。
GSEA的R实现,我推荐用clusterProfiler::gseGO或fgsea包。fgsea在大基因集数量下跑得极快,而且内存占用小,强烈推荐。
library(fgsea) library(msigdbr) # 构建排序向量:命名为Entrez ID,值用log2FC或S2N,从高到低排序 ranks <- setNames(merged$stat, merged$ENTREZID) # stat可以是log2FC或t值 ranks <- sort(ranks, decreasing = TRUE) # 用msigdbr获取基因集,比如hallmark hallmark <- msigdbr(species = "Homo sapiens", category = "H") hallmark_list <- split(hallmark$entrez_gene, hallmark$gs_name) # 跑fgsea set.seed(42) fgsea_res <- fgsea(pathways = hallmark_list, stats = ranks, minSize = 15, maxSize = 500, nPerm = 10000) # 按padj排序,top结果 head(fgsea_res[order(padj), .(pathway, pval, padj, ES, NES)])看到这里你会发现,GSEA里你必须一开始就给每个基因算好一个排序得分。这个得分可以是log2FC、limma的 moderated t 值、DESeq2 的 stat 值等,不管用哪个,命名保持一致,下游分析才不混乱。
4.3 可视化:如何让富集结果“讲得清楚”
结果再好,图不好看也白搭。我按使用频率排几类可视化形式。
气泡图是最通用的展示方式:横轴是GeneRatio或Count,纵轴是通路名,点的大小代表命中基因数,颜色代表P值或Q值。clusterProfiler里的dotplot()一行代码就能出图。这张图我主要用于方法部分展示整体富集谱。
通路图适合KEGG结果:pathview可以把差异基因的表达值映射到KEGG通路图上,颜色深浅直接显示基因上下调。审稿人看到这种图会觉得你确实把机制落在了通路上,而不是停留在条目列表层面。
GO有向无环图(DAG)用于展示GO条目之间的父子关系,适合表现BP过程的层级结构,enrichplot::goplot()可以直接出。注意DAG图信息密集,放文章大图里有点考验排版,我一般在补充材料里用。
GSEA的核心图是running score图:横轴是排序后的基因位置,纵轴是累积富集分数,峰值位置对应leading edge基因集。enrichplot::gseaplot2()可以一次展示running score、基因排序位置以及底部热图三个面板,是GSEA文章里最标配的一幅图。
5. 常见问题与排查技巧实录
表格直接给结论,这几类是我在实操中被问过最多、或自己踩过的坑。
| 现象 | 可能原因 | 排查与解决 |
|---|---|---|
| 富集结果为空或极少条目 | 差异基因太少;背景基因集错误过大;ID匹配失败 | 检查匹配率和差异基因数;把背景换成实际检测基因集;放宽pvalueCutoff暂时看趋势 |
| 富集到的条目全都是细胞组分(CC) | GO三个本体混跑;BP信号弱 | 分开跑BP/CC/MF,重点看BP;如果BP确实无显著条目,接受结果并如实报告 |
| P值接近1,只有很少条目显著 | 背景基因集设置过小;目标基因列表太小;缺失负对照 | 核对universe参数;确认目标基因列表是否真来自这个背景的检测集 |
| GSEA结果显著条目过多(几十上百个) | 基因集重叠度高;排序指标噪声大 | 换hallmark集合先看全局;用NES排序取前几个代表性结果;考虑聚类合并同源通路 |
| GSEA全员不显著,ES也很低 | 排序指标里信号太弱;基因集版本与物种不匹配 | 尝试用limma的t值代替log2FC;检查基因集数据库是否对应物种 |
| 同一通路在不同数据库结论相反 | 数据库注释标准不一致;阈值选择不同 | 不要强求统一;在文章中指出不同层级的证据,注释清楚各数据库版本 |
再补一个我个人的经验性技巧:拿到富集结果后,先看富集到的条目里有没有“冗余扎堆”现象。如果看到一堆条目讲的是同一件事(比如免疫应答的各个子过程全部显著),这通常说明生物学信号真实且强,处理时可以用REVIGO或clusterProfiler的simplify()去掉冗余,挑代表性子集讲。相反,如果显著条目东一个西一个互不关联,那更可能是数据分析流程里的某一步引入了噪声,这时我不急着写结论,而是回头检查差异基因质量、样本关系和批次效应。
另外一个非常实用但容易被忽略的习惯是保存完整session信息。我会在每次富集分析结束后记录R版本、各包版本、数据库下载时间,甚至随机种子。GSEA的排列检验有随机性,设置种子才能让结果严格可复现。sessionInfo()的输出虽不起眼,但真要回复审稿人或补实验,它的价值就体现出来了。
最后说一句掏心窝的话。基因的富集分析不是终点,它本质上是帮你把海量基因列表压缩成几个可验证的生物学假说的过滤器。把ORA和GSEA的原理吃透,把背景基因集和数据库版本这两个最容易被忽略的参数管好,你的富集分析结果就会扎实很多。做分析时少一点“跑完就出图”的惯性,多一点“这个统计量到底在算什么”的追问,这个习惯会让你少走太多弯路。