Failing to recreate Xenium Exlorer-style multi-stain segmentation panel plot

Hi, I am attempting to use Spatiadata's `pl.render_images` to recreate a Xenium Explorer-style plot for Xenium with the multi-stain segmentation panel. See the example here:
<img width="1500" height="842" alt="Image" src="https://github.com/user-attachments/assets/5f647837-5e59-4050-9bab-7c14de084ccb" />

I've processed the data so far as:
```python
import spatialdata as sd
import spatialdata_plot
from spatialdata_io import xenium
from spatialdata.models import get_channel_names
import scanpy as sc
from pathlib import Path
import shutil
import os
import numpy as np
import pandas as pd
import matplotlib
import matplotlib.patches as patches
from matplotlib import pyplot as plt
path = Path().resolve()

# Set root directory
root_dir = ... # Skipping
# Set figure output directory
figdir = os.path.join(
 root_dir,
 'figures'
)
# Set directory to raw Xenium output bundles
raw_data_dir = os.path.join(
 root_dir,
 'data/raw/xenium'
)
# Set path to reference NicheCompass Anndata file
adata_file = os.path.join(
 root_dir,
 'output/nichecompass',
 'FB_SPI_mouse_filt_cells_NicheCompass_annotated__30012026.h5ad'
)
# Mapping Xenium slide names to Xenium output directories
xenium_output_dir = {
 '7dpi-C': 'output-XETG00155__0043039__7dpi-C__20240725__115033',
 'Sham-A': 'output-XETG00155__0043039__Sham-A__20240725__115033',
 '1dpi-C': 'output-XETG00155__0043039__1dpi-C__20240725__115033',
 '3dpi-A': 'output-XETG00155__0043042__3dpi-A__20240725__115033',
 '7dpi-B': 'output-XETG00335__0027406__7dpi-B__20240725__114550',
 '14dpi-B': 'output-XETG00335__0027281__14dpi-B__20240725__114550',
 '28dpi-A': 'output-XETG00155__0043039__28dpi-A__20240725__115033',
}
slide_xenium_data_dir = {}
for slide in xenium_output_dir.keys():
 slide_xenium_data_dir[slide] = os.path.join(
 raw_data_dir,
 xenium_output_dir[slide]
 )

# Load NicheCompass results
adata_nc = sc.read_h5ad(adata_file)

# Load Xenium data as SpatialData
slide = '7dpi-C'
sdata = xenium(
 path=slide_xenium_data_dir[slide],
 n_jobs=8,
 cell_boundaries=True,
 nucleus_boundaries=False,
 morphology_focus=True,
 cells_as_circles=False,
)

# Subset NicheCompass to slide
adata_nc_subset = adata_nc[adata_nc.obs['regionName']==slide,:].copy()

# Fix 8-digit hex-codes to 6-digit
adata_nc_subset.uns['Niches_colors'] = np.array(
 [c[:-2] for c in adata_nc_subset.uns['Niches_colors']]
)
# Copy over to Spatialdata object
sdata.tables['table'].obs_names = sdata.tables['table'].obs['cell_id']
sdata.tables['table'].obs['Niches'] = adata_nc_subset.obs['Niches']
sdata.tables['table'].obs['cell_type_coarse'] = adata_nc_subset.obs['cell_type_coarse']
sdata.tables['table'].obs['cell_type_granular'] = adata_nc_subset.obs['cell_type_granular']
sdata.tables['table'].uns['Niches_colors'] = adata_nc_subset.uns['Niches_colors']
sdata.tables['table'].uns['cell_type_coarse_colors'] = adata_nc_subset.uns['cell_type_coarse_colors']
sdata.tables['table'].uns['cell_type_granular_colors'] = adata_nc_subset.uns['cell_type_granular_colors']

# Reset index of .obs
sdata.tables['table'].obs.reset_index(drop=True, inplace=True)

# Reset table annotation to 'cell_boundaries'
sdata["table"].obs["region"] = "cell_boundaries"
sdata.set_table_annotates_spatialelement("table", region="cell_boundaries")

# Replace NA values
sdata.tables['table'].obs['Niches'] = (
 sdata.tables['table'].obs['Niches'].cat.add_categories(['Unknown'])
 .fillna('Unknown')
)

# Process colour palette
niches_g = sdata.tables['table'].obs['Niches'].cat.categories.tolist()
niches_p = sdata.tables['table'].uns['Niches_colors'].tolist()
# Add a separate value for cells lacking 'Niche' annotation
niches_p.append('#ffffff')
```

