56 lines
2.0 KiB
Python
56 lines
2.0 KiB
Python
#!/usr/bin/env python
|
||
# -*- coding: utf-8 -*-
|
||
"""
|
||
03_integrate_cluster.py — 归一化、Harmony 整合、聚类（第一轮 P0, Step 2）
|
||
|
||
决策记录 #4：Harmony（按 sample）仅用于聚类/注释；定量比较后续用原始归一化表达。
|
||
产出：data/02_clustered.h5ad；output/03_cluster/*.png
|
||
"""
|
||
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")
|
||
OUT.mkdir(parents=True, exist_ok=True)
|
||
sc.settings.figdir = OUT
|
||
plt.rcParams.update({"figure.dpi": 300, "savefig.dpi": 300, "font.size": 9})
|
||
|
||
adata = sc.read_h5ad("data/01_filtered.h5ad")
|
||
print(f"载入: {adata.n_obs:,} 核 × {adata.n_vars:,} 基因")
|
||
|
||
adata.layers["counts"] = adata.X.copy()
|
||
sc.pp.normalize_total(adata, target_sum=1e4)
|
||
sc.pp.log1p(adata)
|
||
adata.raw = adata # 保留归一化表达供下游定量
|
||
|
||
sc.pp.highly_variable_genes(adata, n_top_genes=3000, flavor="seurat_v3", layer="counts")
|
||
print(f"HVG: {adata.var['highly_variable'].sum()}")
|
||
|
||
sc.pp.scale(adata, max_value=10)
|
||
sc.pp.pca(adata, n_comps=50, mask_var="highly_variable")
|
||
|
||
# Harmony 按 sample 整合（harmonypy 2.x 返回 (n_cells, n_pcs)，scanpy 封装不兼容，直接调用）
|
||
import harmonypy as hm
|
||
_ho = hm.run_harmony(adata.obsm["X_pca"], adata.obs, "sample", verbose=False)
|
||
adata.obsm["X_pca_harmony"] = _ho.Z_corr
|
||
|
||
sc.pp.neighbors(adata, use_rep="X_pca_harmony", n_neighbors=15)
|
||
for res in (0.4, 0.8, 1.2):
|
||
sc.tl.leiden(adata, resolution=res, key_added=f"leiden_{res}", flavor="igraph", n_iterations=2)
|
||
print(f"leiden_{res}: {adata.obs[f'leiden_{res}'].nunique()} clusters")
|
||
sc.tl.umap(adata)
|
||
|
||
adata.write_h5ad("data/02_clustered.h5ad")
|
||
|
||
# UMAP 图
|
||
for key in ["sample", "genotype", "leiden_0.4", "leiden_0.8", "leiden_1.2"]:
|
||
fig = sc.pl.umap(adata, color=key, show=False, return_fig=True, size=3,
|
||
title=f"UMAP — {key}")
|
||
fig.savefig(OUT / f"umap_{key.replace('.', 'p')}.png", bbox_inches="tight")
|
||
plt.close(fig)
|
||
print("Done -> data/02_clustered.h5ad, output/03_cluster/")
|