Q: How to render xenium morphology image with original colors?

[tsvvas]

Hello @LucaMarconato!

Thank you very much for this tool. I want to plot a part of morphology image acquired with xenium + multimodal segmentation kit. It is loaded into spatialdata and plotted without errors, but the colors are dim and they are different from what I'd expected (DAPI - #0F73E6, Boundary - #F300A5, Interior RNA - #A4A400, interior protein - #008A00, reference values from Xenium Explorer). I tried passing a list of colors like this:

crop = lambda sdata: spatialdata.bounding_box_query(
    sdata,
    min_coordinate=[17_500, 55_000],
    max_coordinate=[19_500, 57_000],
    axes=("x", "y"),
    target_coordinate_system="global",
)

crop(sdata).pl.render_images("morphology_focus", ["#0F73E6", "#F300A5", "#A4A400", "#008A00"]).pl.show(
    ax=axes[0], title="Morphology image", coordinate_systems="global"
)

The code does not produce an error, it works twice as long, and it results in the same image with colors unchanged.

The key question is how to adjust the colors and saturation of the plot, so that it looked closer to what I see in Xenium Explorer?

Secondary question is how to translate physical coordinates into spatialdata's global coordinate system? Am I supposed to get transformation and then apply transformation to a single or a couple points to get min and max for cropping above? I can make a separate issue with this question, if necessary.

Best,
Vasily

[LucaMarconato]

Hi, for coloring the channels please use the palette argument.

Regarding coordinates yes, once you identify the transformation from the physical coordinate system to the global one, you would need to transform the coordinates and use these in the bounding box query.

We are going to improve the ergonomics of this as described in these issues (except for 338, which is unrelated): https://github.com/scverse/spatialdata/issues?q=sort%3Aupdated-desc%20is%3Aissue%20is%3Aopen%20physical

[LucaMarconato]

Also please check my comment here (which answers almost a duplicate of question 2): scverse/spatialdata-io#205

[timtreis]

Hey @tsvvas,

I assume you have an image with 4 channels that you're trying to map the individual colors to, right? As Luca said, you can use the palette argument for this, see here:

Please let us know if that solves your issue :)

[tsvvas]

Hi @timtreis,

Thank you for checkin in on the issue. I managed to produce an image, but not without issues.

First, the morphology_focus image after reading with spatiadata_io.xenium() has 5 channels for some reason, and I don't know what is the 5th channel. To my best knowledge there should be only 4 images.

Second, there is a small issue where providing only 4 colors causes an error, even if I select only 4 channels. So to produce the image I had to pass five colors like this:

cropped.pl.render_images(
    "morphology_focus",
    palette=["#0F73E6", "#F300A5", "#A4A400", "#008A00", "#FFFFFF"],
    channels=[0, 1, 2, 3],
).pl.show(
    title="Morphology image", 
    coordinate_systems="global",
)

At first I thought that the rendered image is too dim, but it appears that one of the layers is not rendered properly this way. Here is the image I get after executing the code above:

Here is the image of the same region from Xenium Explorer:

And here is the image of the same region from Xenium Explorer without Interior RNA:

So, I am still looking for the way to render the image the same way it is shown in Xenium Explorer.

Best,
Vasily

[timtreis]

The argument is actually called channel without the s, could you try that? It requires 5 colors since you're not limiting the channels and therefore it uses all 5 you have in the image.

[tsvvas]

With channel I get the following error:

cropped.pl.render_images(
    "morphology_focus",
    palette=["#0F73E6", "#F300A5", "#A4A400", "#008A00"],
    channel=[0, 1, 2, 3],
).pl.show(
    title="Morphology image", 
    coordinate_systems="global"
)

Error:

File [spatialdata_plot/pl/basic.py:862]
    856     if wanted_images_on_this_cs:
    857         rasterize = (params_copy.scale is None) or (
    858             isinstance(params_copy.scale, str)
    859             and params_copy.scale != "full"
    860             and (dpi is not None or figsize is not None)
    861         )
--> 862         _render_images(
    863             sdata=sdata,
    864             render_params=params_copy,
    865             coordinate_system=cs,
    866             ax=ax,
    867             fig_params=fig_params,
    868             scalebar_params=scalebar_params,
    869             legend_params=legend_params,
    870             rasterize=rasterize,
    871         )
File [spatialdata_plot/pl/render.py:634]
    632 layers = {}
    633 for ch_index, c in enumerate(channels):
--> 634     layers[c] = img.sel(c=c).copy(deep=True).squeeze()
    636     if not isinstance(render_params.cmap_params, list):
    637         if render_params.cmap_params.norm is not None:
File [xarray/core/dataarray.py:1670]
   1554 def sel(
   1555     self,
   1556     indexers: Mapping[Any, Any] | None = None,
   (...)
   1560     **indexers_kwargs: Any,
   1561 ) -> Self:
   1562     """Return a new DataArray whose data is given by selecting index
   1563     labels along the specified dimension(s).
   1564 
   (...)
   1668     Dimensions without coordinates: points
   1669     """