```
INFO reading 
 [/lustre/scratch124/cellgen/bayraktar/kr23/projects/FB_SPI/data/raw/xenium/output-XETG00155__0043039__7dpi-](https://vscode-remote+ssh-002dremote-002bfarm22.vscode-resource.vscode-cdn.net/lustre/scratch124/cellgen/bayraktar/kr23/projects/FB_SPI/data/raw/xenium/output-XETG00155__0043039__7dpi-)
 C__20240725__115033/cell_feature_matrix.h5 
<timed exec>:3: DeprecationWarning: `cell_boundaries` is being deprecated as an argument to `xenium.xenium` in SpatialData version 0.1, switch to `cells_boundaries` instead.
[/software/conda/users/kr23/spatialdata/lib/python3.10/site-packages/spatialdata_io/_utils.py:48](https://vscode-remote+ssh-002dremote-002bfarm22.vscode-resource.vscode-cdn.net/software/conda/users/kr23/spatialdata/lib/python3.10/site-packages/spatialdata_io/_utils.py:48): UserWarning: The "cell_id" column in the cells metadata table does not match the "cell_id" column in the annotation table. This could be due to trying to read a new version that is not supported yet. Please report this issue.
 return f(*args, **kwargs)
[/software/conda/users/kr23/spatialdata/lib/python3.10/site-packages/spatialdata_io/_utils.py:48](https://vscode-remote+ssh-002dremote-002bfarm22.vscode-resource.vscode-cdn.net/software/conda/users/kr23/spatialdata/lib/python3.10/site-packages/spatialdata_io/_utils.py:48): UserWarning: The cell_id column in the cell_labels_table does not match the cell_id column derived from the cell labels data. This could be due to trying to read a new version that is not supported yet. Please report this issue.
 return f(*args, **kwargs)
CPU times: user 23.3 s, sys: 4.31 s, total: 27.6 s
Wall time: 25.1 s
[/software/conda/users/kr23/spatialdata/lib/python3.10/site-packages/spatialdata/_core/spatialdata.py:511](https://vscode-remote+ssh-002dremote-002bfarm22.vscode-resource.vscode-cdn.net/software/conda/users/kr23/spatialdata/lib/python3.10/site-packages/spatialdata/_core/spatialdata.py:511): UserWarning: Converting `region_key: region` to categorical dtype.
 convert_region_column_to_categorical(table)
```


```python
# Bounding box to region of interest
bb_xmin = 11500
bb_ymin = 4000
bb_w = 4000
bb_h = 4000
bb_xmax = bb_xmin + bb_w
bb_ymax = bb_ymin + bb_h

# Highlight region of interest
f, ax = plt.subplots(figsize=(12,24))
sdata.pl.render_shapes(
 'cell_boundaries',
 method = 'matplotlib',
 color = 'Niches',
 groups = niches_g,
 palette = niches_p,
 alpha_fill = 1
 ).pl.show(
 ax = ax,
 title = f'{slide} (zoom)',
 )
rect = patches.Rectangle((bb_xmin, bb_ymin), bb_w, bb_h, linewidth=1, edgecolor="red", facecolor="none")
ax.add_patch(rect)

# Query the bounding box
sdata_xenium_subset = sdata.query.bounding_box(
 axes=["x", "y"],
 min_coordinate=[bb_xmin, bb_ymin],
 max_coordinate=[bb_xmax, bb_ymax],
 target_coordinate_system="global",
)

sdata_xenium_subset.pl.render_shapes(
 'cell_boundaries',
 method = 'matplotlib',
 color = 'Niches',
 groups = niches_g,
 palette = niches_p,
 fill_alpha = 0.8
 ).pl.show(
 title = f'{slide} (zoom) (niches)',
 figsize = (12,12)
 )
```

The full tissue slide:

<img width="1255" height="517" alt="Image" src="https://github.com/user-attachments/assets/eae9eec2-27ac-4f3c-aee6-8c677756cafb" />

And the queried bounding box:

<img width="1211" height="945" alt="Image" src="https://github.com/user-attachments/assets/d1c91497-08c6-4ace-894f-3d6a638b9457" />


