pathway-enrichment

通路富集

概述

富集分析回答的是“我的基因中哪些生物学过程出现过度代表？”这是差异表达分析、筛选或聚类之后的标准最后一步。主要有两种方法，正确选择方法是最重要的一项决策：

• ORA（过度代表分析）——使用经过阈值筛选的基因列表（例如 padj < 0.05），通过 Fisher 精确检验／超几何检验，测试哪些基因集的重叠程度高于偶然水平。工具：Enrichr、g:Profiler。
• GSEA（基因集富集分析）——使用完整的排序基因列表（不设阈值），测试每个基因集是否集中在列表顶部或底部。预排序 GSEA 使用每个基因的分数（例如 DESeq2 的 stat）。当效应广泛但细微时，这种方法效果更好。

此技能负责协调这些分析、背后的基因集数据库，以及可能导致结果错误或无法发表的解读陷阱。

何时使用此技能

当用户希望进行以下操作时，请使用此技能：

• 在基因列表中查找富集的 GO 术语／KEGG／Reactome／WikiPathways／MSigDB Hallmark 基因集。
• 对 DESeq2、edgeR、limma 或 Scanpy rank_genes_groups 输出运行 GSEA／预排序 GSEA。
• 按样本／细胞计算通路活性分数（ssGSEA、GSVA）。
• 解读、去重和可视化富集结果，或制作论文表格／图。
• 在 ORA 和 GSEA 之间进行选择、挑选基因集库、选择背景集，或解决基因 ID 问题。

对于快速、一次性的 Enrichr 查询，gget 技能（gget enrichr）更轻量；对于原始通路／相互作用 API（Reactome、KEGG、STRING），请参阅 database-lookup 技能。需要完整且经得起检验的富集工作流时，请使用此技能。

选择合适的方法

情况	方法	工具／入口点
你有一个离散的命中列表（差异表达基因、筛选命中基因、簇标记基因）	ORA	gp.enrichr(...) 或 g:Profiler
你有一个完整的排序列表（每个受检基因及其分数）	预排序 GSEA	gp.prerank(...)
你有表达矩阵和类别标签	GSEA	gp.gsea(...)
你想按样本／细胞计算通路分数	ssGSEA／GSVA	gp.ssgsea(...)、gp.gsva(...)
你需要自定义背景集或 500+ 种生物	使用自定义范围的 ORA	g:Profiler（domain_scope='custom'）
你想推断 TF／信号传导的活性（PROGENy、DoRothEA）	活性推断	参阅 references/databases-and-gene-sets.md（decoupler）

不确定时：经过阈值筛选的列表 → ORA；带分数的排序表 → GSEA。切勿对列表设置阈值后再将其输入 GSEA——这样会丢弃 GSEA 所依赖的排序信息。

设置

uv pip install gseapy gprofiler-official
# gseapy pulls pandas, numpy, scipy, matplotlib. Network access is needed for
# Enrichr, g:Profiler, and MSigDB downloads. For fully offline ORA, use a local
# GMT file with gp.enrich() (see references/gseapy.md).

验证并列出可用的基因集库（名称会随时间变化，切勿盲目硬编码）：

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

快速入门

对命中列表运行 ORA（gseapy + Enrichr）

import gseapy as gp

# Enrichr libraries expect HGNC gene SYMBOLS (human: UPPERCASE). Map IDs first if needed.
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,            # in-memory; set a path to also write tables/plots
)
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))

根据 DESeq2 结果运行预排序 GSEA

import gseapy as gp
import pandas as pd

res = pd.read_csv("deseq2_results.csv", index_col=0)   # index = gene symbols
# Rank by the test statistic (sign = direction, magnitude = evidence). This is
# more stable than ranking by log2FoldChange, which is noisy for low-count genes.
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,        # drop tiny/huge sets (noisy or generic)
    permutation_num=1000, seed=123,   # seed = reproducible p-values
    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) 构建排序分数。

核心工作流

要进行经得起检验的分析，请逐步完成以下步骤。中间几步（ID 类型、背景集）最容易在不知不觉中出错。

步骤 1 — 明确输入并选择方法

确认：涉及哪些基因、所属物种、是否有每个基因的分数（→ GSEA）还是只有列表（→ ORA），以及这些基因代表什么比较结果（方向会影响解读）。

步骤 2 — 将基因 ID 转换为正确的命名空间

Enrichr／MSigDB 库以基因符号为键（人类为全大写，鼠类为首字母大写）。如果使用 Ensembl／Entrez ID，请先进行转换。gp.Biomart、g:Profiler g:Convert 和 mygene 的用法请参阅 references/databases-and-gene-sets.md。ID 不匹配却未被察觉，是“没有任何结果显著”的首要原因。

步骤 3 — 选择与问题匹配的基因集文库

Hallmark（广泛主题）→ GO:BP（机制）→ KEGG/Reactome/WikiPathways（整理好的通路）→ C7（免疫）等。不要运行 50 个文库；选择 2–4 个符合生物学背景的文库。文库目录和选择指南：references/databases-and-gene-sets.md。

步骤 4 — 设置背景范围（仅限 ORA）

背景必须是检测实验中可能被检测到的基因（例如所有表达或测试过的基因），而不是整个基因组。错误的背景会夸大显著性。Enrichr 使用固定背景；背景设置很重要时，请使用带有 domain_scope='custom' 和 background 的 g:Profiler，或使用显式指定背景的 gp.enrich()。原理见 references/interpretation.md。

