基于Python的Xenium空转分析流程

Xenium是由 10x Genomics 推出的一种空间转录组学技术，能够在单细胞分辨率的基础上，精确地检测和定位组织切片中的基因表达情况。这项技术结合了分子检测与空间定位的优势，为研究细胞在组织中的位置及其功能关系提供了重要工具。上期推文【跟着Seurat官网学Xenium空转分析】介绍了Xenium的Seurat分析流程，本期推文介绍一下Xenium的Python流程。这里推荐使用Squidpy进行分析，官网教程见https://squidpy.readthedocs.io/en/stable/notebooks/tutorials/tutorial_xenium.html

首先是环境部署：

mamba create -n squid python=3.10
conda activate squid
pip install squidpy -i https://mirrors.aliyun.com/pypi/simple

## 其他安装方式
#mamba install -c conda-forge squidpy
#pip install 'squidpy[interactive]'
#pip install git+https://github.com/scverse/squidpy@main

如果需要衔接Jupyter notebook使用的话，需要在squid这个环境里安装以下几个插件，然后就能在Jupyter notebook里选择squid这个环境：

mamba install -y nb_conda_kernels ipykernel
python -m ipykernel install --user --name squid --display-name squid

一. 读入数据

import spatialdata as sd
from spatialdata_io import xenium

import matplotlib.pyplot as plt
import seaborn as sns

import scanpy as sc
import squidpy as sq

在下载并将示例数据集（https://www.10xgenomics.com/datasets/preview-data-ffpe-human-lung-cancer-with-xenium-multimodal-cell-segmentation-1-standard）提取到名为Xenium的目录后，我们指定数据集的路径以及我们想要存储Zarr文件的位置：

xenium_path = "./Xenium/"
zarr_path = "./Xenium.zarr"
#读得比较慢
sdata = xenium(xenium_path)
sdata.write(zarr_path)

## 读入数据，sd.read_zarr读得比较快
sdata = sd.read_zarr(zarr_path)
sdata

adata = sdata.tables["table"]
adata
#AnnData object with n_obs × n_vars = 161000 × 480
#    obs: 'cell_id', 'transcript_counts', 'control_probe_counts', 'control_codeword_counts', 'unassigned_codeword_counts', #'deprecated_codeword_counts', 'total_counts', 'cell_area', 'nucleus_area', 'region', 'z_level', 'nucleus_count', 'cell_labels'
#    var: 'gene_ids', 'feature_types', 'genome'



    
#    uns: 'spatialdata_attrs'
#    obsm: 'spatial'

adata.obs
adata.obsm["spatial"]
#array([[ 111.34860229,   28.03543472],
#       [ 120.43323517,   75.61804199],
#       [ 126.72792816,   56.90501404],
#       ...,
#       [7191.46337891,  398.10839844],
#       [7185.16210938,  400.09918213],
#       [7179.02294922,  394.43218994]])

Squidpy 会在 .obsm["spatial"] 中查找细胞坐标。通过使用来自 spatialdata-io 的 Xenium 读取器，这些坐标会自动设置。对于更复杂的数据，需要手动设置细胞坐标。有关更多详细信息，可以参考 spatialdata 教程（https://spatialdata.scverse.org/en/latest/tutorials/notebooks/notebooks/examples/squidpy_integration.html）。

二. 质量控制

使用scanpy.pp.calculate_qc_metrics计算质量控制指标：

sc.pp.calculate_qc_metrics(adata, percent_top=(10, 20, 50, 150), inplace=True)

#控制探针和控制代码词的百分比可以从 adata.obs 中计算得出。
cprobes = (
    adata.obs["control_probe_counts"].sum() / adata.obs["total_counts"].sum() * 100
)
cwords = (
    adata.obs["control_codeword_counts"].sum() / adata.obs["total_counts"].sum() * 100
)
print(f"Negative DNA probe count % : {cprobes}")
print(f"Negative decoding count % : {cwords}")

接下来，我们绘制了每个细胞的总转录本的分布情况，每个细胞独特的转录本，以及分割细胞的面积和细胞核面积与细胞的比例：

fig, axs = plt.subplots(1, 4, figsize=(15, 4))

axs[0].set_title("Total transcripts per cell")
sns.histplot(
    adata.obs["total_counts"],
    kde=False,
    ax=axs[0],
)

axs[1].set_title("Unique transcripts per cell")
sns.histplot(
    adata.obs["n_genes_by_counts"],
    kde=False,
    ax=axs[1],
)


