Visium hd

[ArneDefauw]

This adds support for annotating the resulting tables with corresponding labels layers (the bins).

Annotating the table by a labels layer has the advantage that the bins, and their expression levels can be viewed in napari-spatialdata, which is not possible if you annotate via the shapes layers ( it is possible when setting bins_as_squares to False, but this will give you circles, not squares, which is not ideal ).

Viewing and loading expression levels for 016um and 008um bins is quite fast, viewing the 002um bins is also still possible, althoug loading expression levels associated with the bins takes some seconds for each gene, but once loaded is quite fast and convenient to view.

See screenshot for 008um:

It is not as snappy as https://github.com/scverse/spatialdata/blob/8879aff17b8169c6bb6ff3537e4ddcee0204f168/src/spatialdata/_core/operations/rasterize_bins.py#L28, but maybe a little bit more convenient in use.

If this PR would be merged, then this function spatialdata.rasterize_bins would also need to be updated. I would remove the bins parameter, and estimate the transformation only using the tables, i.e. using table.obsm[ 'spatial' ] and the col_key and row_key in table.obs; or get the transformation from the labels layer.

[codecov-commenter]

Codecov Report

Attention: Patch coverage is 11.42857% with 31 lines in your changes missing coverage. Please review.

Project coverage is 46.08%. Comparing base (c4f6cf7) to head (c9616e3).
Report is 21 commits behind head on main.

Files with missing lines	Patch %	Lines
src/spatialdata_io/readers/visium_hd.py	11.42%	31 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #211      +/-   ##
==========================================
+ Coverage   43.68%   46.08%   +2.39%     
==========================================
  Files          23       23              
  Lines        2353     2415      +62     
==========================================
+ Hits         1028     1113      +85     
+ Misses       1325     1302      -23     

Files with missing lines	Coverage Δ
src/spatialdata_io/readers/visium_hd.py	16.83% <11.42%> (-1.46%)	⬇️

... and 1 file with indirect coverage changes

[LucaMarconato]

Thanks for the PR, it looks great, very curios to experiment with it! I have some tasks on my todo list but I hope to get to this soon.

[t-a-m-i]

Hi, great PR, it looks very interesting! I´m trying to use it for my project, but ran into an error. I´m relatively new to using PRs on github, so maybe I already did a mistake during setup.
I created a new conda environment and installed your branch via "pip install git+https://github.com/ArneDefauw/spatialdata-io.git@visium_hd", added an ipykernel for jupyter and tried to use the visium_hd reader function in a notebook. I specified my paths, but get a FileNotFoundError for following code:

from spatialdata_io import visium_hd
sdata = visium_hd(
path_read,
load_all_images=True,
annotate_table_by_labels=True,
fullres_image_file=fullres_image_file,
dataset_id=dataset_id
)

FileNotFoundError: [Errno 2] Unable to synchronously open file (unable to open file: name = '/files/03-SpaceRanger-Results/P31904_402/outs/P31904_402_feature_slice.h5', errno = 2, error message = 'No such file or directory', flags = 0, o_flags = 0)

"P31904_402_feature_slice.h5" indeed does not exist, it should simply be "feature_slice.h5" in "outs" of the "P31904_402" folder. I don´t understand why the error is thrown though, since that part of the code was not changed to my understanding?

Here is the full traceback:
FileNotFoundError Traceback (most recent call last)
Cell In[4], line 2
1 from spatialdata_io import visium_hd
----> 2 sdata = visium_hd(
3 path_read,
4 load_all_images=True,
5 annotate_table_by_labels=True,
6 fullres_image_file=fullres_image_file,
7 dataset_id=dataset_id
8 )

File ~/anaconda3/envs/sd_pr_test/lib/python3.12/site-packages/spatialdata_io/readers/visium_hd.py:115, in visium_hd(path, dataset_id, filtered_counts_file, bin_size, bins_as_squares, annotate_table_by_labels, fullres_image_file, load_all_images, imread_kwargs, image_models_kwargs, anndata_kwargs)
104 def load_image(path: Path, suffix: str, scale_factors: list[int] | None = None) -> None:
105 _load_image(
106 path=path,
107 images=images,
(...)
112 scale_factors=scale_factors,
113 )
--> 115 metadata, hd_layout = _parse_metadata(path, filename_prefix)
116 transform_matrices = _get_transform_matrices(metadata, hd_layout)
117 file_format = hd_layout[VisiumHDKeys.FILE_FORMAT]

File ~/anaconda3/envs/sd_pr_test/lib/python3.12/site-packages/spatialdata_io/readers/visium_hd.py:471, in _parse_metadata(path, filename_prefix)
470 def _parse_metadata(path: Path, filename_prefix: str) -> tuple[dict[str, Any], dict[str, Any]]:
--> 471 with h5py.File(path / f"{filename_prefix}{VisiumHDKeys.FEATURE_SLICE_FILE.value}", "r") as f5:
472 metadata = json.loads(dict(f5.attrs)[VisiumHDKeys.METADATA_JSON])
473 hd_layout = json.loads(metadata[VisiumHDKeys.HD_LAYOUT_JSON])

