147 lines
7.3 KiB
Python
147 lines
7.3 KiB
Python
#!/usr/bin/env python
|
||
# -*- coding: utf-8 -*-
|
||
"""
|
||
12_inflammation_modules.py — P1-④ 炎症模块打分（胶质为主，决策记录 第二轮裁决 #3）
|
||
|
||
模块（小鼠基因符号，手工策展，以 Hallmark/Reactome 标准集为骨干）：
|
||
IFN_alpha_ISG / IFN_gamma / IL6_JAK_STAT3 / TNFA_NFkB / Complement /
|
||
cGAS_STING / NLRP3_inflammasome / UPRmt / ISR /
|
||
Muller_reactive_gliosis / Microglia_DAM / Microglia_homeostatic
|
||
重点组：Muller、Microglia、RGC1-like、RGC2-like；Rod/Cone 作参照。
|
||
Microglia WT 仅 8 核 → 仅描述性（每样本中位数 + 检出率），不做推断。
|
||
|
||
产出：output/12_inflammation/ + data/06_inflammation.h5ad
|
||
"""
|
||
from pathlib import Path
|
||
|
||
import matplotlib
|
||
matplotlib.use("Agg")
|
||
import matplotlib.pyplot as plt
|
||
import numpy as np
|
||
import pandas as pd
|
||
import scanpy as sc
|
||
|
||
OUT = Path("output/12_inflammation")
|
||
OUT.mkdir(parents=True, exist_ok=True)
|
||
plt.rcParams.update({"figure.dpi": 300, "savefig.dpi": 300, "font.size": 8})
|
||
|
||
GENESETS = {
|
||
"IFN_alpha_ISG": "Isg15/Ifit1/Ifit2/Ifit3/Ifit3b/Ifi27/Ifi27l2a/Ifi44/Ifi47/Ifitm1/Ifitm2/Ifitm3/Mx1/Mx2/Oas1a/Oas1b/Oas2/Oas3/Oasl1/Oasl2/Oas1g/Rsad2/Bst2/Usp18/Irf7/Irf9/Stat1/Stat2/Cxcl10/Cxcl9/Cxcl11/Ddx58/Ifih1/Ddx60/Cmpk2/Nlrc5/Psmb8/Psmb9/Psmb10/Tap1/Tap2/Tapbp/B2m/Gbp2/Gbp3/Gbp4/Gbp5/Gbp7/Iigp1/Irgm1/Trim21/Xaf1/Ifi35/Parp9/Parp14/Samd9l/Eif2ak2/Zbp1/Herc6/Rtp4/Lgals3bp/Ccl5",
|
||
"IFN_gamma": "Jak1/Jak2/Stat1/Irf1/Irf8/Ciita/Cd74/H2-Aa/H2-Ab1/H2-Eb1/H2-K1/H2-D1/Icam1/Gbp2/Gbp3/Gbp4/Gbp5/Gbp7/Ifi47/Socs1/Cxcl9/Cxcl10/Cxcl11/Psmb9/Tap1/Tap2",
|
||
"IL6_JAK_STAT3": "Il6/Il6ra/Il6st/Jak1/Jak2/Tyk2/Stat3/Socs3/Nfkbia/Nfkbiz/Tnfaip3/Cxcl1/Ccl2/Icam1/Timp1/Serpina3n/Lcn2/Osmr/Saa1/Saa2/Saa3",
|
||
"TNFA_NFkB": "Tnf/Tnfrsf1a/Tnfrsf1b/Nfkb1/Nfkb2/Rela/Relb/Nfkbia/Nfkbib/Nfkbiz/Tnfaip3/Tnfaip2/Cxcl1/Cxcl2/Cxcl10/Ccl2/Ccl5/Ccl20/Il1b/Il1a/Ptgs2/Icam1/Vcam1/Sod2/Jun/Fos",
|
||
"Complement": "C1qa/C1qb/C1qc/C2/C3/C4b/Cfb/Cfd/Cfh/Cfi/Serping1/Cd55/Cd59a/Cd59b/Cr1l/C3ar1/C5ar1/Itgam/Itgax/Itgb2/Masp1/Masp2/Mbl2/C1ra/C1s1/Clu/Vtn",
|
||
"cGAS_STING": "Mb21d1/Tmem173/Tbk1/Ikbke/Irf3/Ddx41/Ifnb1",
|
||
"NLRP3_inflammasome": "Nlrp3/Pycard/Casp1/Il1b/Il18/Gsdmd/Aim2/Nek7/Txnip",
|
||
"UPRmt": "Atf5/Hspa9/Hspe1/Hspd1/Dnaja3/Lonp1/Clpp/Yme1l1/Afg3l2",
|
||
"ISR": "Atf4/Ddit3/Atf3/Asns/Trib3/Ppp1r15a/Chac1",
|
||
"Muller_reactive_gliosis": "Gfap/Vim/Serpina3n/Apoe/Lcn2/S100a6/Cd44/Timp1/Nes/Clu/Osmr/Ccl2",
|
||
"Microglia_DAM": "Apoe/Spp1/Lpl/Trem2/Tyrobp/Axl/Clec7a/Cst7/Itgax/Cd9/Lgals3/Lilrb4a/Ccl6/Fabp5/Ctsb/Ctsd/Gpnmb",
|
||
"Microglia_homeostatic": "P2ry12/P2ry13/Tmem119/Hexb/Sall1/Cx3cr1/Siglech/Mafb/Selplg/Fcrls/Olfml3/Tgfbr1/Gpr34/Slc2a5",
|
||
}
|
||
|
||
adata = sc.read_h5ad("data/05_rgc2sig.h5ad")
|
||
raw = adata.raw.to_adata()
|
||
raw.obs = adata.obs.copy()
|
||
|
||
gs = {k: [g for g in v.split("/") if g in raw.var_names] for k, v in GENESETS.items()}
|
||
for k, v in gs.items():
|
||
print(f"{k}: {len(v)}/{len(GENESETS[k].split('/'))} genes in data")
|
||
miss = set(GENESETS[k].split("/")) - set(v)
|
||
if miss:
|
||
print(f" missing: {sorted(miss)}")
|
||
|
||
for k, v in gs.items():
|
||
sc.tl.score_genes(raw, gene_list=v, score_name=f"score_{k}", use_raw=False, random_state=0)
|
||
adata.obs[f"score_{k}"] = raw.obs[f"score_{k}"]
|
||
|
||
# ---- 分组定义 ----
|
||
obs = adata.obs
|
||
obs["cmp_group"] = obs["cell_type"].astype(str)
|
||
is_rgc = obs["cell_type"] == "RGC"
|
||
obs.loc[is_rgc, "cmp_group"] = obs.loc[is_rgc, "rgc2_sig_group"].astype(str)
|
||
FOCUS = ["RGC2-like", "RGC1-like", "Muller", "Microglia", "Astrocyte", "Rod", "Cone"]
|
||
|
||
# ---- 每样本中位数 + delta 表 ----
|
||
rows = []
|
||
for grp in FOCUS:
|
||
m = obs["cmp_group"] == grp
|
||
for k in gs:
|
||
col = f"score_{k}"
|
||
per_s = obs.loc[m].groupby("sample")[col].median()
|
||
wt = per_s.get("WT", np.nan)
|
||
rows.append({
|
||
"group": grp, "module": k,
|
||
"n_WT": int((m & (obs["genotype"] == "WT")).sum()),
|
||
"median_WT": round(wt, 4),
|
||
"median_S1": round(per_s.get("Opa1V291D_S1", np.nan), 4),
|
||
"median_S2": round(per_s.get("Opa1V291D_S2", np.nan), 4),
|
||
"delta_S1": round(per_s.get("Opa1V291D_S1", np.nan) - wt, 4),
|
||
"delta_S2": round(per_s.get("Opa1V291D_S2", np.nan) - wt, 4),
|
||
})
|
||
res = pd.DataFrame(rows)
|
||
res.to_csv(OUT / "inflammation_module_scores.csv", index=False)
|
||
print("\n", res[res["group"].isin(["Muller", "Microglia", "RGC2-like"])].to_string(index=False))
|
||
|
||
# ---- 关键基因检出率（每 组×基因型 检出比例） ----
|
||
KEY = ["Ifnb1", "Ifnar1", "Ifnar2", "Jak1", "Jak2", "Stat1", "Stat2", "Stat3", "Irf7", "Irf9",
|
||
"Isg15", "Ifit1", "Ifit2", "Ifit3", "Mx1", "Cxcl10", "Il6", "Tnf", "Il1b", "C3", "C1qa",
|
||
"Mb21d1", "Tmem173", "Nlrp3", "Gfap", "Serpina3n", "Apoe", "Lcn2", "Spp1", "Trem2"]
|
||
rawX = raw.X
|
||
vidx = {g: i for i, g in enumerate(raw.var_names)}
|
||
det_rows = []
|
||
for grp in FOCUS:
|
||
m = obs["cmp_group"] == grp
|
||
for g in KEY:
|
||
if g not in vidx:
|
||
continue
|
||
col = rawX[:, vidx[g]]
|
||
r = {"group": grp, "gene": g}
|
||
for geno in ["WT", "V291D"]:
|
||
mm = (m & (obs["genotype"] == geno)).values
|
||
r[f"det_{geno}"] = round((col[mm].toarray().ravel() > 0).mean(), 4) if mm.sum() else np.nan
|
||
r[f"n_{geno}"] = int(mm.sum())
|
||
det_rows.append(r)
|
||
det = pd.DataFrame(det_rows)
|
||
det.to_csv(OUT / "key_gene_detection_rates.csv", index=False)
|
||
|
||
# ---- 图 1：模块 delta 热图 ----
|
||
piv = res.pivot(index="module", columns="group", values="delta_S1")
|
||
piv2 = res.pivot(index="module", columns="group", values="delta_S2")
|
||
mean_delta = (piv + piv2) / 2
|
||
mean_delta = mean_delta[[c for c in FOCUS if c in mean_delta.columns]]
|
||
fig, ax = plt.subplots(figsize=(6.5, 4.2))
|
||
vmax = np.nanmax(np.abs(mean_delta.values))
|
||
im = ax.imshow(mean_delta.values, cmap="RdBu_r", vmin=-vmax, vmax=vmax, aspect="auto")
|
||
ax.set_xticks(range(len(mean_delta.columns)), mean_delta.columns, rotation=45, ha="right")
|
||
ax.set_yticks(range(len(mean_delta.index)), mean_delta.index, fontsize=7)
|
||
for i in range(mean_delta.shape[0]):
|
||
for j in range(mean_delta.shape[1]):
|
||
v = mean_delta.values[i, j]
|
||
ax.text(j, i, f"{v:+.2f}", ha="center", va="center", fontsize=6,
|
||
color="white" if abs(v) > 0.6 * vmax else "black")
|
||
fig.colorbar(im, label="mean delta module score (V291D - WT)")
|
||
ax.set_title("Inflammation module scores: mutant - WT (mean of S1,S2 deltas)")
|
||
fig.tight_layout()
|
||
fig.savefig(OUT / "inflammation_heatmap.png", bbox_inches="tight")
|
||
|
||
# ---- 图 2：Muller 反应性胶质化 + 补体 分布（每样本） ----
|
||
fig, axes = plt.subplots(1, 3, figsize=(10, 3.2))
|
||
for ax, k in zip(axes, ["Muller_reactive_gliosis", "Complement", "IFN_alpha_ISG"]):
|
||
col = f"score_{k}"
|
||
for grp, color in [("Muller", "#1f77b4"), ("Microglia", "#2ca02c")]:
|
||
for s, ls in [("WT", "-"), ("Opa1V291D_S1", "--"), ("Opa1V291D_S2", ":")]:
|
||
mm = (obs["cmp_group"] == grp) & (obs["sample"] == s)
|
||
v = obs.loc[mm, col]
|
||
if len(v) > 5:
|
||
ax.hist(v, bins=40, histtype="step", density=True, color=color, ls=ls, lw=1.2,
|
||
label=f"{grp} {s.replace('Opa1V291D_', 'V291D_')}")
|
||
ax.set_title(k, fontsize=9); ax.set_xlabel("module score")
|
||
axes[0].set_ylabel("density")
|
||
axes[-1].legend(fontsize=5.5, frameon=False)
|
||
fig.tight_layout()
|
||
fig.savefig(OUT / "glia_module_distributions.png", bbox_inches="tight")
|
||
|
||
adata.write_h5ad("data/06_inflammation.h5ad")
|
||
print("Done ->", OUT, "+ data/06_inflammation.h5ad")
|