Incompatible coordinate transformation with cellpose segmentation result

[mamyAndrianteranagna]

Hello Quentin,

I have this problem that I posted in spatialdata github as an issue but we are wondering if it is not related to sopa:

Here is the issue:
scverse/spatialdata-io#166

Let me know if you need more details.

Regards,
Mamy

[quentinblampey]

Hello @mamyAndrianteranagna, thanks for referring to this issue.

Can you list which versions of sopa, spatialdata, and spatialdata-io do you use? And can you confirm that you read the data with the MERSCOPE reader (i.e. sopa.io.merscope)?

It's weird that the "c" dimension being used in the transformations, because the MERSCOPE reader is only using the "x"/"y" dimensions in the transformations. Then, when using cellpose, we use the same transformation as the image, which shouldn't include "c" either.

Last question: I saw you have an element called 'cellpose_nuclei_6diam_0minArea', which doesn't come from sopa, so what is this object exactly? And how did you compute it?

[mamyAndrianteranagna]

Hello Quentin,

I actually used spatialdata.io to open the MERSCOPE output data before creating the spatial data object:

import spatialdata_io as sdio
sdata = sdio.merscope("/path/to/MERSCOPE/output/region")
sdata.write("/path/to/zarr")

Sorry about that, I didn't notice that there is a sopa module for it. However it works well with baysor.
Therefore, I will change this first then let you know the result.

Many thanks!

Best,
Mamy

[quentinblampey]

Actually, the reader should be the same (sopa.io.merscope is a simple wrapper around spatialdata_io.merscope).

But, still, it's good to re-run it with the latest versions, since it appears that your versions are not up to date!

[mamyAndrianteranagna]

Hi Quentin,

Actually I have the same probleme even with sopa.io.merscope.
I will update everything then re-run it again.
I will let you know.

Best,
Mamy

[mamyAndrianteranagna]

Sorry I forget, here are my current versions of
sopa: v1.1.0
spatialdata: v0.1.2
spatialdata-io: v0.1.2

Best,
Mamy

[mamyAndrianteranagna]

Hello Quentin,

I have now these versions:
sopa: v1.1.2
spatialdata: v0.2.1
spatialdata-io: v0.1.3.post0

but, unfortunately, I still get the same issue.

Here is my code:

# Import raw data and create spatialdata

import sopa.io

regionPath = '/path/to/merscope/output/region'

sdata = sopa.io.merscope(path=regionPath, backend=None, z_layers=3, region_name="region_0")


# Create patches

import spatialdata as sd
import sopa.segmentation

points_key = 'output_region_0_transcripts'

patches = sopa.segmentation.Patches2D(sdata=sdata, element_name=points_key, patch_width=4000, patch_overlap=30)
patches.write()


# Run cellpose nuclei segmentation

method = sopa.segmentation.methods.cellpose_patch(diameter=50, channels="DAPI", flow_threshold=2, cellprob_threshold=-6)

for patch_index in range(len(sdata['sopa_patches'])):
    segmentation.write_patch_cells(tmpDir, patch_index)


# Resolve conflicts

nuclei = sopa.segmentation.StainingSegmentation.read_patches_cells(tmpDir)
nuclei = sopa.segmentation.shapes.solve_conflicts(nuclei)

shapes_key = 'cellpose_DAPI_diam_50_minarea_0_flow_2_cellprob_minus6'
image_key = 'output_region_0_z3'

sopa.segmentation.StainingSegmentation.add_shapes(sdta, nuclei, image_key, shapes_key)


# Aggregate

aggregator = sopa.segmentation.Aggregator(sdata=sdata, image_key=image_key, shapes_key=shapes_key)

aggregator.compute_table(gene_column='gene', average_intensities=False)


# Check spatialdata object

sdata

SpatialData object with:
├── Images
│     └── 'output_region_0_z3': MultiscaleSpatialImage[cyx] (11, 96363, 78210), (11, 48181, 39105), (11, 2224090, 19552), (11, 12045, 9776), (11, 6022, 4888)
├── Points
│     └── 'output_region_0_transcripts': DataFrame with shape: (19797064, 9) (2D points)
├── Shapes
│     ├── 'cellpose_DAPI_diam_50_minarea_0_flow_2_cellprob_minus6': GeoDataFrame shape: (46448, 1) (2D shapes)
│     └── 'sopa_patches': GeoDataFrame shape: (9, 1) (2D shapes)
└── Tables
      └── 'table': AnnData (46448, 315)
with coordinate systems:
▸ 'microns', with elements:
        output_region_0_z3 (Images), output_region_0_transcripts (Points), cellpose_DAPI_diam_50_minarea_0_flow_2_cellprob_minus6 (Shapes), sopa_patches (Shapes)

# Check coordinate system transformations for each layers

