Module: Clustering Stability#
When analyzing (spatial) transcriptomics data, cells are typically clustered into cell types based on their RNA expression. If the segmentation of a spatial transcriptomics dataset went well, we would assume that this clustering is somewhat stable, even if we only cluster on parts of the data. The clustering stability (cs) module includes some functions to apply clustering and assess its robustness. The figure below shows a bad clustering (lots of overlap) vs. a good clustering.
To follow along with this tutorial, you can download the data from here.
[1]:
%load_ext autoreload
%autoreload 2
import warnings
warnings.simplefilter(action="ignore", category=FutureWarning)
[2]:
import importlib
from pathlib import Path
import anndata as ad
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scanpy as sc
import scipy.sparse as sp
import seaborn as sns
import spatialdata as sd
import spatialdata_plot # noqa: F401
import segtraq
segtraq.settings.n_jobs = 8 # Use 8 CPU cores (-1 for all available cores)
[3]:
def filter_zero_count_cells(adata: ad.AnnData) -> ad.AnnData:
"""
Return a view of adata excluding cells with zero total counts.
Does NOT modify the original object.
"""
if sp.issparse(adata.X):
total_counts = np.array(adata.X.sum(axis=1)).flatten()
else:
total_counts = adata.X.sum(axis=1)
mask = total_counts > 0
return adata[mask, :]
[4]:
# example datasets on which we can compare performance
st_bidcell = segtraq.SegTraQ(
sd.read_zarr("../../data/xenium_5K_data/bidcell.zarr"),
images_key="image",
points_background_id=0,
tables_area_key=None,
tables_centroid_x_key="cell_centroid_x",
tables_centroid_y_key="cell_centroid_y",
)
st_segger = segtraq.SegTraQ(
sd.read_zarr("../../data/xenium_5K_data/segger.zarr"),
images_key="image",
tables_centroid_x_key="cell_centroid_x",
tables_centroid_y_key="cell_centroid_y",
)
st_xenium = segtraq.SegTraQ(
sd.read_zarr("../../data/xenium_5K_data/xenium.zarr"),
images_key="image",
tables_centroid_x_key="x_centroid",
tables_centroid_y_key="x_centroid",
)
st_proseg2 = segtraq.SegTraQ(
sd.read_zarr("../../data/xenium_5K_data/proseg2.zarr"),
images_key="image",
points_cell_id_key="assignment",
points_background_id=2**32 - 1,
points_gene_key="gene",
tables_area_key=None,
tables_cell_id_key="cell",
shapes_cell_id_key="cell",
tables_centroid_x_key="centroid_x",
tables_centroid_y_key="centroid_y",
)
st_dict = {"bidcell": st_bidcell, "segger": st_segger, "xenium": st_xenium, "proseg2": st_proseg2}
/g/huber/users/lazic/src/SegTraQ/src/segtraq/SegTraQ.py:160: RuntimeWarning: No area column specified for tables. Area will be automatically computed from shapes.
validate_spatialdata(
/g/huber/users/lazic/src/SegTraQ/src/segtraq/SegTraQ.py:188: UserWarning: Filtering control probes matching prefixes ('NegControlProbe_', 'antisense_', 'NegControlCodeword', 'BLANK_', 'Blank-', 'NegPrb', 'DeprecatedCodeword_', 'UnassignedCodeword_', 'Intergenic_Region_') and transcripts with qv < 20 from transcript points in-place; control probes are also removed from the expression table if present. Configure `min_qv`, `control_prefixes`, and `inplace` via `filter_kwargs`.
sdata_new = _filter_control_and_low_quality_transcripts(
/g/huber/users/lazic/src/SegTraQ/src/segtraq/SegTraQ.py:188: RuntimeWarning: Cell IDs differ between assigned transcript points and the expression table. 0 cells occur only in points (e.g. []), and 84 occur only in the table (e.g. [np.int32(1539), np.int32(15875), np.int32(17923), np.int32(17924), np.int32(17927)]). Cells only occuring in the table can reflect zero-count cells, but differences may also be due to the filtering used to generate the transcript data and expression matrix. Please check that `filter_kwargs`, particularly `min_qv` and `control_prefixes`, are consistent with the filtering used to generate the expression table.
sdata_new = _filter_control_and_low_quality_transcripts(
no parent found for <ome_zarr.reader.Label object at 0x7fff3018ba10>: None
no parent found for <ome_zarr.reader.Label object at 0x7fff12706ad0>: None
/g/huber/users/lazic/src/SegTraQ/src/segtraq/SegTraQ.py:160: UserWarning: Missing 23 cell IDs in shapes: ['idcbpiel-1', 'ipnaeici-1', 'afnajpkp-1', 'ibnnajcn-1', 'ipnffmib-1']... These cells are present in tables, but not in shapes. This might lead to inconsistencies in the spatialdata object.
validate_spatialdata(
/g/huber/users/lazic/src/SegTraQ/src/segtraq/utils.py:2139: UserWarning: Duplicate IDs detected in index 'cell_id' for shapes 'nucleus_boundaries'. Resetting and renaming index to `segtraq_id` to ensure uniqueness.
nucleus_shapes = _ensure_index(
/g/huber/users/lazic/src/SegTraQ/src/segtraq/SegTraQ.py:188: UserWarning: Filtering control probes matching prefixes ('NegControlProbe_', 'antisense_', 'NegControlCodeword', 'BLANK_', 'Blank-', 'NegPrb', 'DeprecatedCodeword_', 'UnassignedCodeword_', 'Intergenic_Region_') and transcripts with qv < 20 from transcript points in-place; control probes are also removed from the expression table if present. Configure `min_qv`, `control_prefixes`, and `inplace` via `filter_kwargs`.
sdata_new = _filter_control_and_low_quality_transcripts(
/g/huber/users/lazic/src/SegTraQ/src/segtraq/SegTraQ.py:188: RuntimeWarning: Cell IDs differ between assigned transcript points and the expression table. 0 cells occur only in points (e.g. []), and 23 occur only in the table (e.g. ['idcbpiel-1', 'ipnaeici-1', 'ibnnajcn-1', 'ipleocce-1', 'japjollp-1']). Cells only occuring in the table can reflect zero-count cells, but differences may also be due to the filtering used to generate the transcript data and expression matrix. Please check that `filter_kwargs`, particularly `min_qv` and `control_prefixes`, are consistent with the filtering used to generate the expression table.
sdata_new = _filter_control_and_low_quality_transcripts(
no parent found for <ome_zarr.reader.Label object at 0x7fff12707110>: None
no parent found for <ome_zarr.reader.Label object at 0x7fff127ec9d0>: None
/g/huber/users/lazic/src/SegTraQ/src/segtraq/utils.py:2139: UserWarning: Duplicate IDs detected in index 'cell_id' for shapes 'nucleus_boundaries'. Resetting and renaming index to `segtraq_id` to ensure uniqueness.
nucleus_shapes = _ensure_index(
/g/huber/users/lazic/src/SegTraQ/src/segtraq/SegTraQ.py:188: UserWarning: Filtering control probes matching prefixes ('NegControlProbe_', 'antisense_', 'NegControlCodeword', 'BLANK_', 'Blank-', 'NegPrb', 'DeprecatedCodeword_', 'UnassignedCodeword_', 'Intergenic_Region_') and transcripts with qv < 20 from transcript points in-place; control probes are also removed from the expression table if present. Configure `min_qv`, `control_prefixes`, and `inplace` via `filter_kwargs`.
sdata_new = _filter_control_and_low_quality_transcripts(
/g/huber/users/lazic/src/SegTraQ/src/segtraq/SegTraQ.py:188: RuntimeWarning: Gene sets differ between assigned transcript points and the expression table. 0 genes occur only in points (e.g. []), and 6 occur only in the table (e.g. ['HPV18-E6', 'F2', 'Human_FGFR3_R248H_ALT:A', 'Human_CTNNB1_S33C_ALT:G', 'CYP2A6']). Genes only occuring in the table can be introduced by cropping (`SpatialData`), but differences may also be due to the filtering used to generate the transcript data and expression matrix. Please check that `filter_kwargs`, particularly `min_qv` and `control_prefixes`, are consistent with the filtering used to generate the expression table.
sdata_new = _filter_control_and_low_quality_transcripts(
/g/huber/users/lazic/src/SegTraQ/src/segtraq/SegTraQ.py:188: RuntimeWarning: Cell IDs differ between assigned transcript points and the expression table. 0 cells occur only in points (e.g. []), and 50 occur only in the table (e.g. ['hlgldogf-1', 'jaeofina-1', 'hjbaimjn-1', 'ibhhmanh-1', 'obnnbpck-1']). Cells only occuring in the table can reflect zero-count cells, but differences may also be due to the filtering used to generate the transcript data and expression matrix. Please check that `filter_kwargs`, particularly `min_qv` and `control_prefixes`, are consistent with the filtering used to generate the expression table.
sdata_new = _filter_control_and_low_quality_transcripts(
/g/huber/users/lazic/src/SegTraQ/src/segtraq/SegTraQ.py:160: RuntimeWarning: No area column specified for tables. Area will be automatically computed from shapes.
validate_spatialdata(
/g/huber/users/lazic/src/SegTraQ/src/segtraq/utils.py:2139: UserWarning: Duplicate IDs detected in index 'cell_id' for shapes 'nucleus_boundaries'. Resetting and renaming index to `segtraq_id` to ensure uniqueness.
nucleus_shapes = _ensure_index(
/g/huber/users/lazic/src/SegTraQ/src/segtraq/SegTraQ.py:188: UserWarning: Filtering control probes matching prefixes ('NegControlProbe_', 'antisense_', 'NegControlCodeword', 'BLANK_', 'Blank-', 'NegPrb', 'DeprecatedCodeword_', 'UnassignedCodeword_', 'Intergenic_Region_') and transcripts with qv < 20 from transcript points in-place; control probes are also removed from the expression table if present. Configure `min_qv`, `control_prefixes`, and `inplace` via `filter_kwargs`.
sdata_new = _filter_control_and_low_quality_transcripts(
/g/huber/users/lazic/src/SegTraQ/src/segtraq/SegTraQ.py:188: RuntimeWarning: Cell IDs differ between assigned transcript points and the expression table. 0 cells occur only in points (e.g. []), and 971 occur only in the table (e.g. [np.int64(16386), np.int64(16390), np.int64(14347), np.int64(10253), np.int64(10254)]). Cells only occuring in the table can reflect zero-count cells, but differences may also be due to the filtering used to generate the transcript data and expression matrix. Please check that `filter_kwargs`, particularly `min_qv` and `control_prefixes`, are consistent with the filtering used to generate the expression table.
sdata_new = _filter_control_and_low_quality_transcripts(
The problem#
Segmentation errors affect downstream analysis in spatial transcriptomics datasets. For example, let’s consider the UMAPs of four different segmentations. You can see below that they differ substantially.
[5]:
fig, axs = plt.subplots(1, 4, figsize=(8, 2))
# Flatten axs in case of single row/column
axs = axs.flatten()
for i, (method, st) in enumerate(st_dict.items()):
# here, we do not want to alter the original data, so we create a deep copy
sdata = st.sdata
adata = sd.deepcopy(sdata.tables["table"])
# this is all done internally by SegTraQ as well before computing the cs metrics,
# so you are not required to run this manually
# normalizing and log-transforming the counts
sc.pp.normalize_total(adata, inplace=True)
sc.pp.log1p(adata)
# computing a PCA and neighbors
sc.pp.pca(adata)
sc.pp.neighbors(adata)
adata = filter_zero_count_cells(adata)
sc.tl.umap(adata)
sc.tl.leiden(adata, flavor="igraph", n_iterations=2)
sc.pl.umap(
adata,
color="leiden",
ax=axs[i],
show=False,
title=method,
legend_loc=None,
)
plt.tight_layout()
plt.show()
/tmp/ipykernel_182991/4287800528.py:14: UserWarning: Some cells have zero counts
sc.pp.normalize_total(adata, inplace=True)
/g/huber/users/lazic/src/SegTraQ/.venv/lib/python3.13/site-packages/scanpy/tools/_umap.py:189: ImplicitModificationWarning: Trying to modify attribute `._uns` of view, initializing view as actual.
adata.uns[key_uns] = dict(params=dict(a=a, b=b))
/g/huber/users/lazic/src/SegTraQ/.venv/lib/python3.13/site-packages/scanpy/tools/_umap.py:189: ImplicitModificationWarning: Trying to modify attribute `._uns` of view, initializing view as actual.
adata.uns[key_uns] = dict(params=dict(a=a, b=b))
/tmp/ipykernel_182991/4287800528.py:14: UserWarning: Some cells have zero counts
sc.pp.normalize_total(adata, inplace=True)
/g/huber/users/lazic/src/SegTraQ/.venv/lib/python3.13/site-packages/scanpy/tools/_umap.py:189: ImplicitModificationWarning: Trying to modify attribute `._uns` of view, initializing view as actual.
adata.uns[key_uns] = dict(params=dict(a=a, b=b))
/tmp/ipykernel_182991/4287800528.py:14: UserWarning: Some cells have zero counts
sc.pp.normalize_total(adata, inplace=True)
/g/huber/users/lazic/src/SegTraQ/.venv/lib/python3.13/site-packages/scanpy/tools/_umap.py:189: ImplicitModificationWarning: Trying to modify attribute `._uns` of view, initializing view as actual.
adata.uns[key_uns] = dict(params=dict(a=a, b=b))
Cluster Connectedness#
Cluster connectedness (CC) measures how well cells with the same cluster assignment are connected in the expression-based neighborhood graph. It computes the fraction of neighbors sharing the same cluster label. In general, the higher this value is, the better.
By default, we compute this metric using Leiden clustering at resolution 0.2. You can adjust this with the resolution parameter. Alternatively, connectedness can be computed using existing cell-type annotations by specifying cell_type_key.
[6]:
ccs = {}
for method, st in st_dict.items():
ccs[method] = st.cs.connectedness(leiden_kwargs={"flavor": "igraph"})
ccs
[6]:
{'bidcell': 0.9747831283931955,
'segger': 0.9079173523451948,
'xenium': 0.9222843180568768,
'proseg2': 0.9320630464896202}
Silhouette Score#
A slightly more elaborate metric is the silhouette score. It measures how similar an object is to its own cluster (cohesion) compared to other clusters (separation). Values range from −1 to +1, where a high value indicates that the object is well matched to its own cluster and poorly matched to neighboring clusters.
By default, we compute this metric using Leiden clustering at resolution 0.2. You can adjust this with the resolution parameter. As for connectedness, the silhouette score can also be computed using existing cell-type annotations by specifying cell_type_key.
[7]:
silhouette_scores = {}
for method, st in st_dict.items():
silhouette_scores[method] = st.cs.silhouette_score(leiden_kwargs={"flavor": "igraph"})
silhouette_scores
[7]:
{'bidcell': -0.10931064188480377,
'segger': 0.05958401411771774,
'xenium': -0.08852291107177734,
'proseg2': -0.06871508061885834}
Purity#
Another way to assess cluster stability is to cluster only on a subset of all cells. For example, if we randomly select 63% (\(1 - e^{-1}\)) of cells and then perform Leiden clustering on those, will our cells typically get assigned to the same cluster or to different ones?
One way to assess this is by comparing the purity between two clusterings. For every cluster in clustering 1, we check how many other clusters it contains in clustering 2.
A purity value of 1 means that the clusters are completely pure, whereas a value closer to 0 means that they are more mixed.
We randomly select 63% of cells, perform clustering on them, and do this five times, to obtain five different clusterings. Then we compare them using the purity score.
Note: this method will recompute the PCA based on a subset of features. If you will be using the PCA in the anndata later, you should recompute it on the whole data.
[8]:
purities = {}
for method, st in st_dict.items():
purities[method] = st.cs.purity(leiden_kwargs={"flavor": "igraph"})
purities
[8]:
{'bidcell': 0.8860854806570154,
'segger': 0.9573111723394222,
'xenium': 0.9138666233610309,
'proseg2': 0.9149728839036125}
Adjusted Rand Index (ARI)#
The adjusted rand index (ARI) is another metric to determine the similarity of different clusterings (again, we create five clusterings based on random subsets of 63% of cells each). Just like the purity, its values range from 0 (no similarity, not a robust clustering) to 1 (exactly the same clusters, high robustness).
[9]:
aris = {}
for method, st in st_dict.items():
aris[method] = st.cs.adjusted_rand_index(leiden_kwargs={"flavor": "igraph"})
aris
[9]:
{'bidcell': 0.6351438532687562,
'segger': 0.7973927932620155,
'xenium': 0.7486167404769389,
'proseg2': 0.7940773354007}
At the end, all of our metrics are stored in the spatialdata object.
[10]:
st_dict["proseg2"].sdata.tables["table"]
[10]:
AnnData object with n_obs × n_vars = 17885 × 5094
obs: 'cell', 'original_cell_id', 'centroid_x', 'centroid_y', 'centroid_z', 'fov', 'cluster', 'volume', 'scale', 'population', 'region', 'cell_area', 'leiden_subset_cells16914_res0.2_seed42', 'leiden_subset_cells10655_res0.2_seed0', 'leiden_subset_cells10655_res0.2_seed1', 'leiden_subset_cells10655_res0.2_seed2', 'leiden_subset_cells10655_res0.2_seed3', 'leiden_subset_cells10655_res0.2_seed4'
uns: 'spatialdata_attrs', 'connectedness_leiden', 'silhouette_score_leiden', 'mean_purity_leiden', 'mean_ari_leiden'
Visualization#
Looking at numbers is one thing, but interpretation will be a lot easier if we visualize our results. The following couple of codeblocks demonstrate how the four methods compare.
[11]:
# putting everything into a dataframe for easier plotting
results_df = pd.DataFrame(
{
"Method": list(ccs.keys()),
"Connectedness": list(ccs.values()),
"Silhouette Score": list(silhouette_scores.values()),
"Purity": list(purities.values()),
"ARI": list(aris.values()),
}
)
[12]:
# plotting the connectedness vs. silhouette score
plt.figure(figsize=(4, 3))
sns.scatterplot(data=results_df, x="Connectedness", y="Silhouette Score", hue="Method", style="Method", s=100)
plt.xlabel("Connectedness (↑)")
plt.ylabel("Silhouette Score (↑)")
plt.title("Connectedness vs. Silhouette Score")
plt.xlim(0, 1)
plt.ylim(-1, 1)
plt.legend(title="Method")
plt.show()
[13]:
# plotting the purity vs. ARI
plt.figure(figsize=(4, 3))
sns.scatterplot(data=results_df, x="Purity", y="ARI", hue="Method", style="Method", s=100)
plt.xlabel("Purity (↑)")
plt.ylabel("ARI (↑)")
plt.title("Purity vs. ARI")
plt.legend(title="Method")
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.show()
Clustering quality using transferred cell-type labels#
Besides Leiden clusters, connectedness and silhouette score can also be computed using cell-type annotations. This allows us to assess how well biologically defined cell types are separated in expression space.
Here, we first perform label transfer from an scRNA-seq reference dataset and then compute both metrics using the transferred cell-type labels.
[14]:
# Running label transfer from an scRNA-seq dataset
adata_ref = ad.read_h5ad("../../data/xenium_5K_data/BC_scRNAseq_Janesick.h5ad")
for _method, st in st_dict.items():
st.run_label_transfer(
adata_ref,
ref_cell_type="celltype_major",
ref_raw_counts_layer="raw",
inplace=True,
)
WARNING: adata.X seems to be already log-transformed.
WARNING: adata.X seems to be already log-transformed.
WARNING: adata.X seems to be already log-transformed.
WARNING: adata.X seems to be already log-transformed.
We can now compute connectedness and silhouette score using the transferred cell-type annotations instead of Leiden clusters by specifying cell_type_key.
[15]:
ccs_labels = {}
silhouette_scores_labels = {}
for method, st in st_dict.items():
ccs_labels[method] = st.cs.connectedness(cell_type_key="transferred_cell_type")
silhouette_scores_labels[method] = st.cs.silhouette_score(cell_type_key="transferred_cell_type")
[16]:
results_labels_df = pd.DataFrame(
{
"Method": list(ccs_labels.keys()),
"Connectedness": list(ccs_labels.values()),
"Silhouette Score": list(silhouette_scores_labels.values()),
}
)
results_labels_df
[16]:
| Method | Connectedness | Silhouette Score | |
|---|---|---|---|
| 0 | bidcell | 0.577303 | -0.004466 |
| 1 | segger | 0.460573 | -0.028275 |
| 2 | xenium | 0.623533 | 0.006978 |
| 3 | proseg2 | 0.703557 | 0.001204 |
Session Info#
[17]:
print(importlib.metadata.version("spatialdata"))
print(importlib.metadata.version("spatialdata_plot"))
print(importlib.metadata.version("scanpy"))
0.8.0
0.4.0
1.12.3