82 lines
3.3 KiB
Python
82 lines
3.3 KiB
Python
#!/usr/bin/env python
|
||
# -*- coding: utf-8 -*-
|
||
"""
|
||
05_composition.py — 细胞组成分析（第一轮 P0, Step 4）
|
||
排除 LowConf；n=1 vs n=2，仅描述性展示。
|
||
产出：output/05_composition/
|
||
"""
|
||
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/05_composition")
|
||
OUT.mkdir(parents=True, exist_ok=True)
|
||
plt.rcParams.update({"figure.dpi": 300, "savefig.dpi": 300, "font.size": 9})
|
||
|
||
adata = sc.read_h5ad("data/03_annotated.h5ad")
|
||
ad = adata[adata.obs["cell_type"] != "LowConf"].copy()
|
||
|
||
# 每样本 × 细胞类型 计数与比例
|
||
cnt = ad.obs.groupby(["sample", "cell_type"], observed=True).size().unstack(fill_value=0)
|
||
frac = cnt.div(cnt.sum(axis=1), axis=0)
|
||
cnt.to_csv(OUT / "counts_per_sample.csv")
|
||
frac.to_csv(OUT / "fractions_per_sample.csv")
|
||
print("各样本细胞数:", cnt.sum(axis=1).to_dict())
|
||
print("\n比例（%）:")
|
||
print((frac * 100).round(2).to_string())
|
||
|
||
# 基因型层面比例（突变=两样本合并）
|
||
cntg = ad.obs.groupby(["genotype", "cell_type"], observed=True).size().unstack(fill_value=0)
|
||
fracg = cntg.div(cntg.sum(axis=1), axis=0)
|
||
fracg.to_csv(OUT / "fractions_per_genotype.csv")
|
||
|
||
# ---- 图 1：堆叠条形图（按样本） ----
|
||
order = frac.loc[["WT", "Opa1V291D_S1", "Opa1V291D_S2"]].mean().sort_values(ascending=False).index
|
||
fig, ax = plt.subplots(figsize=(4.5, 4))
|
||
bottom = np.zeros(3)
|
||
colors = plt.cm.tab20(np.linspace(0, 1, len(order)))
|
||
for ct, c in zip(order, colors):
|
||
vals = frac.loc[["WT", "Opa1V291D_S1", "Opa1V291D_S2"], ct].values
|
||
ax.bar(["WT", "V291D_S1", "V291D_S2"], vals, bottom=bottom, label=ct, color=c, width=0.65)
|
||
bottom += vals
|
||
ax.set_ylabel("Fraction of nuclei")
|
||
ax.legend(bbox_to_anchor=(1.02, 1), loc="upper left", fontsize=7, frameon=False)
|
||
fig.tight_layout()
|
||
fig.savefig(OUT / "composition_stacked_bar.png", bbox_inches="tight")
|
||
plt.close(fig)
|
||
|
||
# ---- 图 2：WT vs V291D 比例变化（点图，每样本一点） ----
|
||
fig, ax = plt.subplots(figsize=(6.5, 3.5))
|
||
xs = np.arange(len(order))
|
||
w = 0.35
|
||
for i, s in enumerate(["WT", "Opa1V291D_S1", "Opa1V291D_S2"]):
|
||
off = -w if s == "WT" else w / 2 * (1 if s.endswith("S1") else 3)
|
||
ax.scatter(xs + off, frac.loc[s, order] * 100,
|
||
c="black" if s == "WT" else "firebrick",
|
||
marker="o" if s == "WT" else "s", s=28, label=s, zorder=3)
|
||
for j, ct in enumerate(order):
|
||
v = frac.loc["WT", ct] * 100
|
||
m = frac.loc[["Opa1V291D_S1", "Opa1V291D_S2"], ct].mean() * 100
|
||
ax.plot([j - w, j + w], [v, m], color="grey", lw=0.8, zorder=1)
|
||
ax.set_xticks(xs, order, rotation=35, ha="right")
|
||
ax.set_ylabel("Fraction (%)")
|
||
ax.legend(frameon=False, fontsize=8)
|
||
fig.tight_layout()
|
||
fig.savefig(OUT / "composition_change_dotplot.png", bbox_inches="tight")
|
||
plt.close(fig)
|
||
|
||
# RGC 相对比例变化
|
||
rgc_wt = frac.loc["WT", "RGC"] * 100
|
||
rgc_mut = frac.loc[["Opa1V291D_S1", "Opa1V291D_S2"], "RGC"] * 100
|
||
print(f"\nRGC 占比: WT {rgc_wt:.2f}% vs V291D {rgc_mut.values.round(2)}% "
|
||
f"(变化 {(rgc_mut.mean()/rgc_wt-1)*100:+.1f}%)")
|
||
mg_wt = frac.loc["WT"].get("Microglia", 0) * 100
|
||
mg_mut = frac.loc[["Opa1V291D_S1", "Opa1V291D_S2"], "Microglia"] * 100
|
||
print(f"Microglia 占比: WT {mg_wt:.2f}% vs V291D {mg_mut.values.round(2)}%")
|
||
print("Done -> output/05_composition/")
|