156 lines
7.0 KiB
Python
156 lines
7.0 KiB
Python
#!/usr/bin/env python
|
||
# -*- coding: utf-8 -*-
|
||
"""
|
||
11_pseudobulk_de.py — P1-③ pseudobulk DE：方向一致性策略（决策记录 Step 5 核心层）
|
||
|
||
策略（WT n=1 无法做正式推断）：
|
||
1. 每 样本×细胞组 对 counts 求和 → CPM + log2
|
||
2. LFC_S1 = S1−WT，LFC_S2 = S2−WT，LFC_null = S2−S1（同基因型阴性校准）
|
||
3. 一致基因：sign(LFC_S1)==sign(LFC_S2) 且 |LFC|>=0.5 且 |LFC_null| < 两次生物 LFC 的较小者
|
||
4. 表达过滤：至少一个样本 CPM>=10
|
||
5. 与作者 S2 各细胞类型表做 Spearman 相关 + 方向一致率
|
||
|
||
细胞组：所有 cell_type（WT>=20 核）+ RGC 按签名分 RGC1-like/RGC2-like（RGC2-like WT 仅 7 核，标记低功效）。
|
||
产出：output/11_pseudobulk_de/
|
||
"""
|
||
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
|
||
from scipy.stats import spearmanr
|
||
|
||
OUT = Path("output/11_pseudobulk_de")
|
||
OUT.mkdir(parents=True, exist_ok=True)
|
||
plt.rcParams.update({"figure.dpi": 300, "savefig.dpi": 300, "font.size": 8})
|
||
MIN_NUCLEI = 20
|
||
LFC_THR = 0.5
|
||
CPM_THR = 10
|
||
|
||
adata = sc.read_h5ad("data/05_rgc2sig.h5ad")
|
||
counts = sp.csr_matrix(adata.layers["counts"])
|
||
|
||
# ---- 细胞组定义：cell_type + RGC 签名亚组 ----
|
||
obs = adata.obs.copy()
|
||
obs["pb_group"] = obs["cell_type"].astype(str)
|
||
is_rgc = obs["cell_type"] == "RGC"
|
||
obs.loc[is_rgc, "pb_group"] = obs.loc[is_rgc, "rgc2_sig_group"].astype(str) # RGC1-like / RGC2-like
|
||
obs.loc[obs["pb_group"] == "LowConf", "pb_group"] = "LowConf"
|
||
|
||
sizes = obs.groupby(["pb_group", "sample"]).size().unstack(fill_value=0)
|
||
sizes.to_csv(OUT / "pseudobulk_group_sizes.csv")
|
||
groups = [g for g in sizes.index if (sizes.loc[g] >= MIN_NUCLEI).all()]
|
||
lowpower = [g for g in sizes.index if g not in groups and (sizes.loc[g] >= 5).all() and g not in ("LowConf",)]
|
||
print("正式组:", groups)
|
||
print("低功效组（描述性）:", lowpower)
|
||
|
||
samples = ["WT", "Opa1V291D_S1", "Opa1V291D_S2"]
|
||
|
||
def pseudobulk(group):
|
||
"""返回 3×genes 的 raw count 矩阵。"""
|
||
mat = np.zeros((3, adata.n_vars))
|
||
for i, s in enumerate(samples):
|
||
idx = np.where((obs["pb_group"] == group) & (obs["sample"] == s))[0]
|
||
mat[i] = np.asarray(counts[idx].sum(axis=0)).ravel()
|
||
return mat
|
||
|
||
def median_of_ratios(mat):
|
||
"""DESeq2 式 size factor：基因几何均值跨样本，样本 size factor = count/geomean 的中位数。
|
||
校正文库组成与全局检出率差异（CPM 在深度悬殊时产生全局 LFC 偏倚）。"""
|
||
pos = (mat > 0).all(axis=0)
|
||
geomean = np.exp(np.log(mat[:, pos]).mean(axis=0))
|
||
ratios = mat[:, pos] / geomean
|
||
sf = np.median(ratios, axis=1)
|
||
return mat / sf[:, None]
|
||
|
||
def compute_de(group):
|
||
mat = pseudobulk(group)
|
||
norm = median_of_ratios(mat) # 归一化后的"每 size-factor"计数
|
||
keep = (mat >= 1).any(axis=0) & ((mat / mat.sum(axis=1, keepdims=True) * 1e6) >= CPM_THR).any(axis=0)
|
||
l2 = np.log2(norm + 0.5)
|
||
df = pd.DataFrame({
|
||
"gene": adata.var_names,
|
||
"cpm_WT": mat[0] / mat[0].sum() * 1e6, "cpm_S1": mat[1] / mat[1].sum() * 1e6,
|
||
"cpm_S2": mat[2] / mat[2].sum() * 1e6,
|
||
"lfc_S1": l2[1] - l2[0], "lfc_S2": l2[2] - l2[0], "lfc_null": l2[2] - l2[1],
|
||
})[keep].copy()
|
||
# 零检出守卫：UP 要求 WT 有基线检出（CPM>=1），DOWN 要求突变有检出，
|
||
# 防止深度不对称造成的"零 vs 有"伪差异（突变文库深，假 UP 偏多）
|
||
detect_ok = np.where(
|
||
np.sign(df["lfc_S1"]) > 0, df["cpm_WT"] >= 1.0,
|
||
np.minimum(df["cpm_S1"], df["cpm_S2"]) >= 1.0)
|
||
consistent = (
|
||
(np.sign(df["lfc_S1"]) == np.sign(df["lfc_S2"]))
|
||
& (df["lfc_S1"].abs() >= LFC_THR) & (df["lfc_S2"].abs() >= LFC_THR)
|
||
& (df["lfc_null"].abs() < np.minimum(df["lfc_S1"].abs(), df["lfc_S2"].abs()))
|
||
& detect_ok
|
||
)
|
||
df["consistent"] = consistent
|
||
df["direction"] = np.where(consistent, np.sign(df["lfc_S1"]).map({1.0: "UP", -1.0: "DOWN"}), "ns")
|
||
return df
|
||
|
||
# ---- 作者 S2 对照 ----
|
||
author_map = {"Rod": "Rod", "Cone": "Cone", "RGC1-like": "RGC-1", "RGC2-like": "RGC-2",
|
||
"RGC": "RGC-2", "Amacrine": "Amacrine", "Muller": "Muller",
|
||
"Uveal_Melanocyte": "Uveal", "Bipolar": "Bipolar",
|
||
"Horizontal": "Horizontal", "Pericyte": "Pericyte"}
|
||
xl = "ref/Data files/S2. snRNAseq-DEG_All_Celltype.xlsx"
|
||
|
||
all_groups = groups + [g for g in lowpower if g in author_map or g in ("RGC2-like",)]
|
||
summary_rows = []
|
||
de_tables = {}
|
||
for g in all_groups:
|
||
df = compute_de(g)
|
||
df["n_WT_nuclei"] = int(sizes.loc[g, "WT"])
|
||
df["low_power"] = g not in groups
|
||
df.to_csv(OUT / f"de_{g.replace('/', '_')}.csv", index=False)
|
||
de_tables[g] = df
|
||
n_up = (df["direction"] == "UP").sum()
|
||
n_dn = (df["direction"] == "DOWN").sum()
|
||
row = {"group": g, "n_genes_tested": len(df), "consistent_UP": int(n_up),
|
||
"consistent_DOWN": int(n_dn), "UP_DOWN_ratio": round(n_up / max(n_dn, 1), 3),
|
||
"n_WT_nuclei": int(sizes.loc[g, "WT"]), "low_power": g not in groups}
|
||
if g in author_map:
|
||
a = pd.read_excel(xl, sheet_name=author_map[g])
|
||
a.columns = [c.lower() for c in a.columns]
|
||
mg = df.merge(a[["gene", "avg_log2fc"]], on="gene", how="inner")
|
||
mg["our_lfc_mean"] = (mg["lfc_S1"] + mg["lfc_S2"]) / 2
|
||
rho, p = spearmanr(mg["our_lfc_mean"], mg["avg_log2fc"])
|
||
sig_a = a[a["p_val"] < 0.001].merge(df, on="gene", how="inner")
|
||
if len(sig_a):
|
||
sig_a = sig_a.assign(our_mean=(sig_a["lfc_S1"] + sig_a["lfc_S2"]) / 2)
|
||
conc = (np.sign(sig_a["avg_log2fc"]) == np.sign(sig_a["our_mean"])).mean()
|
||
else:
|
||
conc = np.nan
|
||
row.update({"author_sheet": author_map[g], "spearman_rho": round(rho, 3),
|
||
"author_sig_gene_concordance": round(conc, 3) if conc == conc else None})
|
||
de_tables[g + ".__author_merge"] = mg
|
||
summary_rows.append(row)
|
||
|
||
summary = pd.DataFrame(summary_rows)
|
||
summary.to_csv(OUT / "de_summary.csv", index=False)
|
||
print("\n", summary.to_string(index=False))
|
||
|
||
# ---- 图：与作者 log2FC 散点对照（Muller / RGC2-like / Rod / RGC1-like） ----
|
||
show = [g for g in ["Muller", "RGC2-like", "RGC1-like", "Rod"] if g + ".__author_merge" in de_tables]
|
||
fig, axes = plt.subplots(1, len(show), figsize=(3.4 * len(show), 3.4))
|
||
if len(show) == 1:
|
||
axes = [axes]
|
||
for ax, g in zip(axes, show):
|
||
mg = de_tables[g + ".__author_merge"]
|
||
ax.scatter(mg["avg_log2fc"], mg["our_lfc_mean"], s=1, alpha=0.15, color="grey", rasterized=True)
|
||
cons = mg[mg["consistent"]]
|
||
ax.scatter(cons["avg_log2fc"], (cons["lfc_S1"] + cons["lfc_S2"]) / 2, s=2, alpha=0.5,
|
||
color="firebrick", rasterized=True)
|
||
rho = summary.loc[summary["group"] == g, "spearman_rho"].iloc[0]
|
||
ax.axhline(0, color="k", lw=0.4); ax.axvline(0, color="k", lw=0.4)
|
||
ax.set_xlabel("authors avg_log2FC"); ax.set_ylabel("our pseudobulk LFC (mean)")
|
||
ax.set_title(f"{g} (Spearman {rho})")
|
||
fig.tight_layout()
|
||
fig.savefig(OUT / "author_comparison_scatter.png", bbox_inches="tight")
|
||
print("Done ->", OUT)
|