148 lines
6.6 KiB
Python
148 lines
6.6 KiB
Python
#!/usr/bin/env python
|
||
# -*- coding: utf-8 -*-
|
||
"""
|
||
04_annotate.py — 细胞类型注释 + microglia 捕捞 + RGC 亚群（第一轮 P0, Step 3）
|
||
|
||
- 基准 marker：原文表 S3（9 类）+ 补充 microglia/astrocyte/血管/少突 marker
|
||
- 按 leiden_0.8 cluster 打分注释（z-scored mean expression）
|
||
- microglia：多 marker 共表达门控（Aif1/C1qa/Tmem119/P2ry12/Hexb）
|
||
- RGC 亚群重聚类，复现 RGC-1/RGC-2
|
||
|
||
产出：data/03_annotated.h5ad；output/04_annotation/
|
||
"""
|
||
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/04_annotation")
|
||
OUT.mkdir(parents=True, exist_ok=True)
|
||
plt.rcParams.update({"figure.dpi": 300, "savefig.dpi": 300, "font.size": 9})
|
||
|
||
MARKERS = {
|
||
"Rod": ["Rho", "Gnat1", "Sag", "Cnga1"],
|
||
"Cone": ["Arr3", "Gnat2", "Opn1mw"],
|
||
"Bipolar": ["Grm6", "Prkca", "Cabp5", "Vsx2"],
|
||
"Horizontal": ["Lhx1", "Onecut2", "Prox1"],
|
||
"Amacrine": ["Gad1", "Gad2", "Tfap2b"],
|
||
"RGC": ["Rbpms", "Slc17a6", "Thy1"],
|
||
"Muller": ["Rlbp1", "Slc1a3", "Glul"],
|
||
"Microglia": ["Aif1", "C1qa", "C1qb", "Tmem119", "Hexb"],
|
||
"Astrocyte": ["Gfap", "S100b", "Aqp4"],
|
||
"Pericyte": ["Pdgfrb", "Rgs5"],
|
||
"Endothelial": ["Pecam1", "Cldn5", "Kdr"],
|
||
"Uveal_Melanocyte": ["Gpnmb", "Tyr", "Mlana"],
|
||
"Oligodendrocyte": ["Mog", "Plp1", "Mbp"],
|
||
}
|
||
|
||
adata = sc.read_h5ad("data/02_clustered.h5ad")
|
||
print(f"载入: {adata.n_obs:,} 核")
|
||
CLU = "leiden_0.8"
|
||
|
||
# ---------- 1. cluster 级 marker 打分 ----------
|
||
genes_all = [g for gs in MARKERS.values() for g in gs if g in adata.raw.var_names]
|
||
# 用 raw（log-normalized）按 cluster 求均值（one-hot × X 手动聚合，稳健）
|
||
import scipy.sparse as sp
|
||
raw_ad = adata.raw.to_adata()
|
||
cats = pd.Categorical(adata.obs[CLU].astype(str))
|
||
onehot = sp.csr_matrix(
|
||
(np.ones(len(cats)), (cats.codes, np.arange(len(cats)))),
|
||
shape=(len(cats.categories), len(cats)),
|
||
)
|
||
sums = onehot @ raw_ad.X
|
||
if sp.issparse(sums):
|
||
sums = sums.toarray()
|
||
counts = np.bincount(cats.codes, minlength=len(cats.categories))
|
||
X_agg = sums / counts[:, None]
|
||
df = pd.DataFrame(X_agg, index=cats.categories.astype(str), columns=raw_ad.var_names)
|
||
df = df[[g for g in genes_all if g in df.columns]]
|
||
z = (df - df.mean()) / df.std(ddof=0)
|
||
z = z.fillna(0.0)
|
||
|
||
ct_score = {}
|
||
for ct, gs in MARKERS.items():
|
||
gs_ok = [g for g in gs if g in z.columns]
|
||
ct_score[ct] = z[gs_ok].mean(axis=1)
|
||
score_df = pd.DataFrame(ct_score).astype(float).fillna(0.0)
|
||
score_df.to_csv(OUT / "cluster_marker_scores.csv")
|
||
|
||
assign = score_df.idxmax(axis=1)
|
||
margin = score_df.apply(lambda r: r.nlargest(2).iloc[0] - r.nlargest(2).iloc[1], axis=1)
|
||
assign_tbl = pd.DataFrame({
|
||
"cluster": score_df.index, "assigned": assign, "margin": margin,
|
||
"n_cells": adata.obs[CLU].value_counts(),
|
||
"top3": score_df.apply(lambda r: ", ".join(f"{k}:{v:.2f}" for k, v in r.nlargest(3).items()), axis=1),
|
||
})
|
||
assign_tbl.to_csv(OUT / "cluster_assignment_auto.csv", index=False)
|
||
print(assign_tbl.sort_values("assigned").to_string(index=False))
|
||
|
||
# ---------- 2. marker dotplot（证据图） ----------
|
||
order = assign_tbl.sort_values(["assigned", "cluster"])["cluster"].tolist()
|
||
fig = sc.pl.dotplot(adata, var_names={k: [g for g in v if g in adata.raw.var_names] for k, v in MARKERS.items()},
|
||
groupby=CLU, categories_order=order, use_raw=True,
|
||
show=False, return_fig=True)
|
||
fig.savefig(OUT / "marker_dotplot.png", bbox_inches="tight")
|
||
plt.close("all")
|
||
|
||
# ---------- 3. 应用注释（margin < 0.5 判为 LowConf，低信度 cluster 不进定量比较） ----------
|
||
mapping = assign.to_dict()
|
||
adata.obs["cell_type"] = adata.obs[CLU].map(mapping).astype(str)
|
||
lowconf = margin[margin < 0.5].index.astype(str)
|
||
adata.obs.loc[adata.obs[CLU].astype(str).isin(lowconf), "cell_type"] = "LowConf"
|
||
adata.obs["cell_type"] = pd.Categorical(adata.obs["cell_type"])
|
||
print(f"\n低信度 cluster（margin<0.5）: {sorted(lowconf)} -> 标记 LowConf")
|
||
print("\n细胞类型分布:")
|
||
print(adata.obs.groupby(["cell_type", "genotype"], observed=True).size().unstack(fill_value=0))
|
||
|
||
# ---------- 4. microglia 捕捞核查 ----------
|
||
adata_raw = adata.raw.to_adata()
|
||
adata_raw.obs = adata.obs
|
||
mg_markers = [g for g in ["Aif1", "C1qa", "C1qb", "Tmem119", "Hexb", "P2ry12", "Cx3cr1"] if g in adata_raw.var_names]
|
||
mg_expr = sc.get.obs_df(adata_raw, keys=mg_markers)
|
||
mg_hits = (mg_expr > 0).sum(axis=1)
|
||
adata.obs["mg_marker_hits"] = mg_hits
|
||
cand = adata.obs[mg_hits >= 3]
|
||
print(f"\nmicroglia 候选（≥3 marker 共表达）: {len(cand)} 核")
|
||
print(cand.groupby(["cell_type", "genotype"], observed=True).size().unstack(fill_value=0))
|
||
if len(cand):
|
||
adata.obs["is_microglia_candidate"] = adata.obs.index.isin(cand.index)
|
||
|
||
# ---------- 5. RGC 亚群重聚类 ----------
|
||
rgc = adata[adata.obs["cell_type"] == "RGC"].copy()
|
||
print(f"\nRGC 亚群重聚类: {rgc.n_obs:,} 核")
|
||
if rgc.n_obs > 200:
|
||
rgc.X = rgc.layers["counts"].copy()
|
||
sc.pp.normalize_total(rgc, target_sum=1e4)
|
||
sc.pp.log1p(rgc)
|
||
sc.pp.highly_variable_genes(rgc, n_top_genes=2000, flavor="seurat_v3", layer="counts")
|
||
sc.pp.scale(rgc, max_value=10)
|
||
sc.pp.pca(rgc, n_comps=30, mask_var="highly_variable")
|
||
import harmonypy as hm
|
||
_ho = hm.run_harmony(rgc.obsm["X_pca"], rgc.obs, "sample", verbose=False)
|
||
rgc.obsm["X_pca_harmony"] = _ho.Z_corr
|
||
sc.pp.neighbors(rgc, use_rep="X_pca_harmony")
|
||
sc.tl.leiden(rgc, resolution=0.2, key_added="rgc_sub", flavor="igraph", n_iterations=2)
|
||
sc.tl.umap(rgc)
|
||
print(rgc.obs.groupby(["rgc_sub", "genotype"], observed=True).size().unstack(fill_value=0))
|
||
# OXPHOS/ETC 模块分（WT 中能量需求最高者对应原文 RGC-2）
|
||
etc_genes = [g for g in rgc.var_names if g.startswith(("mt-Nd", "mt-Co", "mt-Atp", "mt-Cytb", "Nduf", "Cox", "Atp5", "Uqcr", "Sdh"))]
|
||
sc.tl.score_genes(rgc, gene_list=etc_genes, score_name="ETC_score", use_raw=False)
|
||
wt_score = rgc.obs[rgc.obs["genotype"] == "WT"].groupby("rgc_sub", observed=True)["ETC_score"].mean().sort_values(ascending=False)
|
||
print("\nWT 中各 RGC 亚群 ETC 模块分（高者 ≈ 原文 RGC-2）:")
|
||
print(wt_score)
|
||
fig = sc.pl.umap(rgc, color=["rgc_sub", "genotype", "ETC_score"], show=False, return_fig=True, size=8)
|
||
fig.savefig(OUT / "rgc_subcluster_umap.png", bbox_inches="tight")
|
||
plt.close("all")
|
||
# 写回主对象
|
||
adata.obs["rgc_sub"] = np.nan
|
||
adata.obs.loc[rgc.obs_names, "rgc_sub"] = rgc.obs["rgc_sub"].astype(str)
|
||
adata.obs.loc[rgc.obs_names, "ETC_score_rgc"] = rgc.obs["ETC_score"]
|
||
rgc.write_h5ad("data/03b_rgc_subset.h5ad")
|
||
|
||
adata.write_h5ad("data/03_annotated.h5ad")
|
||
print("\nDone -> data/03_annotated.h5ad, output/04_annotation/")
|