95 lines
4.7 KiB
Python
95 lines
4.7 KiB
Python
#!/usr/bin/env python
|
||
# -*- coding: utf-8 -*-
|
||
"""
|
||
13_isg_check.py — P1-⑤ RGC 的 ISG 逐基因核查：通路特异抑制 vs 全局转录抑制
|
||
|
||
做法：
|
||
1. 取 pseudobulk DE 表（11）中 RGC2-like / RGC1-like 的 LFC
|
||
2. ISG 集合（同 12 脚本 IFN_alpha_ISG）逐基因列出：我们 LFC(S1/S2/null) + 作者 RGC-1/RGC-2 log2FC + 检出率
|
||
3. 比较 ISG 的 LFC 分布 vs 同组全基因背景分布：
|
||
ISG 显著负于背景 → 通路特异抑制；与背景一致 → 全局抑制的一部分
|
||
4. 同时给出 RGC 整体 housekeeping 背景作参照
|
||
|
||
产出：output/13_isg_check/
|
||
"""
|
||
from pathlib import Path
|
||
|
||
import matplotlib
|
||
matplotlib.use("Agg")
|
||
import matplotlib.pyplot as plt
|
||
import numpy as np
|
||
import pandas as pd
|
||
|
||
OUT = Path("output/13_isg_check")
|
||
OUT.mkdir(parents=True, exist_ok=True)
|
||
plt.rcParams.update({"figure.dpi": 300, "savefig.dpi": 300, "font.size": 8})
|
||
|
||
ISG = ("Isg15/Ifit1/Ifit2/Ifit3/Ifit3b/Ifi27/Ifi27l2a/Ifi44/Ifi47/Ifitm1/Ifitm2/Ifitm3/Mx1/Mx2/"
|
||
"Oas1a/Oas1b/Oas2/Oasl1/Oasl2/Rsad2/Bst2/Usp18/Irf7/Irf9/Stat1/Stat2/Ddx58/Ifih1/Ddx60/"
|
||
"Cmpk2/Nlrc5/Psmb8/Psmb9/Tap1/Tap2/B2m/Gbp2/Gbp3/Gbp4/Gbp5/Irgm1/Trim21/Xaf1/Ifi35/"
|
||
"Parp9/Parp14/Samd9l/Eif2ak2/Zbp1/Herc6/Rtp4/Lgals3bp/Ccl5").split("/")
|
||
HK = ("Actb/Gapdh/B2m/Hprt/Rpl13a/Rplp0/Tbp/Gusb/Ppia/Rps18/Eef1a1/Sdha/Hmbs/Ywhaz/Tfrc/Pgk1").split("/")
|
||
|
||
xl = "ref/Data files/S2. snRNAseq-DEG_All_Celltype.xlsx"
|
||
ar1 = pd.read_excel(xl, sheet_name="RGC-1"); ar1.columns = [c.lower() for c in ar1.columns]
|
||
ar2 = pd.read_excel(xl, sheet_name="RGC-2"); ar2.columns = [c.lower() for c in ar2.columns]
|
||
ar2 = ar2.set_index("gene"); ar1 = ar1.set_index("gene")
|
||
|
||
rows = []
|
||
for grp, author in [("RGC2-like", ar2), ("RGC1-like", ar1)]:
|
||
de = pd.read_csv(f"output/11_pseudobulk_de/de_{grp}.csv").set_index("gene")
|
||
bg = de["lfc_S1"].add(de["lfc_S2"]) / 2
|
||
for g in ISG:
|
||
if g not in de.index:
|
||
rows.append({"group": grp, "gene": g, "present": False})
|
||
continue
|
||
r = de.loc[g]
|
||
a = author.loc[g] if g in author.index else None
|
||
rows.append({
|
||
"group": grp, "gene": g, "present": True,
|
||
"our_lfc_S1": round(r["lfc_S1"], 3), "our_lfc_S2": round(r["lfc_S2"], 3),
|
||
"our_lfc_null": round(r["lfc_null"], 3), "consistent": r["direction"],
|
||
"cpm_WT": round(r["cpm_WT"], 1),
|
||
"author_log2fc": round(a["avg_log2fc"], 3) if a is not None else np.nan,
|
||
"author_p": a["p_val"] if a is not None else np.nan,
|
||
"bg_median_lfc": round(bg.median(), 3),
|
||
})
|
||
isg_df = pd.DataFrame(rows)
|
||
isg_df.to_csv(OUT / "isg_genelevel_rgc.csv", index=False)
|
||
|
||
# ---- 分布比较：ISG vs 背景（Mann-Whitney，单侧 less） ----
|
||
from scipy.stats import mannwhitneyu
|
||
stat_rows = []
|
||
fig, axes = plt.subplots(1, 2, figsize=(8, 3.4))
|
||
for ax, grp in zip(axes, ["RGC2-like", "RGC1-like"]):
|
||
de = pd.read_csv(f"output/11_pseudobulk_de/de_{grp}.csv")
|
||
de["our_mean"] = (de["lfc_S1"] + de["lfc_S2"]) / 2
|
||
isg_v = de.loc[de["gene"].isin(ISG), "our_mean"].dropna()
|
||
hk_v = de.loc[de["gene"].isin(HK), "our_mean"].dropna()
|
||
bg_v = de.loc[~de["gene"].isin(ISG), "our_mean"].dropna()
|
||
u_isg, p_isg = mannwhitneyu(isg_v, bg_v, alternative="less")
|
||
u_hk, p_hk = mannwhitneyu(hk_v, bg_v, alternative="two-sided")
|
||
stat_rows.append({"group": grp, "n_ISG": len(isg_v), "ISG_median": round(isg_v.median(), 3),
|
||
"bg_median": round(bg_v.median(), 3), "HK_median": round(hk_v.median(), 3),
|
||
"ISG_less_than_bg_p": f"{p_isg:.2e}", "HK_vs_bg_p": f"{p_hk:.2e}"})
|
||
parts = ax.violinplot([bg_v.sample(min(3000, len(bg_v)), random_state=0), isg_v, hk_v],
|
||
showmedians=True)
|
||
for b, c in zip(parts["bodies"], ["lightgrey", "firebrick", "steelblue"]):
|
||
b.set_facecolor(c); b.set_alpha(0.6)
|
||
ax.set_xticks([1, 2, 3], ["background", f"ISG (n={len(isg_v)})", f"HK (n={len(hk_v)})"])
|
||
ax.axhline(0, color="grey", lw=0.5, ls="--")
|
||
ax.set_title(f"{grp}: ISG vs background LFC\np(ISG<bg)={p_isg:.1e}")
|
||
ax.set_ylabel("pseudobulk LFC (mean of S1,S2)")
|
||
pd.DataFrame(stat_rows).to_csv(OUT / "isg_vs_background_stats.csv", index=False)
|
||
print(pd.DataFrame(stat_rows).to_string(index=False))
|
||
|
||
fig.tight_layout()
|
||
fig.savefig(OUT / "isg_vs_background.png", bbox_inches="tight")
|
||
|
||
# 打印 ISG 明细（两样本一致下调的）
|
||
print("\nISG 明细（RGC2-like，|LFC|>=0.3 且两样本同向）:")
|
||
sub = isg_df[(isg_df["group"] == "RGC2-like") & isg_df["present"].fillna(False)].copy()
|
||
sub = sub[(sub["our_lfc_S1"].abs() >= 0.3) & (np.sign(sub["our_lfc_S1"]) == np.sign(sub["our_lfc_S2"]))]
|
||
print(sub[["gene", "our_lfc_S1", "our_lfc_S2", "our_lfc_null", "cpm_WT", "author_log2fc", "bg_median_lfc"]].to_string(index=False))
|
||
print("Done ->", OUT)
|