File ~/anaconda3/envs/sd_pr_test/lib/python3.12/site-packages/h5py/_hl/files.py:562, in File.init(self, name, mode, driver, libver, userblock_size, swmr, rdcc_nslots, rdcc_nbytes, rdcc_w0, track_order, fs_strategy, fs_persist, fs_threshold, fs_page_size, page_buf_size, min_meta_keep, min_raw_keep, locking, alignment_threshold, alignment_interval, meta_block_size, **kwds)
553 fapl = make_fapl(driver, libver, rdcc_nslots, rdcc_nbytes, rdcc_w0,
554 locking, page_buf_size, min_meta_keep, min_raw_keep,
555 alignment_threshold=alignment_threshold,
556 alignment_interval=alignment_interval,
557 meta_block_size=meta_block_size,
558 **kwds)
559 fcpl = make_fcpl(track_order=track_order, fs_strategy=fs_strategy,
560 fs_persist=fs_persist, fs_threshold=fs_threshold,
561 fs_page_size=fs_page_size)
--> 562 fid = make_fid(name, mode, userblock_size, fapl, fcpl, swmr=swmr)
564 if isinstance(libver, tuple):
565 self._libver = libver

File ~/anaconda3/envs/sd_pr_test/lib/python3.12/site-packages/h5py/_hl/files.py:235, in make_fid(name, mode, userblock_size, fapl, fcpl, swmr)
233 if swmr and swmr_support:
234 flags |= h5f.ACC_SWMR_READ
--> 235 fid = h5f.open(name, flags, fapl=fapl)
236 elif mode == 'r+':
237 fid = h5f.open(name, h5f.ACC_RDWR, fapl=fapl)

File h5py/_objects.pyx:54, in h5py._objects.with_phil.wrapper()

File h5py/_objects.pyx:55, in h5py._objects.with_phil.wrapper()

File h5py/h5f.pyx:102, in h5py.h5f.open()

[ArneDefauw]

Hi, great PR, it looks very interesting! I´m trying to use it for my project, but ran into an error. I´m relatively new to using PRs on github, so maybe I already did a mistake during setup. I created a new conda environment and installed your branch via "pip install git+https://github.com/ArneDefauw/spatialdata-io.git@visium_hd", added an ipykernel for jupyter and tried to use the visium_hd reader function in a notebook. I specified my paths, but get a FileNotFoundError for following code:

from spatialdata_io import visium_hd sdata = visium_hd( path_read, load_all_images=True, annotate_table_by_labels=True, fullres_image_file=fullres_image_file, dataset_id=dataset_id )

FileNotFoundError: [Errno 2] Unable to synchronously open file (unable to open file: name = '/files/03-SpaceRanger-Results/P31904_402/outs/P31904_402_feature_slice.h5', errno = 2, error message = 'No such file or directory', flags = 0, o_flags = 0)

"P31904_402_feature_slice.h5" indeed does not exist, it should simply be "feature_slice.h5" in "outs" of the "P31904_402" folder. I don´t understand why the error is thrown though, since that part of the code was not changed to my understanding?

Here is the full traceback: FileNotFoundError Traceback (most recent call last) Cell In[4], line 2 1 from spatialdata_io import visium_hd ----> 2 sdata = visium_hd( 3 path_read, 4 load_all_images=True, 5 annotate_table_by_labels=True, 6 fullres_image_file=fullres_image_file, 7 dataset_id=dataset_id 8 )

File ~/anaconda3/envs/sd_pr_test/lib/python3.12/site-packages/spatialdata_io/readers/visium_hd.py:115, in visium_hd(path, dataset_id, filtered_counts_file, bin_size, bins_as_squares, annotate_table_by_labels, fullres_image_file, load_all_images, imread_kwargs, image_models_kwargs, anndata_kwargs) 104 def load_image(path: Path, suffix: str, scale_factors: list[int] | None = None) -> None: 105 _load_image( 106 path=path, 107 images=images, (...) 112 scale_factors=scale_factors, 113 ) --> 115 metadata, hd_layout = _parse_metadata(path, filename_prefix) 116 transform_matrices = _get_transform_matrices(metadata, hd_layout) 117 file_format = hd_layout[VisiumHDKeys.FILE_FORMAT]

File ~/anaconda3/envs/sd_pr_test/lib/python3.12/site-packages/spatialdata_io/readers/visium_hd.py:471, in _parse_metadata(path, filename_prefix) 470 def _parse_metadata(path: Path, filename_prefix: str) -> tuple[dict[str, Any], dict[str, Any]]: --> 471 with h5py.File(path / f"{filename_prefix}{VisiumHDKeys.FEATURE_SLICE_FILE.value}", "r") as f5: 472 metadata = json.loads(dict(f5.attrs)[VisiumHDKeys.METADATA_JSON]) 473 hd_layout = json.loads(metadata[VisiumHDKeys.HD_LAYOUT_JSON])

