Gene set enrichment analysis

何时使用

当用户已有一组基因或带打分的排序基因表，想知道「哪些生物学被富集」时使用。这是差异表达、筛选实验或聚类后的标准收尾步骤。典型场景：

• 在基因列表中查找富集的 GO / KEGG / Reactome / WikiPathways / MSigDB Hallmark 条目。
• 对 DESeq2、edgeR、limma 或 Scanpy rank_genes_groups 的输出跑 GSEA / 预排序 GSEA。
• 按样本/单细胞给通路活性打分（ssGSEA、GSVA）。
• 对富集结果做解读、去冗余、可视化，或产出发表级表格/图。
• 在 ORA 与 GSEA 之间抉择、选基因集库、定背景、修 gene-ID 问题。

不该用的边界：

• 上游差异表达本身的计算（属于 pydeseq2 / scanpy 等技能）。
• 只想做一次性 Enrichr 快查 —— gget enrichr 更轻量。
• 只想抓取原始通路/互作 API（Reactome、KEGG、STRING）—— 用 database-lookup 类技能。
• 本技能聚焦完整、可辩护、可发表的富集工作流。

步骤

1. 锁定输入并选方法（最关键的一步）

确认：哪些基因、什么物种、有没有每基因打分、代表什么比较（方向影响解读）。两种核心方法二选一：

• ORA（过表征分析）：拿已阈值过滤的基因列表（如 padj < 0.05），用 Fisher 精确/超几何检验看它与哪些基因集重叠超过随机。工具：Enrichr、g:Profiler。
• GSEA（基因集富集分析）：拿完整排序列表（不设阈值），检验每个基因集是否集中在排序的顶端或底端。预排序 GSEA 用每基因打分（如 DESeq2 的 stat）。当效应广泛而微弱时更优。

选型对照：

情形	方法	入口
离散命中列表（DE 基因、筛选命中、聚类标志）	ORA	gp.enrichr(...) 或 g:Profiler
完整排序列表（每个检测基因 + 打分）	预排序 GSEA	gp.prerank(...)
表达矩阵 + 类别标签	GSEA	gp.gsea(...)
每样本/单细胞通路打分	ssGSEA / GSVA	gp.ssgsea(...)、gp.gsva(...)
自定义背景或 500+ 物种	带自定义域的 ORA	g:Profiler（domain_scope='custom'）

经验法则：阈值列表 → ORA；带打分的排序表 → GSEA。绝不要先阈值过滤再喂给 GSEA，那会丢掉 GSEA 赖以工作的排序。

2. 把 gene-ID 转到正确命名空间

Enrichr/MSigDB 库以基因符号为键（人类大写，小鼠首字母大写）。若手头是 Ensembl/Entrez ID，先转换（gp.Biomart、g:Profiler g:Convert、mygene）。ID 静默不匹配是「什么都不显著」的头号原因。

3. 选匹配问题的基因集库

Hallmark（宽主题）→ GO:BP（机制）→ KEGG/Reactome/WikiPathways（精选通路）→ C7（免疫）等。别一口气跑 50 个库，挑 2–4 个贴合生物学的即可。库名会随时间变化，别盲目硬编码，先列出再选。

4. 设定背景全集（仅 ORA）

背景必须是你的实验本可检测到的基因（如所有表达/检测到的基因），而非整个基因组。背景错了会虚高显著性。Enrichr 用固定背景；当背景要紧时，用 g:Profiler 的 domain_scope='custom' + 你的 background，或 gp.enrich() 显式给背景。

5. 运行分析

用下方示例或封装好的 scripts/run_enrichment.py。GSEA 务必设 seed 并报告 permutation_num。

6. 按校正后 p 值过滤

用 Adjusted P-value（ORA，Benjamini–Hochberg）或 FDR q-val（GSEA），不是原始 p 值。常用阈值 0.05；同时看重叠/基因数，避免「命中」只是 2000 基因集里的 1 个。

7. 可视化

点图、条形图、富集图（enrichment map）、GSEA running-score 图，gseapy 内置：gp.dotplot、gp.barplot、gp.enrichment_map、gp.gseaplot。

8. 去冗余与解读

GO 尤其会返回大量近重复条目。用富集图（条目–条目相似度）、leading-edge 重叠或父条目折叠，只报代表性条目。

指令

安装：

uv pip install gseapy gprofiler-official
# gseapy 会带入 pandas/numpy/scipy/matplotlib。Enrichr、g:Profiler、MSigDB 下载需要联网。
# 完全离线 ORA 用本地 GMT 文件配 gp.enrich()。

核对并列出可用基因集库（库名会漂移，别硬编码）：

import gseapy as gp
names = gp.get_library_name(organism="human")   # 200+ Enrichr 库
print([n for n in names if "Reactome" in n or "KEGG" in n or "Hallmark" in n])

封装脚本（自动处理符号清洗、去重、去 NA、从 DESeq2 表构造排序、逐库 FDR 过滤）：

