空间聚类分析之cellcharter

生信菜鸟团

发布于 2025-11-19 19:29:29

870

上期推文介绍了IRIS空转邻域算法【空转邻域分析之R包IRIS】。本期将继续带大家了解另一款发表在 Nature Genetics 上的高分空间聚类算法CellCharter。CellCharter于2023年发表在Nature Genetics上，题为《CellCharter reveals spatial cell niches associated with tissue remodeling and cell plasticity》。作为一个基于 Python 的空间组学分析框架，CellCharter能够在多种空间组学平台生成的数据中，对细胞微环境进行识别、特征刻画与跨样本比较。在多项测试中，CellCharter 在识别大规模空间数据中的组织结构方面表现优异，可处理包含数百万个细胞、数百个样本的复杂数据集。

环境部署

官方github在https://github.com/CSOgroup/cellcharter

mamba create -n cellcharter-env -c conda-forge python=3.12
conda activate cellcharter-env
pip install scvi-tools cellcharter

TORCH和CUDA的版本需要根据每个人的电脑的PyTorch和CUDA版本进行调整，例如：

for PyTorch 2.71 and CUDA 12.6, type:

export TORCH=2.7.1
export CUDA=cu126

pip install torch==${TORCH}+${CUDA} --extra-index-url https://download.pytorch.org/whl/${CUDA}
pip install pyg_lib torch_scatter torch_sparse -f https://data.pyg.org/whl/torch-${TORCH}+${CUDA}.html
# 示例：按需替换具体版本号
#pip install -U jax==0.4.27 jaxlib==0.4.27+cuda12.cudnn89 -f https://storage.googleapis.com/jax-releases/jax_cuda_releases.html

# 安装与 CUDA 12 兼容且 <0.7.0 的 JAX 版本（示例：0.6.*）
pip install "jax[cuda12]<0.7.0""jaxlib<0.7.0" \
  -f https://storage.googleapis.com/jax-releases/jax_cuda_releases.html

pip install --force-reinstall \
nvidia-cudnn-cu12==9.5.1.17 \
nvidia-cublas-cu12==12.6.4.1

Jupyter:

mamba install -y nb_conda_kernels ipykernel
python -m ipykernel install --user --name cellcharter-env --display-name cellcharter

代码实战

本教程（https://cellcharter.readthedocs.io/en/latest/notebooks/codex_mouse_spleen.html#）使用 CellCharter 对小鼠脾脏 CODEX 数据集进行空间聚类的示例，该数据集来自 Goltsev et al., 2018。该数据集包含 3 个健康样本（BALBc-1 至 BALBc-3）和 6 个系统性红斑狼疮小鼠样本（MRL-4 至 MRL-9），总计超过 70 万个细胞，检测约 30 个蛋白标记物。

我们将使用 trVAE 进行降维，并用 CellCharter 来对所有样本一起计算空间聚类。

加载所需Python包

import squidpy as sq
import cellcharter as cc
import scanpy as sc
from lightning.pytorch import seed_everything
seed_everything(0)

1. 数据预处理

数据下载：

adata = cc.datasets.codex_mouse_spleen('./data/codex_mouse_spleen.h5ad')
adata

AnnData object with n_obs × n_vars = 707474 × 29
    obs: 'cell_type', 'i-niche', 'tile', 'area', 'dataset', 'stage', 'sample'
    uns: 'spatial', 'spatial_cluster_colors'
    obsm: 'blanks', 'spatial'

我们现在将使用 scanpy 的 pp.scale 函数对每个 marker 的强度值分别进行标准化。

首先，我们将未标准化的数据保存到 anndata 对象的 raw 属性中，以便在需要时可以随时调用。接着，我们会对每个样本单独进行标准化。不过，由于 scanpy 并没有实现这一功能，所以我们需要手动完成。

adata.raw = adata.copy()
for sample in adata.obs['sample'].cat.categories:
    adata.X[adata.obs['sample'] == sample, :] = sc.pp.scale(adata[adata.obs['sample'] == sample], copy=True).X

2. 降维

鉴于当前空间蛋白质组学技术中的 marker 数量有限，严格来说并不一定需要降维。然而，我们注意到降维有助于降低噪音并加快聚类步骤，尤其是在聚类稳定性分析中更为关键。因为在稳定性分析中，为了推荐最佳的聚类数候选，需要多次重复聚类过程。

在这里，我们将使用一个在该数据集上已预训练好的 trVAE 模型。首先加载该模型。在 map_location 参数中可以指定要加载模型的设备。例如，可以在 CPU 上加载一个在 GPU 上训练的模型。事实上，这个模型最初是在 GPU 上训练的，但我们会在 CPU 上加载它。

