ELIXIR Spatial course
Domain identification and niche reconstruction
Day 3 - Practical 8
8th October 2026
Benedetta Manzato (LUMC) & Emma Vonk (LUMC)
- Graph construction and analysis of spatial transcriptomics data using Squidpy
- Spatial Domain detection with CellCharter
- Spatial Clustering with BANKSY
import os
import warnings
import anndata as ad
import pandas as pd
import scanpy as sc
import squidpy as sq
warnings.filterwarnings("ignore")
/app/.pixi/envs/default/lib/python3.12/site-packages/docrep/decorators.py:43: SyntaxWarning: 'n_jobs' is not a valid key! doc = func(self, args[0].__doc__, *args[1:], **kwargs) /app/.pixi/envs/default/lib/python3.12/site-packages/docrep/decorators.py:43: SyntaxWarning: 'show_progress_bar' is not a valid key! doc = func(self, args[0].__doc__, *args[1:], **kwargs)
Dataset: Xenium In Situ Analysis of Alzheimer's Disease Mouse Model Brain Coronal Sections from One Hemisphere Over a Time Course¶
The dataset consists of 6 coronal mouse brain slices from 2 different conditions (wildtype - ctrl vs TgCRND8 - AD) across 3 timepoints. In this practical, we have information about cell types available in adata.obs['cell_types']. Please note that these annotation are not perfect. For example, there are quite some cells that could not be assigned to a cell type ("Unknown"). These annotations have been made with on leiden clustering and marker genes reported in this document. (adata_loading_and_annotation.ipynb)
In this practical we aim to understand the differences of the mouse brain between the two conditions and across the timepoints using niches and spatial domains.
# load adata
base_dir = '../data/Practical6/'
adata = sc.read_h5ad(f'{base_dir}xenium_mouse_ad_annotated_rotated.h5ad')
adata
AnnData object with n_obs × n_vars = 351714 × 347
obs: 'cell_id', 'transcript_counts', 'control_probe_counts', 'control_codeword_counts', 'unassigned_codeword_counts', 'total_counts', 'cell_area', 'nucleus_area', 'region', 'cell_labels', 'condition', 'time', 'batch_key', 'leiden_res_01', 'celltype', 'leiden_res_02', 'sample'
var: 'gene_ids', 'feature_types', 'genome'
uns: 'celltype_colors', 'dea_leiden_res_01', 'dendrogram_leiden_res_01', 'leiden_res_01', 'leiden_res_01_colors', 'leiden_res_02', 'log1p', 'neighbors', 'pca', 'umap'
obsm: 'X_pca', 'X_umap', 'spatial'
varm: 'PCs'
obsp: 'connectivities', 'distances'
layers: None (.X)
# Creating a DataFrame from 'split', 'fov', and 'condition'
df = adata.obs[['condition', 'time', 'batch_key']]
value_counts = pd.DataFrame(df.values, columns=df.columns).value_counts()
value_counts
condition time batch_key
TgCRND8 17_9 2 62268
wildtype 13_4 5 59935
5_7 4 58685
TgCRND8 5_7 1 58682
wildtype 2_5 3 58231
TgCRND8 2_5 0 53913
Name: count, dtype: int64
1. Cell neighborhood detection via graph construction¶
Spatial transcriptomics data can be represented as graphs with cells as nodes and edges as relations. Depending on the technology (imaging-based or sequencing-based) different assumptions can be made for the graph structure.
Here we will explore graph construction and cell neighborhood analysis using the squidpy sq.gr.spatial_neighbors module on imaging-based data.
Information for graph reconstruciton for sequencing-based data can be found here.
For image-based technologies we construct a graph with coord_type=generic, meaning that the nodes / cells will preserve their spatial location and will not be re-arranged in a spatial grid. For generic graph approaches we can choose to set a fixed radius radius or number of neighbors to connect to (n_neighs).
Below we try both graph construction methods and give them a different key to add to the adata object.
if 'cell_types_colors' in adata.uns:
del adata.uns['cell_types_colors']
sq.gr.spatial_neighbors(adata, n_neighs=5, coord_type="generic", key_added = 'neighs_based_spatial')
sq.pl.spatial_scatter(
adata,
shape=None,
library_key = 'batch_key',
color=["celltype"],
connectivity_key="neighs_based_spatial_connectivities",
title=adata.obs['batch_key'].cat.categories,
ncols=3,
size = 10
)
INFO Creating graph using `None` transform and `1` libraries.
sq.gr.spatial_neighbors(adata, radius=0.2, coord_type="generic", key_added = 'radius_based_spatial')
sq.pl.spatial_scatter(
adata,
shape=None,
library_key = 'batch_key',
color="celltype",
connectivity_key="radius_based_spatial_connectivities",
title=adata.obs['batch_key'].cat.categories,
ncols=3,
size=10,
)
INFO Creating graph using `None` transform and `1` libraries.
First differences across the graphs can be observed by visual inspection, i.e. the cells that are further away from the clear brain structures are connected in the closest neighbor approach but unconnected in the radius-based approach. Some other differences are harder to observe like dependencies in the more dense connected regions. For this Squidpy provides a number of statistics to better the differences between the connectivities of cells.
Using the cell type information we can formulate some hypothesis from the data.
To perform downstream analysis you require a graph structure or radius that connects the cells of interest. For example, if you are interested in the communication between Microglia and Oligodendrocytes you need to set a radius such that the cell typer connections exists.
Task 1: Compute the interaction matrix to understand differences between neighborhood and radius based graphs. Find an appropriate number of neighbors or radius that include sufficient connections between your cell types of interest.
What are the average number of connections per cell type? Which cell types tend to cluster together?
sq.gr.interaction_matrix(adata, cluster_key="celltype", connectivity_key="neighs_based_spatial")
sq.pl.interaction_matrix(adata, cluster_key="celltype", connectivity_key="neighs_based_spatial")
sq.gr.interaction_matrix(adata, cluster_key="celltype", connectivity_key="radius_based_spatial")
sq.pl.interaction_matrix(adata, cluster_key="celltype", connectivity_key="radius_based_spatial")
Task 2: Try out some analysis from Squidpy to identify distinctions between the graphs e.i. using the centrality score and neighborhood enrichment.
sq.gr.centrality_scores(adata, cluster_key = "celltype", connectivity_key = "neighs_based_spatial")
sq.pl.centrality_scores(adata, cluster_key = "celltype")
sq.gr.nhood_enrichment(adata, cluster_key="celltype", connectivity_key = "radius_based_spatial")
0%| | 0/1000 [00:00<?, ?/s]
sq.pl.nhood_enrichment(
adata, cluster_key="celltype", method="average", figsize=(5, 5)
)
2. Spatial Domain detection with CellCharter¶
import logging
import cellcharter as cc
import matplotlib.pyplot as plt
import numpy as np
import scanpy as sc
import scvi
import squidpy as sq
logger = logging.getLogger('pytorch_lightning.utilities.rank_zero')
logger.setLevel(logging.ERROR)
scvi.settings.seed = 12345
scvi.settings.num_threads = 2
Seed set to 12345
There may be cells with very low counts. We will filter them out.
sc.pp.filter_cells(adata, min_counts=3)
2.1 Dimensionality reduction¶
First, we need to run the dimensionality reduction. We will use scVI for this.
Make sure that adata.X contains count data, as it is required by scVI.
In some tutorials you will see that the count data is stored in an AnnData layer called counts. Here we will use the X layer.
scvi.model.SCVI.setup_anndata(adata)
LOAD_MODEL = True
We set the parameters of the neural network. Here we use 1 layer and an embedding size of 10.
These are very common parameters, but you can play with them to see how they affect the results.
model = scvi.model.SCVI(
adata,
n_layers=1,
n_latent=10,
use_layer_norm="both",
use_batch_norm="none",
)
model.train(early_stopping=True, enable_progress_bar=True, max_epochs=10, accelerator="cpu")
GPU available: False, used: False TPU available: False, using: 0 TPU cores 💡 Tip: For seamless cloud logging and experiment tracking, try installing [litlogger](https://pypi.org/project/litlogger/) to enable LitLogger, which logs metrics and artifacts automatically to the Lightning Experiments platform.
Training: 0%| | 0/10 [00:00<?, ?it/s]
`Trainer.fit` stopped: `max_epochs=10` reached.
To make sure that the training has converged, we can plot the training history.
Here we is important to focus on the (validation) reconstruction loss. This shows how well the model is able to reconstruct the original data.
plt.figure(figsize=(5, 5))
plt.plot(
model.history["reconstruction_loss_train"],
label="train",
color="darkgreen",
linewidth=1.25
)
plt.plot(
model.history["reconstruction_loss_validation"],
label="validation",
color="firebrick",
linewidth=1.25
)
plt.legend()
plt.title("reconstruction_loss")
plt.tight_layout()
We extract the latent representation of the data.
adata.obsm['X_scVI'] = model.get_latent_representation(adata).astype(np.float32)
2.2 Neighborhood aggregation¶
Until now, all the analyses have been done ignoring the spatial information.
Now we will use the spatial information to perform the clustering.
First, we need to create a network where cells are connected if they are close to each other using the Delaunay triangulation.
sq.gr.spatial_neighbors(adata, library_key='sample', coord_type='generic')
INFO Creating graph using `None` transform and `6` libraries.
sq.pl.spatial_scatter(
adata,
shape=None,
library_key='sample',
library_id=adata.obs['sample'].cat.categories[0],
color="sample",
size=1,
figsize=(10,10),
connectivity_key="spatial_connectivities",
ncols=1
)
As you can see, the Delaunay triangulation generates very long edges for a few cells.
The most appropriate approach would be to estimate the most biologically relevant distance between cells and use it to remove the edges longer than that distance.
A quick alternative solution is that, since those long edges are sort of outlines, the measure the 99th percentile of the edge lengths and remove the edges longer than that distance.
This process will lead to some isolated cells, but that's not a problem.
sq.gr.spatial_neighbors(adata, library_key='sample', coord_type='generic', percentile=99)
INFO Creating graph using `None` transform and `6` libraries.
sq.pl.spatial_scatter(
adata,
shape=None,
library_key='sample',
library_id=adata.obs['sample'].cat.categories[0],
color="sample",
size=1,
figsize=(10,10),
connectivity_key="spatial_connectivities",
ncols=1
)
We construct the neighborhood aggregated representation by combining every cell's features with the ones of first 3 layers of neighbors.
We will use the scVI latent representation as features and the new features will be stored in adata.obsm['X_cellcharter'].
cc.gr.aggregate_neighbors(adata, n_layers=3, use_rep='X_scVI', out_key='X_cellcharter_temp', sample_key='sample')
0%| | 0/6 [00:00<?, ?it/s]
Task 3: Given that X_scVI contains 10 features, how many features does X_cellcharter contain?
Guess the answer and then write the code to check it.
adata.obsm['X_scVI'].shape
(349217, 10)
#### check 'X_cellcharter_temp'
2.3 Clustering¶
Finally, let's cluster the cells and find their spatial domains.
We will use 18 clusters, we will see later why this is a good choice.
gmm = cc.tl.Cluster(
n_clusters=18,
random_state=12345,
trainer_params={'accelerator': 'cpu'},
)
gmm.fit(adata, use_rep='X_cellcharter_temp')
adata.obs['spatial_domain_temp'] = gmm.predict(adata, use_rep='X_cellcharter_temp')
Output()
GPU available: False, used: False TPU available: False, using: 0 TPU cores 💡 Tip: For seamless cloud logging and experiment tracking, try installing [litlogger](https://pypi.org/project/litlogger/) to enable LitLogger, which logs metrics and artifacts automatically to the Lightning Experiments platform.
Output()
adata.write_h5ad(f"{base_dir}adata_with_spatial_domains.h5ad", compression="gzip")
We can now plot domains and cell types back to back
if 'spatial_domain_temp_colors' in adata.uns:
del adata.uns['spatial_domain_temp_colors']
sq.pl.spatial_scatter(
adata,
shape=None,
library_key='sample',
color=["spatial_domain_temp", "celltype"],
size=1,
figsize=(10,10),
title=np.repeat(adata.obs['sample'].cat.categories, 2),
ncols=2
)