print(sdata.images['output_region_0_z3'].attrs)
print(sdata.points['output_region_0_transcripts'].attrs)
print(sdata.shapes['cellpose_DAPI_diam_50_minarea_0_flow_2_cellprob_minus6'].attrs)
print(sdata.shapes['sopa_patches'].attrs)

{}
{'transform': {'microns': Identity }, 'spatialdata_attrs': {'feature_key': 'gene'}}
{'transform': {'microns': Affine (x, y -> c, x, y)
    [0. 0. 0.]
    [1.07998399e-01 0.00000000e+00 4.20769397e+03]
    [0.00000000e+00 1.07999088e-01 1.62009116e+03]
    [0. 0. 1.]}}
{'transform': {'microns': Identity }}

Let me know if you need more informations.

Regards,
Mamy

[quentinblampey]

Hello @mamyAndrianteranagna,

Despite the weird "c" channel in the transformation, there is something weird happening before.
Starting from spatialdata-io==0.1.3, the main coordinate system for MERSCOPE data is "global" (i.e. the pixel coordinate system of the image), not microns. See the release notes here. And, as sopa.io.merscope is a simple wrapper around it, you should also see the "global" coordinate system.

For instance, using the same versions of sopa/spatialdata/spatialdata_io as you, this is what I get when loading some MERSCOPE data:

>>> sopa.io.merscope("/path/to/region_0")

SpatialData object
├── Images
│     └── 'example_region_0_z3': DataTree[cyx] (2, 45030, 28250), (2, 22515, 14125), (2, 11257, 7062), (2, 5628, 3531), (2, 2814, 1765)
└── Points
      └── 'example_region_0_transcripts': DataFrame with shape: (<Delayed>, 9) (2D points)
with coordinate systems:
    ▸ 'global', with elements:
        example_region_0_z3 (Images), example_region_0_transcripts (Points)

In your case, you have 'microns' instead of 'global'.

Can you confirm that you don't see the 'global' coordinate system when running sopa.io.merscope(path=regionPath, backend=None, z_layers=3, region_name="region_0")?

[mamyAndrianteranagna]

Hello Quentin,

I agree, I suspect also something before, because I never get the 'global' coordinate system like seen everywhere in the exemple data from MERSCOPE. But I assumed that this can be aligned, transformed or corrected by the tool. With baysor results, there is no problem.

So I confirm that I don't see coordinate named 'global' from the beginning.

Thanks for your help!

Best,
Mamy

[quentinblampey]

Well, everything should work whatever the main coordinate system, for both Cellpose and Baysor.

My guess is that there is nothing wrong with Sopa or SpatialData, but instead, there is an issue with your environment.
Indeed, since you don't see the "global" coordinate system, you might have something broken. I know you already gave me the versions of the packages, but can you try this:

import spatialdata
import spatialdata_io
import sopa

print(spatialdata.__version__, spatialdata_io.__version__, sopa.__version__)

Somehow your environment seems to use old versions of these packages, and depending on how old they were, maybe you are experiencing a past bug.

[mamyAndrianteranagna]

Hello Quentin,

import spatialdata
import spatialdata_io
import sopa

print(spatialdata.__version__, spatialdata_io.__version__, sopa.__version__)

0.2.1 0.1.3.post0 1.1.2

Maybe, a solution is to re-process the data using the new Vizgen Post-processing Tool? The data was generated last year and we do not have much informations about how it was pre-processed. But, in my understanding, they choose to use microns instead of pixels for the initial segmentation, etc.
What do you think?

Best,
Mamy

[quentinblampey]

I think the MERSCOPE data format has been the same since approximately early 2023, so everything should be good. Processing the data with VPT will not help (we read the raw images/transcripts, so this shouldn't affect your SpatialData object anyway).

Indeed, the coordinate systems are created when reading the data with spatialdata_io, this is why I'm surprised about the "microns" coordinate system: whatever your data, the code here clearly mentions only the "global" coordinate system, and the string "microns" is not present in the source code.

Sorry for not being very helpful here, I don't understand what is happening. Your version indicates a recent version of spatialdata_io, while your experiment indicates old versions of spatialdata_io 🤔

[mamyAndrianteranagna]

Hi Quentin,

Thanks for these informations!
So, I will re-build a new environment using the latest versions of all tools.
I will let you know.
Thanks for your time!

Best,
Mamy

[mamyAndrianteranagna]

Hi Quentin,

Sorry to come back again with this issue.

Just to inform you that I created a new environment and now we have the 'global' coordinate (not micron) as before but, unfortunately, we still have the image and the cellpose result not aligned and look weird.

Let me know if you need details.

Best,
Mamy

[quentinblampey]

Hi @mamyAndrianteranagna, it's good to see that now you can see the "global" coordinate system, but sorry to hear that the alignment is still wrong.
Do you still have the weird transformation Affine (x, y -> c, x, y) (with the "c" channel)?