model = cc.tl.TRVAE.load('./tutorial_models/codex_mouse_spleen_trvae', adata, map_location='cpu')

AnnData object with n_obs × n_vars = 707474 × 29
    obs: 'cell_type', 'i-niche', 'tile', 'area', 'dataset', 'stage', 'sample'
    uns: 'spatial', 'spatial_cluster_colors'
    obsm: 'blanks', 'spatial'

INITIALIZING NEW NETWORK..............
Encoder Architecture:
 Input Layer in, out and cond: 29 128 1
 Hidden Layer 1 in/out: 128 128
 Mean/Var Layer in/out: 128 10
Decoder Architecture:
 First Layer in, out and cond:  10 128 1
 Hidden Layer 1 in/out: 128 128
 Output Layer in/out:  128 29

接下来，我们将为该数据集的所有细胞提取降维后的嵌入表示。

这里我们使用 dataset 标签作为条件，因为该数据集的所有细胞都来自相同的实验条件，因此 无需进行批次效应校正。

adata.obsm['X_trVAE'] = model.get_latent(adata.X, adata.obs['dataset'])

如果你使用的是不同的数据集，则需要训练你自己的 trVAE 模型。你可以使用 cc.tl.TRVAE 来完成训练，这是基于 scArches 中原始 trVAE 模型的一个轻微修改版本（ChatGPT，尚未测试）：

1) 蛋白数据推荐的预处理

import scanpy as sc
import numpy as np

# 保留原始值
adata.layers["counts"] = adata.X.copy()

# (a) 背景/阴性探针去除（若有）
adata = adata[:, [m for m in adata.var_names if"NegPrb"notin m]].copy()

# (b) 变换：arcsinh 或 log1p（二选一即可）
adata.X = np.arcsinh(adata.X / 5.)       # cofactor=5 可按平台调
# 或：
# sc.pp.normalize_total(adata, target_sum=1e6); sc.pp.log1p(adata)

# (c) 按样本做每 marker 的 z-score（减少跨样本偏移）
for s in adata.obs['sample'].astype('category').cat.categories:
    idx = adata.obs['sample'] == s
    Xs = adata[idx].X
    adata.X[idx] = (Xs - Xs.mean(0)) / (Xs.std(0) + 1e-8)

2) 训练你自己的 trVAE（cc.tl.TRVAE）

import cellcharter as cc

# 指定条件键（批次/样本），决定模型如何对齐不同样本
cond_key = "sample"             # 也可用 'batch' 或 'slide' 等你的字段

# 初始化模型（可调隐层/潜变量维度）
model = cc.tl.TRVAE(
    adata,
    condition_key=cond_key,
    n_hidden=128,         # 隐层大小
    z_dim=10,             # 潜空间维度（与教程一致）
    n_layers_encoder=1,
    n_layers_decoder=1,
    dropout_rate=0.1
)

# 训练（可设验证划分/早停）
model.train(
    max_epochs=200,
    early_stopping=True,
    validation_split=0.1,
    batch_size=2048
)

# 导出低维表示
adata.obsm["X_trVAE"] = model.get_latent(adata.X, adata.obs[cond_key])
# 保存模型供复用
model.save("./models/my_protein_trvae")

3）接入 CellCharter（与教程一致）

import squidpy as sq
# 建图：按样本建独立图
sq.gr.spatial_neighbors(adata, library_key='sample', coord_type='generic', delaunay=True, spatial_key='spatial')
cc.gr.remove_long_links(adata)

# 邻域聚合：用刚才的 X_trVAE 作为特征
cc.gr.aggregate_neighbors(
    adata, n_layers=3, use_rep='X_trVAE',
    out_key='X_cellcharter', sample_key='sample'
)

# 自动选 K 并聚类
autok = cc.tl.ClusterAutoK(n_clusters=(2,10), max_runs=10, convergence_tol=1e-3)
autok.fit(adata, use_rep='X_cellcharter')
adata.obs['cluster_cellcharter'] = autok.predict(adata)

4）何时不必训练 trVAE？

样本数极少、批次效应极弱：直接 PCA(n_components=10) 足够：

from sklearn.decomposition import PCA
adata.obsm['X_pca'] = PCA(n_components=10).fit_transform(adata.X)
# 然后用 X_pca 走 aggregate_neighbors → Cluster

3. CellCharter 的空间聚类

