139 lines
6.7 KiB
Python
139 lines
6.7 KiB
Python
#!/usr/bin/env python
|
||
# -*- coding: utf-8 -*-
|
||
"""
|
||
08_replicate_metabolism.py — 阳性对照：复现原文能量代谢结论（第一轮 P0, Step 7）
|
||
|
||
基因集：作者补充表 S10 中各通路的 leading-edge 基因（geneID 列），与我们的数据取交集。
|
||
比较：RGC 整体 / RGC_highETC（≈原文 RGC-2）/ 其他主要细胞类型 × 基因型。
|
||
另加：核糖体蛋白集（作者 RGC-2 第 2 显著条目），验证"全局抑制"猜想。
|
||
统计：细胞级 Wilcoxon 仅作探索（pseudoreplication 警示）；同时报告每样本中位数。
|
||
|
||
产出：output/08_replication/
|
||
"""
|
||
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
|
||
from scipy.stats import mannwhitneyu
|
||
|
||
OUT = Path("output/08_replication")
|
||
OUT.mkdir(parents=True, exist_ok=True)
|
||
plt.rcParams.update({"figure.dpi": 300, "savefig.dpi": 300, "font.size": 9})
|
||
|
||
# 作者 leading-edge 基因集（S10 RGC2）
|
||
GENESETS = {
|
||
"ETC_WP295": "Ndufs3/Ndufb6/Cox17/Ndufa1/Uqcrfs1/Slc25a4/Ndufb4/Ndufa3/Ndufa5/Ndufb7/Ndufv2/Ndufs5/Uqcrb/Ndufa6/Uqcrc2/Cox7b/Uqcr11/Slc25a5/Ndufb9/Ndufv1/Ndufa10/Ndufc1/Uqcr10/Sdhc/Ndufc2/Ndufb10/Ndufb5/Cox6b1/Ndufa2/Cox8a/Slc25a14/Cox5a/Surf1/Ndufs7",
|
||
"CI_biogenesis_RE": "Ndufs3/Ndufb6/Ndufa1/Ndufb4/Ndufa3/Ndufa5/Ndufb7/Ndufaf3/Ndufv2/Ndufs5/Ndufa6/Ndufb8/Nubpl/Ndufb9/Ndufv1/Ndufaf5/Ndufa13/Ndufa10/Ndufc1/Ecsit/Ndufc2/Ndufb10/Ndufb5/Ndufa2/Tmem126b/Ndufs7",
|
||
"Glycolysis_WP157": "Gapdh/Pgam2/Tpi1/Aldoc/Ldhb/Eno2/Aldoa/Pkm/Got2/Pdha1/Mdh1/Ldha/Pgk1/Pgam1",
|
||
"Ribosomal_WP163": None, # 运行时从基因名 Rpl/Rps 前缀取
|
||
}
|
||
MITO_AUTOPHAGY = "Pink1/Prkn/Bnip3/Bnip3l/Fundc1/Map1lc3b/Sqstm1/Ulk1/Atg5/Atg7/Becn1/Gabarap".split("/")
|
||
|
||
adata = sc.read_h5ad("data/03_annotated.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() if v}
|
||
gs["Ribosomal_WP163"] = [g for g in raw.var_names if g.startswith(("Rpl", "Rps"))][:200]
|
||
gs["Mitophagy_Autophagy"] = [g for g in MITO_AUTOPHAGY if g in raw.var_names]
|
||
for k, v in gs.items():
|
||
print(f"{k}: {len(v)} genes")
|
||
|
||
for k, v in gs.items():
|
||
sc.tl.score_genes(raw, gene_list=v, score_name=f"score_{k}", use_raw=False)
|
||
adata.obs[f"score_{k}"] = raw.obs[f"score_{k}"]
|
||
|
||
# RGC_highETC 定义：WT RGC 中 ETC 模块分 top 30% 的阈值，套用到全部 RGC
|
||
is_rgc = adata.obs["cell_type"] == "RGC"
|
||
wt_rgc_thr = adata.obs.loc[is_rgc & (adata.obs["genotype"] == "WT"), "score_ETC_WP295"].quantile(0.70)
|
||
adata.obs["rgc_group"] = np.where(
|
||
~is_rgc, "non-RGC",
|
||
np.where(adata.obs["score_ETC_WP295"] >= wt_rgc_thr, "RGC_highETC(~RGC-2)", "RGC_lowETC(~RGC-1)"))
|
||
print("\nRGC 分组核数:")
|
||
print(adata.obs[is_rgc].groupby(["rgc_group", "genotype"], observed=True).size().unstack(fill_value=0))
|
||
adata.write_h5ad("data/04_scored.h5ad") # 保存模块分与 RGC 分组
|
||
|
||
# ---- 比较表：每细胞类型 × 基因型 模块分（中位数 + 探索性 Wilcoxon） ----
|
||
rows = []
|
||
cts = ["RGC_highETC(~RGC-2)", "RGC_lowETC(~RGC-1)", "Muller", "Rod", "Cone", "Bipolar", "Amacrine", "Microglia"]
|
||
for ct in cts:
|
||
if ct.startswith("RGC_"):
|
||
mask = adata.obs["rgc_group"] == ct
|
||
else:
|
||
mask = adata.obs["cell_type"] == ct
|
||
for k in gs:
|
||
col = f"score_{k}"
|
||
wt = adata.obs.loc[mask & (adata.obs["genotype"] == "WT"), col]
|
||
mu = adata.obs.loc[mask & (adata.obs["genotype"] == "V291D"), col]
|
||
if len(wt) < 10 or len(mu) < 10:
|
||
continue
|
||
u, p = mannwhitneyu(mu, wt, alternative="two-sided")
|
||
# 每样本中位数
|
||
per_sample = adata.obs.loc[mask].groupby("sample")[col].median()
|
||
rows.append({
|
||
"group": ct, "geneset": k, "n_WT": len(wt), "n_V291D": len(mu),
|
||
"median_WT": round(wt.median(), 4), "median_V291D": round(mu.median(), 4),
|
||
"delta": round(mu.median() - wt.median(), 4),
|
||
"wilcoxon_p_exploratory": f"{p:.2e}",
|
||
"per_sample_median": "; ".join(f"{s}:{v:.3f}" for s, v in per_sample.items()),
|
||
})
|
||
res = pd.DataFrame(rows)
|
||
res.to_csv(OUT / "module_score_comparison.csv", index=False)
|
||
print("\n", res[res["group"].str.contains("RGC")].to_string(index=False))
|
||
|
||
# ---- 图：RGC 两组 + Muller 的模块分分布（按样本分面） ----
|
||
fig, axes = plt.subplots(1, len(gs), figsize=(3.2 * len(gs), 3.2), sharey=False)
|
||
for ax, k in zip(axes, gs):
|
||
col = f"score_{k}"
|
||
for grp, color in [("RGC_highETC(~RGC-2)", "#d62728"), ("RGC_lowETC(~RGC-1)", "#ff9896"), ("Muller", "#1f77b4")]:
|
||
for s, ls in [("WT", "-"), ("Opa1V291D_S1", "--"), ("Opa1V291D_S2", ":")]:
|
||
if grp.startswith("RGC"):
|
||
m = (adata.obs["rgc_group"] == grp) & (adata.obs["sample"] == s)
|
||
else:
|
||
m = (adata.obs["cell_type"] == grp) & (adata.obs["sample"] == s)
|
||
v = adata.obs.loc[m, col]
|
||
if len(v) > 5:
|
||
ax.hist(v, bins=40, histtype="step", density=True, color=color, ls=ls, lw=1.2,
|
||
label=f"{grp.split('(')[0]} {s.replace('Opa1V291D_', 'V291D_')}")
|
||
ax.set_title(k, fontsize=9)
|
||
ax.set_xlabel("module score")
|
||
axes[0].set_ylabel("density")
|
||
axes[-1].legend(fontsize=6, frameon=False, loc="upper right")
|
||
fig.tight_layout()
|
||
fig.savefig(OUT / "module_score_distributions.png", bbox_inches="tight")
|
||
plt.close(fig)
|
||
|
||
# ---- 基因级抽查：作者 RGC-2 下调 top 基因在我们 RGC_highETC 中的方向 ----
|
||
s2 = pd.read_excel("ref/Data files/S2. snRNAseq-DEG_All_Celltype.xlsx", sheet_name="RGC-2")
|
||
s2.columns = [c.lower() for c in s2.columns]
|
||
s2_sig = s2[(s2["p_val"] < 0.001) & (s2["avg_log2fc"] < -0.5)]
|
||
rawX = raw.X
|
||
vidx = {g: i for i, g in enumerate(raw.var_names)}
|
||
gene_rows = []
|
||
for g in s2_sig["gene"]:
|
||
if g not in vidx:
|
||
continue
|
||
i = vidx[g]
|
||
m_hi = adata.obs["rgc_group"] == "RGC_highETC(~RGC-2)"
|
||
wt_m = m_hi & (adata.obs["genotype"] == "WT")
|
||
mu_m = m_hi & (adata.obs["genotype"] == "V291D")
|
||
if wt_m.sum() < 10 or mu_m.sum() < 10:
|
||
continue
|
||
col_wt = rawX[wt_m.values, i].mean()
|
||
col_mu = rawX[mu_m.values, i].mean()
|
||
gene_rows.append({"gene": g, "author_log2fc": s2.loc[s2["gene"] == g, "avg_log2fc"].iloc[0],
|
||
"our_meanWT": col_wt, "our_meanV291D": col_mu,
|
||
"our_log2fc_approx": np.log2((col_mu + 1e-9) / (col_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_RGChighETC.csv", index=False)
|
||
print(f"\n基因级方向一致率（作者 RGC-2 显著下调基因 vs 我们 RGC_highETC）: {concord*100:.1f}% (n={len(gd)})")
|
||
|
||
print("Done -> output/08_replication/")
|