# 从命中列表跑 ORA（每行一个基因符号）
python scripts/run_enrichment.py ora \
  --genes deg_symbols.txt \
  --libraries MSigDB_Hallmark_2020 GO_Biological_Process_2023 KEGG_2021_Human \
  --organism human --outdir results/

# 从 DESeq2 结果 CSV 跑预排序 GSEA（自动用 stat 构造排序）
python scripts/run_enrichment.py gsea \
  --deseq2 deseq2_results.csv \
  --libraries MSigDB_Hallmark_2020 GO_Biological_Process_2023 \
  --organism human --outdir results/ --seed 123

# 从显式两列排序文件（gene,score）跑预排序 GSEA
python scripts/run_enrichment.py gsea --rnk ranked_genes.csv --outdir results/

--help 查看全部选项（背景文件、FDR 阈值、min/max 集合大小、置换次数）。

示例

ORA：对命中列表（gseapy + Enrichr）

import gseapy as gp

# Enrichr 库要 HGNC 基因符号（人类大写）。需要时先映射 ID。
genes = [g.strip() for g in open("deg_symbols.txt") if g.strip()]

enr = gp.enrichr(
    gene_list=genes,
    gene_sets=["MSigDB_Hallmark_2020", "GO_Biological_Process_2023",
               "KEGG_2021_Human", "Reactome_2022"],
    organism="human",
    outdir=None,            # 内存模式；给路径则同时写表/图
)
res = enr.results
sig = res[res["Adjusted P-value"] < 0.05].sort_values("Adjusted P-value")
print(sig[["Gene_set", "Term", "Overlap", "Adjusted P-value", "Combined Score", "Genes"]].head(20))

预排序 GSEA：从 DESeq2 结果

import gseapy as gp
import pandas as pd

res = pd.read_csv("deseq2_results.csv", index_col=0)   # index = 基因符号
# 按检验统计量排序（符号=方向，幅度=证据）。比按 log2FoldChange 更稳，
# 后者对低 count 基因噪声大。
rnk = res["stat"].dropna().sort_values(ascending=False)
rnk.index = rnk.index.str.upper()
rnk = rnk[~rnk.index.duplicated(keep="first")]

pre = gp.prerank(
    rnk=rnk,
    gene_sets=["MSigDB_Hallmark_2020", "GO_Biological_Process_2023"],
    min_size=15, max_size=500,        # 丢掉过小/过大的集合（噪声或太泛）
    permutation_num=1000, seed=123,   # seed = 可复现的 p 值
    threads=4, outdir=None,
)
out = pre.res2d.sort_values("FDR q-val")
print(out[["Term", "ES", "NES", "NOM p-val", "FDR q-val", "Lead_genes"]].head(20))

若没有 stat 列，用 sign(log2FoldChange) * -log10(pvalue) 构造排序。

注意事项

以下问题导致大多数错误或不可复现的结果：

1. gene-ID / 物种不匹配 —— 符号 vs Ensembl、人类 vs 小鼠大小写。映射好 ID 并正确设 organism，否则匹配会静默掉到约零。
2. 背景错误（ORA） —— 用整个基因组而非检测到/表达的基因集，会虚高 p 值。必要时设自定义背景。
3. GSEA 前做了阈值过滤 —— GSEA 需要完整排序列表；只有 ORA 用切过的列表。
4. 仅按 log2FoldChange 给 GSEA 排序 —— 对低 count 基因不稳；优先 stat 或 sign(LFC) * -log10(p)。
5. 跨库多重检验 —— FDR 是在单库内计算的；跑多个库会成倍增加检验。报告逐库 FDR 并保守对待。
6. 冗余 GO 条目 —— 别报同一条目的 40 个变体；折叠并展示代表。
7. 显著 ≠ 相关 —— 检查重叠数和基因集大小；小集合轻易就显著。
8. ORA 列表过短/过长 —— <10 基因功效不足；>2000 失去特异性（改用 GSEA）。
9. 缺复现元数据 —— Enrichr/GO 库有版本且会随时间漂移。记录库名+日期，并给 GSEA 设 seed。

互见

• 上游（基因来源）：差异表达（pydeseq2 的 DE 基因 + GSEA 用的 stat）、单细胞标志基因（scanpy rank_genes_groups）、筛选命中、蛋白质组命中。
• 数据库 / ID：Reactome、KEGG、STRING、Gene Ontology API；gget enrichr 快查、ID 映射。
• 下游：科学可视化（自定义图）、网络图（富集图）、科学写作/文献综述（解读与引用）、统计分析（多重检验细节）。

资源：gseapy 文档 https://gseapy.readthedocs.io/ ·仓库 https://github.com/zqfang/GSEApy ；g:Profiler https://biit.cs.ut.ee/gprofiler/ ；Enrichr https://maayanlab.cloud/Enrichr/ ；MSigDB https://www.gsea-msigdb.org/gsea/msigdb/ ；GSEA 方法 Subramanian et al. (2005) PNAS, DOI: 10.1073/pnas.0506580102。

---

本条采编自 K-Dense-AI/scientific-agent-skills（MIT 许可）。