113 lines
5.1 KiB
Python
113 lines
5.1 KiB
Python
#!/usr/bin/env python
|
||
# -*- coding: utf-8 -*-
|
||
"""
|
||
14_global_shift.py — P1-⑥ 全局转录抑制鉴别（计划 Step 5.5）
|
||
|
||
回答：各细胞类型在突变后是否存在"每核转录总量/检出基因数"的整体下降（全局塌陷），
|
||
还是仅特定通路改变？三条独立证据：
|
||
A. 每核 total_counts / n_genes_by_counts 分布（细胞级，不依赖 pseudobulk 样本量）
|
||
B. housekeeping 模块分（16 个经典 HK 基因，细胞级 score_genes）
|
||
C. pseudobulk DE 的 UP/DOWN 比（script 11 汇总，指示极端基因的方向偏斜）
|
||
|
||
产出：output/14_global_shift/
|
||
"""
|
||
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/14_global_shift")
|
||
OUT.mkdir(parents=True, exist_ok=True)
|
||
plt.rcParams.update({"figure.dpi": 300, "savefig.dpi": 300, "font.size": 8})
|
||
|
||
HK = ("Actb/Gapdh/B2m/Hprt/Rpl13a/Rplp0/Tbp/Gusb/Ppia/Rps18/Eef1a1/Sdha/Hmbs/Ywhaz/Tfrc/Pgk1").split("/")
|
||
|
||
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)
|
||
GROUPS = ["Rod", "Cone", "Bipolar", "Amacrine", "Muller", "RGC1-like", "RGC2-like",
|
||
"Microglia", "Pericyte", "Oligodendrocyte"]
|
||
|
||
# ---- A. 每核转录本指标 ----
|
||
rows = []
|
||
for grp in GROUPS:
|
||
m = obs["cmp_group"] == grp
|
||
for geno in ["WT", "V291D"]:
|
||
mm = m & (obs["genotype"] == geno)
|
||
if mm.sum() < 5:
|
||
continue
|
||
rows.append({"group": grp, "genotype": geno, "n": int(mm.sum()),
|
||
"median_total_counts": obs.loc[mm, "total_counts"].median(),
|
||
"median_n_genes": obs.loc[mm, "n_genes_by_counts"].median()})
|
||
a_tab = pd.DataFrame(rows)
|
||
a_piv = a_tab.pivot(index="group", columns="genotype", values="median_n_genes")
|
||
a_piv["ratio_V291D_WT"] = (a_piv["V291D"] / a_piv["WT"]).round(3)
|
||
a_cnt = a_tab.pivot(index="group", columns="genotype", values="median_total_counts")
|
||
a_cnt.columns = [f"counts_{c}" for c in a_cnt.columns]
|
||
a_tab_out = a_piv.join(a_cnt)
|
||
a_tab_out.to_csv(OUT / "pernucleus_transcript_load.csv")
|
||
print("A. 每核中位检出基因数（比值 V291D/WT）:")
|
||
print(a_tab_out.round(2).to_string())
|
||
|
||
# ---- B. HK 模块分（细胞级） ----
|
||
raw = adata.raw.to_adata()
|
||
raw.obs = obs
|
||
hk_in = [g for g in HK if g in raw.var_names]
|
||
print(f"\nB. HK 基因在数据中 {len(hk_in)}/16")
|
||
sc.tl.score_genes(raw, gene_list=hk_in, score_name="score_HK", use_raw=False, random_state=0)
|
||
obs["score_HK"] = raw.obs["score_HK"]
|
||
hk_rows = []
|
||
for grp in GROUPS:
|
||
m = obs["cmp_group"] == grp
|
||
per_g = obs.loc[m].groupby("genotype")["score_HK"].median()
|
||
per_s = obs.loc[m].groupby("sample")["score_HK"].median()
|
||
hk_rows.append({"group": grp, "n_WT": int((m & (obs["genotype"] == "WT")).sum()),
|
||
"HK_median_WT": round(per_g.get("WT", np.nan), 4),
|
||
"delta_S1": round(per_s.get("Opa1V291D_S1", np.nan) - per_g.get("WT", np.nan), 4),
|
||
"delta_S2": round(per_s.get("Opa1V291D_S2", np.nan) - per_g.get("WT", np.nan), 4)})
|
||
hk_tab = pd.DataFrame(hk_rows)
|
||
hk_tab.to_csv(OUT / "housekeeping_module_scores.csv", index=False)
|
||
print(hk_tab.to_string(index=False))
|
||
|
||
# ---- C. 汇总 UP/DOWN 比（来自 script 11） ----
|
||
summ = pd.read_csv("output/11_pseudobulk_de/de_summary.csv")
|
||
print("\nC. pseudobulk UP/DOWN 比:")
|
||
print(summ[["group", "consistent_UP", "consistent_DOWN", "UP_DOWN_ratio", "low_power"]].to_string(index=False))
|
||
|
||
# ---- 图 ----
|
||
fig, axes = plt.subplots(1, 3, figsize=(13, 3.6))
|
||
order = [g for g in GROUPS if g in a_tab_out.index]
|
||
x = np.arange(len(order))
|
||
axes[0].bar(x, a_tab_out.loc[order, "ratio_V291D_WT"], color="steelblue")
|
||
axes[0].axhline(1, color="grey", ls="--", lw=0.7)
|
||
axes[0].set_xticks(x, order, rotation=45, ha="right", fontsize=7)
|
||
axes[0].set_ylabel("median n_genes ratio (V291D / WT)")
|
||
axes[0].set_title("A. Per-nucleus detected genes")
|
||
hk_plot = hk_tab.set_index("group").loc[[g for g in order if g in hk_tab["group"].values]]
|
||
hk_order = hk_plot.index
|
||
xh = np.arange(len(hk_order))
|
||
axes[1].bar(xh - 0.2, hk_plot["delta_S1"], width=0.4, color="firebrick", alpha=0.7, label="S1-WT")
|
||
axes[1].bar(xh + 0.2, hk_plot["delta_S2"], width=0.4, color="darkorange", alpha=0.9, label="S2-WT")
|
||
axes[1].axhline(0, color="grey", ls="--", lw=0.7)
|
||
axes[1].set_xticks(xh, hk_order, rotation=45, ha="right", fontsize=7)
|
||
axes[1].set_ylabel("delta HK module score"); axes[1].legend(fontsize=7, frameon=False)
|
||
axes[1].set_title("B. Housekeeping module (16 genes)")
|
||
s2 = summ.set_index("group")
|
||
s2_order = [g for g in order if g in s2.index]
|
||
xs = np.arange(len(s2_order))
|
||
axes[2].bar(xs, np.log10(s2.loc[s2_order, "UP_DOWN_ratio"]),
|
||
color=["grey" if s2.loc[g, "low_power"] else "teal" for g in s2_order])
|
||
axes[2].axhline(0, color="grey", ls="--", lw=0.7)
|
||
axes[2].set_xticks(xs, s2_order, rotation=45, ha="right", fontsize=7)
|
||
axes[2].set_ylabel("log10(UP/DOWN) consistent genes")
|
||
axes[2].set_title("C. Pseudobulk DE direction skew (grey=low power)")
|
||
fig.tight_layout()
|
||
fig.savefig(OUT / "global_shift.png", bbox_inches="tight")
|
||
print("Done ->", OUT)
|