NG文章的空间转录组没给h5怎么办:用Python手动构建AnnData,让空间坐标、图片、尺度信息都对上

import scanpy as scimport matplotlib.pyplot as pltimport seaborn as snsimport gzip, ioimport numpy as np, pandas as pdimport scipy.io, scipy.sparse as spimport h5pyimport os, jsonfrom PIL import Image

MTX_DIR = r"E:\xy\2026\2.19\GSE267904_RAW\\matrix\\"                  # 含 barcodes.tsv.gz / features.tsv.gz / matrix.mtx.gzSPATIAL_DIR = r"E:\xy\2026\2.19\GSE267904_RAW\\spatial"             # 含 tissue_positions_list.csv / scalefactors_json.json / tissue_*_image.pngLIBRARY_ID = "072-A1"

adata = sc.read_10x_mtx(    MTX_DIR,    var_names="gene_symbols",   # 你也可以换成 gene_ids    cache=True)adata.var_names_make_unique() 

sf_path = os.path.join(SPATIAL_DIR,                        "scalefactors_json.json")withopen(sf_path, "r", encoding="utf-8") as f:    scalefactors = json.load(f)

pos_path = os.path.join(SPATIAL_DIR, "tissue_positions_list.csv")pos = pd.read_csv(pos_path, header=None)# 10x Visium 标准（无表头）6列：# barcode, in_tissue, array_row, array_col, pxl_row_in_fullres, pxl_col_in_fullrespos.columns = ["barcode", "in_tissue", "array_row",               "array_col", "pxl_row_in_fullres",                "pxl_col_in_fullres"]pos = pos.set_index("barcode")pos.index.name = None

obs_names = adata.obs_names.astype(str)if obs_names[0].endswith("-1") and (not str(pos.index[0]).endswith("-1")):    pos.index = pos.index.astype(str) + "-1"elif (not obs_names[0].endswith("-1")) and str(pos.index[0]).endswith("-1"):    pos.index = pos.index.astype(str).str.replace("-1$", "", regex=True)# 只保留 adata 里有的条码顺序pos = pos.reindex(obs_names)

adata.obs = adata.obs.join(pos[["in_tissue", "array_row",                                "array_col", "pxl_row_in_fullres",                                "pxl_col_in_fullres"]])# 重点：obsm["spatial"] 放 fullres 坐标；不要乘 tissue_hires_scalef / lowres_scalef# 坐标顺序：x=col, y=row（更符合图像坐标：横轴列、纵轴行）adata.obsm["spatial"] = pos[["pxl_col_in_fullres",                              "pxl_row_in_fullres"]].to_numpy()

img_candidates = [    ("hires", os.path.join(SPATIAL_DIR, "tissue_hires_image.png")),    ("lowres", os.path.join(SPATIAL_DIR, "tissue_lowres_image.png")),]images = {}for k, p in img_candidates:    if os.path.exists(p):        im = Image.open(p).convert("RGB")        images[k] = (np.asarray(im) / 255.0).astype(np.float32)iflen(images) == 0:    print("WARNING: 没找到 tissue_hires_image.png 或 tissue_lowres_image.png，后续 sc.pl.spatial 只能画点图不叠图。")

adata.uns["spatial"] = {LIBRARY_ID: {}}adata.uns["spatial"][LIBRARY_ID]["scalefactors"] = scalefactorsadata.uns["spatial"][LIBRARY_ID]["images"] = images# 如果同时有 hires/lowres，默认用 hires；没有就用 lowresif "hires" in images:    adata.uns["spatial"][LIBRARY_ID]["use_quality"] = "hires"elif "lowres" in images:    adata.uns["spatial"][LIBRARY_ID]["use_quality"] = "lowres"

print(adata)print("scalefactors:", adata.uns["spatial"][LIBRARY_ID]["scalefactors"])# 叠加组织图（优先 hires）img_key = "hires" if "hires" in images else ("lowres" if "lowres" in images else None)if img_key is not None:    sc.pl.spatial(adata, library_id=LIBRARY_ID, img_key=img_key, spot_size=1.2)else:    # 没有图片也能画（只是没有底图）    sc.pl.embedding(adata, basis="spatial")

# 标记线粒体基因adata.var["mt"] = adata.var_names.str.startswith("Mt-")sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True)print(f"#cells after MT filter: {adata.n_obs}")sc.pp.filter_genes(adata, min_cells=10)sc.pp.normalize_total(adata, inplace=True)sc.pp.log1p(adata)sc.pp.highly_variable_genes(adata, flavor="seurat", n_top_genes=2000)sc.pp.pca(adata)sc.pp.neighbors(adata)sc.tl.umap(adata)sc.tl.leiden(adata, key_added="clusters" )plt.rcParams["figure.figsize"] = (4, 4)sc.pl.umap(adata,            color=["total_counts", "n_genes_by_counts", "clusters"],             wspace=0.4)sc.pl.umap(adata, color=["total_counts" ], wspace=0.4)sc.pl.umap(adata, color=[ "n_genes_by_counts"], wspace=0.4)sc.pl.umap(adata, color=[ "clusters"], wspace=0.4)plt.rcParams["figure.figsize"] = (8, 8)sc.pl.spatial(adata, img_key="hires",               color=["total_counts", "n_genes_by_counts"])sc.pl.spatial(adata, img_key="hires",               color="clusters", size=1.5)sc.tl.rank_genes_groups(adata, "clusters", method="t-test")sc.pl.rank_genes_groups_heatmap(adata, groups="2", n_genes=10, groupby="clusters")sc.pl.spatial(adata, img_key="hires",               color=["clusters" ])

sc.pl.spatial(adata, img_key="hires",               color=["Gapdh"], alpha=0.9)

adata.write_h5ad(r"E:\xy\2026\2.19\A1_visium_fixed.h5ad")