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:

I've processed the data so far as:

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)

# 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:

And the queried bounding box:

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:

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.

Individual channels look fine:

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)
    )

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)
    )

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)
    )

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

###
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),
    )

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:

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:

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:

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

首先，通过对查询到的 morphology_focus 图像执行 render_images 调用来复现该问题，并将单个通道与多通道 stack 渲染进行比较。了解 palette、channel 和 norm 参数在多通道路径中的应用方式。当组合的 Xenium 通道能够对归一化作出响应，并使用所请求的颜色和对比度进行渲染时，即表示完成。