现在我们开始计算空间聚类。 CellCharter 将每个样本中的所有细胞编码为一个网络。每个细胞是一个节点，当两个细胞在组织中物理位置上足够接近时，它们之间会被连上一条边。我们可以使用 squidpy 的 gr.spatial_neighbors 函数来构建这个网络。

sq.gr.spatial_neighbors(adata, library_key='sample', coord_type='generic', delaunay=True)

然而，squidpy 使用的 Delaunay 三角剖分方法存在一个缺点，即它可能会在相距较远的细胞之间生成连接。通过放大观察 BALBc-3 样本的某一区域，我们可以清楚地看到这一点。

sq.pl.spatial_scatter(
    adata, 
    color='sample', 
    library_key='sample', 
    img=None, 
    title=['BALBc-3'],
    size=2500,
    connectivity_key='spatial_connectivities',
    edges_width=0.3,
    legend_loc=None,
    crop_coord=(900000, 900000, 1700000, 1700000 ),
    library_id=['BALBc-3']
    )

../_images/6f6bbaff9257aa5bc365317873dc8719c182c56cd95da99800daec6bd851122d.png

我们可以选择一个最大半径，并将其作为参数传递给函数，但这个半径会因数据集的空间分辨率不同而有所差异。为此，我们开发了一种名为 gr.remove_longLinks 的方法，它会移除那些距离超过某个百分位数的连接（默认是第 99 个百分位数）。可以看到，这种方法能显著减少过长的连接，而不会影响其余的正常连接。

cc.gr.remove_long_links(adata)

sq.pl.spatial_scatter(
    adata, 
    color='sample', 
    library_key='sample', 
    img=None, 
    title=['BALBc-3'],
    size=2500,
    connectivity_key='spatial_connectivities',
    edges_width=0.3,
    legend_loc=None,
    crop_coord=(900000, 900000, 1700000, 1700000 ),
    library_id=['BALBc-3']
    )

../_images/95113322ca51e091df58a32d459725198860fc4f1c7816d3a3136467b166c4d1.png

下一步是邻域聚合（neighborhood aggregation）。其过程是：将每个细胞自身的特征，与从其相邻细胞中逐层聚合得到的特征拼接在一起，直到设定的层数 n_layers。聚合函数用于将多个邻居细胞的向量整合为一个单一的特征向量，默认采用 mean 平均函数。

在本例中，我们使用 3 层邻居，因此每个细胞会得到一个长度为 40 的特征向量。其中包含：来自 trVAE 的降维向量（10 维），以及从 3 层邻居中分别聚合得到的 3 个 10 维向量。

cc.gr.aggregate_neighbors(adata, n_layers=3, use_rep='X_trVAE')

我们接下来构建一个高斯混合模型（Gaussian Mixture Model, GMM），并将聚类数设定为 11。这个数值来源于聚类稳定性分析的结果。

gmm = cc.tl.Cluster(
    n_clusters=11, 
    random_state=12345,
    # If running on GPU
    #trainer_params=dict(accelerator='gpu', devices=1)
)

然后，我们将模型拟合到数据上，并预测所有细胞的聚类标签。

gmm.fit(adata, use_rep='X_cellcharter')
adata.obs['spatial_cluster'] = gmm.predict(adata, use_rep='X_cellcharter')

我们可以对部分样本的空间聚类结果进行可视化。在本例中，我们选择了一个健康样本，以及分别来自早期、中期和晚期的样本。

sq.pl.spatial_scatter(
    adata, 
    color='spatial_cluster', 
    library_key='sample', 
    img=None, 
    title=['BALBc-1', 'MRL-4 (early)', 'MRL-8 (intermediate)', 'MRL-9 (late)'],
    size=2500,
    ncols=1,
    library_id=['BALBc-1', 'MRL-4', 'MRL-8', 'MRL-9'],
)

../_images/cceb49546f2d1110c5047aea4208b9b130bc14089859e854a7eafd406d8c369b.png

cc.gr.enrichment(adata, group_key='spatial_cluster', label_key='cell_type')
cc.pl.enrichment(adata, group_key='spatial_cluster', label_key='cell_type', figsize=(4,4), fontsize=6, dot_scale=1)

../_images/3157ba2632601f8eae15caa460665e68bb038ff86453fa4ef9744d8117b1d695.png

4. 邻近性分析（Proximity analysis）

下一步是描述样本中空间簇（spatial clusters）的相对排列情况，并观察在健康状态与疾病状态下，不同空间簇之间的相对邻近关系是否发生变化。我们首先从样本标签中提取条件标签（condition labels）。

