← 生物信息学入门 · 四个实验

实验 3：富集分析（ORA 超几何检验 · BH 校正 · GSEA）

预计用时：6–8 小时。
实验 2 给出了约 1000 个差异基因，但一张基因列表本身讲不出生物学故事。这个实验把基因列表翻译成"通路"。
输入统一用实验 2 参考解的结果（course fetch 3 会自动算好并缓存），所以实验 2 没做完也可以先做这个。

一、问题与数据

问题：地塞米松上调和下调的基因，分别集中在哪些生物学过程里？这些过程和糖皮质激素"抗炎"的作用对得上吗？

数据内容来源
data/cache/airway_de_results.csv实验 2 的差异表达结果，约 19000 个基因由实验 2 参考解生成
data/genesets/h.all.v2024.1.Hs.symbols.gmtMSigDB Hallmark 基因集：50 个精选通路，每个约 30–200 个基因Broad Institute MSigDB v2024.1（CC BY 4.0）

Hallmark 是 MSigDB 从几千个原始基因集里提炼出的 50 个"标志性"生物学过程，冗余少，很适合入门。

二、目标

  1. 读懂 GMT 这种基因集格式
  2. 手写 BH 多重检验校正，理解 FDR 控制的是什么
  3. 把"通路富集"写成一个抽球问题，用超几何分布算 p 值
  4. 手写 ORA，并想清楚背景基因集该怎么选
  5. 用 GSEA 做不设阈值的富集分析，比较它和 ORA 的结果有什么不同
  6. 把结果放回生物学里解释

三、前置

概念一句话
基因集 / 通路一组功能相关的基因，比如"TNF-α 通过 NF-κB 的信号转导"通路里的 200 个基因
ORA先用阈值挑出差异基因，再问：这个通路里的基因，在差异基因中占的比例是否高于随机预期？
超几何分布袋子里有 M 个球，其中 n 个红球。不放回地抽 N 个，抽到 k 个红球的概率
背景（universe）袋子里原本有哪些球。应该是实验里检测得到的基因，不是全基因组的两万个
fold enrichment实际命中数 ÷ 期望命中数，大于 1 表示富集
BH / FDR50 个通路就是 50 次检验。BH 校正后的 padj < 0.05，意思是在"显著"的通路里，假阳性所占比例约 5%
GSEA不设阈值，把所有基因按变化方向和强度排好序，看某个通路的基因是否集中在列表的顶端或底端
NESGSEA 的标准化富集分数。正值表示富集在上调端，负值表示富集在下调端

两个关键的生物学背景：

NFKBIA 也属于这类基因，它在本数据里同样上调（log2FC ≈ 0.8，padj 显著），只是没过 |log2FC| ≥ 1 的阈值，所以 ORA 看不到它。GSEA 不设阈值，能把它算进去。

四、任务

先运行 course fetch 3 和 course start 3，然后在 交付/enrich.py 里完成。

任务 1：读 GMT read_gmt

任务 2：BH 校正 bh_adjust

任务 3：超几何检验 hypergeom_pvalue

任务 4：ORA ora

任务 5：GSEA rank_metric / run_gsea

任务 6：course run 3

交付/ 里会生成：ora_up.csv、ora_down.csv、ora_dotplot.png、gsea.csv、gsea_nes.png、结果摘要.md。

思考题（写进 交付/思考题.md）

  1. 打开 ora_up.csv，看 TNFA_SIGNALING_VIA_NFKB 那一行命中的 26 个基因。查其中 5 个，比如 DUSP1、TSC22D1、NFIL3、KLF9、PER1，它们是促进炎症还是抑制炎症？这和"激素抗炎"矛盾吗？
  2. 把背景改成全部 Hallmark 基因（4384 个），或者全基因组（约 2 万个），重新跑 ORA。p 值会变大还是变小？哪种背景是对的，为什么？
  3. ORA 和 GSEA 找到的通路有哪些一样、哪些不一样？各举一个例子，说明一种方法能发现、另一种发现不了的情况。
  4. P53_PATHWAY 同时出现在上调和下调的列表里，这说明什么？"通路被激活"和"通路里有基因发生了变化"是一回事吗？

选做

五、验收标准

检查点通过标准对应测试
GMT50 个基因集，P53 通路 200 个基因任务1_*
BH手算例子正确；和 scipy 一致；结果 ≤ 1任务2_*
超几何和 Fisher 单侧检验一致任务3_*
ORA背景 3483，查询 164；上调第一名 TNFA/NF-κB（26 个），下调第一名 P53（20 个）任务4_*
GSEATNFA/NF-κB 的 NES > 1.5，FDR < 0.05任务5_*
交付course check 3 全部通过；6 个文件齐全；思考题已写任务6_*

六、常见错误