-> 1670     ds = self._to_temp_dataset().sel(
   1671         indexers=indexers,
   1672         drop=drop,
   1673         method=method,
   1674         tolerance=tolerance,
   1675         **indexers_kwargs,
   1676     )
   1677     return self._from_temp_dataset(ds)
File [xarray/core/dataset.py:3184]
   3116 """Returns a new dataset with each array indexed by tick labels
   3117 along the specified dimension(s).
   3118 
   (...)
   3181 
   3182 """
   3183 indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "sel")
-> 3184 query_results = map_index_queries(
   3185     self, indexers=indexers, method=method, tolerance=tolerance
   3186 )
   3188 if drop:
   3189     no_scalar_variables = {}
File [xarray/core/indexing.py:193]
    191         results.append(IndexSelResult(labels))
    192     else:
--> 193         results.append(index.sel(labels, **options))
    195 merged = merge_sel_results(results)
    197 # drop dimension coordinates found in dimension indexers
    198 # (also drop multi-index if any)
    199 # (.sel() already ensures alignment)
File [xarray/core/indexes.py:791]
    789                 indexer = self.index.get_loc(label_value)
    790             except KeyError as e:
--> 791                 raise KeyError(
    792                     f"not all values found in index {coord_name!r}. "
    793                     "Try setting the `method` keyword argument (example: method='nearest')."
    794                 ) from e
    796 elif label_array.dtype.kind == "b":
    797     indexer = label_array

KeyError: "not all values found in index 'c'. Try setting the `method` keyword argument (example: method='nearest')."

[tsvvas]

A fully reproducible example:

! curl -O https://cf.10xgenomics.com/samples/xenium/2.0.0/Xenium_V1_humanLung_Cancer_FFPE/Xenium_V1_humanLung_Cancer_FFPE_outs.zip --output-dir ../data
! mkdir ../data/Xenium_V1_humanLung_Cancer_FFPE_outs
! unzip ../data/Xenium_V1_humanLung_Cancer_FFPE_outs.zip -d ../data/Xenium_V1_humanLung_Cancer_FFPE_outs

DATA_PATH = "../data/Xenium_V1_humanLung_Cancer_FFPE_outs"
sdata = spio.xenium(DATA_PATH, n_jobs=-1)

OUT_PATH = f"../output/{date.today()}"
FILE_PATH = os.path.join(OUT_PATH, "data.zarr")
os.mkdir(OUT_PATH)
sdata.write(FILE_PATH)

sdata = spatialdata.read_zarr(FILE_PATH)

sdata.pl.render_images(
    "morphology_focus",
    palette=["#0F73E6", "#F300A5", "#A4A400", "#008A00"],
    channel=[0, 1, 2, 3],
).pl.show(
    title="Morphology image", 
    coordinate_systems="global"
)

Interestingly, if I provide five colors, the error becomes:
AttributeError: 'DataTree' object has no attribute 'c'

>>> spatialdata.__version__
'0.2.3'

[timtreis]

Hm, I think the true underlying issue here is that the image is stored as a datatree (multi-scale image) and we're not respecting that when sub-selecting for channels. Will check your MRE later, thank you!

If you want, you could try to overwrite the multi-scale image with one of the scales of the datatree. It probably looks like the second one here, right?

[tsvvas]

Yes, this is how it looks like in the example above

What is the 5th channel here? Or should I ask this in the spatialdata_io issues?

[timtreis]

I'll probably create a new issue for the issue you probably identified, that we cannot pick channels out of datatrees robustly. So we might close this afterwards.

Could you plot the individual channels of your image separately before opening an issue in spatialdata-io? I'd be surprised if we create a non-sense channel during reading. The separate plot of the 5 channels might help troubleshoot here

[tsvvas]

Hi @timtreis,

Seems like the 5th channel is empty:

>>> sdata.images["morphology_focus"]["scale0"]["image"].isel(c=4).sum(dim=["y", "x"]).values
array(0, dtype=uint64)

[tsvvas]

I tried the overwriting with single-scale image approach and still got the same error (continuing the example with the data above):

img = sdata.images["morphology_focus"].copy()
img = img.drop_nodes(["scale1", "scale2", "scale3", "scale4"])
sdata.images["morpology_single_scale"] = img
sdata.pl.render_images(
    "morpology_single_scale",
    palette=["#0F73E6", "#F300A5", "#A4A400", "#008A00"],
    channels=[0, 1, 2, 3],
).pl.show(
    title="Morphology image", 
    coordinate_systems="global"
)

The error is still ValueError: If 'palette' is provided, its length must match the number of channels.

Providing 5 colors in the palette still results in the absence of Interior RNA channel in the rendered image.

[timtreis]

@LucaMarconato have we seen a 5th channel popping up before?

[tsvvas]

@timtreis The data from spatialdata tutorial on xenium analysis has 5 channels in "morphology_focus" image. But for me the key problem now is that I cannot render properly the other 4 channels.

[LucaMarconato]

@LucaMarconato have we seen a 5th channel popping up before?

Yes, adding a dummy fifth channel was a workaround when the support for the rgb parameter in spatial-image (and napari-spatialdata) was not available. We can remove that from the xenium() reader now, I'll track this in a separate issue.