Now I am aware of this issue: https://github.com/scverse/spatialdata-plot/issues/370,
but don't even seem to get a similar result:
```python
sdata_xenium_subset.pl.render_images(
 "morphology_focus",
 palette = ["#0F73E6", "#F300A5", "#A4A400", "#008A00"],
 channel = ['DAPI', 'ATP1A1/CD45/E-Cadherin', '18S', 'AlphaSMA/Vimentin']
).pl.show(
 title="Morphology image", 
 coordinate_systems="global",
 colorbar = False,
 figsize = (6,6)
)

[/software/conda/users/kr23/spatialdata/lib/python3.10/site-packages/tifffile/tifffile.py:8797](https://vscode-remote+ssh-002dremote-002bfarm22.vscode-resource.vscode-cdn.net/software/conda/users/kr23/spatialdata/lib/python3.10/site-packages/tifffile/tifffile.py:8797): UserWarning: <tifffile.TiffPage 0 @16> reading array from closed file
 warnings.warn(
INFO Your image has 4 channels. Sampling categorical colors and using multichannel strategy 'stack' to render.
```

<img width="611" height="590" alt="Image" src="https://github.com/user-attachments/assets/3f2db5b9-0b28-487f-9c13-18440c84efe3" /> 



Individual channels look fine:
```python
sdata_xenium_subset.pl.render_images(
 'morphology_focus',
 palette = ['#0f73e6'],
 channel = ['DAPI'],
 alpha = 1
 # norm = matplotlib.colors.Normalize(0, 5000, clip=True)
 ).pl.show(
 title = f'{slide} (zoom) (DAPI)',
 colorbar = False,
 figsize = (6,6)
 )
```

<img width="611" height="590" alt="Image" src="https://github.com/user-attachments/assets/a024e528-2916-48d9-ae0a-19c21ee37e86" />

```python
sdata_xenium_subset.pl.render_images(
 'morphology_focus',
 palette = ['#a4a400'],
 channel = ['18S'],
 # norm = matplotlib.colors.Normalize(0, 5000, clip=True)
 ).pl.show(
 title = f'{slide} (zoom) (18S)',
 colorbar = False,
 figsize = (6,6)
 )
```

<img width="611" height="590" alt="Image" src="https://github.com/user-attachments/assets/4d73c264-c2ee-4606-8d7e-5783c71da7f2" />


```python
sdata_xenium_subset.pl.render_images(
 'morphology_focus',
 palette = ['#f300a5'],
 channel = ['ATP1A1/CD45/E-Cadherin'],
 # norm = matplotlib.colors.Normalize(0, 8000, clip=True)
 ).pl.show(
 title = f'{slide} (zoom) (ATP1A1/CD45/E-Cadherin)',
 colorbar = False,
 figsize = (6,6)
 )
```

<img width="611" height="590" alt="Image" src="https://github.com/user-attachments/assets/7af699c7-69a3-4ad0-8de9-c42526eb5ae1" />


However none of the options to do normalisation clipping seem to affect any of the combined plots, neither when I do this for instance:

```python
###
sdata_xenium_subset.pl.render_images(
 'morphology_focus',
 palette = [
 '#0f73e6',
 '#a4a400',
 '#f300a5'
 ],
 channel = [
 'DAPI',
 '18S',
 'ATP1A1/CD45/E-Cadherin'
 ],
 # norm = matplotlib.colors.Normalize(0, 8000, clip=True)
 ).pl.show(
 title = f'{slide} (zoom) (DAPI - 18S - E-Cadherin)',
 colorbar = False,
 figsize = (6,6),
 )
```

<img width="611" height="590" alt="Image" src="https://github.com/user-attachments/assets/3f1ba35b-75d3-47a0-8993-f02849793822" />

Is exactly the same look regardless of setting e.g. `norm = matplotlib.colors.Normalize(0, 10000, clip=True)`, `norm = matplotlib.colors.Normalize(0, 1.5, clip=True)`, or `norm = matplotlib.colors.Normalize(0, 2000, clip=True)`.


This is what is found in the queried object:
```python
sdata_xenium_subset

SpatialData object
├── Images
│ └── 'morphology_focus': DataTree[cyx] (4, 4000, 4000), (4, 2000, 2000), (4, 1000, 1000), (4, 500, 500), (4, 250, 250)
├── Labels
│ ├── 'cell_labels': DataTree[yx] (4000, 4000), (2000, 2000), (1000, 1000), (500, 500), (250, 250)
│ └── 'nucleus_labels': DataTree[yx] (4000, 4000), (2000, 2000), (1000, 1000), (500, 500), (250, 250)
├── Points
│ └── 'transcripts': DataFrame with shape: (<Delayed>, 13) (3D points)
├── Shapes
│ └── 'cell_boundaries': GeoDataFrame shape: (4500, 1) (2D shapes)
└── Tables
 └── 'table': AnnData (4500, 474)
with coordinate systems:
 ▸ 'global', with elements:
 morphology_focus (Images), cell_labels (Labels), nucleus_labels (Labels), transcripts (Points), cell_boundaries (Shapes)
```

