convert: Visium HD bins (spatial extra)

AGENTS.md

@@ -10,7 +10,9 @@ point of the tool is that memory is set by `--chunk`, not by the size of the
 input, so code that loads a whole matrix or column is a regression unless it is
 behind an explicit `--in-memory`. pandas and scipy are not dependencies either:
 they come only with the `pickle` extra, for pycisTopic pickles
-(`conversion/cistopic.py`), and are imported there alone.
+(`conversion/cistopic.py`), and are imported there alone. pyarrow likewise
+comes only with the `spatial` extra, for Visium HD's Parquet positions, and
+is imported only in `formats/parquet.py`.
 
 Naming:
 - The PyPI distribution is `pyadata-cli`.

Dockerfile

@@ -61,9 +61,10 @@ COPY . .
 # them; that reclaims about 9 MB, the cache being mostly hardlinks into the
 # venv rather than separate copies.
 #
-# The `pickle` extra (pandas, scipy) is in: a pipeline cannot add it to the
-# image later, and without it pycisTopic pickles do not convert.
-RUN uv sync --locked --extra pickle \
+# The `pickle` (pandas, scipy) and `spatial` (pyarrow) extras are in: a
+# pipeline cannot add them to the image later, and without them pycisTopic
+# pickles and Visium HD do not convert.
+RUN uv sync --locked --extra pickle --extra spatial \
     && rm -rf /tmp/.cache /tmp/uv-*.lock \
     && mkdir -p /tmp/.cache \
     && chmod 1777 /tmp/.cache

README.md

@@ -21,8 +21,9 @@ adata convert filtered_feature_bc_matrix.h5 -o counts.h5ad
 ```
 
 The command is `adata`; the distribution is `pyadata-cli`, because `adata-cli` was
-taken on PyPI. `pip install 'pyadata-cli[pickle]'` adds pycisTopic pickles. The
-container image is `quay.io/cellgeni/adata-cli`:
+taken on PyPI. `pip install 'pyadata-cli[pickle]'` adds pycisTopic pickles, and
+`'pyadata-cli[spatial]'` Visium HD. The container image, with both, is
+`quay.io/cellgeni/adata-cli`:
 
 ```bash
 docker run --rm -v "$PWD":/data quay.io/cellgeni/adata-cli:latest adata view /data/data.h5ad

STATUS.md

@@ -20,13 +20,29 @@ The conversion epic's sub-issues (#21–#30) are all done and ticked in #20.
 
 ## In progress
 
-Nothing. Spatial (standard Visium and Xenium) is out in 0.17.0; Visium HD
-is first among the follow-ups below.
+### Visium HD (`feat/visium-hd`, for 0.17.1)
+
+Decisions:
+- One bin size per output, `--bin-size` (8 µm by default, as 10x
+  recommends); a bin folder given directly reads as it is.
+- pyarrow comes with a new `spatial` extra, in the image too, and is
+  imported only by `formats/parquet.py`.
+- The Parquet positions are read as a buffered stream, a set of columns per
+  pass; the join holds barcodes as bytes; images are written in strips as
+  they are read. 2 µm on 10x's tiny sample: 3.5 s, 470 MB without images.
+
+Steps:
+- [x] `formats/parquet.py`, HD in `visium.py`, `--bin-size`
+- [x] Tests: HD fixture, flags and refusals, the missing extra, a performance
+  guard, a scanpy and pandas comparison
+- [x] Real data: 10x's SR 4.0.1 tiny HD sample at 2, 8 and 16 µm
+- [x] Docs
+- [ ] PR, review, merge; release 0.17.1
 
 ## Follow-ups noticed
 
-- Visium HD bins and Space Ranger 4's cell segmentation, behind a `spatial`
-  extra (pyarrow) for their Parquet positions.
+- Space Ranger 4's cell segmentation (`segmented_outputs/`): cells as
+  GeoJSON polygons, mapped to bins by `barcode_mappings.parquet`.
 - Spatial bundles whose matrix is a MEX folder only (refused for now, with a
   pointer to `--from 10x-mtx`).
 - Space Ranger's own `analysis/` (clusters, UMAP) for Visium, as Xenium's.

docs/developers.md