File ~/anaconda3/envs/sd_pr_test/lib/python3.12/site-packages/h5py/_hl/files.py:562, in File.init(self, name, mode, driver, libver, userblock_size, swmr, rdcc_nslots, rdcc_nbytes, rdcc_w0, track_order, fs_strategy, fs_persist, fs_threshold, fs_page_size, page_buf_size, min_meta_keep, min_raw_keep, locking, alignment_threshold, alignment_interval, meta_block_size, **kwds) 553 fapl = make_fapl(driver, libver, rdcc_nslots, rdcc_nbytes, rdcc_w0, 554 locking, page_buf_size, min_meta_keep, min_raw_keep, 555 alignment_threshold=alignment_threshold, 556 alignment_interval=alignment_interval, 557 meta_block_size=meta_block_size, 558 **kwds) 559 fcpl = make_fcpl(track_order=track_order, fs_strategy=fs_strategy, 560 fs_persist=fs_persist, fs_threshold=fs_threshold, 561 fs_page_size=fs_page_size) --> 562 fid = make_fid(name, mode, userblock_size, fapl, fcpl, swmr=swmr) 564 if isinstance(libver, tuple): 565 self._libver = libver

File ~/anaconda3/envs/sd_pr_test/lib/python3.12/site-packages/h5py/_hl/files.py:235, in make_fid(name, mode, userblock_size, fapl, fcpl, swmr) 233 if swmr and swmr_support: 234 flags |= h5f.ACC_SWMR_READ --> 235 fid = h5f.open(name, flags, fapl=fapl) 236 elif mode == 'r+': 237 fid = h5f.open(name, h5f.ACC_RDWR, fapl=fapl)

File h5py/_objects.pyx:54, in h5py._objects.with_phil.wrapper()

File h5py/_objects.pyx:55, in h5py._objects.with_phil.wrapper()

File h5py/h5f.pyx:102, in h5py.h5f.open()

Hi @t-a-m-i ,

Does this work using spatialdata-io (last release), installed via pip install spatialdata-io? And then do:

from spatialdata_io import visium_hd
sdata = visium_hd(
path_read,
load_all_images=True,
fullres_image_file=fullres_image_file,
dataset_id=dataset_id
)

I don't think the error you encounter is related to this PR.

[t-a-m-i]

yes, I was also confused about the error since it didn´t seem related to the PR. I didn´t realise though that there was a new release. I get the same error with the latest version of spatialdata-io (v0.1.5), with v0.1.4 everything works fine, I´ll do a bug report then, thanks!

[LucaMarconato]

Thanks @t-a-m-i for opening #216, I'll follow up there.

[LucaMarconato]

Notebook relevant for the PR review: https://scverse.zulipchat.com/#narrow/channel/315824-spatial/topic/SpatialData/near/479568220

[LucaMarconato]

Thanks again @ArneDefauw and @t-a-m-i. Today I studied your code and reviewed and tested this PR. Here are my comments.

1. The PR works great and adds an important use case, because unlike the on-demand rasterize_bins(), here we can save the data to disk (as observed the performance is slower than the on-demand rasterize_bins(), but faster than when plotting shapes).

2. The example from @t-a-m-i shows that sometimes one may want to convert the shapes from an existing SpatialData .zarr Visium HD object to labels. The approach that she took was to extend the rasterize_bins() functions to replicate the code introduced in this PR. For this reason I believe that the code in this PR would better suited directly in spatialdata rather than in spatialdata-io. Some further arguments for this are that in this way this function would be accessible also for other technologies, such as Stereo-seq and Open-ST, and the code would be cleaner since right now some parts duplicate code already present in rasterize_bins().

   I therefore kindly ask you @ArneDefauw if you could open a new PR in spatialdata and move the current code there; in the current PR you would then simply call the new implementation of rasterize_bins(). To distinguish between the on-demand rasterize_bins() (which returns images) and the rasterize_bins() returning labels, I would add a new argument return_regions_as_labels: bool = False (as we do for rasterize()). Also, in spatialdata it would be important to include tests for this new functionality.

3. Last comment, there is only one logic in the code that I'd change: I would not shift the instance_id by 1 because this could break compatibility (or create confusion) with previously written datasets; rather I would add a new column to the table to be used only when mapping the table to the newly created labels object.

[ArneDefauw]

Hi @LucaMarconato , thanks for having a look!
It makes sense to move the code for this PR to spatialdata, and I probably have some time next week to do it. I will then also update this PR.

[ArneDefauw]

I think this PR is ready for review again, together with related PR scverse/spatialdata#811

Note, we probably need to update this notebook https://spatialdata.scverse.org/en/stable/tutorials/notebooks/notebooks/examples/technology_visium_hd.html after the PR's are merged. I would be happy to help.