步骤 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 — 可视化

gseapy 内置点图、条形图、富集图和 GSEA 运行分数图（gp.dotplot、gp.barplot、gp.enrichment_map、gp.gseaplot）。详见 references/gseapy.md。

步骤 8 — 消除冗余并解读

GO 尤其容易返回许多近乎重复的术语。使用富集图（术语之间的相似性）、前沿基因重叠或父级术语进行合并，并报告有代表性的术语。解读框架和论文表格格式见 references/interpretation.md。

辅助脚本

scripts/run_enrichment.py 可端到端运行 ORA 或 GSEA，并生成结果表和点图，同时处理样板流程（符号清理、去重、移除 NA、根据 DESeq2 表构建排名、按文库进行 FDR 筛选）。

# ORA from a hit list (one gene symbol per line)
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/

# Preranked GSEA from a DESeq2 results CSV (auto-builds the rank from `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

# Preranked GSEA from an explicit 2-column rank file (gene,score)
python scripts/run_enrichment.py gsea --rnk ranked_genes.csv --outdir results/

运行 python scripts/run_enrichment.py --help 查看所有选项（背景文件、FDR 阈值、基因集大小最小值／最大值、置换次数）。

常见陷阱

以下问题会导致大多数错误或无法复现的结果：

1. 基因 ID／物种不匹配 — 符号与 Ensembl 不匹配，或人类与小鼠的大小写不匹配。请映射 ID 并正确设置 organism，否则匹配项会在不知不觉中降至约 0。
2. 背景错误（ORA） — 使用整个基因组，而非经过测试／表达的基因集，会夸大 p 值。在背景设置很重要时，请设置自定义背景。
3. GSEA 之前先进行阈值筛选 — GSEA 需要完整的排序列表；只有 ORA 使用经过筛选的列表。
4. 仅按 log2FoldChange 对 GSEA 排名 — 对低计数基因而言不稳定；优先使用 stat 或 sign(LFC) * -log10(p)。
5. 跨文库进行多重检验 — FDR 是在每个文库内计算的；运行许多文库会增加检验次数。请报告每个文库的 FDR，并保持谨慎。
6. GO 术语冗余 — 不要报告同一术语的 40 种变体；请合并并展示有代表性的术语。
7. 显著性 ≠ 相关性 — 检查重叠数量和基因集大小；很小的基因集很容易达到显著性。
8. ORA 列表过短／过长 — 少于 10 个基因时统计效能不足；超过 2000 个基因则特异性不足（可考虑改用 GSEA）。
9. 缺少可复现性元数据 — Enrichr／GO 文库带有版本信息，并会随时间变化。请记录文库名称和日期，并为 GSEA 设置 seed。

与其他技能集成

• 上游（基因来源）： pydeseq2（差异表达基因及用于 GSEA 的 stat）、scanpy（rank_genes_groups 标记基因／分数）、depmap／pytdc（筛选命中项）、蛋白质组学技能（pyopenms、matchms）。
• 数据库／ID： database-lookup（Reactome、KEGG、STRING、Gene Ontology API）、gget（gget enrichr 快速路径、gget info 用于 ID 映射）、bioservices。
• 下游： scientific-visualization（自定义图表）、networkx（富集图网络）、scientific-writing／literature-review（解读和引用）、statistical-analysis（多重检验详情）。

参考文件

需要深入了解时，请阅读相关文件：

• references/gseapy.md — 完整的 gseapy API：enrichr、离线 enrich、prerank、gsea、ssgsea、gsva、Msigdb、Biomart、get_library_name/read_gmt、所有绘图、结果列含义、GMT/离线用法，以及故障排查（速率限制、结果为空）。
• references/databases-and-gene-sets.md — GO、KEGG、Reactome、WikiPathways、MSigDB 集合、Enrichr 库命名、g:Profiler 数据源、物种处理、基因 ID 转换、根据问题选择库，以及 Reactome/STRING API 和 decoupler 活性推断的参考信息。
• references/interpretation.md — ORA 与 GSEA 统计、背景基因集选择、多重检验方法（BH、g:SCS 与 Bonferroni）、核心富集基因、冗余项精简、效应与显著性的区别、发表用表格模板，以及可复现性检查清单。

资源

• gseapy 文档：https://gseapy.readthedocs.io/ · 仓库：https://github.com/zqfang/GSEApy
• g:Profiler：https://biit.cs.ut.ee/gprofiler/ · Python 客户端：https://pypi.org/project/gprofiler-official/
• Enrichr：https://maayanlab.cloud/Enrichr/ · MSigDB：https://www.gsea-msigdb.org/gsea/msigdb/
• GSEA 方法：Subramanian 等（2005），PNAS，DOI：10.1073/pnas.0506580102

相关技能

• bulk-rnaseq：端到端 bulk RNA-seq 编排工具 — 从原始 FASTQ 读段开始，经过 QC 和剪切（FastQC、fastp/Trim Galore）、比对和定量…
• bioservices：统一的 Python 接口，可访问 40 多种生物信息学服务。
• gget：通过 gget CLI 和 Python 包查询 20 多个生物信息学数据库和分析服务，涵盖 Ensembl 基因搜索、信息查询和…