139 lines
7.2 KiB
Python
139 lines
7.2 KiB
Python
#!/usr/bin/env python
|
||
# -*- coding: utf-8 -*-
|
||
"""
|
||
10_rgc2_signature.py — P1-② RGC-2 重定义：作者 DEG 签名打分法（替代"WT ETC top30%"循环定义）
|
||
|
||
签名构建（数据文件 S2，仅用 WT 检出率，不用疾病 DE 信息）：
|
||
RGC-2 身份基因 = pct.WT(RGC-2) >= 0.5 且 pct.WT(RGC-2) - pct.WT(RGC-1) >= 0.25，按差值取 top 100。
|
||
RGC-1 签名对称构建。分类器 = score(RGC2签名) - score(RGC1签名)，>0 判为 RGC-2-like。
|
||
⚠️ 已知局限：RGC-2 的身份基因（微管/ETC/核糖体）在突变 RGC-2 中被强烈下调，
|
||
无法剔除疾病 DE 基因（剔了签名只剩 2 个基因）——突变细胞的分类因此偏向"保留身份"的较健康细胞，
|
||
会低估效应量，解读时注意。
|
||
"""
|
||
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/10_rgc2_signature")
|
||
OUT.mkdir(parents=True, exist_ok=True)
|
||
plt.rcParams.update({"figure.dpi": 300, "savefig.dpi": 300, "font.size": 8})
|
||
|
||
# ---- 1. 从作者 S2 构建身份签名（pct.WT only，不用 DE 过滤） ----
|
||
r1 = pd.read_excel("ref/Data files/S2. snRNAseq-DEG_All_Celltype.xlsx", sheet_name="RGC-1")
|
||
r2 = pd.read_excel("ref/Data files/S2. snRNAseq-DEG_All_Celltype.xlsx", sheet_name="RGC-2")
|
||
for df in (r1, r2):
|
||
df.columns = [c.lower() for c in df.columns]
|
||
m = r1.merge(r2, on="gene", suffixes=("_r1", "_r2"))
|
||
m["pct_diff_r2"] = m["pct.wt_r2"] - m["pct.wt_r1"]
|
||
|
||
sig2_df = m[(m["pct.wt_r2"] >= 0.5) & (m["pct_diff_r2"] >= 0.25)].nlargest(100, "pct_diff_r2")
|
||
sig1_df = m[(m["pct.wt_r1"] >= 0.5) & (m["pct_diff_r2"] <= -0.25)].nsmallest(100, "pct_diff_r2")
|
||
|
||
adata = sc.read_h5ad("data/04_scored.h5ad")
|
||
sig2 = [g for g in sig2_df["gene"] if g in adata.raw.var_names]
|
||
sig1 = [g for g in sig1_df["gene"] if g in adata.raw.var_names]
|
||
print(f"RGC-2 签名: {len(sig2_df)} -> 数据中 {len(sig2)} 基因")
|
||
print(f"RGC-1 签名: {len(sig1_df)} -> 数据中 {len(sig1)} 基因")
|
||
sig2_df.to_csv(OUT / "authors_RG2_signature.csv", index=False)
|
||
sig1_df.to_csv(OUT / "authors_RG1_signature.csv", index=False)
|
||
|
||
# ---- 2. 在全部 RGC 上打分 ----
|
||
raw = adata.raw.to_adata()
|
||
raw.obs = adata.obs.copy()
|
||
rgc_mask = adata.obs["cell_type"] == "RGC"
|
||
rgc = raw[rgc_mask].copy()
|
||
sc.tl.score_genes(rgc, gene_list=sig2, score_name="sig2", use_raw=False, random_state=0)
|
||
sc.tl.score_genes(rgc, gene_list=sig1, score_name="sig1", use_raw=False, random_state=0)
|
||
rgc.obs["sig_diff"] = rgc.obs["sig2"] - rgc.obs["sig1"]
|
||
rgc.obs["rgc2_sig_group"] = np.where(rgc.obs["sig_diff"] > 0, "RGC2-like", "RGC1-like")
|
||
|
||
adata.obs["score_RG2_signature"] = np.nan
|
||
adata.obs["score_RG1_signature"] = np.nan
|
||
adata.obs["rgc2_sig_group"] = "non-RGC"
|
||
adata.obs.loc[rgc.obs_names, "score_RG2_signature"] = rgc.obs["sig2"]
|
||
adata.obs.loc[rgc.obs_names, "score_RG1_signature"] = rgc.obs["sig1"]
|
||
adata.obs.loc[rgc.obs_names, "rgc2_sig_group"] = rgc.obs["rgc2_sig_group"].values
|
||
|
||
print("\n签名分组核数:")
|
||
print(rgc.obs.groupby(["rgc2_sig_group", "genotype"], observed=True).size().unstack(fill_value=0))
|
||
print("\n与 P0 ETC 分组重叠:")
|
||
print(pd.crosstab(rgc.obs["rgc2_sig_group"], adata.obs.loc[rgc.obs_names, "rgc_group"]))
|
||
|
||
# ---- 3. 新分组下的模块分比较（复用 04_scored 中的 score_ 列） ----
|
||
mods = ["score_ETC_WP295", "score_CI_biogenesis_RE", "score_Glycolysis_WP157",
|
||
"score_Ribosomal_WP163", "score_Mitophagy_Autophagy"]
|
||
rows = []
|
||
for grp in ["RGC2-like", "RGC1-like"]:
|
||
gm = adata.obs["rgc2_sig_group"] == grp
|
||
for col in mods:
|
||
wt = adata.obs.loc[gm & (adata.obs["genotype"] == "WT"), col]
|
||
per_s = adata.obs.loc[gm].groupby("sample")[col].median()
|
||
s1, s2 = per_s.get("Opa1V291D_S1", np.nan), per_s.get("Opa1V291D_S2", np.nan)
|
||
rows.append({"group": grp, "module": col.replace("score_", ""),
|
||
"n_WT": int((gm & (adata.obs["genotype"] == "WT")).sum()),
|
||
"median_WT": round(wt.median(), 4),
|
||
"delta_S1": round(s1 - wt.median(), 4), "delta_S2": round(s2 - wt.median(), 4)})
|
||
mod_res = pd.DataFrame(rows)
|
||
mod_res.to_csv(OUT / "module_scores_newgrouping.csv", index=False)
|
||
print("\n", mod_res.to_string(index=False))
|
||
|
||
# ---- 4. 基因级一致率（作者 RGC-2 显著下调基因 vs 我们 RGC2-like） ----
|
||
s2_sig = r2[(r2["p_val"] < 0.001) & (r2["avg_log2fc"] < -0.5)]
|
||
rawX = raw.X
|
||
vidx = {g: i for i, g in enumerate(raw.var_names)}
|
||
gm = adata.obs["rgc2_sig_group"] == "RGC2-like"
|
||
wt_m = (gm & (adata.obs["genotype"] == "WT")).values
|
||
mu_m = (gm & (adata.obs["genotype"] == "V291D")).values
|
||
print(f"\nRGC2-like: WT n={wt_m.sum()}, V291D n={mu_m.sum()}")
|
||
gene_rows = []
|
||
for g in s2_sig["gene"]:
|
||
if g not in vidx:
|
||
continue
|
||
i = vidx[g]
|
||
c_wt = rawX[wt_m, i].mean()
|
||
c_mu = rawX[mu_m, i].mean()
|
||
gene_rows.append({"gene": g, "author_log2fc": r2.loc[r2["gene"] == g, "avg_log2fc"].iloc[0],
|
||
"our_log2fc_approx": np.log2((c_mu + 1e-9) / (c_wt + 1e-9))})
|
||
gd = pd.DataFrame(gene_rows)
|
||
if len(gd):
|
||
concord = (np.sign(gd["author_log2fc"]) == np.sign(gd["our_log2fc_approx"])).mean()
|
||
gd.to_csv(OUT / "genelevel_concordance_RGC2like.csv", index=False)
|
||
print(f"基因级方向一致率（作者 RGC-2 下调基因 vs 我们 RGC2-like）: {concord*100:.1f}% (n={len(gd)})")
|
||
|
||
# ---- 5. 图 ----
|
||
fig, axes = plt.subplots(1, 3, figsize=(12, 3.6))
|
||
# (a) sig_diff 分布按基因型
|
||
for geno, c in [("WT", "black"), ("V291D", "firebrick")]:
|
||
v = rgc.obs.loc[rgc.obs["genotype"] == geno, "sig_diff"]
|
||
axes[0].hist(v, bins=50, histtype="step", density=True, color=c, lw=1.4, label=f"{geno} (n={len(v)})")
|
||
axes[0].axvline(0, color="grey", ls="--", lw=0.8)
|
||
axes[0].set_xlabel("score(RGC-2 sig) - score(RGC-1 sig)")
|
||
axes[0].set_ylabel("density"); axes[0].legend(fontsize=7, frameon=False)
|
||
axes[0].set_title("Signature score distribution in our RGCs")
|
||
# (b) 新旧分组重叠（百分比堆叠）
|
||
ct = pd.crosstab(rgc.obs["rgc2_sig_group"], adata.obs.loc[rgc.obs_names, "rgc_group"], normalize="index")
|
||
ct.plot(kind="bar", stacked=True, ax=axes[1], color=["#d62728", "#ff9896", "lightgrey"])
|
||
axes[1].set_ylabel("fraction"); axes[1].set_title("New signature group vs P0 ETC group")
|
||
axes[1].legend(fontsize=6, frameon=False); axes[1].set_xticklabels(axes[1].get_xticklabels(), rotation=0)
|
||
# (c) 新分组模块分 delta
|
||
x = np.arange(len(mods)); w = 0.35
|
||
for off, grp, c in [(-w/2, "RGC2-like", "#d62728"), (w/2, "RGC1-like", "#ff9896")]:
|
||
sub = mod_res[mod_res["group"] == grp].set_index("module")
|
||
d1 = sub["delta_S1"].reindex([m.replace("score_", "") for m in mods])
|
||
d2 = sub["delta_S2"].reindex([m.replace("score_", "") for m in mods])
|
||
axes[2].bar(x + off - w/4, d1, width=w/2, color=c, alpha=0.6, label=f"{grp} S1")
|
||
axes[2].bar(x + off + w/4, d2, width=w/2, color=c, alpha=1.0, label=f"{grp} S2")
|
||
axes[2].set_xticks(x, [m.replace("score_", "").replace("_", "\n") for m in mods], fontsize=6)
|
||
axes[2].axhline(0, color="grey", lw=0.6); axes[2].set_ylabel("delta median score vs WT")
|
||
axes[2].legend(fontsize=5, frameon=False); axes[2].set_title("Metabolism modules in new groups")
|
||
fig.tight_layout()
|
||
fig.savefig(OUT / "rgc2_signature_redefinition.png", bbox_inches="tight")
|
||
|
||
adata.write_h5ad("data/05_rgc2sig.h5ad")
|
||
print("Done ->", OUT, "+ data/05_rgc2sig.h5ad")
|