"""从零实现一个教学版 GSEA，并对应项目中的两种基因排序方法。 运行：python Learn/gsea_from_scratch.py 依赖：numpy、pandas、matplotlib 输出：Learn/gsea_output/gsea_results.tsv 和 enrichment_curves.png 说明：本代码用于理解算法。正式分析应使用 clusterProfiler、fgsea 或 gseapy， 因为成熟工具对置换策略、并列值、NES 和多重检验有更完整的处理。 """ from pathlib import Path from math import erfc import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np import pandas as pd SEED = 20260719 def make_demo_data(n_genes=500, seed=SEED): """构造一份可重复的差异表达表和三个容易观察的基因集。""" rng = np.random.default_rng(seed) genes = np.array([f"Gene_{i:03d}" for i in range(1, n_genes + 1)]) log2fc = rng.normal(0, 0.7, n_genes) # 人为制造两个生物学信号：前 25 个上调，随后 25 个下调。 log2fc[:25] += 2.5 log2fc[25:50] -= 2.5 z = log2fc / 0.7 + rng.normal(0, 0.5, n_genes) pvalue = np.fromiter( (erfc(value / np.sqrt(2)) for value in np.abs(z)), dtype=float, count=n_genes ).clip(1e-300, 1.0) de = pd.DataFrame({"gene_symbol": genes, "log2fc": log2fc, "pvalue": pvalue}) gene_sets = { "UP_PATHWAY": set(genes[:25]), "DOWN_PATHWAY": set(genes[25:50]), "RANDOM_PATHWAY": set(rng.choice(genes[50:], size=25, replace=False)), } return de, gene_sets def rank_genes(de, method="signed_logp"): """生成降序排列的基因列表，对应 04_enrichment.R 的 rank 计算。""" data = de.dropna(subset=["gene_symbol", "log2fc", "pvalue"]).copy() if method == "signed_logp": # 上调为正、下调为负；P 越小，分数绝对值越大。 data["rank_score"] = np.sign(data["log2fc"]) * -np.log10( data["pvalue"].clip(lower=1e-300) ) elif method == "log2fc": data["rank_score"] = data["log2fc"] else: raise ValueError("method 必须是 'signed_logp' 或 'log2fc'") # 同名基因保留绝对排名分数最大的一条，与 R 脚本的处理一致。 data["abs_score"] = data["rank_score"].abs() data = data.sort_values("abs_score", ascending=False).drop_duplicates("gene_symbol") return data.sort_values("rank_score", ascending=False)[["gene_symbol", "rank_score"]] def enrichment_score(ranked, gene_set, weight=1.0): """计算加权运行和及 ES；ES 是运行和偏离 0 最远的位置。""" genes = ranked["gene_symbol"].to_numpy() scores = ranked["rank_score"].to_numpy() hits = np.isin(genes, list(gene_set)) n_hit, n_total = int(hits.sum()), len(genes) if n_hit == 0 or n_hit == n_total: raise ValueError("基因集与排序列表的交集必须非空且不能覆盖全部基因") hit_weights = np.abs(scores[hits]) ** weight # 命中时按排名权重上升；未命中时均匀下降。最终运行和应回到 0。 increments = np.full(n_total, -1.0 / (n_total - n_hit)) increments[hits] = hit_weights / hit_weights.sum() running = np.cumsum(increments) max_i, min_i = int(np.argmax(running)), int(np.argmin(running)) es = running[max_i] if abs(running[max_i]) >= abs(running[min_i]) else running[min_i] peak_i = max_i if es >= 0 else min_i # 正 ES 的 leading edge 位于峰值之前；负 ES 位于谷值之后。 leading = genes[: peak_i + 1][hits[: peak_i + 1]] if es >= 0 else genes[peak_i:][hits[peak_i:]] return float(es), running, hits, leading def permutation_test(ranked, gene_set, n_perm=500, seed=SEED): """随机抽取相同大小的基因集，得到零分布、名义 P 值和 NES。""" rng = np.random.default_rng(seed) observed, running, hits, leading = enrichment_score(ranked, gene_set) genes = ranked["gene_symbol"].to_numpy() null_es = np.empty(n_perm) for i in range(n_perm): random_set = set(rng.choice(genes, size=len(gene_set), replace=False)) null_es[i] = enrichment_score(ranked, random_set)[0] same_sign = null_es[null_es >= 0] if observed >= 0 else -null_es[null_es < 0] observed_abs = observed if observed >= 0 else -observed # 加 1 避免有限次置换产生 P=0。 pvalue = (1 + np.sum(same_sign >= observed_abs)) / (1 + len(same_sign)) nes = observed / same_sign.mean() return observed, nes, pvalue, running, hits, leading, null_es def benjamini_hochberg(pvalues): """用 BH 方法把多个基因集的名义 P 值校正为 FDR。""" pvalues = np.asarray(pvalues, dtype=float) order = np.argsort(pvalues) ranked_p = pvalues[order] adjusted = ranked_p * len(pvalues) / np.arange(1, len(pvalues) + 1) adjusted = np.minimum.accumulate(adjusted[::-1])[::-1].clip(0, 1) result = np.empty_like(adjusted) result[order] = adjusted return result def run_demo(rank_method="signed_logp", n_perm=500): de, gene_sets = make_demo_data() ranked = rank_genes(de, rank_method) details, rows = {}, [] for i, (name, genes) in enumerate(gene_sets.items()): result = permutation_test(ranked, genes, n_perm=n_perm, seed=SEED + i) es, nes, pvalue, running, hits, leading, null_es = result details[name] = result rows.append({ "pathway": name, "size": int(hits.sum()), "ES": es, "NES": nes, "pvalue": pvalue, "leading_edge": ";".join(leading), }) results = pd.DataFrame(rows) results["padj_BH"] = benjamini_hochberg(results["pvalue"]) results = results.sort_values("padj_BH") return ranked, results, details def plot_curves(ranked, results, details, output_file): """画出经典 GSEA 图：运行富集分数、命中刻线和排名统计量。""" fig, axes = plt.subplots(3, 1, figsize=(10, 8), sharex=True, gridspec_kw={"height_ratios": [3, 0.7, 1.4]}) colors = {"UP_PATHWAY": "#C43D3D", "DOWN_PATHWAY": "#2878A5", "RANDOM_PATHWAY": "#666666"} for pathway in results["pathway"]: es, nes, _, running, hits, *_ = details[pathway] axes[0].plot(running, color=colors[pathway], label=f"{pathway} NES={nes:.2f}") axes[1].vlines(np.flatnonzero(hits), 0, 1, color=colors[pathway], alpha=0.65, linewidth=0.7) axes[0].axhline(0, color="black", linewidth=0.7) axes[0].set_ylabel("Running ES") axes[0].legend(frameon=False) axes[1].set_yticks([]) axes[1].set_ylabel("Hits") axes[2].plot(ranked["rank_score"].to_numpy(), color="#333333", linewidth=0.9) axes[2].axhline(0, color="#999999", linewidth=0.7) axes[2].set_ylabel("Rank score") axes[2].set_xlabel("Genes ordered from high to low rank score") fig.suptitle("GSEA from scratch: enrichment along a ranked gene list") fig.tight_layout() fig.savefig(output_file, dpi=160, bbox_inches="tight") plt.close(fig) def main(): output_dir = Path(__file__).resolve().parent / "gsea_output" output_dir.mkdir(exist_ok=True) ranked, results, details = run_demo(rank_method="signed_logp", n_perm=500) results.to_csv(output_dir / "gsea_results.tsv", sep="\t", index=False) ranked.to_csv(output_dir / "ranked_genes.tsv", sep="\t", index=False) plot_curves(ranked, results, details, output_dir / "enrichment_curves.png") print("\nGSEA 的逻辑：") print("1. 用 signed_logp 或 log2fc 将所有基因从高到低排序。") print("2. 遇到通路基因时运行和上升，否则下降；最大偏离值就是 ES。") print("3. 随机基因集产生零分布，计算 nominal P；ES 除以零分布均值为 NES。") print("4. 对多个通路的 P 值做 BH 校正，得到教学版 padj_BH。\n") print(results[["pathway", "size", "ES", "NES", "pvalue", "padj_BH"]].to_string(index=False)) print(f"\n结果已写入：{output_dir}") if __name__ == "__main__": main()