spatialdata-io icon indicating copy to clipboard operation
spatialdata-io copied to clipboard

cosmx - ValueError: Affine matrix must be homogeneous.

Open jamesboot opened this issue 1 year ago • 7 comments

Hi,

I'm trying to read in some CosMx data - I have a dataset formatted into a folder with the CellLabels and CellComposite directories and then the relevant CSV files. When I execute the code below:

adata3 = sp.cosmx(sample_dir,
                  dataset_id = 'GCIN10692_TMA2')

I get the following error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[21], line 1
----> 1 adata3 = sp.cosmx(sample_dir,
      2                   dataset_id = 'GCIN10692_TMA2')

File [/path/to/in839/analysis/init-rocker-renv/venv/lib/python3.12/site-packages/spatialdata_io/readers/cosmx.py:148](http://localhost:8888/lab/tree/venv/lib/python3.12/site-packages/spatialdata_io/readers/cosmx.py#line=147), in cosmx(path, dataset_id, transcripts, imread_kwargs, image_models_kwargs)
    146     glob = table[idx, :].obs[[CosmxKeys.X_GLOBAL_CELL, CosmxKeys.Y_GLOBAL_CELL]].values
    147     out = estimate_transform(ttype="affine", src=loc, dst=glob)
--> 148     affine_transforms_to_global[fov] = Affine(
    149         # out.params, input_coordinate_system=input_cs, output_coordinate_system=output_cs
    150         out.params,
    151         input_axes=("x", "y"),
    152         output_axes=("x", "y"),
    153     )
    155 table.obsm["global"] = table.obs[[CosmxKeys.X_GLOBAL_CELL, CosmxKeys.Y_GLOBAL_CELL]].to_numpy()
    156 table.obsm["spatial"] = table.obs[[CosmxKeys.X_LOCAL_CELL, CosmxKeys.Y_LOCAL_CELL]].to_numpy()

File [/path/to/in839/analysis/init-rocker-renv/venv/lib/python3.12/site-packages/spatialdata/transformations/transformations.py:523](http://localhost:8888/lab/tree/venv/lib/python3.12/site-packages/spatialdata/transformations/transformations.py#line=522), in Affine.__init__(self, matrix, input_axes, output_axes)
    521     raise ValueError("Invalid shape for affine matrix.")
    522 if not np.array_equal(self.matrix[-1, :-1], np.zeros(len(input_axes))):
--> 523     raise ValueError("Affine matrix must be homogeneous.")
    524 assert self.matrix[-1, -1] == 1.0

ValueError: Affine matrix must be homogeneous.

Any recommendations on how to address this?

Cheers

jamesboot avatar Oct 02 '24 09:10 jamesboot

Hi, can you set a breakpoint/add a print and copy-paste the matrix here please?

LucaMarconato avatar Oct 02 '24 12:10 LucaMarconato

The expression matrix?

jamesboot avatar Oct 03 '24 15:10 jamesboot

No I mean the affine matrix, since the error is coming from there. In this case the out.params variable. Thanks.

LucaMarconato avatar Oct 03 '24 15:10 LucaMarconato

Sorry for the slow response. Took me a while to try and debug and get a sensible output that might make sense for you.

In the end I modified the spatialdata_io/readers/cosmx.py script to include some print statements in the for loop to print the out variable and the out.params variable around line 148.

Output:

printing out
<AffineTransform(matrix=
    [[ 8.10368448e-01, -5.52281315e-01,  1.23965615e+04],
     [ 3.19759226e-01, -6.79466884e-02,  1.53052481e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 8.10368448e-01 -5.52281315e-01  1.23965615e+04]
 [ 3.19759226e-01 -6.79466884e-02  1.53052481e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00,  1.52916069e-09,  1.45160000e+04],
     [ 6.02308243e-06, -9.99996320e-01,  1.56674966e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00  1.52916069e-09  1.45160000e+04]
 [ 6.02308243e-06 -9.99996320e-01  1.56674966e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00, -4.03841338e-11,  1.87750000e+04],
     [ 1.40818796e-06, -9.99999676e-01,  1.56674993e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00 -4.03841338e-11  1.87750000e+04]
 [ 1.40818796e-06 -9.99999676e-01  1.56674993e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00, -1.25917955e-09,  2.30410000e+04],
     [ 2.84186251e-07, -9.99987522e-01,  1.56674951e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00 -1.25917955e-09  2.30410000e+04]
 [ 2.84186251e-07 -9.99987522e-01  1.56674951e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00,  1.92296269e-16,  1.02500000e+04],
     [-2.88444403e-16, -1.00000000e+00,  1.52408000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00  1.92296269e-16  1.02500000e+04]
 [-2.88444403e-16 -1.00000000e+00  1.52408000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00, -1.92296269e-16,  1.45160000e+04],
     [ 4.50694380e-18, -1.00000000e+00,  1.52408000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00 -1.92296269e-16  1.45160000e+04]
 [ 4.50694380e-18 -1.00000000e+00  1.52408000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00, -1.92296269e-16,  1.87750000e+04],
     [-1.08166651e-16, -1.00000000e+00,  1.52408000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00 -1.92296269e-16  1.87750000e+04]
 [-1.08166651e-16 -1.00000000e+00  1.52408000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.0000000e+00,  0.0000000e+00,  2.3041000e+04],
     [ 1.1417591e-16, -1.0000000e+00,  1.5240800e+05],
     [ 0.0000000e+00,  0.0000000e+00,  1.0000000e+00]])>
printing out.params
[[ 1.0000000e+00  0.0000000e+00  2.3041000e+04]
 [ 1.1417591e-16 -1.0000000e+00  1.5240800e+05]
 [ 0.0000000e+00  0.0000000e+00  1.0000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00,  0.00000000e+00,  2.73080000e+04],
     [ 5.76888806e-16, -1.00000000e+00,  1.52408000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00  0.00000000e+00  2.73080000e+04]
 [ 5.76888806e-16 -1.00000000e+00  1.52408000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00,  0.00000000e+00,  1.02500000e+04],
     [ 3.61306661e-16, -1.00000000e+00,  1.48150000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00  0.00000000e+00  1.02500000e+04]
 [ 3.61306661e-16 -1.00000000e+00  1.48150000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00,  1.92296269e-16,  1.45160000e+04],
     [ 2.40370336e-17, -1.00000000e+00,  1.48150000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00  1.92296269e-16  1.45160000e+04]
 [ 2.40370336e-17 -1.00000000e+00  1.48150000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00,  1.92296269e-16,  1.87750000e+04],
     [ 4.20648088e-17, -1.00000000e+00,  1.48150000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00  1.92296269e-16  1.87750000e+04]
 [ 4.20648088e-17 -1.00000000e+00  1.48150000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00,  1.92296269e-16,  2.30410000e+04],
     [ 2.13328673e-16, -1.00000000e+00,  1.48150000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00  1.92296269e-16  2.30410000e+04]
 [ 2.13328673e-16 -1.00000000e+00  1.48150000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00,  5.76888806e-16,  2.73080000e+04],
     [-1.41217572e-16, -1.00000000e+00,  1.48150000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00  5.76888806e-16  2.73080000e+04]
 [-1.41217572e-16 -1.00000000e+00  1.48150000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00, -2.88444403e-16,  1.02500000e+04],
     [ 3.36518470e-16, -1.00000000e+00,  1.43883000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00 -2.88444403e-16  1.02500000e+04]
 [ 3.36518470e-16 -1.00000000e+00  1.43883000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00,  0.00000000e+00,  1.45160000e+04],
     [-4.20648088e-17, -1.00000000e+00,  1.43883000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00  0.00000000e+00  1.45160000e+04]
 [-4.20648088e-17 -1.00000000e+00  1.43883000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00,  1.92296269e-16,  1.87750000e+04],
     [-1.44222201e-16, -1.00000000e+00,  1.43883000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00  1.92296269e-16  1.87750000e+04]
 [-1.44222201e-16 -1.00000000e+00  1.43883000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00, -1.92296269e-16,  2.30410000e+04],
     [-1.20185168e-17, -1.00000000e+00,  1.43883000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00 -1.92296269e-16  2.30410000e+04]
 [-1.20185168e-17 -1.00000000e+00  1.43883000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00,  3.84592537e-16,  2.73080000e+04],
     [ 2.64407369e-16, -1.00000000e+00,  1.43883000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00  3.84592537e-16  2.73080000e+04]
 [ 2.64407369e-16 -1.00000000e+00  1.43883000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00,  0.00000000e+00,  1.45160000e+04],
     [ 9.61481343e-17, -1.00000000e+00,  1.39616000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00  0.00000000e+00  1.45160000e+04]
 [ 9.61481343e-17 -1.00000000e+00  1.39616000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00,  0.00000000e+00,  1.87750000e+04],
     [ 2.93702504e-16, -1.00000000e+00,  1.39616000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00  0.00000000e+00  1.87750000e+04]
 [ 2.93702504e-16 -1.00000000e+00  1.39616000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00, -5.76888806e-16,  2.30410000e+04],
     [-9.61481343e-17, -1.00000000e+00,  1.39616000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00 -5.76888806e-16  2.30410000e+04]
 [-9.61481343e-17 -1.00000000e+00  1.39616000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00,  2.30755522e-15,  2.73080000e+04],
     [ 2.74022183e-15, -1.00000000e+00,  1.39616000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00  2.30755522e-15  2.73080000e+04]
 [ 2.74022183e-15 -1.00000000e+00  1.39616000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00, -2.94213291e-14,  3.56500000e+04],
     [-3.65362910e-15, -1.00000000e+00,  1.57066000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00 -2.94213291e-14  3.56500000e+04]
 [-3.65362910e-15 -1.00000000e+00  1.57066000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00,  0.00000000e+00,  3.99160000e+04],
     [-9.61481343e-17, -1.00000000e+00,  1.57066000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00  0.00000000e+00  3.99160000e+04]
 [-9.61481343e-17 -1.00000000e+00  1.57066000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00,  1.92296269e-16,  4.41830000e+04],
     [-2.40370336e-17, -1.00000000e+00,  1.57066000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00  1.92296269e-16  4.41830000e+04]
 [-2.40370336e-17 -1.00000000e+00  1.57066000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[ 1.00000000e+00,  1.92296269e-16,  4.84500000e+04],
     [ 1.92296269e-16, -1.00000000e+00,  1.57066000e+05],
     [ 0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])>
printing out.params
[[ 1.00000000e+00  1.92296269e-16  4.84500000e+04]
 [ 1.92296269e-16 -1.00000000e+00  1.57066000e+05]
 [ 0.00000000e+00  0.00000000e+00  1.00000000e+00]]
printing out
<AffineTransform(matrix=
    [[nan, nan, nan],
     [nan, nan, nan],
     [nan, nan, nan]])>
printing out.params
[[nan nan nan]
 [nan nan nan]
 [nan nan nan]]

Followed by error:

[/nemo/stp/babs/working/bootj/projects/swantonc/imran.noorani/in839/analysis/init-rocker-renv/venv/lib/python3.12/site-packages/numpy/core/numeric.py:407](http://localhost:8888/venv/lib/python3.12/site-packages/numpy/core/numeric.py#line=406): RuntimeWarning: invalid value encountered in cast
  multiarray.copyto(res, fill_value, casting='unsafe')
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[9], line 2
      1 # Trial the spatialdata_io way of loading data - FAILS - GitHub issue raised.
----> 2 adata1 = sp.cosmx(sample_dir, dataset_id = 'GCIN10692_TMA2')

File [/nemo/stp/babs/working/bootj/projects/swantonc/imran.noorani/in839/analysis/init-rocker-renv/venv/lib/python3.12/site-packages/spatialdata_io/readers/cosmx.py:152](http://localhost:8888/venv/lib/python3.12/site-packages/spatialdata_io/readers/cosmx.py#line=151), in cosmx(path, dataset_id, transcripts, imread_kwargs, image_models_kwargs)
    150     print('printing out.params')
    151     print(out.params)
--> 152     affine_transforms_to_global[fov] = Affine(
    153         # out.params, input_coordinate_system=input_cs, output_coordinate_system=output_cs
    154         out.params,
    155         input_axes=("x", "y"),
    156         output_axes=("x", "y"),
    157     )
    159 table.obsm["global"] = table.obs[[CosmxKeys.X_GLOBAL_CELL, CosmxKeys.Y_GLOBAL_CELL]].to_numpy()
    160 table.obsm["spatial"] = table.obs[[CosmxKeys.X_LOCAL_CELL, CosmxKeys.Y_LOCAL_CELL]].to_numpy()

File [/nemo/stp/babs/working/bootj/projects/swantonc/imran.noorani/in839/analysis/init-rocker-renv/venv/lib/python3.12/site-packages/spatialdata/transformations/transformations.py:523](http://localhost:8888/venv/lib/python3.12/site-packages/spatialdata/transformations/transformations.py#line=522), in Affine.__init__(self, matrix, input_axes, output_axes)
    521     raise ValueError("Invalid shape for affine matrix.")
    522 if not np.array_equal(self.matrix[-1, :-1], np.zeros(len(input_axes))):
--> 523     raise ValueError("Affine matrix must be homogeneous.")
    524 assert self.matrix[-1, -1] == 1.0

ValueError: Affine matrix must be homogeneous.

jamesboot avatar Oct 07 '24 09:10 jamesboot

I believe this is caused by your source data including FOVs with just a single cell. If these FOVs are kept, you run into a problem when estimating the affine transformation from local to global coordinates. Attempting to do that with just a pair of a single coordinate leads to estimate_transform() returning a np.nan matrix, which in turn throws the error you're seeing. I've fixed this locally with the following patch:

diff --git a/src/spatialdata_io/readers/cosmx.py b/src/spatialdata_io/readers/cosmx.py
index 57e9f5d..343b940 100644
--- a/src/spatialdata_io/readers/cosmx.py
+++ b/src/spatialdata_io/readers/cosmx.py
@@ -129,6 +129,10 @@ def cosmx(
     )
     adata.var_names = counts.columns

+    # Filter out single-cell FOVs since we cannot define a transform to global from a single cell
+    num_cells = adata.obs[['fov']].groupby('fov').size()
+    adata = adata[adata.obs['fov'].isin(num_cells[num_cells > 2].index)]
+
     table = TableModel.parse(
         adata,
         region=list(set(adata.obs[CosmxKeys.REGION_KEY].astype(str).tolist())),

happy to open a PR if that's helpful.

clwgg avatar Oct 07 '24 21:10 clwgg

Thanks, that patch works for me! Might be useful to have this incorporated. I suspect this is because we have round tissue sections, but the FOVs are square grids therefore inevitably there are FOVs on the edge of the tissue which aren't capturing much tissue if any tissue at all.

Thanks again!

jamesboot avatar Oct 09 '24 09:10 jamesboot

Thanks for the details @jamesboot and @clwgg for the solution. @clwgg could you please make a PR with the fix? Thanks 😊

LucaMarconato avatar Oct 09 '24 13:10 LucaMarconato