@@ -19,9 +19,10 @@ uv run adata --help
 uv run pytest -m "not integration"     # the unit suite CI gates on, a few minutes
 ```
 
-The tool itself depends on h5py, zarr, numpy, Pillow, Typer and Rich.
-anndata, pandas and scipy come with the `dev` extra, as the reference the
-tests check against; scanpy, loompy, mudata and older anndata releases are
+The tool itself depends on h5py, zarr, numpy, Pillow, Typer and Rich; the
+`pickle` extra adds pandas and scipy, and the `spatial` extra pyarrow.
+anndata, pandas, scipy and pyarrow come with the `dev` extra, as the
+reference the tests check against; scanpy, loompy, mudata and older anndata releases are
 installed by the `integration` tests in environments of their own.
 
 ## Repository layout

docs/guides/convert.md

@@ -41,6 +41,7 @@ no extension.
 |---|---|---|
 | Cell Ranger's `*_feature_bc_matrix.h5` | `adata convert in.h5 -o out.h5ad` | [10x HDF5](../reference/formats/10x-h5.md) |
 | A Space Ranger `outs/` folder (standard Visium) | `adata convert outs -o out.h5ad` | [Visium](../reference/formats/visium.md) |
+| A Visium HD `outs/` folder | `adata convert outs -o out.h5ad --bin-size 8` (needs `pyadata-cli[spatial]`) | [Visium HD](../reference/formats/visium.md#visium-hd) |
 | A Xenium output bundle | `adata convert bundle -o out.h5ad` | [Xenium](../reference/formats/xenium.md) |
 | CellBender `output.h5` | `adata convert output.h5 -o out.h5ad` | [CellBender](../reference/formats/cellbender.md) |
 | A `matrix.mtx` folder from Cell Ranger, STARsolo or GEO | `adata convert folder -o out.h5ad` | [MTX folders](../reference/formats/10x-mtx.md) |

docs/guides/pipelines.md

@@ -13,7 +13,7 @@ nav_order: 6
 ## The container image
 
 Every release is published as `quay.io/cellgeni/adata-cli:<version>`, and the
-latest as `:latest`. The image has the `pickle` extra and
+latest as `:latest`. The image has the `pickle` and `spatial` extras and
 [duckdb](https://duckdb.org) in it, and sets no entrypoint, as Nextflow
 needs, so name the command:
 

docs/index.md

@@ -30,13 +30,14 @@ container or a pipeline step where the scientific stack would be overkill.
 ```bash
 pip install pyadata-cli
 pip install 'pyadata-cli[pickle]'   # also read pycisTopic pickles (adds pandas, scipy)
+pip install 'pyadata-cli[spatial]'  # also read Visium HD (adds pyarrow)
 ```
 
 The command is `adata`. The distribution is named `pyadata-cli` because
 `adata-cli` was already taken on PyPI by an unrelated project. Python 3.12 or
 later.
 
-Or use the container, which has the `pickle` extra and duckdb in it:
+Or use the container, which has both extras and duckdb in it:
 
 ```bash
 docker run --rm -v "$PWD":/data quay.io/cellgeni/adata-cli:latest \

docs/reference/commands.md

@@ -120,6 +120,7 @@ adata convert counts.h5ad -o for_cellbender.h5 --layer counts
 adata convert Solo.out/Gene/filtered -o counts.h5ad
 adata convert counts.h5ad -o counts_mtx --to 10x-mtx
 adata convert spaceranger_run/outs -o section.h5ad
+adata convert visium_hd_run/outs -o bins.h5ad --bin-size 8
 adata convert xenium_bundle -o cells.h5ad
 adata convert splitpipe_out -o parse.h5ad
 adata convert velocyto.loom -o velocyto.h5ad
@@ -146,7 +147,7 @@ adata convert multiome.h5ad -o multiome.h5mu
 | `h5ad` | `.h5ad` | yes | yes |
 | `zarr` | `.zarr` (v2 or v3) | yes | yes |
 | [`10x-h5`](formats/10x-h5.md) | `.h5` | Cell Ranger 2 and 3+; the matrix alone of Space Ranger and Xenium | Cell Ranger 3+ layout |
-| [`visium`](formats/visium.md) | a folder | Space Ranger `outs/`: counts, spot positions, images (standard Visium; not HD yet) | no |
+| [`visium`](formats/visium.md) | a folder | Space Ranger `outs/`: counts, spot or bin positions, images (standard Visium; Visium HD with `pyadata-cli[spatial]`) | no |
 | [`xenium`](formats/xenium.md) | a folder | a Xenium bundle: counts, cells and centroids, Xenium's analysis | no |
 | [`cellbender`](formats/cellbender.md) | `.h5` | CellBender `remove-background` 0.2 and later | no |
 | [`10x-mtx`](formats/10x-mtx.md) | a folder (`--to 10x-mtx`) | Cell Ranger 2 and 3+, STARsolo (velocyto too) | Cell Ranger 3+ layout, gzipped |

docs/reference/formats/index.md