adata.obs['condition'] = adata.obs['sample'].str.split('-').str[0].astype('category')

我们选择与 健康状态 相关的样本，并计算这些样本中不同空间簇之间的 邻域富集（neighborhood enrichment）。

adata_balbc = adata[adata.obs['condition'] == 'BALBc']
cc.gr.nhood_enrichment(
    adata_balbc,
    cluster_key='spatial_cluster',
)

cc.pl.nhood_enrichment(
    adata_balbc,
    cluster_key='spatial_cluster',
    annotate=True,
    vmin=-1,
    vmax=1,
    figsize=(3,3),
    fontsize=5,
)

../_images/1d8245e880f864369a4397172d23957c1c1bed42a5e8913445dfc65b7d256fcc.png

我们对 红斑狼疮（lupus）状态 的样本重复进行了同样的分析。

adata_mrl = adata[adata.obs['condition'] == 'MRL']
cc.gr.nhood_enrichment(
    adata_mrl,
    cluster_key='spatial_cluster',
)

cc.pl.nhood_enrichment(
    adata_mrl,
    cluster_key='spatial_cluster',
    annotate=True,
    vmin=-1,
    vmax=1,
    figsize=(3,3),
    fontsize=5,
)

../_images/ebe953955a02f1ed553f61a3a50ecb16f5f907be0f2c650ebe8587ad840b2137.png

虽然可以通过将两幅图并排放置来比较邻域富集的差异，但更方便的方法是使用一幅图来直接显示两种状态之间的差异（即差异邻域富集，differential neighborhood enrichment）。这正是 diff_nhood_enrichment 的功能。可选地，我们还可以为每一对簇计算 p 值并设置显著性阈值。这会增加函数的计算时间，因为需要通过置换检验来计算 p 值，在本例中设定为 n_perms=100 次置换。

作者这里实现了一个解析版本（analytical version）的非差异邻域富集方法（即前一步运行的版本），其速度非常快，否则就需要运行“置换的置换”。此外，diff_nhood_enrichment 还支持多进程并行，其核心数可以通过 n_jobs 参数来设定。

cc.gr.diff_nhood_enrichment(
    adata,
    cluster_key='spatial_cluster',
    condition_key='condition',
    library_key='sample',
    pvalues=True,
    n_jobs=15,
    n_perms=100
)

cc.pl.diff_nhood_enrichment(
    adata,
    cluster_key='spatial_cluster',
    condition_key='condition',
    condition_groups=['MRL', 'BALBc'],
    annotate=True,
    figsize=(3,3),
    significance=0.05,
    fontsize=5
)

../_images/137bfd5a3832610b2d1a2217aea810548c8b126af056e87deb608c2bb259aacb.png

与论文中的结果一致，我们观察到在狼疮脾脏中，B 滤泡（C8，粉色）失去了与边缘区（C7，深蓝色）的邻近关系。

参数 condition_groups 的顺序决定了哪一种状态被视为“观察值”，哪一种被视为“期望值”。因此，如果将顺序反转，就会得到完全相反的趋势。

cc.pl.diff_nhood_enrichment(
    adata,
    cluster_key='spatial_cluster',
    condition_key='condition',
    condition_groups=['BALBc', 'MRL'],
    annotate=True,
    figsize=(3,3),
    significance=0.05,
    fontsize=5
)

../_images/f67c552b7a4a2cac6cf5efbbc1d520daf989ad246cfc26e25b7869e21f5464ab.png

5. 形状特征分析（Shape characterization）

另一种描述组织样本空间结构变化的方法，是通过测量在两种条件下共有的空间簇的形状变化。

由于一个空间簇可能在同一样本中的多个区域重复出现，因此我们需要通过计算空间网络的连通分量（connected components）来找到单独的局部簇。

cc.gr.connected_components(adata, cluster_key='spatial_cluster')

然后，我们计算每个局部簇的边界。

cc.tl.boundaries(adata, min_hole_area_ratio=0.1)

在在本教程中实现了论文中描述的四种度量指标：线性度（linearity）、卷曲度（curl）、伸长度（elongation）和纯度（purity）。在本例中，我们将仅比较前两种指标。

cc.tl.linearity(adata)
cc.tl.curl(adata)

cc.pl.shape_metrics(
    adata,
    condition_key='condition', 
    condition_groups=['BALBc', 'MRL'],
    cluster_key='spatial_cluster', 
    cluster_id=[10], 
    title='C4 - B-PALS shape',
    figsize=(4,3),
    fontsize=6
)