10X空转visium HD多样本分析策略二【scanpy+spatialdata】

上一节我们介绍了python版visium HD数据的一种分析策略（python版10X空转visium HD分析策略一【scanpy+squidpy】）。其实是将visium HD的bins当作visium的spots进行的分析，虽然没有太大问题，但毕竟不一样，基础分析没有任何问题，后续继续其他分析可能不太行。这里我们介绍visium HD分析的第二种策略，也是官方推荐的主流方式，即【scanpy+spatialdata】。SpatialData的优势主要体现在：多尺度空间表达 + 几何对象统一管理 + 下游空间操作能力。此外spatialdata能够支持更大数据量的分析。在实践过程中发现，scanpy+squidpy分析visium HD数据可能更符合大家先前的分析习惯及流程，且对于文件的要求与visium spots一致，对于单样本或者少量样本分析是不错的选择；但是SpatialData对于更多样本的分析具有能力，此外其对input文件的需求有额外的，接下来会介绍到。

1-SpatialData简介及安装相关分析库

SpatialData是一个数据框架，包含一个符合FAIR原则（一套针对科研数据与软件的国际通用规范）的存储格式，以及一组Python库，用于对单模态和多模态空间组学数据进行高性能访问、对齐与处理。SpatialData数据格式与之前的anndata是不同的。

库的论文发表在nature methods：
Marconato, L., Palla, G., Yamauchi, K.A. et al. SpatialData: an open and universal data framework for spatial omics. Nat Methods (2024). https:///10.1038/s41592-024-02212-x

Github链接：
https://github.com/scverse/spatialdata

#请安装相关的库#最好新建一个环境进行分析，避免包冲突#pip install squidpy -i https://pypi.tuna./simple#pip install spatialdata -i https://pypi.tuna./simple#pip install geosketch -i https://pypi.tuna./simple#pip install spatialdata_plot -i https://pypi.tuna./simple#pip install spatialdata-io -i https://pypi.tuna./simple

import scanpy as scimport numpy as npimport pandas as pdimport squidpy as sqimport matplotlib.pyplot as pltimport spatialdata as spdimport spatialdata_plot as spltimport spatialdata_io as soimport geosketch as sketchimport jsonimport gcfrom PIL import Image

首先再认识一下spaceranger visium HD数据output文件，如果使用scanpy+squidpy读取，只需要binned_outputs文件夹的数据即可，加载需要分辨率的数据分析即可；spatialdata对于input数据的要求需要有binned_outputs文件夹，以及feature_slice.h5【这是必须的】，否则其读取函数会报错！这里我们使用的数据还是10X的数据，文章在https://www./articles/s41588-025-02193-3,我们下载了两例数据用于演示多样本分析，结直肠癌数据以对照数据https://www./platforms/visium/product-family/dataset-human-crc。

cd P5_CRC├── binned_outputs│   ├── square_002um│   ├── square_008um│   └── square_016um└── feature_slice.h5

cd P5_NAT.├── binned_outputs│   ├── square_002um│   │   ├── filtered_feature_bc_matrix│   │   ├── filtered_feature_bc_matrix.h5│   │   ├── raw_feature_bc_matrix│   │   ├── raw_feature_bc_matrix.h5│   │   ├── raw_probe_bc_matrix.h5│   │   └── spatial│   ├── square_008um│   │   ├── analysis│   │   ├── cloupe.cloupe│   │   ├── filtered_feature_bc_matrix│   │   ├── filtered_feature_bc_matrix.h5│   │   ├── raw_feature_bc_matrix│   │   ├── raw_feature_bc_matrix.h5│   │   └── spatial│   └── square_016um│       ├── analysis│       ├── filtered_feature_bc_matrix│       ├── filtered_feature_bc_matrix.h5│       ├── raw_feature_bc_matrix│       ├── raw_feature_bc_matrix.h5│       └── spatial└── feature_slice.h5

2-数据读取

首先我们单独读取一个样本演示读入函数用法,以及spatial object的保存与再次读入，并探究一下SpatialData数据格式：数据的读取使用so.visium_hd函数。函数使用简单，提供数据路径及分辨率即可，有点像seurat读取visium HD的赶脚！

