170 lines
8.0 KiB
Python
170 lines
8.0 KiB
Python
#!/usr/bin/env python
|
||
# -*- coding: utf-8 -*-
|
||
"""
|
||
09_muller_downsample.py — P1-① Müller 全面上调的"真激活 vs 深度假象"裁决
|
||
|
||
背景：突变 Müller 的 ETC/CI/核糖体/糖酵解模块分全面上调（P0 意外发现）。
|
||
但突变文库核数多（1070/1319 vs WT 183）且测序深（中位 UMI 高 ~35%），
|
||
模块分可能受检出率膨胀驱动。本脚本做匹配敏感性对照（决策记录 第二轮裁决 #1）：
|
||
|
||
1. 突变 Müller 随机抽至 183 核（=WT 核数）
|
||
2. 每核按 WT Müller total_counts 经验分布逐个下采样 UMI（深度也匹配）
|
||
3. 重算 5 个代谢模块分，比较 WT vs 下采样突变
|
||
4. 5 个随机种子重复，看方向是否稳定
|
||
|
||
阳性对照：RGC_highETC 的糖酵解/线粒体自噬下调（P0 已复现的真信号）做同样处理，
|
||
验证"下采样后仍能检出真实差异"的功效。WT RGC_highETC 仅 21 核，结果仅供参考。
|
||
|
||
产出：output/09_downsample/
|
||
"""
|
||
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/09_downsample")
|
||
OUT.mkdir(parents=True, exist_ok=True)
|
||
plt.rcParams.update({"figure.dpi": 300, "savefig.dpi": 300, "font.size": 8})
|
||
N_SEEDS = 5
|
||
|
||
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",
|
||
"Mitophagy_Autophagy": "Pink1/Prkn/Bnip3/Bnip3l/Fundc1/Map1lc3b/Sqstm1/Ulk1/Atg5/Atg7/Becn1/Gabarap",
|
||
}
|
||
|
||
adata = sc.read_h5ad("data/04_scored.h5ad") # 含 rgc_group 与 counts 层
|
||
counts = adata.layers["counts"]
|
||
if not hasattr(counts, "tocsr"):
|
||
counts = counts # 已是稀疏
|
||
import scipy.sparse as sp
|
||
counts = sp.csr_matrix(counts)
|
||
|
||
|
||
def downsample_group(obs_mask, wt_ref_counts_per_cell, rng, label):
|
||
"""对 mask 选中的核：先抽至 len(wt_ref) 个核，再逐核下采样到 WT 经验深度。返回下采样后的 counts 子矩阵与细胞名。"""
|
||
idx = np.where(obs_mask.values)[0]
|
||
n_target = len(wt_ref_counts_per_cell)
|
||
if len(idx) <= n_target:
|
||
chosen = idx
|
||
else:
|
||
chosen = rng.choice(idx, size=n_target, replace=False)
|
||
sub = counts[chosen].copy()
|
||
targets = rng.choice(wt_ref_counts_per_cell, size=len(chosen), replace=True)
|
||
rows, cols, data = [], [], []
|
||
for i in range(sub.shape[0]):
|
||
row = sub.getrow(i)
|
||
tot = row.sum()
|
||
tgt = targets[i]
|
||
if tot <= tgt or tot == 0:
|
||
continue
|
||
# 逐 UMI 不放回抽样
|
||
gene_idx = row.indices
|
||
gene_cnt = row.data.astype(int)
|
||
expanded = np.repeat(gene_idx, gene_cnt)
|
||
keep = rng.choice(expanded, size=int(tgt), replace=False)
|
||
uc, cc = np.unique(keep, return_counts=True)
|
||
rows.append(np.full(len(uc), i)); cols.append(uc); data.append(cc)
|
||
if rows:
|
||
r = np.concatenate(rows); c = np.concatenate(cols); d = np.concatenate(data)
|
||
sub = sp.csr_matrix((d, (r, c)), shape=sub.shape)
|
||
names = adata.obs_names[chosen]
|
||
return names, sub
|
||
|
||
|
||
def score_subset(wt_names, wt_mat, mu_names, mu_mat, gs):
|
||
"""WT + 下采样突变拼一个对象，normalize+log1p+score，返回两组中位数差。"""
|
||
import anndata as ad
|
||
X = sp.vstack([wt_mat, mu_mat])
|
||
obs = pd.DataFrame(
|
||
{"group": ["WT"] * len(wt_names) + ["MU"] * len(mu_names)},
|
||
index=list(wt_names) + list(mu_names))
|
||
a = ad.AnnData(X=X, obs=obs, var=pd.DataFrame(index=adata.var_names))
|
||
sc.pp.normalize_total(a, target_sum=1e4)
|
||
sc.pp.log1p(a)
|
||
out = {}
|
||
for k, v in gs.items():
|
||
sc.tl.score_genes(a, gene_list=v, score_name="s", use_raw=False, random_state=1)
|
||
med = a.obs.groupby("group")["s"].median()
|
||
out[k] = med["MU"] - med["WT"]
|
||
return out
|
||
|
||
|
||
# 基因集与数据取交集
|
||
gs = {k: [g for g in v.split("/") if g in adata.var_names] for k, v in GENESETS.items()}
|
||
gs["Ribosomal_WP163"] = [g for g in adata.var_names if g.startswith(("Rpl", "Rps"))][:200]
|
||
print({k: len(v) for k, v in gs.items()})
|
||
|
||
# ---- 任务组：Muller（裁决对象） + RGC_highETC（阳性对照） ----
|
||
rgc_hi = adata.obs["rgc_group"] == "RGC_highETC(~RGC-2)"
|
||
tasks = {
|
||
"Muller": adata.obs["cell_type"] == "Muller",
|
||
"RGC_highETC": rgc_hi,
|
||
}
|
||
results = []
|
||
for task_name, mask in tasks.items():
|
||
wt_mask = mask & (adata.obs["genotype"] == "WT")
|
||
wt_names = adata.obs_names[wt_mask]
|
||
wt_mat = counts[np.where(wt_mask.values)[0]]
|
||
wt_depths = np.asarray(wt_mat.sum(axis=1)).ravel()
|
||
print(f"\n== {task_name}: WT n={wt_mask.sum()}, WT 深度中位 {np.median(wt_depths):.0f}")
|
||
for sample in ["Opa1V291D_S1", "Opa1V291D_S2"]:
|
||
mu_mask = mask & (adata.obs["sample"] == sample)
|
||
for seed in range(N_SEEDS):
|
||
rng = np.random.default_rng(1000 + seed)
|
||
mu_names, mu_mat = downsample_group(mu_mask, wt_depths, rng, task_name)
|
||
deltas = score_subset(wt_names, wt_mat, mu_names, mu_mat, gs)
|
||
for k, d in deltas.items():
|
||
results.append({"group": task_name, "sample": sample, "seed": seed,
|
||
"geneset": k, "delta_median_score": round(d, 4)})
|
||
print(f" {sample}: 完成 {N_SEEDS} 种子")
|
||
|
||
res = pd.DataFrame(results)
|
||
res.to_csv(OUT / "downsample_replicates.csv", index=False)
|
||
summary = res.groupby(["group", "sample", "geneset"])["delta_median_score"].agg(["mean", "std", "count"]).round(4)
|
||
summary.to_csv(OUT / "downsample_summary.csv")
|
||
print("\n", summary.to_string())
|
||
|
||
# ---- 对照：原始全量数据的 delta（从 04_scored.h5ad 读模块分中位数） ----
|
||
full = sc.read_h5ad("data/04_scored.h5ad")
|
||
full_rows = []
|
||
for task_name, mask in tasks.items():
|
||
for k in gs:
|
||
col = f"score_{k}"
|
||
wt_med = full.obs.loc[mask & (full.obs["genotype"] == "WT"), col].median()
|
||
for s in ["Opa1V291D_S1", "Opa1V291D_S2"]:
|
||
mu_med = full.obs.loc[mask & (full.obs["sample"] == s), col].median()
|
||
full_rows.append({"group": task_name, "sample": s, "geneset": k,
|
||
"delta_full": round(mu_med - wt_med, 4)})
|
||
fulldf = pd.DataFrame(full_rows)
|
||
fulldf.to_csv(OUT / "fulldata_deltas.csv", index=False)
|
||
|
||
# ---- 图：全量 vs 下采样 delta 对照 ----
|
||
gs_keys = list(gs.keys())
|
||
fig, axes = plt.subplots(1, 2, figsize=(11, 4), sharey=False)
|
||
for ax, task_name in zip(axes, ["Muller", "RGC_highETC"]):
|
||
x = np.arange(len(gs_keys))
|
||
w = 0.35
|
||
for off, s, c in [(-w/2, "Opa1V291D_S1", "firebrick"), (w/2, "Opa1V291D_S2", "darkorange")]:
|
||
full_d = fulldf[(fulldf["group"] == task_name) & (fulldf["sample"] == s)].set_index("geneset")["delta_full"].reindex(gs_keys)
|
||
ax.bar(x + off, full_d, width=w, color=c, alpha=0.35, label=f"{s.split('_')[-1]} full", hatch="//")
|
||
sub = res[(res["group"] == task_name) & (res["sample"] == s)]
|
||
m = sub.groupby("geneset")["delta_median_score"].mean().reindex(gs_keys)
|
||
e = sub.groupby("geneset")["delta_median_score"].std().reindex(gs_keys)
|
||
ax.bar(x + off, m, width=w * 0.9, yerr=e, color=c, alpha=0.9, capsize=2,
|
||
label=f"{s.split('_')[-1]} downsampled (5 seeds)")
|
||
ax.set_xticks(x, [g.replace("_WP295", "").replace("_RE", "").replace("_WP157", "").replace("_WP163", "").replace("_", "\n") for g in gs_keys], fontsize=7)
|
||
ax.axhline(0, color="grey", lw=0.6)
|
||
ax.set_title(task_name)
|
||
ax.legend(fontsize=6, frameon=False)
|
||
axes[0].set_ylabel("Δ median module score (MU − WT)")
|
||
fig.suptitle("Does Muller upregulation survive downsampling (nuclei + depth matched to WT)?", y=1.02)
|
||
fig.tight_layout()
|
||
fig.savefig(OUT / "downsample_adjudication.png", bbox_inches="tight")
|
||
print("Done ->", OUT)
|