Is there a way to subset cell_boundaries using an AnnData table?

[Pancreas-Pratik]

This is what I have done with cell_circles, I was wondering if I could do the same with cell_boundaries?

Initially, I more or less did this to create a subset AnnData:

#keep cells that have 10 counts or greater (remove cells that have less than 10 counts)
sc.pp.filter_cells(sdata.tables["table"], min_counts=10)

#keep genes that are in 5 cells or greater (remove genes that are expressed in less than 5 cells)
sc.pp.filter_genes(sdata.tables["table"], min_cells=5)

and then I did this to subset cell_circles within the sdata:

sdata_filtered=spatialdata.match_sdata_to_table(sdata=sdata, table_name='table', table=sdata.tables["table"], how='right')

sdata_filtered

SpatialData object
├── Shapes
│     └── 'cell_circles': GeoDataFrame shape: (154472, 2) (2D shapes)
└── Tables
      └── 'table': AnnData (154472, 377)
with coordinate systems:
    ▸ 'global', with elements:
        cell_circles (Shapes)

sdata.shapes['cell_circles']=sdata_filtered.shapes['cell_circles']

Now, I was wondering if I could do something similar to this type of subset that was done on cell_circles, but if I could do it on cell_boundaries also. I tried match_element_to_table(), however it did not work, I could post the error if needed, but I suspect that this is known and expected.

I believe this has to do with the dimensions of GeoDataFrame shape, and how cell_circles appears to have "something" that cell_boundaries does not have:

'cell_boundaries': GeoDataFrame shape: (162254, 1) (2D shapes)
'cell_circles': GeoDataFrame shape: (154472, 2) (2D shapes)

Thank you very much in advance.

[LucaMarconato]

Hi @Pancreas-Pratik please find the explanation and solution to the problem reported in this code here

##
import numpy as np
import os
import scanpy as sc
import spatialdata as sd

# please adjust the path to your dataset
f = "../spatialdata-sandbox/xenium_rep1_io/data.zarr"
assert os.path.exists(f)

sdata = sd.read_zarr(f)
sdata

# keep cells that have 10 counts or greater (remove cells that have less than 10 counts)
sc.pp.filter_cells(sdata.tables["table"], min_counts=10)

# keep genes that are in 5 cells or greater (remove genes that are expressed in less than 5 cells)
sc.pp.filter_genes(sdata.tables["table"], min_cells=5)

sdata_filtered = sd.match_sdata_to_table(
    sdata=sdata, table_name="table", table=sdata.tables["table"], how="right"
)


# this gives:
#     UserWarning: The element `cell_boundaries` is not annotated by the table. Skipping
filtered_cell_boundaries = sd.match_element_to_table(
    sdata, element_name="cell_boundaries", table_name="table"
)[0]["cell_boundaries"]
assert len(filtered_cell_boundaries) != len(sdata["table"])

# solution: make the table annotate the cell_boundaries instead of the cell_circles

# the indices are mapped correctly (here we do not need that the columns are the
# same, but that the circles i corresponds to the boundaries i, for all i. This is true
# for Xenium data by construction)
assert np.array_equal(sdata["cell_circles"].index, sdata["cell_boundaries"].index)

sdata["table"].obs["region"] = "cell_boundaries"
sdata.set_table_annotates_spatialelement(table_name="table", region="cell_boundaries")

# now it works
filtered_cell_boundaries = sd.match_element_to_table(
    sdata, element_name="cell_boundaries", table_name="table"
)[0]["cell_boundaries"]
assert len(filtered_cell_boundaries) == len(sdata["table"])

Two comments:

• we are currently improving the ergonomics around these type of operations with a new API called match_sdata_to_table(), merged here allow filtering by ids #627 and with a work-in-progress PR called for a new API called filter_table_by_query(), discussed here Added filter_table_by_query #894. Also, squidpy will offer APIs similar to the scanpy ones. The implementation for the moment will be separate because first we need to enable the join APIs (used in the functions above), to return a view and not always a copy. This is being worked out in this spatialdata PR here refactor filter_table_by_elements #701. The squidpy PR is this one here: Add filter function that subsets the entire object squidpy#967
• Please when reporting bug attach a fully copy-pastable reproducible example. In general, we tend to prioritize issues where this is available. Thanks!

[Pancreas-Pratik]

Thank you very much @LucaMarconato
This worked!

My apologies for not providing a fully copy-pastable reproducible example.

Respectfully,
Pratik

[StevenWijnen]

I have a similar question but than for a CosMx dataset.
Since I want to subset the spatialdata obejct based on the filtered cells I'ce converted my labels to shapes using to_polygon as suggested. I got the same warning

UserWarning: The element 
`2_shapes` is not annotated by the table. Skipping
  return self.value(*args). 

But when I try to set 2_shapes as suggested here I get the following error:

`in check_target_region_column_symmetry(table, region_key, target)
74 symmetric_difference = found_regions.symmetric_difference(target_element_set)
75 if symmetric_difference:
---> 76 raise ValueError(
77 f"Mismatch(es) found between regions in region column in obs and target element: "
78 f"{', '.join(diff for diff in symmetric_difference)}"
79 )

ValueError: Mismatch(es) found between regions in region column in obs and target element: 2_labels, 2_shapes`

Full code:

sdata = sd.read_zarr(path + "/zarr_store")
for e in sdata.labels.keys():
    id = e.replace("_labels", "")
    label = sdata.labels[e]
    s = sd.to_polygons(label)
    sdata.shapes[f"{id}_shapes"] = s

fov_2 = sdata.subset(["2_labels", "2_shapes", "2_points"])
fov_2.table.obs["fov_labels"] = fov_2.table.obs["fov_labels"].cat.remove_unused_categories()

fov_2["table"].obs["region"] = "2_shapes"
fov_2.set_table_annotates_spatialelement(table_name="table", region="2_shapes")

[Pancreas-Pratik]

@StevenWijnen I just tried this pull request forfilter_table_by_query(): #894
It is amazing for filtering, and almost instant with subsetting/filtering. You may be able to bypass the error you have by using it.

If it saves you any time, this how I installed the dev pull request commit:
pip install git+https://github.com/scverse/spatialdata.git@3fb202eb3f4a6482d8e20778c30e0790503a2209

You can double check on the bottom of https://github.com/scverse/spatialdata/pull/894/commits for the full SHA.

EDIT: @StevenWijnen I apologize. Upon reading your comment more carefully, I realize that perhaps filter_table_by_query may not be the answer you need.