@@ -18,7 +18,7 @@ flags every conversion takes (`--output`, `--from`, `--to`, `--compression`,
 | [10x Genomics HDF5](10x-h5.md) | `10x-h5` | Cell Ranger 2 and 3+; the matrix alone of Space Ranger and Xenium | Cell Ranger 3+ layout |
 | [CellBender](cellbender.md) | `cellbender` | `remove-background` 0.2 and later, with its latents | no |
 | [10x and STARsolo MTX folders](10x-mtx.md) | `10x-mtx` | Cell Ranger 2 and 3+, STARsolo (velocyto too), GEO prefixes | Cell Ranger 3+ layout, gzipped |
-| [Visium (Space Ranger)](visium.md) | `visium` | an `outs/` folder: counts, spot positions, tissue images, scale factors (Space Ranger 1–4; not HD yet) | no |
+| [Visium (Space Ranger)](visium.md) | `visium` | an `outs/` folder: counts, spot or bin positions, tissue images, scale factors (Space Ranger 1–4; Visium HD needs `pyadata-cli[spatial]`) | no |
 | [Xenium](xenium.md) | `xenium` | an output bundle: counts, cells and centroids, Xenium's clusters, UMAP and PCA | no |
 | [Parse Biosciences](parse.md) | `parse` | split-pipe: a run, a sample or a `DGE_*` folder | no |
 | [loom and SCope](loom.md) | `loom` | loompy, velocyto, kallisto\|bustools, SCope; spec 2 and 3 | spec 3.0.0; a SCope loom with `--scope` |

docs/reference/formats/visium.md

@@ -7,11 +7,13 @@ nav_order: 4
 
 # Visium (Space Ranger)
 
-A Space Ranger output folder for standard Visium (55 µm spots), from Space
-Ranger 1 to 4: the counts, where each spot sits on the slide, the tissue
+A Space Ranger output folder, for standard Visium (55 µm spots, Space Ranger
+1 to 4) and for [Visium HD](#visium-hd) (2, 8 and 16 µm bins, Space Ranger 3
+and 4): the counts, where each spot or bin sits on the slide, the tissue
 images and the scale factors between them. The result is laid out as
 `squidpy.read.visium` lays it out, so scanpy's and squidpy's spatial plots
-work on it, and CI checks it against squidpy.
+work on it. CI checks standard Visium against squidpy, and Visium HD, which
+squidpy does not read, against scanpy with pandas reading the positions.
 
 Give the `outs/` folder, or the run folder holding it. It is recognised by
 `spatial/` (the positions and `scalefactors_json.json`) beside a
@@ -67,7 +69,42 @@ longest side.
 Refused, with the reason:
 
 - output of `spaceranger aggr`, whose matrix holds several libraries;
-- Visium HD (`binned_outputs/`), whose positions are Parquet and are not
-  read yet. Each bin's `filtered_feature_bc_matrix.h5` converts as a
-  [10x HDF5 file](10x-h5.md) meanwhile;
 - a folder with the matrix as Matrix Market only, with no `.h5`.
+
+## Visium HD
+
+Space Ranger 3 and 4 write a Visium HD run's counts at three bin sizes, each
+in its own folder under `binned_outputs/` (`square_002um`, `square_008um`,
+`square_016um`), laid out as a standard Visium `outs/` but with the
+positions as `spatial/tissue_positions.parquet`. Give the `outs/` folder and
+`--bin-size` picks the bins, 8 µm by default as 10x recommends for
+analysis; or give a bin folder itself.
+
+```bash
+pip install 'pyadata-cli[spatial]'    # pyarrow, which reads the Parquet positions
+
+adata convert outs -o bins_8um.h5ad
+adata convert outs -o bins_2um.h5ad --bin-size 2
+adata convert outs/binned_outputs/square_016um -o bins_16um.h5ad
+```
+
+| Flag | Meaning |
+|---|---|
+| `--bin-size` | The bins to read, in µm (default 8; refused for standard Visium) |
+
+The result is laid out as above: `in_tissue`, `array_row` and `array_col` in
+obs, `obsm["spatial"]` in full-resolution pixels (float64, as Space Ranger
+records them), and the bin folder's own images and scale factors, which
+carry `bin_size_um` and `microns_per_pixel`, in `uns["spatial"]`.
+
+The positions list every bin on the slide, about 11 million at 2 µm, where
+the filtered matrix holds the bins under tissue. They are read twice, a
+buffered stream of only the columns needed: the barcodes, matched to the
+matrix's as they come, then the numbers, written in the matrix's order. On
+10x's Visium HD tiny sample, the 2 µm bins (330,035 in the matrix) convert in
+3.5 seconds with a peak of 470 MB with `--no-images`, and 610 MB with the
+100 MB hires image; pandas alone takes about 1 GB to read the positions. `--no-images` leaves the
+images out.
+
+The cell segmentation Space Ranger 4 adds (`segmented_outputs/`) is not read
+yet; its matrix converts as a [10x HDF5 file](10x-h5.md).

docs/troubleshooting.md

@@ -152,12 +152,11 @@ A coordinate, centroid or projection column holds text. Missing values are
 fine (they become NaN); a word is not. The file has been edited, or is not
 what its name says.
 
-### `This is Visium HD output, which is not read yet: its positions are Parquet.`
+### `This Visium HD output has no 2 µm bins; it has 8 µm, 16 µm. Pick one with --bin-size.`
 
-Visium HD's bin positions exist only as `tissue_positions.parquet`, which a
-later release will read. Meanwhile, convert one bin size's matrix as a plain
-10x matrix, without positions:
-`adata convert binned_outputs/square_008um/filtered_feature_bc_matrix.h5 -o bins.h5ad`.
+`--bin-size` names a folder under `binned_outputs/` (`square_002um` for 2).
+Space Ranger writes 2, 8 and 16 µm bins, but a copied or trimmed output may
+hold fewer; the message lists those there are.
 
 ### `This matrix holds 2 libraries (A1, B1), as spaceranger aggr writes`
 
@@ -252,6 +251,11 @@ wins alongside those where it loses.
 pandas and scipy are not dependencies of `adata`; the `pickle` extra adds
 them, and only pycisTopic pickles need them. The container image has them.
 
+### `Reading Visium HD positions needs pyarrow, which reads its Parquet tables: pip install 'pyadata-cli[spatial]'.`
+
+Visium HD records where each bin sits only as Parquet. The `spatial` extra
+adds pyarrow, which only Visium HD needs; the container image has it.
+
 ## Reading the output elsewhere
 
 ### anndata warns `OldFormatWarning`

pyproject.toml

@@ -59,9 +59,15 @@ pickle = [
     "pandas>=2.2",
     "scipy>=1.13",
 ]
+# Reading Visium HD: its bin positions are Parquet only. Imported only by
+# the Visium converter, and only for an HD input.
+spatial = [
+    "pyarrow>=17",
+]
 dev = [
     "pandas>=2.2",
     "scipy>=1.13",
+    "pyarrow>=17",
     "pytest>=8.3.4",
     "pytest-cov>=6.0.0",
     # Turns a hang into a failure with a traceback. Without it a stuck test

src/adata/cli.py

@@ -530,6 +530,7 @@ def concat(
     "unfiltered": "--unfiltered",
     "library_id": "--library-id",
     "no_images": "--no-images",
+    "bin_size": "--bin-size",
     "dense": "--dense",
     "obs_names_attr": "--obs-names-attr",
     "var_names_attr": "--var-names-attr",
@@ -712,6 +713,12 @@ def convert(
         help="Visium input: leave the tissue images out of uns['spatial'] "
         "(the scale factors stay)",
     ),
+    bin_size: Optional[int] = typer.Option(
+        None,
+        "--bin-size",
+        help="Visium HD input: the bin size to read, in µm (default 8; "
+        "Space Ranger writes 2, 8 and 16)",
+    ),
     dense: bool = typer.Option(
         False,
         "--dense",
@@ -809,8 +816,9 @@ def convert(
     CellBender's own loader names them. A 10x matrix is read as scanpy's
     read_10x_h5 or read_10x_mtx (gex_only=False) reads it, keeping integer
     counts, and written as Cell Ranger writes it. A Space Ranger outs/
-    folder (standard Visium) keeps its spot positions, images and scale
-    factors where squidpy puts them, and a Xenium bundle its cell table,
+    folder keeps its spot positions, images and scale factors where squidpy
+    puts them; a Visium HD one is read at the bin size --bin-size names
+    (needs pyadata-cli[spatial]); and a Xenium bundle its cell table,
     centroids and Xenium's own clusters, UMAP and PCA. MTX folders from Cell
     Ranger and STARsolo, velocyto's included, and Parse split-pipe output
     are read in two streaming passes, as are text tables (.csv, .tsv, .txt,
@@ -840,6 +848,7 @@ def convert(
         adata convert Solo.out/Gene/filtered -o counts.h5ad
         adata convert Solo.out/Velocyto/raw -o velocyto.h5ad
         adata convert spaceranger_run/outs -o section.h5ad
+        adata convert visium_hd_run/outs -o bins.h5ad --bin-size 8
         adata convert xenium_bundle -o cells.h5ad
         adata convert counts.h5ad -o counts_mtx --to 10x-mtx --layer counts
         adata convert cellbender_output.h5 -o cb.h5ad --raw raw_feature_bc_matrix.h5
@@ -946,6 +955,7 @@ def convert(
                 unfiltered=unfiltered,
                 library_id=library_id,
                 no_images=no_images,
+                bin_size=bin_size,
                 dense=dense,
                 obs_names_attr=obs_names_attr,
                 var_names_attr=var_names_attr,

src/adata/conversion/base.py

@@ -47,6 +47,7 @@ class Options:
     # Reading Space Ranger output.
     library_id: Optional[str] = None
     no_images: bool = False
+    bin_size: Optional[int] = None
 
     # Reading a loom.
     dense: bool = False