151 lines
6.7 KiB
Python
151 lines
6.7 KiB
Python
#!/usr/bin/env python
|
||
# -*- coding: utf-8 -*-
|
||
"""
|
||
17_ifn_heatmap_ma.py — 补齐 Step 5/6 缺失的可视化
|
||
|
||
1. IFN–JAK–STAT 逐基因热图：Müller/RGC2-like/RGC1-like 用 pseudobulk LFC，
|
||
microglia 用检出率差（Δdet = det_V291D − det_WT，深度已说明仅方向性）。
|
||
2. MA 图：Müller 与 RGC2-like 的 LFC 均值 vs 平均表达，一致基因标红、关键基因标名。
|
||
3. 细胞级 Wilcoxon（探索层）对照：展示细胞级 DE 的 UP/DOWN 偏斜（被深度污染，
|
||
类作者 artifact）vs pseudobulk 的结构化改变。
|
||
产出：output/17_visualization/
|
||
"""
|
||
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
|
||
import scipy.sparse as sp
|
||
|
||
OUT = Path("output/17_visualization")
|
||
OUT.mkdir(parents=True, exist_ok=True)
|
||
plt.rcParams.update({"figure.dpi": 300, "savefig.dpi": 300, "font.size": 8})
|
||
|
||
IFN_GENES = ["Ifnar1", "Ifnar2", "Jak1", "Jak2", "Stat1", "Stat2", "Stat3", "Irf7", "Irf9",
|
||
"Isg15", "Ifit1", "Ifit2", "Ifit3", "Ifitm3", "Mx1", "Oasl2", "Bst2", "Usp18",
|
||
"Ccl5", "Gbp2", "Gbp4", "Icam1", "Cxcl10", "Il6", "C1qa", "C3"]
|
||
|
||
# ---- 1. IFN 逐基因热图 ----
|
||
groups_pb = ["Muller", "RGC2-like", "RGC1-like"]
|
||
lfc_mat = {}
|
||
for grp in groups_pb:
|
||
de = pd.read_csv(f"output/11_pseudobulk_de/de_{grp}.csv").set_index("gene")
|
||
lfc_mat[grp] = de["lfc_S1"].add(de["lfc_S2"]) / 2 # 均值 LFC
|
||
|
||
det = pd.read_csv("output/12_inflammation/key_gene_detection_rates.csv")
|
||
det_mg = det[det["group"] == "Microglia"].set_index("gene")
|
||
det_delta = det_mg["det_V291D"] - det_mg["det_WT"] # 检出率差
|
||
|
||
genes = [g for g in IFN_GENES if g in lfc_mat["Muller"].index]
|
||
mat = pd.DataFrame({grp: lfc_mat[grp].reindex(genes) for grp in groups_pb})
|
||
mat["Microglia_Δdet"] = det_delta.reindex(genes)
|
||
|
||
fig, ax = plt.subplots(figsize=(5.5, 7))
|
||
vmax = max(2.5, np.nanmax(np.abs(mat[groups_pb].values)))
|
||
im = ax.imshow(mat[groups_pb].values, cmap="RdBu_r", vmin=-vmax, vmax=vmax, aspect="auto")
|
||
ax.set_xticks(range(len(groups_pb)), [g.replace("-like", "") for g in groups_pb])
|
||
ax.set_yticks(range(len(genes)), genes, fontsize=7)
|
||
for i in range(len(genes)):
|
||
for j in range(len(groups_pb)):
|
||
v = mat[groups_pb].values[i, j]
|
||
if np.isnan(v):
|
||
ax.text(j, i, "n.d.", ha="center", va="center", fontsize=5, color="grey")
|
||
continue
|
||
ax.text(j, i, f"{v:+.1f}", ha="center", va="center", fontsize=6,
|
||
color="white" if abs(v) > 0.7 * vmax else "black")
|
||
ax.set_title("IFN-JAK-STAT gene-level LFC (V291D - WT, pseudobulk)")
|
||
fig.colorbar(im, ax=ax, label="mean LFC (S1,S2)")
|
||
fig.tight_layout()
|
||
fig.savefig(OUT / "ifn_gene_heatmap.png", bbox_inches="tight")
|
||
plt.close(fig)
|
||
|
||
# microglia 检出率差单独一列
|
||
fig, ax = plt.subplots(figsize=(2.4, 7))
|
||
mg_genes = [g for g in IFN_GENES if g in det_delta.index]
|
||
v = det_delta.reindex(mg_genes).values
|
||
im = ax.imshow(v[:, None], cmap="RdBu_r", vmin=-0.6, vmax=0.6, aspect="auto")
|
||
ax.set_yticks(range(len(mg_genes)), mg_genes, fontsize=7)
|
||
ax.set_xticks([])
|
||
for i, g in enumerate(mg_genes):
|
||
ax.text(0, i, f"{det_delta[g]:+.2f}", ha="center", va="center", fontsize=6)
|
||
ax.set_title("Microglia detection-rate diff\n(WT n=8 nuclei, direction only)")
|
||
fig.colorbar(im, ax=ax, label="Δdet")
|
||
fig.tight_layout()
|
||
fig.savefig(OUT / "ifn_gene_heatmap_microglia.png", bbox_inches="tight")
|
||
plt.close(fig)
|
||
|
||
# ---- 2. MA 图 ----
|
||
def ma_plot(grp, ax):
|
||
de = pd.read_csv(f"output/11_pseudobulk_de/de_{grp}.csv")
|
||
de["lfc_mean"] = (de["lfc_S1"] + de["lfc_S2"]) / 2
|
||
de["log_cpm"] = np.log10((de["cpm_WT"] + de["cpm_S1"] + de["cpm_S2"]) / 3 + 1)
|
||
ax.scatter(de["log_cpm"], de["lfc_mean"], s=1.5, alpha=0.2, color="grey", rasterized=True)
|
||
cons = de[de["direction"].isin(["UP", "DOWN"])]
|
||
ax.scatter(cons["log_cpm"], cons["lfc_mean"], s=3, alpha=0.5,
|
||
c=["firebrick" if d == "UP" else "steelblue" for d in cons["direction"]], rasterized=True)
|
||
# 标关键基因
|
||
key = ["Apoe", "Clu", "Gfap", "Serpina3n", "Lcn2", "Vim", "Spp1", "Timp1", "Cd44",
|
||
"Ndufs3", "Cox7b", "Uqcrc2", "Gapdh", "Pgk1", "Ifit2", "Cmpk2", "Stat1"]
|
||
for g in key:
|
||
if g in de.set_index("gene").index:
|
||
r = de[de["gene"] == g].iloc[0]
|
||
ax.annotate(g, (r["log_cpm"], r["lfc_mean"]), fontsize=6, alpha=0.85)
|
||
ax.axhline(0, color="k", lw=0.5)
|
||
ax.set_xlabel("log10(mean CPM)")
|
||
ax.set_ylabel("mean LFC (V291D − WT)")
|
||
ax.set_title(f"{grp} MA plot (red=consistent UP, blue=DOWN)")
|
||
|
||
fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))
|
||
ma_plot("Muller", axes[0])
|
||
ma_plot("RGC2-like", axes[1])
|
||
fig.tight_layout()
|
||
fig.savefig(OUT / "ma_plot.png", bbox_inches="tight")
|
||
plt.close(fig)
|
||
|
||
# ---- 3. 细胞级 Wilcoxon（探索层）vs pseudobulk 对照 ----
|
||
adata = sc.read_h5ad("data/06_inflammation.h5ad")
|
||
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)
|
||
|
||
rows = []
|
||
for grp in ["Muller", "RGC2-like", "RGC1-like", "Rod"]:
|
||
sub = adata[obs["cmp_group"] == grp].copy()
|
||
# 细胞级 Wilcoxon（探索层，明确 pseudoreplication + 深度污染局限）
|
||
sc.tl.rank_genes_groups(sub, groupby="genotype", groups=["V291D"], reference="WT",
|
||
method="wilcoxon", use_raw=True)
|
||
res = sc.get.rank_genes_groups_df(sub, group="V291D")
|
||
res = res[res["pvals_adj"] < 0.05]
|
||
n_up = (res["logfoldchanges"] > 0.25).sum()
|
||
n_dn = (res["logfoldchanges"] < -0.25).sum()
|
||
# pseudobulk 对照
|
||
de = pd.read_csv(f"output/11_pseudobulk_de/de_{grp}.csv")
|
||
pb_up = (de["direction"] == "UP").sum()
|
||
pb_dn = (de["direction"] == "DOWN").sum()
|
||
rows.append({"group": grp, "cell_UP": int(n_up), "cell_DOWN": int(n_dn),
|
||
"cell_ratio": round(n_up / max(n_dn, 1), 2),
|
||
"pseudo_UP": int(pb_up), "pseudo_DOWN": int(pb_dn),
|
||
"pseudo_ratio": round(pb_up / max(pb_dn, 1), 2)})
|
||
cw = pd.DataFrame(rows)
|
||
cw.to_csv(OUT / "celllevel_vs_pseudobulk.csv", index=False)
|
||
print("\n细胞级 Wilcoxon（探索层）vs pseudobulk UP/DOWN 对照：")
|
||
print(cw.to_string(index=False))
|
||
|
||
fig, ax = plt.subplots(figsize=(6, 4))
|
||
x = np.arange(len(cw))
|
||
w = 0.35
|
||
ax.bar(x - w/2, np.log10(cw["cell_ratio"]), width=w, color="grey", alpha=0.7, label="cell-level Wilcoxon (depth-biased)")
|
||
ax.bar(x + w/2, np.log10(cw["pseudo_ratio"]), width=w, color="teal", label="pseudobulk median-of-ratios")
|
||
ax.axhline(0, color="k", lw=0.5)
|
||
ax.set_xticks(x, cw["group"])
|
||
ax.set_ylabel("log10(UP/DOWN)")
|
||
ax.set_title("Cell-level vs pseudobulk DE direction bias (grey = depth-biased, author-like artifact)")
|
||
ax.legend(fontsize=7, frameon=False)
|
||
fig.tight_layout()
|
||
fig.savefig(OUT / "celllevel_vs_pseudobulk.png", bbox_inches="tight")
|
||
print("Done ->", OUT)
|