axs[2].set_title("Area of segmented cells")
sns.histplot(
    adata.obs["cell_area"],
    kde=False,
    ax=axs[2],
)

axs[3].set_title("Nucleus ratio")
sns.histplot(
    adata.obs["nucleus_area"] / adata.obs["cell_area"],
    kde=False,
    ax=axs[3],
)

然后进行过滤：

细胞水平的过滤：使用scanpy.pp.filter_cells根据所需的最小count数筛选细胞。
基因水平的过滤：使用scanpy.pp.filter_genes根据所需的最小细胞数过滤基因。

两者的参数都是根据上面的图指定的。它们被设置为分别过滤出最小count数和最小细胞的细胞和基因。其他筛选标准还可以是细胞面积、DAPI信号或最少的唯一转录本。

sc.pp.filter_cells(adata, min_counts=10)
sc.pp.filter_genes(adata, min_cells=5)

三. 标准分析流程

这里和scanpy单细胞/空转流程差不多，具体为：

使用scanpy.pp.normalize_total进行标准化，
使用sc.pp.log1p函数取对数，
sc.pp.pca做主成分分析，
sc.pp.neighbors计算邻域图，sc.tl.umap降维，sc.tl.leiden聚类。

adata.layers["counts"] = adata.X.copy()
sc.pp.normalize_total(adata, inplace=True)
sc.pp.log1p(adata)
sc.pp.pca(adata)
sc.pp.neighbors(adata)
sc.tl.umap(adata)
sc.tl.leiden(adata)

可视化检查结果：

sc.pl.umap(
    adata,
    color=[
        "total_counts",
        "n_genes_by_counts",
        "leiden",
    ],
    wspace=0.4,
)

sq.pl.spatial_scatter(
    adata,
    library_id="spatial",
    shape=None,
    color=[
        


    
"leiden",
    ],
    wspace=0.4,
)

四. 空间统计分析

此外，Squidpy 官方教程展示了如何基于空间领域图和细胞类型注释计算中心性评分。

计算的评分都是空间图的描述性统计信息，包括以下几种：

接近中心性（closeness centrality）：衡量某一群体与其他节点之间的接近程度。
聚类系数（clustering coefficient）：衡量节点之间聚集的程度。
度中心性（degree centrality）：群体成员与非群体成员之间连接的比例。

该数据集包含存储在 adata.obs 中的 Leiden 聚类组注释，这些注释用于计算中心性评分。

4.1 构建空间邻域图

首先，基于空间坐标信息，我们可以使用squidpy.gr.spatial_neighbors函数计算连通性矩阵（connectivity matrix）来计算中心性得分。然后使用 coord_type="generic"函数并通过指定delaunay=True使用Delaunay三角剖分对邻居进行分类。

sq.gr.spatial_neighbors(adata, coord_type="generic", delaunay=True)

4.2 计算中心性得分

基于Leiden识别到的聚类，使用 squidpy.gr.centrality_scores函数计算中心性得分：

sq.gr.centrality_scores(adata, cluster_key="leiden")

结果通过使用squidpy.pl.centrality_scores函数进行可视化，主要绘制平均中心性、接近中心性和度中心性等。




    
sq.pl.centrality_scores(adata, cluster_key="leiden", figsize=(16, 5))

4.3 计算共现概率

本示例展示了如何计算共现概率。

我用AI查了一下什么是共现概率，以下仅供参考：
共现概率用于衡量某一簇（或细胞类型）在空间上是否倾向于出现在另一簇（或另一种细胞类型）附近。它可以帮助研究特定细胞类型之间的相互作用及空间分布规律。
共现概率 >1：目标簇倾向于与条件簇空间上共同出现（存在正关联）。
共现概率 ≈1：两簇的空间分布无显著关联。
共现概率 <1：目标簇倾向于与条件簇空间上不相邻（存在负关联）。

共现得分定义为：

我们可以使用 squidpy.gr.co_occurrence 函数计算共现得分。共现概率比值的结果可以通过 squidpy.pl.co_occurrence 可视化。

此外，我们可以通过 squidpy.pl.spatial_scatter 在空间坐标中进一步可视化组织结构，并叠加与 Leiden 聚类一致的表达基因。

首先，我们创建一个表，这是原始AnnData对象的子集，并将其存储在表格中：

sdata.tables["subsample"] = sc.pp.subsample(adata, fraction=0.5, copy=True)
adata_subsample = sdata.tables["subsample"]

