富集分析对比:GSEA/ORA/SPIA
一句话概述:ORA(过代表分析)用差异基因列表做富集,GSEA(基因集富集分析)用所有基因的排序做富集不需要人为截断,SPIA(信号通路影响分析)还额外考虑通路拓扑结构——三者回答不同层次的生物学问题。
核心知识点速查表
| 概念 | 说明 |
|---|
| ORA | 过代表分析(白话:我的差异基因在某个通路里是不是特别多?) |
| GSEA | 基因集富集分析(白话:某个通路的基因是不是整体偏高/偏低?) |
| SPIA | 信号通路影响分析(白话:考虑基因在通路中的位置和关系) |
| ssGSEA | 单样本GSEA(给每个样本算一个通路活性分数) |
| fgsea | 快速GSEA的R实现(比原版GSEA快100倍) |
| 基因集 | 一组功能相关的基因(如GO term、KEGG pathway) |
一、三种方法对比
| 特性 | ORA | GSEA | SPIA |
|---|
| 输入 | 差异基因列表 | 所有基因+排序指标 | 差异基因+logFC |
| 需要阈值? | 是(要先筛DEG) | 否(★优势) | 是 |
| 统计方法 | 超几何/Fisher | 排列检验 | 拓扑+ORA |
| 考虑通路结构 | 否 | 否 | 是(★优势) |
| 考虑基因数值 | 否(只用有/无) | 是(用logFC排序) | 是(用logFC) |
| 速度 | 最快 | 中等 | 慢 |
| 常用工具 | clusterProfiler | fgsea/GSEA | SPIA包 |
| 适用场景 | 快速初筛 | 全面分析(★推荐) | 通路层面深入 |
二、R语言实操
2.1 ORA(过代表分析)
# === ORA分析 ===
library(clusterProfiler)
library(org.Hs.eg.db) # 人类注释包
# 准备差异基因列表(需要Entrez ID)
deg_genes <- bitr(
rownames(sig_genes), # 差异基因名
fromType = "SYMBOL", # 从基因符号
toType = "ENTREZID", # 转为Entrez ID
OrgDb = org.Hs.eg.db # 注释数据库
)
# GO富集(ORA)
go_ora <- enrichGO(
gene = deg_genes$ENTREZID, # 差异基因
OrgDb = org.Hs.eg.db, # 注释数据库
ont = "BP", # BP/CC/MF
pAdjustMethod = "BH", # 多重检验校正
pvalueCutoff = 0.05, # p值阈值
qvalueCutoff = 0.05, # q值阈值
readable = TRUE # 基因ID转为Symbol
)
# KEGG富集(ORA)
kegg_ora <- enrichKEGG(
gene = deg_genes$ENTREZID, # 差异基因
organism = "hsa", # 物种(人)
pvalueCutoff = 0.05
)
# 可视化
dotplot(go_ora, showCategory=15) # 气泡图
barplot(go_ora, showCategory=15) # 柱状图
cnetplot(go_ora, showCategory=5) # 基因-通路网络图
2.2 GSEA(基因集富集分析)
# === GSEA分析 ===
library(clusterProfiler)
library(fgsea)
# 准备排序基因列表(★不需要筛选差异基因)
# 用所有基因的logFC或-log10(p)*sign(logFC)排序
gene_list <- res$log2FoldChange # DESeq2的logFC
names(gene_list) <- rownames(res) # 基因名
gene_list <- sort(gene_list, decreasing=TRUE) # 降序排列
gene_list <- gene_list[!is.na(gene_list)] # 去掉NA
# 转换为Entrez ID
gene_df <- bitr(names(gene_list), "SYMBOL", "ENTREZID", org.Hs.eg.db)
gene_list_entrez <- gene_list[gene_df$SYMBOL]
names(gene_list_entrez) <- gene_df$ENTREZID
# GO GSEA
gsea_go <- gseGO(
geneList = gene_list_entrez, # 排序的基因列表
OrgDb = org.Hs.eg.db, # 注释数据库
ont = "BP", # GO类别
minGSSize = 10, # 最小基因集大小
maxGSSize = 500, # 最大基因集大小
pvalueCutoff = 0.05, # p值阈值
verbose = FALSE
)
# KEGG GSEA
gsea_kegg <- gseKEGG(
geneList = gene_list_entrez,
organism = "hsa",
minGSSize = 10,
pvalueCutoff = 0.05
)
# GSEA可视化
gseaplot2(gsea_go, geneSetID=1:3, # 富集得分图
pvalue_table=TRUE) # 显示p值表
ridgeplot(gsea_go, showCategory=15) # 山脊图
dotplot(gsea_go, showCategory=15) # 气泡图
# === 使用fgsea(更快) ===
library(fgsea)
library(msigdbr)
# 获取基因集(MSigDB)
pathways <- msigdbr(species="Homo sapiens", category="H") # Hallmark基因集
pathways_list <- split(pathways$entrez_gene, pathways$gs_name)
# 运行fgsea
fgsea_res <- fgsea(
pathways = pathways_list, # 基因集列表
stats = gene_list_entrez, # 排序的基因列表
minSize = 15, # 最小基因集大小
maxSize = 500 # 最大基因集大小
)
# 查看结果
fgsea_res[order(pval), ][1:10, ] # 前10个显著通路
2.3 SPIA(通路影响分析)
# === SPIA分析 ===
library(SPIA)
# 准备输入
sig_entrez <- deg_genes$ENTREZID # 差异基因Entrez ID
all_entrez <- bitr(rownames(res), "SYMBOL", "ENTREZID", org.Hs.eg.db)$ENTREZID
# 需要logFC向量,名字为Entrez ID
de_logfc <- res$log2FoldChange[res$padj < 0.05 & !is.na(res$padj)]
names(de_logfc) <- bitr(
rownames(res)[res$padj < 0.05 & !is.na(res$padj)],
"SYMBOL", "ENTREZID", org.Hs.eg.db
)$ENTREZID
# 运行SPIA
spia_res <- spia(
de = de_logfc, # 差异基因logFC(名字为Entrez ID)
all = all_entrez, # 所有检测到的基因
organism = "hsa", # 物种
nB = 2000 # 排列次数
)
# 查看结果
head(spia_res)
# pSize: 通路中基因数
# NDE: 通路中差异基因数
# pNDE: ORA p值(基因数量)
# pPERT: 扰动p值(通路拓扑)
# pG: 全局p值(综合)
# Status: Activated/Inhibited
plotP(spia_res) # 二维证据图
三、ssGSEA(单样本GSEA)
# === ssGSEA:每个样本一个通路活性分数 ===
library(GSVA)
# 准备表达矩阵(需要标准化后的)
expr_matrix <- assay(vst(dds)) # VST标准化
# 获取基因集
pathways <- msigdbr(species="Homo sapiens", category="H")
pathways_list <- split(pathways$gene_symbol, pathways$gs_name)
# 运行ssGSEA
ssgsea_scores <- gsva(
expr_matrix, # 表达矩阵(基因×样本)
pathways_list, # 基因集列表
method = "ssgsea", # ssGSEA方法
kcdf = "Gaussian" # 核密度函数
)
# 输出: 通路×样本的活性分数矩阵
# 可视化
pheatmap(ssgsea_scores,
annotation_col = sample_info["condition"],
scale = "row")
四、面试高频考点
Q1: ORA和GSEA的本质区别?
- ORA:你先画一条线(如padj<0.05)把基因分成"差异"和"非差异"两堆,然后看差异基因里某个通路的基因是不是比随机多
- GSEA:不画线,把所有基因按表达变化排队,看某个通路的基因是集中在队伍前面(上调)还是后面(下调)
- 白话:ORA问"考上清华的人里学数学的多不多?",GSEA问"学数学的人在全校排名中整体偏前还是偏后?"
Q2: GSEA的NES是什么?
- NES = Normalized Enrichment Score(标准化富集分数)
- NES > 0:通路在处理组中上调
- NES < 0:通路在处理组中下调
- |NES| > 1.5 且 FDR < 0.25:通常认为显著
Q3: 什么时候用SPIA?
- 当你关心通路的拓扑结构时(如上游基因vs下游基因的影响不同)
- 信号通路中基因有明确的上下游关系
- 白话:ORA只看你家有多少人感冒了,SPIA还看谁先传染给谁的
五、常见报错
| 报错 | 原因 | 解决 |
|---|
No gene can be mapped | 基因ID格式不对 | 用bitr转换ID格式 |
GSEA: gene list has no valid entries | 排序列表有问题 | 检查是否有NA/重复 |
fgsea: too few genes | 基因集太小 | 降低minSize |
SPIA: requires KEGG data | 缺少通路数据 | 先安装KEGG数据库 |
速查表
# === 富集分析速查 ===
library(clusterProfiler)
# ORA(差异基因列表→富集)
enrichGO(gene=entrez_ids, OrgDb=org.Hs.eg.db, ont="BP")
enrichKEGG(gene=entrez_ids, organism="hsa")
# GSEA(全部基因排序→富集,★推荐)
gene_list <- sort(logFC, decreasing=TRUE)
gseGO(geneList=gene_list, OrgDb=org.Hs.eg.db, ont="BP")
gseKEGG(geneList=gene_list, organism="hsa")
# fgsea(快速GSEA)
fgsea(pathways=pathway_list, stats=gene_list)
# SPIA(通路拓扑)
spia(de=de_logfc, all=all_genes, organism="hsa")
# ssGSEA(单样本打分)
gsva(expr, pathways, method="ssgsea")
# 选择: 快速初筛→ORA | 全面分析→GSEA | 通路拓扑→SPIA | 样本打分→ssGSEA