ora demo

从零实现一个教学版 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_genes500, seedSEED): 构造一份可重复的差异表达表和三个容易观察的基因集。 rng np.random.default_rng(seed) genes np.array([fGene_{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)), dtypefloat, countn_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:], size25, replaceFalse)), } return de, gene_sets def rank_genes(de, methodsigned_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(lower1e-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, ascendingFalse).drop_duplicates(gene_symbol) return data.sort_values(rank_score, ascendingFalse)[[gene_symbol, rank_score]] def enrichment_score(ranked, gene_set, weight1.0): 计算加权运行和及 ESES 是运行和偏离 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_perm500, seedSEED): 随机抽取相同大小的基因集得到零分布、名义 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, sizelen(gene_set), replaceFalse)) 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 避免有限次置换产生 P0。 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, dtypefloat) 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_methodsigned_logp, n_perm500): 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_permn_perm, seedSEED 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), sharexTrue, 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, colorcolors[pathway], labelf{pathway} NES{nes:.2f}) axes[1].vlines(np.flatnonzero(hits), 0, 1, colorcolors[pathway], alpha0.65, linewidth0.7) axes[0].axhline(0, colorblack, linewidth0.7) axes[0].set_ylabel(Running ES) axes[0].legend(frameonFalse) axes[1].set_yticks([]) axes[1].set_ylabel(Hits) axes[2].plot(ranked[rank_score].to_numpy(), color#333333, linewidth0.9) axes[2].axhline(0, color#999999, linewidth0.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, dpi160, bbox_inchestight) plt.close(fig) def main(): output_dir Path(__file__).resolve().parent / gsea_output output_dir.mkdir(exist_okTrue) ranked, results, details run_demo(rank_methodsigned_logp, n_perm500) results.to_csv(output_dir / gsea_results.tsv, sep\t, indexFalse) ranked.to_csv(output_dir / ranked_genes.tsv, sep\t, indexFalse) plot_curves(ranked, results, details, output_dir / enrichment_curves.png) print(\nGSEA 的逻辑) print(1. 用 signed_logp 或 log2fc 将所有基因从高到低排序。) print(2. 遇到通路基因时运行和上升否则下降最大偏离值就是 ES。) print(3. 随机基因集产生零分布计算 nominal PES 除以零分布均值为 NES。) print(4. 对多个通路的 P 值做 BH 校正得到教学版 padj_BH。\n) print(results[[pathway, size, ES, NES, pvalue, padj_BH]].to_string(indexFalse)) print(f\n结果已写入{output_dir}) if __name__ __main__: main()