and `images`:
```python
sdata_xenium_subset.images
{'morphology_focus': <xarray.DataTree>
Group: [/](https://vscode-remote+ssh-002dremote-002bfarm22.vscode-resource.vscode-cdn.net/)
├── Group: [/scale0](https://vscode-remote+ssh-002dremote-002bfarm22.vscode-resource.vscode-cdn.net/scale0)
│ Dimensions: (c: 4, y: 4000, x: 4000)
│ Coordinates:
│ * c (c) <U22 352B 'DAPI' ... 'AlphaSMA/Vimentin'
│ * y (y) float64 32kB 0.5 1.5 2.5 3.5 ... 3.998e+03 3.998e+03 4e+03
│ * x (x) float64 32kB 0.5 1.5 2.5 3.5 ... 3.998e+03 3.998e+03 4e+03
│ Data variables:
│ image (c, y, x) uint16 128MB dask.array<chunksize=(2, 4000, 4000), meta=np.ndarray>
├── Group: [/scale1](https://vscode-remote+ssh-002dremote-002bfarm22.vscode-resource.vscode-cdn.net/scale1)
│ Dimensions: (c: 4, y: 2000, x: 2000)
│ Coordinates:
│ * c (c) <U22 352B 'DAPI' ... 'AlphaSMA/Vimentin'
│ * y (y) float64 16kB 1.0 3.0 5.0 7.0 ... 3.995e+03 3.997e+03 3.999e+03
│ * x (x) float64 16kB 1.0 3.0 5.0 7.0 ... 3.995e+03 3.997e+03 3.999e+03
│ Data variables:
│ image (c, y, x) uint16 32MB dask.array<chunksize=(2, 2000, 2000), meta=np.ndarray>
├── Group: [/scale2](https://vscode-remote+ssh-002dremote-002bfarm22.vscode-resource.vscode-cdn.net/scale2)
│ Dimensions: (c: 4, y: 1000, x: 1000)
│ Coordinates:
│ * c (c) <U22 352B 'DAPI' ... 'AlphaSMA/Vimentin'
│ * y (y) float64 8kB 2.0 6.0 10.0 14.0 ... 3.99e+03 3.994e+03 3.998e+03
│ * x (x) float64 8kB 2.0 6.0 10.0 14.0 ... 3.99e+03 3.994e+03 3.998e+03
│ Data variables:
...
 * y (y) float64 2kB 8.0 24.0 40.0 56.0 ... 3.96e+03 3.976e+03 3.992e+03
 * x (x) float64 2kB 8.0 24.0 40.0 56.0 ... 3.96e+03 3.976e+03 3.992e+03
 Data variables:
 image (c, y, x) uint16 500kB dask.array<chunksize=(4, 250, 250), meta=np.ndarray>}
Output is truncated. View as a [scrollable element](command:cellOutput.enableScrolling?a43a3130-bed1-4d79-94bd-0ac7729b722a) or open in a [text editor](command:workbench.action.openLargeOutput?a43a3130-bed1-4d79-94bd-0ac7729b722a). Adjust cell output [settings](command:workbench.action.openSettings?%5B%22%40tag%3AnotebookOutputLayout%22%5D)...
```

Are you familiar with any ways to improve this plotting? I would love to have the custom control that Spatialdata offers to generate high-quality figures for publications, but so far the different scales of stains doesn't seem to work well.

Finally, I have also tried this issue: https://github.com/scverse/spatialdata-plot/issues/460, but got exactly the same issue in the end as the author mentioned in their last comment.

All my code was run with:
```python
print(sd.__version__)
print(spatialdata_io.__version__)
print(spatialdata_plot.__version__)
0.5.0
0.3.0
0.2.12
```

I hope you can help me out!

Best,
Koen

Navigation Menu

Uh oh!

Failing to recreate Xenium Exlorer-style multi-stain segmentation panel plot #534

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Uh oh!

Failing to recreate Xenium Exlorer-style multi-stain segmentation panel plot #534

Description

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions