65 lines
2.7 KiB
Python
65 lines
2.7 KiB
Python
#!/usr/bin/env python
|
||
# -*- coding: utf-8 -*-
|
||
"""
|
||
16_gsea_prerank.py — 补齐 Step 6 的 GSEA preranked 富集分析（Hallmark）
|
||
|
||
排序统计量：pseudobulk median-of-ratios LFC 的 S1/S2 均值（= 方向一致性 LFC）。
|
||
基因集：MSigDB Hallmark（小鼠，50 条）——作者未用 Hallmark（仅 Reactome/WikiPathways），
|
||
故本步是相对作者的增量；Reactome 富集直接用作者 S3–S12 交叉验证。
|
||
组：Muller / RGC2-like / RGC1-like（RGC2-like 标低功效）。
|
||
输出：output/16_gsea/ prerank 结果表 + NES 条形图。
|
||
"""
|
||
from pathlib import Path
|
||
|
||
import matplotlib
|
||
matplotlib.use("Agg")
|
||
import matplotlib.pyplot as plt
|
||
import numpy as np
|
||
import pandas as pd
|
||
import gseapy as gp
|
||
|
||
OUT = Path("output/16_gsea")
|
||
OUT.mkdir(parents=True, exist_ok=True)
|
||
plt.rcParams.update({"figure.dpi": 300, "savefig.dpi": 300, "font.size": 8})
|
||
|
||
lib = gp.get_library(name="MSigDB_Hallmark_2020", organism="Mouse")
|
||
print(f"Hallmark: {len(lib)} sets")
|
||
|
||
GROUPS = ["Muller", "RGC2-like", "RGC1-like"]
|
||
all_results = {}
|
||
for grp in GROUPS:
|
||
de = pd.read_csv(f"output/11_pseudobulk_de/de_{grp}.csv")
|
||
rnk = de[["gene", "lfc_S1", "lfc_S2"]].copy()
|
||
rnk["score"] = (rnk["lfc_S1"] + rnk["lfc_S2"]) / 2
|
||
rnk = rnk[["gene", "score"]].dropna().sort_values("score", ascending=False)
|
||
print(f"== {grp}: {len(rnk)} 基因")
|
||
res = gp.prerank(rnk=rnk, gene_sets=lib, min_size=15, max_size=500,
|
||
permutation_num=1000, seed=0, threads=4, no_plot=True)
|
||
out = res.res2d.copy()
|
||
out.to_csv(OUT / f"gsea_{grp}_Hallmark.csv", index=False)
|
||
all_results[grp] = out
|
||
sig = out[out["FDR q-val"].astype(float) < 0.25]
|
||
print(f" FDR<0.25: {len(sig)} 条")
|
||
if len(sig):
|
||
print(sig[["Term", "NES", "FDR q-val"]].sort_values("NES").to_string(index=False))
|
||
|
||
# ---- NES 条形图（三组并排，取所有 Hallmark 条目的 NES） ----
|
||
terms = all_results["Muller"]["Term"].tolist()
|
||
nes = pd.DataFrame({g: all_results[g].set_index("Term")["NES"].astype(float).reindex(terms)
|
||
for g in GROUPS})
|
||
nes = nes.loc[nes.abs().max(axis=1) > 0.5]
|
||
order = nes.abs().max(axis=1).sort_values().index
|
||
fig, ax = plt.subplots(figsize=(7, 9))
|
||
y = np.arange(len(order))
|
||
w = 0.28
|
||
for off, g, c in [(-w, "Muller", "#1f77b4"), (0, "RGC2-like", "#d62728"), (w, "RGC1-like", "#ff9896")]:
|
||
ax.barh(y + off, nes.loc[order, g], height=w, color=c, label=g, alpha=0.85)
|
||
ax.set_yticks(y, [t.replace("HALLMARK_", "").replace("_", " ").title() for t in order], fontsize=7)
|
||
ax.axvline(0, color="grey", lw=0.6)
|
||
ax.set_xlabel("NES")
|
||
ax.legend(fontsize=8, frameon=False)
|
||
ax.set_title("Hallmark GSEA prerank (NES, |NES|>0.5 in any group)")
|
||
fig.tight_layout()
|
||
fig.savefig(OUT / "gsea_hallmark_nes.png", bbox_inches="tight")
|
||
print("Done ->", OUT)
|