sdata_crc = so.visium_hd(    path="./P5_CRC/",#理想情况下是spaceranger outs文件夹路径，这里就是包含了binned_outputs和feature_slice.h5文件路径    dataset_id = 'P5_CRC',#sample识别ID，多样本可用于区分sample    filtered_counts_file=True,#使用filtered_feature_bc_matrix.h5，过滤后的表达矩阵    bin_size='008',#选择square-bin的分辨率，一般有2，8，16三个分辨率数据。如果加载多个分辨率数据，传入一个列表['008','016']    fullres_image_file=None,#全分辨率图像的路径，这里选择不加载，数据太大    load_all_images=True,#如果设置为False，则只加载hires/lowres image，不加载cytassist_image    var_names_make_unique=True #确保基因名称唯一性)sdata_crc

SpatialData object├── Images│     ├── 'P5_CRC_cytassist_image': DataArray[cyx] (3, 3000, 3199)│     ├── 'P5_CRC_hires_image': DataArray[cyx] (3, 5298, 6000)│     └── 'P5_CRC_lowres_image': DataArray[cyx] (3, 530, 600)├── Shapes│     └── 'P5_CRC_square_008um': GeoDataFrame shape: (541968, 1) (2D shapes)└── Tables      └── 'square_008um': AnnData (541968, 18085)with coordinate systems:    ▸ 'P5_CRC', with elements:        P5_CRC_cytassist_image (Images), P5_CRC_hires_image (Images), P5_CRC_lowres_image (Images), P5_CRC_square_008um (Shapes)    ▸ 'P5_CRC_downscaled_hires', with elements:        P5_CRC_hires_image (Images), P5_CRC_square_008um (Shapes)    ▸ 'P5_CRC_downscaled_lowres', with elements:        P5_CRC_lowres_image (Images), P5_CRC_square_008um (Shapes)

可以看到，SpatialData object是由几部分组成的，Iamge存储了组织切片的显微镜图像，通常用于可视化或与空间点对齐。Shapes是quare-bin 的几何边界。Tables是anndata数据，是表达矩阵，用于后续数据分析。

多样本数据处理：

#字典传入不同sample信息#这里字典的key最好设置为样本名称samples = {    "P5_CRC":["./P5_CRC/"],    "P5_NAT":["./P5_NAT/"]}
sdatas = []#循环读取for key, inputs in samples.items():
    sdata = so.visium_hd(path=inputs[0],                 dataset_id = key,                 filtered_counts_file=True,                 bin_size='008',                 fullres_image_file=None,                 load_all_images=True,                 var_names_make_unique=True)
    for table in sdata.tables.values():        table.var_names_make_unique()        table.obs["sample"] = key
    sdatas.append(sdata)

合并数据：这里的合并不是真正意义上的合并（不是图片合并哦），可以说是拼接不同样本SpatialData对象！

comb_sdata = spd.concatenate(sdatas, concatenate_tables=True)#concatenate_tables为True表示anndata也合并comb_sdata.write("comb_sdata", overwrite=True)#做好保存数据

3-数据质控及降维聚类

这一部分内容就是相当于进行了scRNA的步骤【还是一样的，数据质控标准需要结合自己实际的数据】数据的质控当然是针对表达矩阵的数据，需要使用scanpy，和其流程一样的处理！去处一些低质量的bins。

首先提取adata，上面的数据结构介绍中以及知晓，表达矩阵以及metadata储存在table。
Tables
└── 'square_008um': AnnData (977741, 18085)

adata = comb_sdata["square_008um"] adata

AnnData object with n_obs × n_vars = 977741 × 18085    obs: 'in_tissue', 'array_row', 'array_col', 'location_id', 'region', 'sample'    uns: 'spatialdata_attrs'    obsm: 'spatial'

#计算线粒体基因比例adata.var["mt"] = adata.var_names.str.startswith(("MT-"))sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], inplace=True, percent_top=None)#plot QC前sc.pl.violin(adata=adata, keys=["total_counts"], stripplot=False, inner="box",groupby="sample",)sc.pl.violin(adata=adata, keys=["n_genes_by_counts"], groupby="sample", stripplot=False, inner="box")sc.pl.violin(adata=adata, keys=["pct_counts_mt"], groupby="sample", stripplot=False, inner="box")

#过滤基因细胞，线粒体基因比例不做过滤了，仁者见仁，智者见智，也需要结合生物学意义#有些地方认为，空转线粒体基因比例高低可能具有生物学意义，不建议过滤。sc.pp.filter_genes(adata, min_cells=50)sc.pp.filter_cells(adata, min_counts=50)sc.pp.filter_cells(adata, max_counts=3000)

#plot QC后sc.pl.violin(adata=adata, keys=["total_counts"], stripplot=False, inner="box",groupby="sample",)sc.pl.violin(adata=adata, keys=["n_genes_by_counts"], groupby="sample", stripplot=False, inner="box")

#标准化sc.pp.normalize_total(adata, target_sum = None)sc.pp.log1p(adata)sc.tl.pca(adata)# Elbow plotsc.pl.pca_variance_ratio(adata, log=True,n_pcs=50)

运行到这里，就涉及一个去批次及整合的问题了，因为这里演示的是肿瘤和正常对照样本，所以不适合去批次，会让一些生物学差异抹除。事实上，对于空转数据整合矩阵是否去批次还是存在争议的，也需要看具体的数据。如果是连续切片，或者同类组织技术重复整合去批次没有问题，如果需要进行空间组织结构比较，不太适合去批次。如果需要去批次，按照scanpy流程介绍的方法选择运行。

因为visium HD bins多，数据很大，看这里我们整合的两个样本包含949126个bins，那么对于下游的处理能力又很大的要求，如果样本更多的数据整合进行分析，负担会非常重。在python中，也可以像seurat中提到的方式一样，对数据进行下采样，然后映射回整体，可以有效节约空间时间。如果不需要下采样，直接按照流程run即可。这里测试一下服务器能力，不采样直接运行，看看能不能行：在python中，完全可以处理，R中不敢想象！

sc.pp.neighbors(adata, n_neighbors=30, use_rep="X_pca",metric="correlation")sc.tl.leiden(adata, flavor="igraph", key_added="clusters", resolution=0.5,random_state=0)sc.tl.umap(adata,min_dist=0.5, spread=2, random_state=0)

# Plot UMAPsc.pl.umap(adata, color=["clusters"], title="UMAP by Clusters")sc.pl.umap(adata, color=["sample"], title="UMAP by Sample")

4-空间组织可视化聚类

这里可能会有一个关心的问题，我们做降维聚类的时候是单独提取了adata出来进行的，那么这个结果需要返回spatialdata obj吗？其实您关注一下就会知道，结果是直接返回的，之前的comb_sdata中tabels已经发生了改变！所以不需要额外的操作！

comb_sdata

SpatialData object, with associated Zarr store: /home/data_analysis/10X_space/CRC_visiumHD/comb_sdata├── Images│     ├── 'P5_CRC_cytassist_image': DataArray[cyx] (3, 3000, 3199)│     ├── 'P5_CRC_hires_image': DataArray[cyx] (3, 5298, 6000)│     ├── 'P5_CRC_lowres_image': DataArray[cyx] (3, 530, 600)│     ├── 'P5_NAT_cytassist_image': DataArray[cyx] (3, 3004, 3200)│     ├── 'P5_NAT_hires_image': DataArray[cyx] (3, 4153, 6000)│     └── 'P5_NAT_lowres_image': DataArray[cyx] (3, 415, 600)├── Shapes│     ├── 'P5_CRC_square_008um': GeoDataFrame shape: (541968, 1) (2D shapes)│     └── 'P5_NAT_square_008um': GeoDataFrame shape: (435773, 1) (2D shapes)└── Tables      └── 'square_008um': AnnData (858806, 17860)with coordinate systems:    ▸ 'P5_CRC', with elements:        P5_CRC_cytassist_image (Images), P5_CRC_hires_image (Images), P5_CRC_lowres_image (Images), P5_CRC_square_008um (Shapes)    ▸ 'P5_CRC_downscaled_hires', with elements:        P5_CRC_hires_image (Images), P5_CRC_square_008um (Shapes)    ▸ 'P5_CRC_downscaled_lowres', with elements:        P5_CRC_lowres_image (Images), P5_CRC_square_008um (Shapes)    ▸ 'P5_NAT', with elements:        P5_NAT_cytassist_image (Images), P5_NAT_hires_image (Images), P5_NAT_lowres_image (Images), P5_NAT_square_008um (Shapes)    ▸ 'P5_NAT_downscaled_hires', with elements:        P5_NAT_hires_image (Images), P5_NAT_square_008um (Shapes)    ▸ 'P5_NAT_downscaled_lowres', with elements:        P5_NAT_lowres_image (Images), P5_NAT_square_008um (Shapes)

comb_sdata.pl.render_images("P5_CRC_hires_image").pl.render_shapes("P5_CRC_square_008um", color="clusters").pl.show(coordinate_systems="P5_CRC_downscaled_lowres", title='P5_CRC')

comb_sdata.pl.render_images("P5_NAT_hires_image").pl.render_shapes("P5_NAT_square_008um", color="clusters").pl.show(coordinate_systems="P5_NAT_downscaled_lowres", title='P5_NAT')

还可以可视化吧局部区域，需要设置坐标系范围！为例高清展示，这里使用P5_CRC_hires_image高分辨率图片。需要注意，坐标系相比于lower的需要✖️10.使用spd.bounding_box_query函数截取坐标系，min_coordinate=[3750, 3750],max_coordinate=[5000, 4250],axes=("x", "y")3个参数依次对应x，y坐标范围！

spd.bounding_box_query(        comb_sdata,        min_coordinate=[3750, 3750],        max_coordinate=[5000, 4250],        axes=("x", "y"),        target_coordinate_system='P5_CRC_downscaled_hires').pl.render_images("P5_CRC_hires_image").pl.show(coordinate_systems="P5_CRC_downscaled_hires", title='P5_CRC')

spd.bounding_box_query(        comb_sdata,        min_coordinate=[3750, 3750],        max_coordinate=[5000, 4250],        axes=("x", "y"),        target_coordinate_system='P5_CRC_downscaled_hires').pl.render_images("P5_CRC_hires_image").pl.render_shapes("P5_CRC_square_008um", color="clusters").pl.show(coordinate_systems="P5_CRC_downscaled_hires", title='P5_CRC')

5-marker genes & annotation

当然不能直接按照scRNA那样进行注释，不过cluster marker可以粗略的看看组织结构celltype分布。更精细的需要进行细胞分割或者反卷积！

#cluster marker genessc.tl.rank_genes_groups(adata = adata, groupby="clusters", method="wilcoxon")

# cannonical markers for annotationmarker_genes = {    "Fibroblasts": ["COL1A1", "MMP2", 'VIM'],    "Goblet cells": ["FCGBP", "MUC2", "CLCA1"],    "Enterocyte":["EPCAM","KRT8","KRT20","FABP2"],    "Plasma B cells":['JCHAIN','IGKC','IGHG1'],    "B cell":['MS4A1','CD74','CD19','CD22'],    "Smooth muscle":['TAGLN','DES','MYH11'],    "Tumor":["CEACAM6",'REG1B',"REG1A"]}# Plot dotplot for initial cluster assessmentsc.pl.dotplot(adata = adata, var_names = marker_genes, groupby="clusters", standard_scale="var")

comb_sdata.pl.render_images("P5_CRC_hires_image").pl.render_shapes("P5_CRC_square_008um",color="CEACAM6").pl.show(coordinate_systems="P5_CRC_downscaled_lowres", title='P5_CRC')

之后我们会继续演示细胞分割结果的读入及geosketch下采样聚类分析！

觉得分享有用的点个赞再走呗！