sq.gr.co_occurrence(
    adata_subsample,
    cluster_key="leiden",
)
sq.pl.co_occurrence(
    adata_subsample,
    cluster_key="leiden",
    clusters="12",
    figsize=(10, 10),
)
sq.pl.spatial_scatter(
    adata_subsample,
    color="leiden",
    shape=None,
    size=2,
)

如何解读图中数据？

横轴（distance）：表示目标簇与条件簇之间的距离。距离越小，表示细胞更接近；距离增大时，表示计算更远范围内的共现概率。
纵轴（value）：表示共现概率值，越高的值表示目标簇在距离条件簇更近的区域有更高的出现概率。趋于平坦时，说明目标簇与条件簇的空间关系已无显著差异（趋于全局分布的平均概率）。
不同曲线（颜色/图例）：每条曲线表示一个特定簇（或细胞类型）的共现概率随距离变化的趋势。

从图中得到的结论：

曲线的峰值位置：表示目标簇在条件簇的某一特定距离范围内共现概率最高。例如，距离为 500 时某条曲线达到最大值，说明在距离 500 微米的范围内，两种簇之间共现概率最高。
曲线的下降速度：下降越快，表示目标簇与条件簇的空间关联越局限于短距离。缓慢下降则表示目标簇与条件簇的关联在较大范围内均匀存在。
曲线趋于平坦：当距离较大（如 5000 微米）时，所有簇的共现概率接近 1，说明在远距离下，空间关系趋于随机分布。

4.4 邻域富集分析

这个例子展示了如何运行邻居富集分析。它根据细胞簇连接图的接近度计算富集分数。将观察到的事件数量与排列进行比较，并计算z-score。

该数据集包含位于 adata.obs 中的细胞类型注释，这些注释用于计算邻域富集。我们使用 squidpy.gr.nhood_enrichment 计算邻域富集得分。使用squidpy.pl.nhood_enrichment函数可视化结果：

sq.gr.nhood_enrichment(adata, cluster_key="leiden")

fig, ax = plt.subplots(1, 2, figsize=(13, 7))
sq.pl.nhood_enrichment(
    adata,
    cluster_key="leiden",
    figsize=(8, 8),
    title="Neighborhood enrichment adata",
    ax=ax[0],
)
sq.pl.spatial_scatter(adata_subsample, color="leiden", shape=None, size=2, ax=ax[1])

4.5 计算Moran's I得分

本示例展示了如何计算 Moran's I 全局空间自相关统计量。

Moran's I 全局空间自相关统计量用于评估特征（例如基因）在组织中是否呈现聚集、分散或随机的分布模式。

我们可以使用 squidpy.gr.spatial_autocorr 并设置 mode='moran' 来计算 Moran' I 得分。在此之前，需要通过 squidpy.gr.spatial_neighbors 构建一个空间图。此外，我们还会对需要评估的基因进行子集筛选。

Moran's I 是一种常用的空间自相关统计指标，用于评估某种特征（如基因表达）的空间分布模式是聚集（clustered）、分散（dispersed）还是随机（random）。它是描述空间依赖性的统计量。

sq.gr.spatial_neighbors(adata_subsample, coord_type="generic", delaunay=True)
sq.gr.spatial_autocorr(
    adata_subsample,
    mode="moran",
    n_perms=100,
    n_jobs=1,
)
adata_subsample.uns["moranI"].head(10)

我们可以使用 squidpy.pl.spatial_scatter 可视化其中的一些基因。我们还可以通过设置 mode='geary' 来计算一个与之密切相关的自相关统计量，即 Geary's C。有关更多信息，请参阅 squidpy.gr.spatial_autocorr 的文档。

sq.pl.spatial_scatter(
    adata_subsample,
    library_id="spatial",
    color=[
        "AREG",
        "MET",
    ],
    shape=None,
    size=2,
    img=False,
)

对于绘图，我们还可以使用spatialdata-plot:

import spatialdata_plot

gene_name = ["AREG", "MET"]
for name in gene_name:
    sdata.pl.render_images("morphology_focus").pl.render_shapes(
        "cell_circles",
        color=name,
        table_name="table",
        use_raw=False,
    ).pl.show(
        title=f"{name} expression over Morphology image",
        coordinate_systems="global",
        figsize=(10, 5),
    )

4.6 使用napri-spatialdata进行交互分析和可视化

此外，可以使用napari-spatialdata通过交互式GUI可视化Xenium数据：

from napari_spatialdata import Interactive

Interactive(sdata)

图中显示了所有细胞中AREG的表达：