scverse
Q: How to render xenium morphology image with original colors?
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
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
Also please check my comment here (which answers almost a duplicate of question 2): https://github.com/scverse/spatialdata-io/issues/205
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 :)
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
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.
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')."
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'
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?
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?
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
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)
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 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 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.