现象原因解决
p 值整体偏大一截sf(k) 少减了 1sf(k - 1, M, n, N)
BH 结果不单调漏了"取累计最小值"这一步np.minimum.accumulate(x[::-1])[::-1]
BH 结果顺序全乱了忘了按原顺序放回out[order] = adjusted
背景变成 19211，p 值偏小把没有注释的基因也算进了背景背景先和所有基因集的并集取交集
GSEA 报 "duplicate index"同名基因没有去重保留 \
GSEA 结果排序怪怪的NES 是字符串，按字母顺序排了先 astype(float)
每次运行 GSEA 结果都略有不同置换检验是随机的固定 seed

参考实现：enrich

卡住超过 20 分钟再看。对照着看自己是哪一步想岔了 —— 抄一遍没用。

  1. 展开参考解
    看答案
    """实验 3 · 富集分析（参考解）

输入：实验 2 的差异表达结果（airway，地塞米松 vs 未处理）。
问题：上调、下调的那几百个基因，集中在哪些生物学通路里？
"""
from __future__ import annotations

import sys
from pathlib import Path

import numpy as np
import pandas as pd
from scipy.stats import hypergeom

HERE = Path(__file__).resolve().parent
ROOT = HERE.parents[2]
sys.path.insert(0, str(ROOT / "tools"))
import datasets  # noqa: E402
import de_cache  # noqa: E402


# ---------------------------------------------------------------- 任务 1：读基因集
def read_gmt(path) -> dict[str, set[str]]:
    """GMT：每行 = 基因集名 \\t 描述/链接 \\t 基因1 \\t 基因2 …"""
    sets: dict[str, set[str]] = {}
    with open(path, encoding="utf-8") as f:
        for line in f:
            parts = line.rstrip("\r\n").split("\t")
            if len(parts) < 3:
                continue
            sets[parts[0]] = {g for g in parts[2:] if g}
    return sets


# ---------------------------------------------------------------- 任务 2：多重检验校正
def bh_adjust(pvalues) -> np.ndarray:
    """Benjamini–Hochberg：p_(i) × m / i，再从大到小取累计最小值，封顶 1；按原顺序返回。"""
    p = np.asarray(pvalues, dtype=float)
    m = len(p)
    if m == 0:
        return p
    order = np.argsort(p)
    ranked = p[order] * m / np.arange(1, m + 1)
    ranked = np.minimum.accumulate(ranked[::-1])[::-1]
    out = np.empty(m)
    out[order] = np.minimum(ranked, 1.0)
    return out


# ---------------------------------------------------------------- 任务 3：超几何检验
def hypergeom_pvalue(k: int, M: int, n: int, N: int) -> float:
    """背景 M 个基因，其中 n 个属于这个通路；随机抽 N 个，抽中 ≥ k 个通路基因的概率。"""
    return float(hypergeom.sf(k - 1, M, n, N))


# ---------------------------------------------------------------- 任务 4：ORA
def ora(genes, universe, gene_sets: dict[str, set[str]], min_size: int = 15, max_size: int = 500) -> pd.DataFrame:
    """过表达分析（over-representation analysis）。

    背景 = universe（检测过的基因）∩ 至少出现在一个基因集里的基因——
    和 clusterProfiler 的默认做法一致：没有任何注释的基因不可能"命中"，不该稀释背景。
    """
    annotated = set().union(*gene_sets.values())
    bg = set(universe) & annotated
    query = set(genes) & bg
    M, N = len(bg), len(query)
    rows = []
    for term, members in gene_sets.items():
        s = members & bg
        n = len(s)
        if n < min_size or n > max_size:
            continue
        hits = sorted(query & s)
        k = len(hits)
        expected = N * n / M if M else 0.0
        rows.append({
            "term": term, "overlap": k, "set_size": n, "query_size": N, "universe_size": M,
            "expected": expected, "fold_enrichment": k / expected if expected else np.nan,
            "pvalue": hypergeom_pvalue(k, M, n, N), "genes": ",".join(hits),
        })
    table = pd.DataFrame(rows)
    if table.empty:
        return table
    table["padj"] = bh_adjust(table["pvalue"])
    return table.sort_values(["pvalue", "term"]).reset_index(drop=True)


# ---------------------------------------------------------------- 任务 5：GSEA（排序式）
def rank_metric(results: pd.DataFrame) -> pd.Series:
    """用 DESeq2 的 Wald 统计量（stat）给所有检测过的基因排序：正=处理后升高，负=降低。
    去掉 stat 为 NaN、没有基因名的；同名基因保留 |stat| 最大的一个；从大到小排。"""
    df = results.loc[results["stat"].notna() & (results["gene_name"] != ""), ["gene_name", "stat"]]
    df = df.assign(abs_stat=df["stat"].abs()).sort_values("abs_stat", ascending=False)
    df = df.drop_duplicates("gene_name")
    return df.set_index("gene_name")["stat"].sort_values(ascending=False)


