Proposal: spatial-transcriptomics adapter + Visium proof-of-concept

[brendancol]

Reason or Problem

Spatial transcriptomics (ST) pairs spatial coordinates with high-dimensional gene-expression vectors, the same shape of data xarray-spatial already models as a many-band, Dask/GPU-backed DataArray. Two ST regimes are becoming raster-native:

• 10x Visium HD produces a continuous 2 µm bin grid: a regular north-up lattice with one value per gene per bin, structurally identical to a multiband raster.
• Imaging platforms (MERFISH, seqFISH, Xenium, CosMX) emit individual transcript molecules or cell centroids as points, at billions-of-points scale.

The existing ST tooling (Squidpy on AnnData/GeoPandas, SpatialData) covers the statistics catalog well, but it is graph/point oriented and handles very large dense grids awkwardly. xarray-spatial's out-of-core Dask + GPU raster engine fits that scale regime, and none of its raster ops are currently reachable from an ST data model.

User persona

"Priya", a computational biologist / bioinformatics engineer on a spatial-omics team. She works in the scanpy/scverse Python stack and is comfortable with AnnData and numpy, but she is not a GIS specialist. Her Visium HD runs are large enough that per-gene spatial statistics in the point/graph tools are slow, and she has GPUs available but no easy way to use them for spatial-domain analysis. She wants to find spatially variable genes and expression hotspots over a tissue section without rewriting her pipeline.

Proposal

A thin adapter that bridges an ST data object into xarray-spatial's raster model, then lets the existing ops run with Dask/GPU scale-out. This is not a reimplementation of Squidpy. It is a bridge plus a documented workflow.

Design:

• A from_anndata(adata, genes=..., like=None, resolution=...) helper in a new small module (for example xrspatial.experimental.transcriptomics) that takes spot/bin coordinates from adata.obsm['spatial'] and expression columns and returns a (gene, y, x) DataArray. Visium HD bins map directly onto their grid; classic Visium spots bin to their lattice. The point-to-grid step is a plain coordinate-binning scatter (see the proof-of-concept below); rasterize stays available for the polygon and line cases.
• Once the data is in the DataArray model, the existing ops apply unchanged: focal.hotspots (Getis-Ord Gi*) per gene for expression hotspots, focal.focal_stats/convolve_2d for smoothing, proximity for cell-type distance fields, zonal for per-domain aggregation, and kde for molecule density from imaging platforms.
• Keep it at the experimental/advanced tier with no new hard dependencies. AnnData is handled as an optional import, following the existing optional-import pattern.

Usage / Visium proof-of-concept (shapely-free):

import numpy as np
import xarray as xr
from xrspatial.focal import hotspots
from xrspatial.convolution import circle_kernel

# --- inputs from a Visium/Visium HD AnnData ---
# xy: (n_spots, 2) tissue-space coordinates from adata.obsm['spatial']
# expr: (n_spots,) counts for one gene of interest (e.g. adata[:, 'EPCAM'].X)
xy   = np.asarray(adata.obsm['spatial'], dtype='float64')
expr = np.asarray(adata[:, gene].X).ravel().astype('float64')

# 1. bin spot coordinates onto a regular grid -> (y, x) DataArray, no GDAL/shapely
x0, y0 = xy.min(0)
ix = np.floor((xy[:, 0] - x0) / bin_size).astype('int64')
iy = np.floor((xy[:, 1] - y0) / bin_size).astype('int64')
h, w = iy.max() + 1, ix.max() + 1

grid = np.full((h, w), np.nan, dtype='float64')
grid[iy, ix] = expr                                    # scatter expression into bins
agg = xr.DataArray(
    grid, dims=('y', 'x'),
    coords={'y': y0 + (np.arange(h) + 0.5) * bin_size,
            'x': x0 + (np.arange(w) + 0.5) * bin_size},
    name=gene,
)

# 2. Getis-Ord Gi* hotspots of expression over a local neighborhood
kernel = circle_kernel(bin_size, bin_size, radius=3 * bin_size)
hot = hotspots(agg, kernel)   # signed z-score classes: hot / cold / not-significant

# 3. `hot` is a north-up raster of statistically significant expression
#    hot/cold spots for that gene, ready to overlay on the H&E image,
#    and it runs on numpy, cupy, or Dask by swapping the input backend
#    (wrap `grid` with dask.array / cupy before building `agg`).

For Visium HD the same code runs against the native 2 µm bin grid; for classic Visium it runs against the ~55 µm spot lattice. Swapping agg for a Dask- or CuPy-backed DataArray scales the identical call out-of-core or onto the GPU.

Value: Reuses shipped, cross-backend ops (rasterize, focal.hotspots, proximity, zonal, kde) to reach a new domain with a small surface area. The differentiator is scale: GPU/Dask focal and hotspot statistics over dense expression grids, where the incumbent point/graph tooling slows down. It also grows the user base beyond GIS and earth observation into spatial-omics.

Stakeholders and Impacts

Stakeholders are the maintainers (a small new experimental module plus an optional AnnData import) and ST practitioners in the scverse ecosystem. The impact is additive because it reuses existing ops; the only new code is the adapter and docs/example. It touches a new experimental module, optional-import handling, and one user-guide example notebook.

Drawbacks

• Squidpy and SpatialData already own the ST analysis workflow and integrate with single-cell downstream steps this library does not touch, so there is a risk of duplicating a served niche if scope creeps past the bridge.
• ST domain conventions (coordinate systems, units, QC) differ from GIS, and getting them subtly wrong erodes trust with the target persona.
• It adds a domain the maintainers may not have expertise to support long-term.

Alternatives

• Do nothing and point ST users at Squidpy (status quo).
• Ship only documentation showing how to hand-roll the DataArray from AnnData, with no adapter code.
• Contribute the Dask/GPU raster ops upstream to SpatialData instead of pulling ST into this repo.

Unresolved Questions

• Scope: bridge-only, or a curated set of ST-flavored convenience wrappers?
• Where the adapter lives (experimental vs a separate contrib package) and how AnnData is depended on.
• Whether to target Visium HD first (the cleanest raster fit) and defer imaging platforms.

Additional Notes or Context

The best first deliverable is the Visium proof-of-concept above, written as a user-guide example against a public Visium HD dataset, to validate the ergonomics before committing to an adapter API.