82 lines
3.7 KiB
Python
82 lines
3.7 KiB
Python
#!/usr/bin/env python
|
||
# -*- coding: utf-8 -*-
|
||
"""
|
||
15_umap_tsne_2x2.py — 2×2 降维可视化对照：UMAP / tSNE × 细胞类型 / 基因型
|
||
|
||
tSNE 在 Harmony 校正的 PCA 空间（X_pca_harmony，n_pcs=50）计算，与 UMAP 同输入，
|
||
仅作可视化对照（两者都不参与聚类与任何定量，见报告说明）。
|
||
着色：细胞类型（14 类，小群置顶以便可见）与基因型（WT vs V291D）。
|
||
产出：output/03_cluster/umap_tsne_2x2.png；tSNE 坐标写回 data/04_scored.h5ad (obsm["X_tsne"])。
|
||
"""
|
||
from pathlib import Path
|
||
|
||
import matplotlib
|
||
matplotlib.use("Agg")
|
||
import matplotlib.pyplot as plt
|
||
import numpy as np
|
||
import scanpy as sc
|
||
|
||
OUT = Path("output/03_cluster")
|
||
plt.rcParams.update({"figure.dpi": 300, "savefig.dpi": 300, "font.size": 9})
|
||
|
||
adata = sc.read_h5ad("data/04_scored.h5ad")
|
||
|
||
# ---- 计算 tSNE（若尚未计算） ----
|
||
if "X_tsne" not in adata.obsm:
|
||
sc.tl.tsne(adata, use_rep="X_pca_harmony", n_pcs=50, random_state=0, perplexity=30)
|
||
adata.write_h5ad("data/04_scored.h5ad")
|
||
print("X_tsne 就绪:", adata.obsm["X_tsne"].shape)
|
||
|
||
# ---- 细胞类型调色板（14 类，高区分度；LowConf 用灰色表"低置信度"） ----
|
||
cts = adata.obs["cell_type"].astype(str)
|
||
ct_order = cts.value_counts().index.tolist() # 大类在前，先画底层
|
||
CT_COLORS = {
|
||
"Rod": "#4e79a7", "LowConf": "#a8a8a8", "Muller": "#f28e2b", "Amacrine": "#e15759",
|
||
"Bipolar": "#76b7b2", "Cone": "#59a14f", "RGC": "#edc948", "Pericyte": "#b07aa1",
|
||
"Uveal_Melanocyte": "#9c755f", "Oligodendrocyte": "#ff9da7", "Endothelial": "#8cd17d",
|
||
"Horizontal": "#86bcb6", "Microglia": "#d37295", "Astrocyte": "#f7b6d2",
|
||
}
|
||
cmap_ct = {ct: CT_COLORS[ct] for ct in ct_order}
|
||
# 小群手工放大，便于在散点中辨认
|
||
SIZE = {"Microglia": 16, "Astrocyte": 16, "RGC": 11, "Horizontal": 11, "Endothelial": 9}
|
||
ALPHA = {"Rod": 0.55, "Amacrine": 0.6, "Muller": 0.6}
|
||
|
||
fig, axes = plt.subplots(2, 2, figsize=(11, 11))
|
||
|
||
def draw(ax, emb, color_by):
|
||
if color_by == "cell_type":
|
||
# 大类先画底层，小群后画顶层
|
||
for ct in reversed(ct_order): # 反过来：最小的最后画 -> 顶层
|
||
m = (cts == ct).values
|
||
ax.scatter(emb[m, 0], emb[m, 1], s=SIZE.get(ct, 5), c=[cmap_ct[ct]],
|
||
alpha=ALPHA.get(ct, 0.75), linewidths=0, rasterized=True, label=ct)
|
||
ax.legend(markerscale=2.2, fontsize=6.5, frameon=False, loc="center left",
|
||
bbox_to_anchor=(1.0, 0.5), ncol=1)
|
||
else:
|
||
# 少数样本（WT）画顶层并放大加白边，避免被多数组（V291D）淹没
|
||
mu = (adata.obs["genotype"] == "V291D").values
|
||
wt = (adata.obs["genotype"] == "WT").values
|
||
ax.scatter(emb[mu, 0], emb[mu, 1], s=4, c="#d62728", alpha=0.30,
|
||
linewidths=0, rasterized=True, label=f"V291D (n={mu.sum()})")
|
||
ax.scatter(emb[wt, 0], emb[wt, 1], s=15, c="#1f77b4", alpha=0.95,
|
||
edgecolors="white", linewidths=0.3, rasterized=True, label=f"WT (n={wt.sum()})")
|
||
ax.legend(markerscale=2.2, fontsize=7, frameon=False, loc="best")
|
||
ax.set_xticks([]); ax.set_yticks([])
|
||
for sp in ax.spines.values():
|
||
sp.set_visible(False)
|
||
|
||
draw(axes[0, 0], adata.obsm["X_umap"], "cell_type")
|
||
draw(axes[0, 1], adata.obsm["X_umap"], "genotype")
|
||
draw(axes[1, 0], adata.obsm["X_tsne"], "cell_type")
|
||
draw(axes[1, 1], adata.obsm["X_tsne"], "genotype")
|
||
|
||
axes[0, 0].set_title("UMAP — cell type")
|
||
axes[0, 1].set_title("UMAP — genotype")
|
||
axes[1, 0].set_title("tSNE (Harmony) — cell type")
|
||
axes[1, 1].set_title("tSNE (Harmony) — genotype")
|
||
|
||
fig.suptitle("UMAP vs tSNE embedding comparison (both on Harmony-corrected PCA)", y=0.98)
|
||
fig.tight_layout(rect=[0, 0, 1, 0.96])
|
||
fig.savefig(OUT / "umap_tsne_2x2.png", bbox_inches="tight")
|
||
print("Done ->", OUT / "umap_tsne_2x2.png")
|