def run_gsea(ranking: pd.Series, gene_sets: dict[str, set[str]], permutations: int = 1000, seed: int = 7) -> pd.DataFrame:
    import gseapy

    pre = gseapy.prerank(
        rnk=ranking, gene_sets={k: sorted(v) for k, v in gene_sets.items()},
        permutation_num=permutations, min_size=15, max_size=500, seed=seed,
        threads=1, outdir=None, verbose=False,
    )
    res = pre.res2d.copy()
    res["NES"] = res["NES"].astype(float)
    res["FDR q-val"] = res["FDR q-val"].astype(float)
    return res.sort_values("NES", ascending=False).reset_index(drop=True)


# ---------------------------------------------------------------- 任务 6：出图与报告
def load_inputs():
    res = de_cache.load()
    gene_sets = read_gmt(datasets.path_of("hallmark_gmt"))
    tested = res[res["padj"].notna() & (res["gene_name"] != "")]
    universe = set(tested["gene_name"])
    up = set(tested.loc[tested["class"] == "up", "gene_name"])
    down = set(tested.loc[tested["class"] == "down", "gene_name"])
    return res, gene_sets, universe, up, down


def main(out_dir) -> None:
    import matplotlib

    matplotlib.use("Agg")
    import matplotlib.pyplot as plt

    out = Path(out_dir)
    out.mkdir(parents=True, exist_ok=True)
    res, gene_sets, universe, up, down = load_inputs()

    tables = {"up": ora(up, universe, gene_sets), "down": ora(down, universe, gene_sets)}
    for direction, t in tables.items():
        t.to_csv(out / f"ora_{direction}.csv", index=False, encoding="utf-8-sig")

    # 点图：两个方向各取前 10
    fig, axes = plt.subplots(1, 2, figsize=(12, 5), sharex=False)
    for ax, (direction, t) in zip(axes, tables.items()):
        top = t.head(10).iloc[::-1]
        sc = ax.scatter(top["fold_enrichment"], top["term"].str.replace("HALLMARK_", ""),
                        s=top["overlap"] * 8, c=-np.log10(top["padj"]), cmap="viridis")
        ax.set(xlabel="fold enrichment", title=f"ORA · {direction}-regulated (Hallmark)")
        fig.colorbar(sc, ax=ax, label="-log10 padj")
    fig.tight_layout()
    fig.savefig(out / "ora_dotplot.png", dpi=150)
    plt.close(fig)

    # GSEA
    ranking = rank_metric(res)
    gsea = run_gsea(ranking, gene_sets)
    gsea.to_csv(out / "gsea.csv", index=False, encoding="utf-8-sig")
    fig, ax = plt.subplots(figsize=(7, 8))
    g = gsea.sort_values("NES")
    colors = np.where(g["FDR q-val"] < 0.05, np.where(g["NES"] > 0, "#d62728", "#1f77b4"), "#bbbbbb")
    ax.barh(g["Term"].str.replace("HALLMARK_", ""), g["NES"], color=colors)
    ax.tick_params(axis="y", labelsize=6)
    ax.set(xlabel="NES (dex vs untreated)", title="GSEA preranked · Hallmark (colored: FDR < 0.05)")
    fig.tight_layout()
    fig.savefig(out / "gsea_nes.png", dpi=150)
    plt.close(fig)

    lines = ["# 实验 3 结果摘要（airway · Hallmark）\n",
             f"- 背景（检测过且有 Hallmark 注释的基因）：{tables['up']['universe_size'].iloc[0]} 个",
             f"- 上调基因 {len(up)} 个，其中有注释的 {tables['up']['query_size'].iloc[0]} 个；"
             f"下调 {len(down)} 个，有注释的 {tables['down']['query_size'].iloc[0]} 个\n"]
    for direction, t in tables.items():
        lines += [f"## ORA · {direction}\n", "| term | overlap | set | fold | padj |", "|---|---:|---:|---:|---:|"]
        lines += [f"| {r.term} | {r.overlap} | {r.set_size} | {r.fold_enrichment:.2f} | {r.padj:.2e} |"
                  for r in t.head(8).itertuples()]
        lines.append("")
    sig = gsea[gsea["FDR q-val"] < 0.05]
    lines += ["## GSEA（FDR < 0.05）\n", "| term | NES | FDR |", "|---|---:|---:|"]
    lines += [f"| {t} | {n:.2f} | {q:.3f} |" for t, n, q in zip(sig["Term"], sig["NES"], sig["FDR q-val"])]
    (out / "结果摘要.md").write_text("\n".join(lines) + "\n", encoding="utf-8")
    print("\n".join(lines))
← 实验 2：批量 RNA-seq 差异表达（真实数据 · PyDESeq2）目录