diff --git a/AGENTS.md b/AGENTS.md
index 1b1d6f6..1d07f4f 100644
--- a/AGENTS.md
+++ b/AGENTS.md
@@ -66,6 +66,8 @@ Each layer calls only the layers below it:
- `csr_build.py` builds a CSR matrix from entries in any order, through
row-range buckets on disk. The transpose uses it, and so does anything
that reads entries out of order;
+ - `scatter.py` writes rows read in one order (a table beside a matrix) in
+ another, through the same kind of scratch buckets;
- `query.py` and `select.py` hold the expression language behind
`--obs-query`.
- **`elements/`**: the AnnData on-disk spec.
@@ -99,7 +101,11 @@ Each layer calls only the layers below it:
an output only with `--to`. Order in `FORMATS` matters: the first sniffer
that matches wins, so a more specific format (`cellbender`, which is also a
10x matrix) goes before the general one, and an output extension picks the
- first format that can be written.
+ first format that can be written. The folder formats `visium` and `xenium`
+ go before `10x-mtx`, which would take a Space Ranger `outs/` for a Cell
+ Ranger 2 genome layout. `conversion/spatial.py` matches a table to a
+ matrix's barcodes while `scan_table` reads it (`on_index`), and writes it
+ in the matrix's order.
- **`formats/`**: the per-file-type halves of `export` and `import`: npy, mtx,
csv, json and images. `mtx.py` is a streaming Matrix Market codec (two
passes in, block writes out) shared by `import`/`export sparse` and the
diff --git a/README.md b/README.md
index e650d41..e6273a1 100644
--- a/README.md
+++ b/README.md
@@ -5,7 +5,8 @@ without loading them into memory and without a Python session. `adata` inspects,
subsets, splits, concatenates, recasts and converts stores, and moves elements in
and out as CSV, JSON, `.npy`, `.mtx` or images. It implements the AnnData on-disk
spec directly against h5py and zarr, so it streams stores larger than RAM, and
-converts to and from 10x, STARsolo, CellBender, Parse, loom, MuData and more.
+converts to and from 10x, STARsolo, Visium, Xenium, CellBender, Parse, loom, MuData
+and more.
**Documentation: **
diff --git a/STATUS.md b/STATUS.md
index 5bd17f8..6a4b476 100644
--- a/STATUS.md
+++ b/STATUS.md
@@ -19,11 +19,40 @@ The conversion epic's sub-issues (#21–#30) are all done and ticked in #20.
## In progress
-Nothing. The documentation site is live (#76, released in 0.16.2); the
-spatial work below can resume.
+### Spatial import: standard Visium and Xenium (`feat/spatial`, for 0.17.0)
+
+Decisions:
+- Standard Visium (Space Ranger 1–4) and Xenium (XOA 1–4) first, with no new
+ dependency: both come as CSV, JSON and PNG. Visium HD, whose positions are
+ Parquet only, waits for a `spatial` extra with pyarrow.
+- Laid out as squidpy (Visium) and spatialdata-io's table (Xenium) lay
+ them out; images stay the PNGs' uint8 pixels. Xenium brings its own
+ clusters, UMAP and PCA, read from `analysis.tar.gz` as it ships.
+- A table beside the matrix is matched to its barcodes while scanned
+ (`conversion/spatial.py`), and written straight through when in the
+ matrix's order (Xenium's cells, on every real bundle checked) or through
+ `core/scatter.py` when not.
+
+Steps:
+- [x] `table.py`: headerless files, `on_index`, reusable columns
+- [x] `core/scatter.py`
+- [x] Visium: `visium.py`, flags `--library-id`, `--no-images`, `--unfiltered`
+- [x] Xenium: `xenium.py`, analysis from the folder or the archive
+- [x] Performance guards; squidpy and 10x-recipe comparisons; `import-xenium` benchmark
+- [x] Real data by hand: Visium SR 1.1 and SR 2.1, Xenium XOA 3.0 and 4.0 tiny
+- [x] Docs: two format pages, tables, guide, Troubleshooting, Developers
+- [ ] PR, review, merge; release 0.17.0
## Follow-ups noticed
+- Visium HD bins and Space Ranger 4's cell segmentation, behind a `spatial`
+ extra (pyarrow) for their Parquet positions.
+- 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.
+- CellBender's `_write_embeddings` onto `core/scatter.py`, which generalises
+ it.
+
- Fixed in 0.16.2 (#77): `export image` without `-o` crashed (`'NoneType'
object has no attribute 'parent'`); it now writes PNG to stdout, as every
other export writes there. Found by the docs review.
@@ -37,7 +66,5 @@ spatial work below can resume.
## Next candidates (#31)
-- Spatial: Visium / Xenium (matrix via 10x h5, coordinates to
- `obsm["spatial"]`)
- alevin-fry USA mode, kallisto | bustools
- Seurat `.h5Seurat`
diff --git a/benchmarks/cases.py b/benchmarks/cases.py
index 8e7bd69..21ab01f 100644
--- a/benchmarks/cases.py
+++ b/benchmarks/cases.py
@@ -489,6 +489,43 @@ def _py(body: str) -> str:
],
tags=["headline"],
),
+ Case(
+ name="import-xenium",
+ question=(
+ "Read a Xenium bundle of the same counts -- the cell matrix, a "
+ "cells.csv.gz of centroids and per-cell counts, and analysis "
+ "clusters -- into an h5ad with obsm['spatial']. adata-cli matches "
+ "the cell table to the matrix as it scans it and streams both; "
+ "10x's scanpy recipe loads the matrix and reads the tables with "
+ "pandas. adata-cli keeps the matrix's own gzip, scanpy writes lzf."
+ ),
+ sidecar="xenium",
+ contenders=[
+ Contender(
+ "adata-cli",
+ argv=["adata", "convert", "{sidecar}", "-o", "{output}"],
+ ),
+ Contender(
+ "scanpy",
+ env="scanpy",
+ script=_py(
+ "import os, pandas as pd, scanpy as sc\n"
+ "from pathlib import Path\n"
+ "bundle = Path(os.environ['BENCH_SIDECAR'])\n"
+ "obj = sc.read_10x_h5(bundle / 'cell_feature_matrix.h5', gex_only=False)\n"
+ "cells = pd.read_csv(bundle / 'cells.csv.gz', index_col='cell_id')\n"
+ "cells = cells.reindex(obj.obs_names)\n"
+ "obj.obsm['spatial'] = cells[['x_centroid', 'y_centroid']].to_numpy()\n"
+ "obj.obs = cells.drop(columns=['x_centroid', 'y_centroid'])\n"
+ "clusters = pd.read_csv(bundle / 'analysis/clustering/"
+ "gene_expression_graphclust/clusters.csv', index_col='Barcode')\n"
+ "obj.obs['graphclust'] = clusters['Cluster'].astype(str)"
+ ".reindex(obj.obs_names).astype('category')\n"
+ "obj.write_h5ad(OUT, compression=COMPRESSION)\n"
+ ),
+ ),
+ ],
+ ),
Case(
name="import-csv",
question=(
diff --git a/benchmarks/page_template.md b/benchmarks/page_template.md
index 46a0ec0..c64a0d2 100644
--- a/benchmarks/page_template.md
+++ b/benchmarks/page_template.md
@@ -20,10 +20,11 @@ for time, and a table that hid the cost would not be worth publishing.
## What is measured
-Nineteen cases on a GitHub-hosted runner: concat (inner, outer, `--merge same`),
+Twenty cases on a GitHub-hosted runner: concat (inner, outer, `--merge same`),
subset by query, h5ad→zarr (`convert`), `view`, `ls`, `create`, `export` dataframe /
-array / sparse, `import dataframe`, MTX folder import (`convert`), `split`, `recast` narrowing and layout, and a
-peak-RSS-against-input-size sweep. Inputs
+array / sparse, `import dataframe`, MTX folder, text table and Xenium bundle imports
+(`convert`), `split`, `recast` narrowing and layout, and a peak-RSS-against-input-size
+sweep. Inputs
are 50,000 × 20,000 at 5% density, plus the 2,000 × 36,601 shape that hung in 0.5.1.
`recast` narrowing starts from a copy of the main input stored as float64 values
with int64 indices, once with X uncompressed (as `write_h5ad` leaves it) and once
diff --git a/benchmarks/run.py b/benchmarks/run.py
index ec6e30f..3c5c9a5 100644
--- a/benchmarks/run.py
+++ b/benchmarks/run.py
@@ -187,6 +187,8 @@ def _build_sidecar(case: Case, inputs: List[Path], workdir: Path) -> Optional[Pa
return _build_mtx_sidecar(case, inputs, workdir)
if case.sidecar == "dense-csv":
return _build_dense_csv_sidecar(case, inputs, workdir)
+ if case.sidecar == "xenium":
+ return _build_xenium_sidecar(case, inputs, workdir)
if case.sidecar != "csv":
return None
# Named for its input too: a work directory reused across tiers must not
@@ -245,6 +247,74 @@ def _build_dense_csv_sidecar(case: Case, inputs: List[Path], workdir: Path) -> P
return path
+def _build_xenium_sidecar(case: Case, inputs: List[Path], workdir: Path) -> Path:
+ """A Xenium bundle holding the first input's counts, as XOA lays one out.
+
+ The matrix is exported by adata-cli (integer counts, then the 10x h5),
+ and the cell table and clusters are written in the matrix's order, as
+ Xenium writes them, with made-up centroids and per-cell counts.
+ """
+ import gzip
+
+ import h5py
+ import numpy as np
+
+ from adata.elements.read import element_len, read_str_chunk, resolve_index
+
+ path = workdir / f"{case.name}-{inputs[0].stem}-sidecar"
+ if path.exists():
+ return path
+ partial = path.with_name(path.name + ".partial")
+ if partial.exists():
+ shutil.rmtree(partial)
+ partial.mkdir()
+ counts = workdir / f"{case.name}-{inputs[0].stem}-counts.h5ad"
+ if counts.exists():
+ counts.unlink()
+ for command in (
+ ["adata", "recast", str(inputs[0]), "X", "-o", str(counts), "--dtype", "int32"],
+ ["adata", "convert", str(counts), "-o", str(partial / "cell_feature_matrix.h5"),
+ "--gene-names-column", "gene_symbol", "--feature-types-column",
+ "feature_type", "--force"],
+ ):
+ done = subprocess.run(command, capture_output=True, text=True)
+ if done.returncode != 0:
+ raise RuntimeError(f"could not build the Xenium sidecar:\n{done.stderr or done.stdout}")
+ counts.unlink()
+ (partial / "experiment.xenium").write_text('{"pixel_size": 0.2125, "major_version": 5}')
+
+ rng = np.random.default_rng(0)
+ clusters_dir = partial / "analysis" / "clustering" / "gene_expression_graphclust"
+ clusters_dir.mkdir(parents=True)
+ with h5py.File(inputs[0], "r") as handle, gzip.open(
+ partial / "cells.csv.gz", "wt", compresslevel=6
+ ) as cells, open(clusters_dir / "clusters.csv", "w") as clusters:
+ index, _ = resolve_index(handle["obs"], "obs")
+ n = element_len(index)
+ cells.write(
+ "cell_id,x_centroid,y_centroid,transcript_counts,control_probe_counts,"
+ "cell_area,nucleus_area,segmentation_method\n"
+ )
+ clusters.write("Barcode,Cluster\n")
+ step = 50_000
+ for lo in range(0, n, step):
+ names = read_str_chunk(index, lo, min(lo + step, n), {}, parent_group=handle["obs"])
+ k = len(names)
+ x, y = rng.uniform(0, 8000, k), rng.uniform(0, 8000, k)
+ counts_, probes = rng.integers(0, 500, k), rng.integers(0, 3, k)
+ area, nucleus = rng.uniform(20, 300, k), rng.uniform(5, 60, k)
+ method = rng.choice(["Segmented by boundary stain", "Segmented by nucleus expansion"], k)
+ cluster = rng.integers(1, 20, k)
+ cells.writelines(
+ f"{c},{x[i]:.4f},{y[i]:.4f},{counts_[i]},{probes[i]},{area[i]:.3f},{nucleus[i]:.3f},{method[i]}\n"
+ for i, c in enumerate(names)
+ )
+ clusters.writelines(f"{c},{cluster[i]}\n" for i, c in enumerate(names))
+ # Renamed only when whole: a failed build must not be reused as a sidecar.
+ partial.rename(path)
+ return path
+
+
def _build_mtx_sidecar(case: Case, inputs: List[Path], workdir: Path) -> Path:
"""A Cell Ranger 3 MTX folder holding the first input's counts.
diff --git a/docs/developers.md b/docs/developers.md
index e920efa..9587174 100644
--- a/docs/developers.md
+++ b/docs/developers.md
@@ -60,7 +60,9 @@ cli.py ──► commands/ ──► core/, conversion/, formats/ ──► elem
- **`core/`** holds the algorithms that stream. `subset.py`'s writer serves
`subset` and `split`, and reads each matrix once for every output.
`csr_build.py` builds CSR from entries in any order through row-range buckets
- on disk; the transpose and every out-of-order reader use it.
+ on disk; the transpose and every out-of-order reader use it. `scatter.py`
+ does the same for the rows of a table that comes in another order than the
+ matrix beside it.
- **`elements/spec.py`** is the single source of truth for `encoding-type` and
`encoding-version`. `read.py` accepts every layout from anndata 0.7 on;
`write.py` always writes the current one; `strings.py` reconciles how HDF5
@@ -76,9 +78,13 @@ a converter to another format.
1. Add a `Format` to `FORMATS` in `conversion/__init__.py`, with a content
sniffer. The first sniffer to match wins, so a specific format
- (`cellbender`, which is also a 10x matrix) goes before a general one.
+ (`cellbender`, which is also a 10x matrix) goes before a general one, and
+ the folder formats (`visium`, `xenium`) go before `10x-mtx`, which would
+ otherwise claim their MEX subfolders.
2. Write the converter module; `tenx.py` is the model, and `features.py` holds
- the feature table every 10x-style format shares.
+ the feature table every 10x-style format shares. A format that pairs a 10x
+ matrix with a table (`visium.py`, `xenium.py`) uses `spatial.py` to match
+ the table to the matrix's barcodes and write it in their order.
3. Add fixtures shaped like the real tool's output to `tests/conftest.py`, and
a comparison against the reader people use in
`tests/test_external_readers.py`.
diff --git a/docs/guides/convert.md b/docs/guides/convert.md
index e882708..41c1153 100644
--- a/docs/guides/convert.md
+++ b/docs/guides/convert.md
@@ -18,8 +18,9 @@ to turn a loom into a 10x matrix, convert it to `.h5ad` first, then from that.
```
10x HDF5, 10x/STARsolo MTX ─┐ ┌─► 10x HDF5, 10x MTX
-CellBender, Parse, text ├─► .h5ad ◄──► .zarr ◄──┤
-loom, HDF5, pycisTopic │ ├─► loom (SCope too)
+Visium, Xenium, CellBender │ │
+Parse, text, HDF5 ├─► .h5ad ◄──► .zarr ◄──┤
+loom, pycisTopic │ ├─► loom (SCope too)
MuData ─────────────────────┘ └─► MuData
```
@@ -38,7 +39,9 @@ no extension.
| You have | Convert with | Read the page |
|---|---|---|
-| Cell Ranger, Space Ranger or Xenium `*_feature_bc_matrix.h5` | `adata convert in.h5 -o out.h5ad` | [10x HDF5](../reference/formats/10x-h5.md) |
+| 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 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) |
| A Parse Biosciences split-pipe run | `adata convert run -o out.h5ad` | [Parse](../reference/formats/parse.md) |
diff --git a/docs/index.md b/docs/index.md
index d16ae0d..f1932a0 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -68,10 +68,11 @@ Status goes to stderr and data to stdout, so every command pipes cleanly.
checks this against stores written by anndata 0.8 to 0.13, and checks that
each of those releases reads what `adata` writes.
- **Other formats, through `convert`**: 10x Genomics HDF5 and MTX folders
- (Cell Ranger, STARsolo, velocyto), CellBender output, Parse Biosciences
- runs, looms (velocyto's and SCope's), text tables, plain HDF5 matrices,
+ (Cell Ranger, STARsolo, velocyto), Visium and Xenium spatial outputs with
+ their coordinates and images, CellBender output, Parse Biosciences runs,
+ looms (velocyto's and SCope's), text tables, plain HDF5 matrices,
pycisTopic pickles and MuData. CI checks them against the readers people
- use for them: scanpy, loompy, anndata, mudata.
+ use for them: scanpy, squidpy, loompy, anndata, mudata.
- **Outputs that keep the source's choices**: Zarr version, index dtypes,
chunking and compression carry over unless a flag changes them.
- **Outputs that appear whole or not at all**: every store is written beside
diff --git a/docs/reference/commands.md b/docs/reference/commands.md
index d4427fd..830a2d3 100644
--- a/docs/reference/commands.md
+++ b/docs/reference/commands.md
@@ -119,6 +119,8 @@ adata convert filtered_feature_bc_matrix.h5 -o counts.h5ad
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 xenium_bundle -o cells.h5ad
adata convert splitpipe_out -o parse.h5ad
adata convert velocyto.loom -o velocyto.h5ad
adata convert counts.h5ad -o counts.loom
@@ -143,7 +145,9 @@ 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+, Space Ranger, Xenium | Cell Ranger 3+ layout |
+| [`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 |
+| [`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 |
| [`parse`](formats/parse.md) | a folder | Parse Biosciences split-pipe: a run, a sample or a `DGE_*` folder | no |
diff --git a/docs/reference/formats/10x-h5.md b/docs/reference/formats/10x-h5.md
index 3fb3ea5..4b7183d 100644
--- a/docs/reference/formats/10x-h5.md
+++ b/docs/reference/formats/10x-h5.md
@@ -9,7 +9,10 @@ nav_order: 1
Cell Ranger's `filtered_feature_bc_matrix.h5` and `raw_feature_bc_matrix.h5`,
and the same layout from Space Ranger, Xenium and CellBender. The file is
-recognised by its content, whatever its name.
+recognised by its content, whatever its name. Given alone, a Space Ranger or
+Xenium matrix converts without its positions: convert the
+[Space Ranger folder](visium.md) or the [Xenium bundle](xenium.md) to keep
+them.
```bash
adata convert filtered_feature_bc_matrix.h5 -o counts.h5ad
diff --git a/docs/reference/formats/cistopic.md b/docs/reference/formats/cistopic.md
index 04e4b0d..bd230db 100644
--- a/docs/reference/formats/cistopic.md
+++ b/docs/reference/formats/cistopic.md
@@ -2,7 +2,7 @@
title: pycisTopic pickles
parent: Formats
grand_parent: Reference
-nav_order: 8
+nav_order: 10
---
# pycisTopic pickles
diff --git a/docs/reference/formats/hdf5.md b/docs/reference/formats/hdf5.md
index 76cbd72..53d00a0 100644
--- a/docs/reference/formats/hdf5.md
+++ b/docs/reference/formats/hdf5.md
@@ -2,7 +2,7 @@
title: HDF5 matrices
parent: Formats
grand_parent: Reference
-nav_order: 7
+nav_order: 9
---
# HDF5 matrices
diff --git a/docs/reference/formats/index.md b/docs/reference/formats/index.md
index ed68071..8abc228 100644
--- a/docs/reference/formats/index.md
+++ b/docs/reference/formats/index.md
@@ -15,9 +15,11 @@ flags every conversion takes (`--output`, `--from`, `--to`, `--compression`,
| Format | `--from` / `--to` | Read | Write |
|---|---|---|---|
-| [10x Genomics HDF5](10x-h5.md) | `10x-h5` | Cell Ranger 2 and 3+, Space Ranger, Xenium | Cell Ranger 3+ layout |
+| [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 |
+| [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` |
| [Text tables](text.md) | `text` | a dense matrix in `.csv`, `.tsv`, `.txt` and the like, gzipped too | no |
diff --git a/docs/reference/formats/loom.md b/docs/reference/formats/loom.md
index 8e7a830..3c21311 100644
--- a/docs/reference/formats/loom.md
+++ b/docs/reference/formats/loom.md
@@ -2,7 +2,7 @@
title: loom and SCope
parent: Formats
grand_parent: Reference
-nav_order: 5
+nav_order: 7
---
# loom and SCope looms
diff --git a/docs/reference/formats/mudata.md b/docs/reference/formats/mudata.md
index 07042b0..60956c0 100644
--- a/docs/reference/formats/mudata.md
+++ b/docs/reference/formats/mudata.md
@@ -2,7 +2,7 @@
title: MuData
parent: Formats
grand_parent: Reference
-nav_order: 9
+nav_order: 11
---
# MuData
diff --git a/docs/reference/formats/parse.md b/docs/reference/formats/parse.md
index ff5b65e..c9405e3 100644
--- a/docs/reference/formats/parse.md
+++ b/docs/reference/formats/parse.md
@@ -2,7 +2,7 @@
title: Parse Biosciences
parent: Formats
grand_parent: Reference
-nav_order: 4
+nav_order: 6
---
# Parse Biosciences split-pipe
diff --git a/docs/reference/formats/text.md b/docs/reference/formats/text.md
index f9a260b..9cb8072 100644
--- a/docs/reference/formats/text.md
+++ b/docs/reference/formats/text.md
@@ -2,7 +2,7 @@
title: Text tables
parent: Formats
grand_parent: Reference
-nav_order: 6
+nav_order: 8
---
# Text tables
diff --git a/docs/reference/formats/visium.md b/docs/reference/formats/visium.md
new file mode 100644
index 0000000..5ab58ea
--- /dev/null
+++ b/docs/reference/formats/visium.md
@@ -0,0 +1,73 @@
+---
+title: Visium (Space Ranger)
+parent: Formats
+grand_parent: Reference
+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
+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.
+
+Give the `outs/` folder, or the run folder holding it. It is recognised by
+`spatial/` (the positions and `scalefactors_json.json`) beside a
+`*_feature_bc_matrix.h5`, including one named with its sample, as 10x's
+downloads are (`V1_Mouse_Brain_filtered_feature_bc_matrix.h5`).
+
+```bash
+adata convert outs -o section.h5ad
+adata convert outs -o every_spot.h5ad --unfiltered
+adata convert outs -o section.h5ad --library-id brain_A1 --no-images
+adata convert outs/filtered_feature_bc_matrix.h5 -o section.h5ad --from visium
+```
+
+| Space Ranger | AnnData |
+|---|---|
+| `filtered_feature_bc_matrix.h5` (`raw_` with `--unfiltered`) | `X`, var, as a [10x HDF5 file](10x-h5.md) is read |
+| `spatial/tissue_positions.csv` (`tissue_positions_list.csv`, no header, in Space Ranger 1): `in_tissue`, `array_row`, `array_col` | obs columns of those names, int64 |
+| the same: `pxl_col_in_fullres`, `pxl_row_in_fullres` | `obsm["spatial"]`, `[x, y]` in full-resolution pixels; int64, or float64 when the file has decimals |
+| `spatial/tissue_hires_image.png`, `tissue_lowres_image.png` | `uns["spatial"][library]["images"]["hires"]`, `["lowres"]`: the PNGs' pixels, uint8 |
+| `spatial/scalefactors_json.json` | `uns["spatial"][library]["scalefactors"]`, as it is |
+| the matrix's `chemistry_description`, `software_version` | `uns["spatial"][library]["metadata"]` |
+| the matrix's other run metadata | `uns["cellranger"]`, so an export to `.h5` keeps it |
+
+| Flag | Meaning |
+|---|---|
+| `--unfiltered` | Read `raw_feature_bc_matrix.h5`: every spot on the slide, not only those under tissue |
+| `--library-id` | The key under `uns["spatial"]` (default: the matrix's library id, as scanpy and squidpy use) |
+| `--no-images` | Leave the images out; the scale factors and metadata stay |
+
+The [10x HDF5 flags](10x-h5.md) (`--genome`, `--var-names`,
+`--make-unique`, `--feature-type`, `--split-feature-types`) apply as well;
+with `--split-feature-types` every output gets the positions and images.
+
+## Differences from squidpy
+
+- Features are named by gene id, with the symbols in `var["gene_symbols"]`,
+ and every feature type is kept; squidpy names them by symbol and keeps
+ Gene Expression only. `--var-names gene_symbols --feature-type "Gene
+ Expression"` gives squidpy's.
+- `X` keeps Space Ranger's integer counts, where squidpy (through scanpy)
+ makes them float32.
+
+## How it is read
+
+obs follows the matrix's barcode order. The positions list every spot on
+the slide, in their own order, so they are matched to the matrix's barcodes
+as they are scanned and written in its order, a block at a time; they are a
+few thousand rows (some 14,000 on an 11 mm slide) however deep the
+sequencing. A barcode with no position, or two, is refused before anything
+is written. The images are read whole: a hires PNG is 2,000 pixels on its
+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`.
diff --git a/docs/reference/formats/xenium.md b/docs/reference/formats/xenium.md
new file mode 100644
index 0000000..d31248d
--- /dev/null
+++ b/docs/reference/formats/xenium.md
@@ -0,0 +1,60 @@
+---
+title: Xenium
+parent: Formats
+grand_parent: Reference
+nav_order: 5
+---
+
+# Xenium
+
+A Xenium Onboard Analysis output bundle (XOA 1 to 4): the cell x feature
+counts, the table of cells with their centroids, and the clustering, UMAP
+and PCA Xenium computed. Give the bundle folder; it is recognised by
+`experiment.xenium` beside `cell_feature_matrix.h5`.
+
+```bash
+adata convert output-XETG00000__0001234__Region_1 -o region1.h5ad
+adata convert xenium_bundle -o region1.h5ad --feature-type "Gene Expression"
+```
+
+| Xenium | AnnData |
+|---|---|
+| `cell_feature_matrix.h5` | `X`, var, as a [10x HDF5 file](10x-h5.md) is read: Gene Expression and the control probes and codewords, each a feature type |
+| `cells.csv.gz`: `cell_id` | the obs names (the matrix's barcodes) |
+| `cells.csv.gz`: `x_centroid`, `y_centroid` | `obsm["spatial"]`, `[x, y]` in microns |
+| `cells.csv.gz`: every other column (`transcript_counts`, `cell_area`, `nucleus_area`, `segmentation_method`, ...) | obs columns, typed as pandas reads them; text as a categorical when it repeats a value, as anndata writes it |
+| `analysis/clustering/gene_expression_graphclust/clusters.csv` | `obs["graphclust"]`, a categorical |
+| `analysis/umap/gene_expression_2_components/projection.csv` | `obsm["X_umap"]` |
+| `analysis/pca/gene_expression_10_components/projection.csv` | `obsm["X_pca"]` |
+| the matrix's run metadata | `uns["cellranger"]` |
+
+Xenium ships `analysis/` packed as `analysis.tar.gz`. When it has not been
+unpacked, the three tables above are copied out of it, by name, into a
+temporary folder beside the output, removed when the conversion ends. A cell an
+analysis table lacks is left empty (no cluster, NaN), with a warning
+giving the count.
+
+The [10x HDF5 flags](10x-h5.md) (`--var-names`, `--make-unique`,
+`--feature-type`, `--split-feature-types`) apply. Cell boundaries,
+transcripts and the morphology images are not read.
+
+## Differences from spatialdata-io
+
+spatialdata-io's `xenium` reader returns a SpatialData, whose table holds
+the same cells and `obsm["spatial"]`. Its obs also has `cell_id` as a
+column, which here is the obs names, and a `region` column naming the
+SpatialData element the table annotates, which an AnnData on its own has
+no use for. CI checks the import against 10x's own recipe for scanpy:
+`read_10x_h5`, the cell table reindexed to it, the analysis tables.
+
+## How it is read
+
+The cell table can hold a million rows. It is matched to the matrix's
+barcodes while it is scanned, and Xenium writes the two in the same order,
+so it is then written straight through, a block at a time. A table in
+another order is written in the matrix's order through scratch space beside
+the output; cells the matrix lacks are left out, with a warning. A cell of
+the matrix with no row in the table is refused before anything is written.
+
+A bundle with `cells.parquet` and no `cells.csv.gz` is refused: Parquet is
+not read yet.
diff --git a/docs/troubleshooting.md b/docs/troubleshooting.md
index 25e3c44..f97e011 100644
--- a/docs/troubleshooting.md
+++ b/docs/troubleshooting.md
@@ -134,6 +134,53 @@ The input belongs to another store, or is transposed.
`subset` with no selection would copy the store. Use `convert` to copy, or add
a selection.
+## Visium and Xenium
+
+These are refused before anything is written too.
+
+### `12 of the matrix's 2,702 barcodes have no row in spatial/tissue_positions.csv (the first is 'AAAC...-1')`
+
+The positions and the matrix come from different runs, or the positions
+were edited. Every spot of the matrix needs a position; spots in the
+positions but not the matrix are fine. A Xenium cell missing from
+`cells.csv.gz` is refused the same way, and so is a table (positions, cells
+or an analysis table) with a row for none of the matrix's barcodes.
+
+### `tissue_positions.csv: pxl_col_in_fullres holds values that are not numbers.`
+
+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.`
+
+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`.
+
+### `This matrix holds 2 libraries (A1, B1), as spaceranger aggr writes`
+
+An aggregated matrix's spots come from several slides, with no one image to
+place them on. Convert each library's own `outs/`, and `adata concat` the
+results if you need them together.
+
+### `... has filtered_feature_bc_matrix/ as a Matrix Market folder but no filtered_feature_bc_matrix.h5`
+
+The positions are read only beside the `.h5` matrix. Convert the folder with
+`--from 10x-mtx` for the counts alone. A Xenium bundle with only
+`cell_feature_matrix/` says the same.
+
+### `... has no cells.csv.gz. It has cells.parquet only, which is not read yet.`
+
+Xenium writes `cells.csv.gz` beside `cells.parquet` in every version so far;
+a bundle trimmed to Parquet cannot be read yet.
+
+### `--unfiltered does not apply when converting xenium to h5ad.`
+
+Xenium has one matrix. `--unfiltered`, `--library-id` and `--no-images` are
+Visium's.
+
## Leftovers from an interrupted run
### `'out.adata-tmp.h5ad' already exists, probably left by an interrupted run. Check what it holds, remove it, and try again.`
diff --git a/src/adata/cli.py b/src/adata/cli.py
index e100595..71e3dbb 100644
--- a/src/adata/cli.py
+++ b/src/adata/cli.py
@@ -528,6 +528,8 @@ def concat(
"sample": "--sample",
"all_samples": "--all-samples",
"unfiltered": "--unfiltered",
+ "library_id": "--library-id",
+ "no_images": "--no-images",
"dense": "--dense",
"obs_names_attr": "--obs-names-attr",
"var_names_attr": "--var-names-attr",
@@ -695,7 +697,20 @@ def convert(
unfiltered: bool = typer.Option(
False,
"--unfiltered",
- help="Parse input: read DGE_unfiltered (every barcode) not DGE_filtered",
+ help="Parse or Visium input: read every barcode (DGE_unfiltered, "
+ "raw_feature_bc_matrix.h5), not only cells",
+ ),
+ library_id: Optional[str] = typer.Option(
+ None,
+ "--library-id",
+ help="Visium input: the key under uns['spatial'] (default: the "
+ "matrix's library id, as scanpy and squidpy use)",
+ ),
+ no_images: bool = typer.Option(
+ False,
+ "--no-images",
+ help="Visium input: leave the tissue images out of uns['spatial'] "
+ "(the scale factors stay)",
),
dense: bool = typer.Option(
False,
@@ -793,7 +808,10 @@ def convert(
target can express them. CellBender output keeps its latent variables, as
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. MTX folders from Cell
+ 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,
+ 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,
gzipped too), read as anndata's read_csv and read_text read them; a plain
@@ -821,6 +839,8 @@ def convert(
adata convert counts.h5ad -o for_cellbender.h5 --layer counts
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 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
adata convert cellbender_output.h5 -o cells.h5ad --min-cell-probability 0.5
@@ -924,6 +944,8 @@ def convert(
sample=sample,
all_samples=all_samples,
unfiltered=unfiltered,
+ library_id=library_id,
+ no_images=no_images,
dense=dense,
obs_names_attr=obs_names_attr,
var_names_attr=var_names_attr,
diff --git a/src/adata/conversion/__init__.py b/src/adata/conversion/__init__.py
index 36e3a19..0b3dee2 100644
--- a/src/adata/conversion/__init__.py
+++ b/src/adata/conversion/__init__.py
@@ -17,7 +17,9 @@
import h5py
-from adata.conversion import cellbender, cistopic, hdf5, loom, mtx, mudata, parse, tenx, text
+from adata.conversion import (
+ cellbender, cistopic, hdf5, loom, mtx, mudata, parse, tenx, text, visium, xenium,
+)
from adata.formats import pickled
from adata.conversion.anndata_copy import AnnDataCopy
from adata.conversion.base import ConversionPlan, Converter, Format, Options
@@ -103,6 +105,23 @@ def _sniff_zarr(path: Path) -> bool:
importer=tenx.TenxImport(),
exporter=tenx.TenxExport(),
),
+ # Folders, and before 10x-mtx: a Space Ranger outs/ or a Xenium bundle
+ # holds its matrix as MEX folders too, which 10x-mtx would take for a
+ # Cell Ranger 2 genome layout. A bare .h5 is no folder, so it stays 10x-h5.
+ "visium": Format(
+ name="visium",
+ description="Space Ranger output for standard Visium: matrix, spot positions, images",
+ extensions=(),
+ sniff=visium.sniff,
+ importer=visium.VisiumImport(),
+ ),
+ "xenium": Format(
+ name="xenium",
+ description="Xenium output bundle: cell matrix, cell table and centroids, Xenium's analysis",
+ extensions=(),
+ sniff=xenium.sniff,
+ importer=xenium.XeniumImport(),
+ ),
"loom": Format(
name="loom",
description="Loom (loompy, velocyto, kallisto|bustools, SCope; spec 2 and 3)",
diff --git a/src/adata/conversion/base.py b/src/adata/conversion/base.py
index 6dec9db..9931fdd 100644
--- a/src/adata/conversion/base.py
+++ b/src/adata/conversion/base.py
@@ -39,11 +39,15 @@ class Options:
raw: Optional[Path] = None
min_cell_probability: Optional[float] = None
- # Reading Parse split-pipe output.
+ # Reading Parse split-pipe output (and, --unfiltered, Space Ranger's).
sample: Optional[str] = None
all_samples: bool = False
unfiltered: bool = False
+ # Reading Space Ranger output.
+ library_id: Optional[str] = None
+ no_images: bool = False
+
# Reading a loom.
dense: bool = False
obs_names_attr: Optional[str] = None
diff --git a/src/adata/conversion/spatial.py b/src/adata/conversion/spatial.py
new file mode 100644
index 0000000..93ceff5
--- /dev/null
+++ b/src/adata/conversion/spatial.py
@@ -0,0 +1,187 @@
+"""What the spatial formats share: tables matched to a matrix, and images.
+
+Space Ranger's spot positions and Xenium's cell table each list their rows
+in their own order, beside a 10x matrix whose barcode order AnnData keeps.
+`BarcodeJoin` matches a table to the matrix while the table is scanned, and
+`write_rows` writes it in the matrix's order: straight through when the two
+already agree, as Xenium's do, and through `core.scatter` when not.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from pathlib import Path
+from typing import Any, Callable, Dict, List, Optional, Sequence
+
+import numpy as np
+
+from adata.conversion.features import STRING_STEP
+from adata.core.scatter import RowScatter
+from adata.elements.read import decode_str_array
+from adata.formats.table import TableScan, column_values, iter_table_blocks
+
+
+@dataclass
+class Join:
+ """Where each row of a table goes among the matrix's barcodes."""
+
+ n_out: int
+ #: The output row of each table row, -1 for one the matrix lacks; None
+ #: when the table lists the matrix's barcodes exactly, in its order.
+ target: Optional[np.ndarray]
+ #: Table rows the matrix has no barcode for.
+ n_extra: int = 0
+ #: Matrix barcodes the table has no row for.
+ n_missing: int = 0
+ #: Matrix barcodes the table has more than one row for.
+ n_repeated: int = 0
+ example_missing: str = ""
+
+ def check(self, what: str, *, allow_missing: bool = False) -> None:
+ """Refuse a join that would leave a barcode with no row, or two."""
+ if self.n_repeated:
+ raise ValueError(
+ f"{what} has more than one row for {self.n_repeated:,} of the "
+ "matrix's barcodes."
+ )
+ if self.n_out and self.n_missing == self.n_out:
+ raise ValueError(
+ f"{what} has a row for none of the matrix's {self.n_out:,} barcodes; "
+ "it must come from the same run as the matrix."
+ )
+ if self.n_missing and not allow_missing:
+ raise ValueError(
+ f"{self.n_missing:,} of the matrix's {self.n_out:,} barcodes have "
+ f"no row in {what} (the first is {self.example_missing!r}); "
+ "it must come from the same run as the matrix."
+ )
+
+
+class BarcodeJoin:
+ """Matches a table's index to a matrix's barcodes, a block at a time.
+
+ Given to `scan_table` as `on_index`, so the match costs no pass of its
+ own. While each block equals the matrix's next barcodes it only compares
+ them; at the first difference it switches to looking each value up in a
+ sorted copy of the barcodes, a numpy array rather than a dict, whose
+ per-string overhead is several times the barcode itself.
+ """
+
+ def __init__(self, barcodes: Any) -> None:
+ self._barcodes = barcodes
+ self.n_out = int(barcodes.shape[0])
+ self._in_order = True
+ self._done = 0
+ self._targets: List[np.ndarray] = []
+ self._ranked: Optional[np.ndarray] = None
+ self._order: Optional[np.ndarray] = None
+
+ def _slice(self, lo: int, hi: int) -> np.ndarray:
+ return np.asarray(decode_str_array(np.asarray(self._barcodes[lo:hi])), dtype=str)
+
+ def __call__(self, values: Sequence[str]) -> None:
+ values = np.asarray(values, dtype=str)
+ lo, hi = self._done, self._done + values.size
+ if self._in_order:
+ if hi <= self.n_out and np.array_equal(self._slice(lo, hi), values):
+ self._done = hi
+ return
+ self._leave_order()
+ self._targets.append(self._lookup(values))
+ self._done = hi
+
+ def _leave_order(self) -> None:
+ # Every row so far matched the barcode at its own position.
+ self._in_order = False
+ self._targets = [np.arange(self._done, dtype=np.int64)]
+ wanted = np.concatenate(
+ [self._slice(lo, min(lo + STRING_STEP, self.n_out)) for lo in range(0, self.n_out, STRING_STEP)]
+ or [np.zeros(0, dtype=str)]
+ )
+ self._order = np.argsort(wanted, kind="stable")
+ self._ranked = wanted[self._order]
+
+ def _lookup(self, values: np.ndarray) -> np.ndarray:
+ out = np.full(values.size, -1, dtype=np.int64)
+ if self._ranked.size == 0:
+ return out
+ at = np.minimum(np.searchsorted(self._ranked, values), self._ranked.size - 1)
+ hit = self._ranked[at] == values
+ out[hit] = self._order[at[hit]]
+ return out
+
+ def finish(self) -> Join:
+ if self._in_order and self._done == self.n_out:
+ return Join(n_out=self.n_out, target=None)
+ if self._in_order:
+ self._leave_order()
+ target = np.concatenate(self._targets) if self._targets else np.zeros(0, dtype=np.int64)
+ counts = np.bincount(target[target >= 0], minlength=self.n_out)
+ missing = np.flatnonzero(counts == 0)
+ return Join(
+ n_out=self.n_out,
+ target=target,
+ n_extra=int((target < 0).sum()),
+ n_missing=int(missing.size),
+ n_repeated=int((counts > 1).sum()),
+ example_missing=str(self._slice(missing[0], missing[0] + 1)[0]) if missing.size else "",
+ )
+
+
+@dataclass
+class Sink:
+ """One array to fill from a table: `values` turns a block's columns
+ into the block's values; `fill` is what rows the table lacks get."""
+
+ dest: Any
+ values: Callable[[Dict[str, List[str]]], np.ndarray]
+ fill: Any = None
+
+
+def write_rows(scan: TableScan, join: Join, sinks: Sequence[Sink], scratch_parent: Any) -> None:
+ """Fill each sink from one read of the table, in the matrix's order."""
+ if join.target is None:
+ lo = 0
+ for names, block in iter_table_blocks(scan):
+ hi = lo + len(names)
+ for sink in sinks:
+ sink.dest[lo:hi] = sink.values(block)
+ lo = hi
+ return
+ with RowScatter(join.target, join.n_out, scratch_parent) as scatter:
+ columns = [scatter.column(sink.dest, sink.fill) for sink in sinks]
+ for _, block in iter_table_blocks(scan):
+ for column, sink in zip(columns, sinks):
+ column.add(sink.values(block))
+
+
+def require_numeric(scan: TableScan, names: Sequence[str]) -> None:
+ """Refuse, at plan time, a coordinate or projection column that is not
+ all numbers (missing values are allowed, and become NaN)."""
+ kinds = {c.name: scan.dtype_of(c) for c in scan.columns}
+ absent = [n for n in names if n not in kinds]
+ if absent:
+ raise ValueError(f"{scan.path.name} has no {', '.join(absent)} column.")
+ bad = [n for n in names if kinds[n] not in ("int64", "float64")]
+ if bad:
+ raise ValueError(
+ f"{scan.path.name}: {', '.join(bad)} {'holds' if len(bad) == 1 else 'hold'} "
+ "values that are not numbers."
+ )
+
+
+def numbers(scan: TableScan, names: Sequence[str], block: Dict[str, List[str]], dtype: Any = np.float64) -> np.ndarray:
+ """A block's values of numeric columns, side by side, NaN where missing."""
+ by_name = {c.name: c for c in scan.columns}
+ return np.stack(
+ [column_values(scan, by_name[n], block[n], {}).astype(dtype) for n in names], axis=1
+ )
+
+
+def read_image(path: Path) -> np.ndarray:
+ """An image's pixels as the file holds them: uint8 for an 8-bit PNG, as
+ squidpy reads them, H x W or H x W x C."""
+ from PIL import Image
+
+ with Image.open(path) as image:
+ return np.asarray(image)
diff --git a/src/adata/conversion/tenx.py b/src/adata/conversion/tenx.py
index 7f370cc..3f0d6b0 100644
--- a/src/adata/conversion/tenx.py
+++ b/src/adata/conversion/tenx.py
@@ -90,6 +90,9 @@
"original_gem_groups",
"software_version",
"version",
+ # Xenium 4's protein values are scaled by this; spatialdata-io divides.
+ "protein_scaling_factor",
+ "protein_metric",
)
diff --git a/src/adata/conversion/visium.py b/src/adata/conversion/visium.py
new file mode 100644
index 0000000..5fcbebb
--- /dev/null
+++ b/src/adata/conversion/visium.py
@@ -0,0 +1,352 @@
+"""Space Ranger output for standard Visium: counts, spot positions, images.
+
+A Space Ranger `outs/` folder holds the 10x matrix and, in `spatial/`, where
+each spot sits on the slide (`tissue_positions.csv`, or Space Ranger 1's
+`tissue_positions_list.csv` with no header), the scale factors between the
+full-resolution image and the two PNGs beside them, and the PNGs.
+
+The AnnData produced is laid out as `squidpy.read.visium` lays it out, which
+scanpy and squidpy plot from:
+
+* obs: `in_tissue`, `array_row`, `array_col`;
+* `obsm["spatial"]`: each spot's full-resolution pixel, `[col, row]` (x, y);
+* `uns["spatial"][library_id]`: `images` (`hires`, `lowres`, the PNGs'
+ uint8 pixels), `scalefactors` (the JSON as it is) and `metadata`
+ (`chemistry_description`, `software_version`).
+
+The matrix is read as a 10x HDF5 file is (`tenx`), so features are named by
+gene id and every feature type is kept, where squidpy names them by symbol
+and keeps Gene Expression only. The positions list every spot on the slide,
+a few thousand however large the data, and are matched to the matrix's
+barcodes as they are scanned.
+
+Visium HD is not read yet: its positions are Parquet.
+"""
+
+from __future__ import annotations
+
+import json
+from contextlib import ExitStack
+from dataclasses import dataclass
+from pathlib import Path
+from typing import Any, Dict, List, Optional, Sequence
+
+import h5py
+import numpy as np
+from rich.console import Console
+
+from adata.conversion import tenx
+from adata.conversion.base import ConversionPlan, Converter, Options
+from adata.conversion.features import STRING_STEP, Target, describe_import, write_obs_index, write_var
+from adata.conversion.spatial import BarcodeJoin, Join, Sink, numbers, read_image, require_numeric, write_rows
+from adata.core.subset import fan_out_sparse_matrix
+from adata.elements.read import decode_str_array
+from adata.elements.write import (
+ ensure_anndata_skeleton,
+ set_column_order,
+ write_dense,
+ write_mapping,
+ write_scalar,
+)
+from adata.formats.json_data import write_json_value
+from adata.formats.table import TableScan, column_values, create_columns, scan_table
+from adata.formats.textio import open_binary
+from adata.storage import create_dataset, is_zarr_path, open_store
+from adata.elements import spec
+
+#: The columns of `tissue_positions.csv`, which Space Ranger 1 wrote without
+#: a header, in this order.
+POSITION_COLUMNS = (
+ "barcode", "in_tissue", "array_row", "array_col",
+ "pxl_row_in_fullres", "pxl_col_in_fullres",
+)
+_POSITION_FILES = ("tissue_positions.csv", "tissue_positions_list.csv")
+_IMAGES = ("hires", "lowres")
+_FILTERED, _RAW = "filtered_feature_bc_matrix", "raw_feature_bc_matrix"
+
+
+def _is_hd(folder: Path) -> bool:
+ return (folder / "binned_outputs").is_dir() or (
+ folder / "spatial" / "tissue_positions.parquet"
+ ).exists()
+
+
+def _matrix_files(folder: Path, name: str) -> List[Path]:
+ """`.h5` in the folder, or one prefixed with the sample name, as
+ 10x's downloads are (`V1_Mouse_Brain_filtered_feature_bc_matrix.h5`)."""
+ return sorted(
+ p for p in folder.glob(f"*{name}.h5")
+ if p.name == f"{name}.h5" or p.name.endswith(f"_{name}.h5")
+ )
+
+
+def _locate(path: Path) -> Optional[Path]:
+ """The Space Ranger output folder at `path` or in its `outs/`."""
+ for folder in (path, path / "outs"):
+ if not folder.is_dir():
+ continue
+ if _is_hd(folder):
+ return folder
+ spatial = folder / "spatial"
+ if (spatial / "scalefactors_json.json").is_file() and any(
+ (spatial / name).is_file() for name in _POSITION_FILES
+ ) and (
+ _matrix_files(folder, _FILTERED) or _matrix_files(folder, _RAW)
+ or (folder / _FILTERED).is_dir() or (folder / _RAW).is_dir()
+ ):
+ return folder
+ return None
+
+
+def sniff(path: Path) -> bool:
+ """A Space Ranger output folder (or the run folder holding its `outs/`).
+
+ Visium HD folders are claimed too, only so that they are refused by name
+ rather than taken for a Cell Ranger 2 MTX layout.
+ """
+ path = Path(path)
+ return path.is_dir() and not is_zarr_path(path) and _locate(path) is not None
+
+
+@dataclass
+class _Plan:
+ matrix: Any
+ matrix_path: Path
+ targets: List[Target]
+ positions: TableScan
+ join: Join
+ library_id: str
+ scalefactors: Dict[str, Any]
+ images: Dict[str, Path]
+ filtered_barcodes: bool
+ zarr_format: Optional[int]
+
+
+def _bundle(source: Path, options: Options) -> tuple:
+ """(the outs folder, the matrix file) for this input."""
+ if source.is_file():
+ if "--unfiltered" in options.given:
+ raise ValueError(
+ "--unfiltered picks between the matrices of a Space Ranger outs "
+ f"folder; {source.name} is a matrix already. Give the folder."
+ )
+ return source.parent, source
+ outs = _locate(source)
+ if outs is None:
+ raise ValueError(
+ f"{source} is not a Space Ranger output folder: it needs spatial/ "
+ "(tissue positions and scalefactors_json.json) beside a "
+ "*_feature_bc_matrix.h5."
+ )
+ if _is_hd(outs):
+ raise ValueError(
+ "This is Visium HD output, which is not read yet: its positions are "
+ "Parquet. Convert one bin's matrix as a 10x matrix in the meantime "
+ "(binned_outputs/square_008um/filtered_feature_bc_matrix.h5)."
+ )
+ name = _RAW if options.unfiltered else _FILTERED
+ found = _matrix_files(outs, name)
+ if len(found) > 1:
+ raise ValueError(
+ f"{outs} holds several {name}.h5 files ({', '.join(p.name for p in found)}); "
+ "give the one to read with --from visium."
+ )
+ if not found:
+ if (outs / name).is_dir():
+ raise ValueError(
+ f"{outs} has {name}/ as a Matrix Market folder but no {name}.h5, "
+ f"which this reads. Convert the folder with --from 10x-mtx, "
+ "without the positions and images."
+ )
+ raise ValueError(f"{outs} has no {name}.h5.")
+ return outs, found[0]
+
+
+def _has_header(path: Path) -> bool:
+ with open_binary(path) as fh:
+ first = fh.readline().decode("utf-8-sig", "replace")
+ return first.split(",", 1)[0].strip() == "barcode"
+
+
+def _library_id(attrs: Dict[str, Any], outs: Path, options: Options) -> str:
+ if options.library_id is not None:
+ if not options.library_id or "/" in options.library_id:
+ raise ValueError(
+ f"--library-id {options.library_id!r} cannot be a key of uns: it "
+ "must be non-empty, with no '/'."
+ )
+ return options.library_id
+ ids = [
+ v.decode("utf-8", "replace") if isinstance(v, bytes) else str(v)
+ for v in np.atleast_1d(attrs.get("library_ids", []))
+ ]
+ if len(ids) > 1:
+ raise ValueError(
+ f"This matrix holds {len(ids)} libraries ({', '.join(ids)}), as "
+ "spaceranger aggr writes; their spots share no one image. Convert "
+ "each library's own output."
+ )
+ if ids:
+ return ids[0]
+ return outs.parent.name if outs.name == "outs" else outs.name
+
+
+class VisiumImport(Converter):
+ options = tenx.TenxImport.options | frozenset({"--unfiltered", "--library-id", "--no-images"})
+
+ def plan(
+ self,
+ source: Path,
+ output: Path,
+ src_format: str,
+ dst_format: str,
+ options: Options,
+ ) -> ConversionPlan:
+ source = Path(source)
+ outs, matrix_path = _bundle(source, options)
+ spatial = outs / "spatial"
+ position_file = next((spatial / n for n in _POSITION_FILES if (spatial / n).is_file()), None)
+ scalefactor_file = spatial / "scalefactors_json.json"
+ if position_file is None or not scalefactor_file.is_file():
+ raise ValueError(
+ f"{outs} has no spatial/tissue_positions.csv and "
+ "spatial/scalefactors_json.json beside the matrix."
+ )
+
+ with h5py.File(matrix_path, "r") as handle:
+ if not tenx.sniff(matrix_path):
+ raise ValueError(f"{matrix_path.name} is not a 10x HDF5 matrix.")
+ matrix, outputs, targets, warnings = tenx.plan_matrix_import(handle, output, options)
+ library_id = _library_id(matrix.attrs, outs, options)
+ group = handle[matrix.group]
+ joiner = BarcodeJoin(group["barcodes"])
+ positions = scan_table(
+ position_file,
+ index="barcode",
+ names=None if _has_header(position_file) else POSITION_COLUMNS,
+ on_index=joiner,
+ )
+ filtered_barcodes = "filtered_barcodes" in group
+ join = joiner.finish()
+ require_numeric(positions, POSITION_COLUMNS[1:])
+ join.check(f"spatial/{position_file.name}")
+
+ scalefactors = json.loads(scalefactor_file.read_text())
+ images: Dict[str, Path] = {}
+ if not options.no_images:
+ for name in _IMAGES:
+ path = spatial / f"tissue_{name}_image.png"
+ if path.is_file():
+ images[name] = path
+ else:
+ warnings.append(f"No spatial/{path.name}; uns['spatial'] will have no {name} image.")
+
+ described = describe_import(matrix.n_barcodes, matrix.n_features, matrix.feature_types, targets)
+ return ConversionPlan(
+ source=source,
+ src_format=src_format,
+ dst_format=dst_format,
+ outputs=outputs,
+ summary=(
+ f"Space Ranger (Visium) → {dst_format}, {described}, positions from "
+ f"{position_file.name}, library {library_id!r}"
+ + (f", images {', '.join(images)}" if images else "")
+ ),
+ warnings=warnings,
+ detail=_Plan(
+ matrix=matrix,
+ matrix_path=matrix_path,
+ targets=targets,
+ positions=positions,
+ join=join,
+ library_id=library_id,
+ scalefactors=scalefactors,
+ images=images,
+ filtered_barcodes=filtered_barcodes,
+ zarr_format=(options.zarr_format or 3) if dst_format == "zarr" else None,
+ ),
+ )
+
+ def write(self, plan: ConversionPlan, staged: Sequence[Path], console: Console) -> None:
+ detail: _Plan = plan.detail
+ matrix = detail.matrix
+ images = {name: read_image(path) for name, path in detail.images.items()}
+
+ with h5py.File(detail.matrix_path, "r") as handle, ExitStack() as stack:
+ group = handle[matrix.group]
+ roots = [
+ stack.enter_context(open_store(path, "w", zarr_format=detail.zarr_format)).root
+ for path in staged
+ ]
+ for root, target in zip(roots, detail.targets):
+ with console.status("[cyan]Writing obs, var and the spot positions...[/]"):
+ _write_obs(root, group, detail)
+ write_var(root, target)
+ tenx._write_uns(root, matrix)
+ _write_spatial_uns(root, detail, images)
+
+ with console.status("[cyan]Copying X...[/]"):
+ fan_out_sparse_matrix(
+ tenx._AsCsr(group, (matrix.n_barcodes, matrix.n_features)),
+ "X",
+ [(root, None, target.keep) for root, target in zip(roots, detail.targets)],
+ )
+
+
+def _write_obs(root: Any, group: h5py.Group, detail: _Plan) -> None:
+ """obs (the barcodes, and the positions' columns), and `obsm["spatial"]`."""
+ barcodes = group["barcodes"]
+ n = barcodes.shape[0]
+ obs = write_obs_index(
+ root, n,
+ (
+ [str(v) for v in decode_str_array(np.asarray(barcodes[lo : lo + STRING_STEP]))]
+ for lo in range(0, n, STRING_STEP)
+ ),
+ )
+ ensure_anndata_skeleton(root)
+
+ scan = detail.positions
+ pixels = ("pxl_col_in_fullres", "pxl_row_in_fullres") # x, y
+ kinds = {c.name: scan.dtype_of(c) for c in scan.columns}
+ columns, lookups = create_columns(obs, scan, n, exclude=pixels)
+ dtype = np.int64 if all(kinds[p] == "int64" for p in pixels) else np.float64
+ spatial = create_dataset(root["obsm"], "spatial", shape=(n, 2), dtype=dtype)
+ spec.set_encoding(spatial, spec.ARRAY)
+
+ by_name = {c.name: c for c in scan.columns}
+ sinks = [
+ Sink(
+ dest,
+ lambda block, column=by_name[name]: column_values(
+ scan, column, block[column.name], lookups.get(column.name, {})
+ ),
+ )
+ for name, dest in columns.items()
+ ]
+ sinks.append(Sink(spatial, lambda block: numbers(scan, pixels, block, dtype)))
+ write_rows(scan, detail.join, sinks, root)
+
+ order = list(columns)
+ if detail.filtered_barcodes:
+ write_dense(obs, "filtered_barcodes", np.asarray(group["filtered_barcodes"][...]).astype(bool))
+ order.append("filtered_barcodes")
+ set_column_order(obs, order)
+
+
+def _write_spatial_uns(root: Any, detail: _Plan, images: Dict[str, np.ndarray]) -> None:
+ """`uns["spatial"][library_id]`, as squidpy writes it."""
+ library = write_mapping(write_mapping(root["uns"], "spatial"), detail.library_id)
+ if images:
+ holder = write_mapping(library, "images")
+ for name, pixels in images.items():
+ write_dense(holder, name, pixels)
+ write_json_value(library, "scalefactors", detail.scalefactors)
+ metadata = write_mapping(library, "metadata")
+ for key in ("chemistry_description", "software_version"):
+ value = detail.matrix.attrs.get(key)
+ if value is not None:
+ write_scalar(
+ metadata, key,
+ value.decode("utf-8", "replace") if isinstance(value, (bytes, np.bytes_)) else str(value),
+ )
diff --git a/src/adata/conversion/xenium.py b/src/adata/conversion/xenium.py
new file mode 100644
index 0000000..83abbd8
--- /dev/null
+++ b/src/adata/conversion/xenium.py
@@ -0,0 +1,352 @@
+"""Xenium Onboard Analysis output: counts, cells, and Xenium's own analysis.
+
+A Xenium output bundle holds the cells x features matrix
+(`cell_feature_matrix.h5`, a 10x matrix), a table of the cells
+(`cells.csv.gz`: centroid, counts per kind of probe, cell and nucleus area),
+and, in `analysis/`, the clustering, UMAP and PCA Xenium Explorer shows. It
+is recognised by `experiment.xenium`.
+
+The AnnData produced:
+
+* obs: the cell table's columns, typed as pandas reads them, indexed by
+ `cell_id` (the matrix's barcodes); `graphclust`, Xenium's graph-based
+ clusters, as a categorical;
+* `obsm["spatial"]`: each cell's centroid, `[x, y]` in microns, as
+ spatialdata-io places it; `obsm["X_umap"]` and `obsm["X_pca"]` from
+ `analysis/`.
+
+Xenium ships `analysis/` packed as `analysis.tar.gz`; the three tables used
+are read from it when it has not been unpacked.
+
+The cell table can hold a million rows. It is matched to the matrix's
+barcodes while it is scanned, and Xenium writes both in the same order, so
+it is then written straight through; a table in another order goes through
+`core.scatter`. Boundaries, transcripts and morphology images are not read.
+"""
+
+from __future__ import annotations
+
+import shutil
+import tarfile
+import tempfile
+from contextlib import ExitStack
+from dataclasses import dataclass, field
+from pathlib import Path
+from typing import Any, Dict, List, Optional, Sequence
+
+import h5py
+import numpy as np
+from rich.console import Console
+
+from adata.conversion import tenx
+from adata.conversion.base import ConversionPlan, Converter, Options
+from adata.conversion.features import STRING_STEP, Target, describe_import, write_obs_index, write_var
+from adata.conversion.spatial import BarcodeJoin, Join, Sink, numbers, require_numeric, write_rows
+from adata.core.subset import fan_out_sparse_matrix
+from adata.elements import spec
+from adata.elements.read import decode_str_array
+from adata.elements.write import create_categorical, ensure_anndata_skeleton, set_column_order
+from adata.formats.table import TableScan, column_values, create_columns, drop_unused, scan_table
+from adata.storage import create_dataset, is_zarr_path, open_store
+from adata.util.names import natural_sorted
+
+_CENTROIDS = ("x_centroid", "y_centroid")
+_CELL_TABLES = ("cells.csv.gz", "cells.csv")
+
+#: Xenium's analysis outputs: (path under analysis/, where it goes, value columns).
+_CLUSTERS = Path("clustering/gene_expression_graphclust/clusters.csv")
+_PROJECTIONS = {
+ "X_umap": Path("umap/gene_expression_2_components/projection.csv"),
+ "X_pca": Path("pca/gene_expression_10_components/projection.csv"),
+}
+
+
+def sniff(path: Path) -> bool:
+ """A Xenium output bundle: `experiment.xenium` beside the cell matrix."""
+ path = Path(path)
+ return (
+ path.is_dir()
+ and not is_zarr_path(path)
+ and (path / "experiment.xenium").is_file()
+ and ((path / "cell_feature_matrix.h5").is_file() or (path / "cell_feature_matrix").is_dir())
+ )
+
+
+@dataclass
+class _Table:
+ """A table beside the matrix, scanned and matched to its barcodes."""
+
+ scan: TableScan
+ join: Join
+
+
+@dataclass
+class _Plan:
+ matrix: Any
+ matrix_path: Path
+ targets: List[Target]
+ cells: _Table
+ clusters: Optional[_Table]
+ projections: Dict[str, _Table] = field(default_factory=dict)
+ zarr_format: Optional[int] = None
+ #: Where tables taken out of analysis.tar.gz wait to be written; removed
+ #: when the plan is.
+ unpacked: Optional[tempfile.TemporaryDirectory] = None
+
+
+def _bundle(source: Path) -> Path:
+ if source.is_file() and source.name != "cell_feature_matrix.h5":
+ raise ValueError(
+ f"--from xenium reads a bundle: give its folder, or its "
+ f"cell_feature_matrix.h5, not {source.name}."
+ )
+ bundle = source.parent if source.is_file() else source
+ if not (bundle / "cell_feature_matrix.h5").is_file():
+ if (bundle / "cell_feature_matrix").is_dir():
+ raise ValueError(
+ f"{bundle} has the cell matrix as a Matrix Market folder only, not "
+ "cell_feature_matrix.h5, which this reads. Convert "
+ "cell_feature_matrix/ with --from 10x-mtx, without the cell table."
+ )
+ raise ValueError(f"{bundle} has no cell_feature_matrix.h5.")
+ return bundle
+
+
+def _unpack_analysis(archive: Path, beside: Path) -> tempfile.TemporaryDirectory:
+ """The analysis tables this reads, out of Xenium's analysis.tar.gz.
+
+ Only those members, each copied to a path built here rather than the
+ one the archive names, so the archive cannot write anywhere else. Into
+ a folder beside the output rather than the system's temporary one,
+ which on a cluster node is often too small for a large run's PCA.
+ """
+ holder = tempfile.TemporaryDirectory(prefix=".adata-xenium-", dir=beside)
+ wanted = {str(Path("analysis") / p): p for p in (_CLUSTERS, *_PROJECTIONS.values())}
+ with tarfile.open(archive, "r:*") as tar:
+ for member in tar:
+ relative = wanted.get(member.name.lstrip("./"))
+ if relative is None or not member.isfile():
+ continue
+ target = Path(holder.name) / relative
+ target.parent.mkdir(parents=True, exist_ok=True)
+ with tar.extractfile(member) as source, open(target, "wb") as out:
+ shutil.copyfileobj(source, out)
+ return holder
+
+
+def _scan(path: Path, index: str, barcodes: Any) -> _Table:
+ joiner = BarcodeJoin(barcodes)
+ scan = scan_table(path, index=index, on_index=joiner)
+ return _Table(scan, joiner.finish())
+
+
+class XeniumImport(Converter):
+ options = tenx.TenxImport.options
+
+ def plan(
+ self,
+ source: Path,
+ output: Path,
+ src_format: str,
+ dst_format: str,
+ options: Options,
+ ) -> ConversionPlan:
+ source = Path(source)
+ bundle = _bundle(source)
+ matrix_path = bundle / "cell_feature_matrix.h5"
+ cell_file = next((bundle / n for n in _CELL_TABLES if (bundle / n).is_file()), None)
+ if cell_file is None:
+ parquet = " It has cells.parquet only, which is not read yet." if (
+ bundle / "cells.parquet"
+ ).is_file() else ""
+ raise ValueError(f"{bundle} has no cells.csv.gz.{parquet}")
+
+ analysis = bundle / "analysis"
+ unpacked = None
+ packed_but_empty = False
+ if not analysis.is_dir() and (bundle / "analysis.tar.gz").is_file():
+ unpacked = _unpack_analysis(bundle / "analysis.tar.gz", Path(output).resolve().parent)
+ analysis = Path(unpacked.name)
+ packed_but_empty = not any(analysis.rglob("*.csv"))
+ warnings: List[str] = []
+ with h5py.File(matrix_path, "r") as handle:
+ if not tenx.sniff(matrix_path):
+ raise ValueError(f"{matrix_path.name} is not a 10x HDF5 matrix.")
+ matrix, outputs, targets, warnings = tenx.plan_matrix_import(handle, output, options)
+ barcodes = handle[matrix.group]["barcodes"]
+ cells = _scan(cell_file, "cell_id", barcodes)
+ clusters = None
+ projections: Dict[str, _Table] = {}
+ if analysis.is_dir():
+ if (analysis / _CLUSTERS).is_file():
+ clusters = _scan(analysis / _CLUSTERS, "Barcode", barcodes)
+ for key, path in _PROJECTIONS.items():
+ if (analysis / path).is_file():
+ projections[key] = _scan(analysis / path, "Barcode", barcodes)
+
+ require_numeric(cells.scan, _CENTROIDS)
+ if clusters is not None and "Cluster" not in clusters.scan.header:
+ raise ValueError(f"{_CLUSTERS.name} in analysis has no Cluster column.")
+ for table in projections.values():
+ require_numeric(table.scan, [c.name for c in table.scan.columns])
+ cells.join.check(cell_file.name)
+ if cells.join.n_extra:
+ warnings.append(
+ f"{cells.join.n_extra:,} cells of {cell_file.name} are not in the "
+ "matrix; they are left out."
+ )
+
+ found = []
+ for label, table in [("clusters", clusters), *projections.items()]:
+ if table is None:
+ continue
+ found.append(label)
+ table.join.check(f"analysis {label}", allow_missing=True)
+ if table.join.n_missing:
+ warnings.append(
+ f"{table.join.n_missing:,} cells have no {label} in analysis/; "
+ "they are left empty."
+ )
+ if packed_but_empty:
+ warnings.append(
+ "analysis.tar.gz holds none of the clustering, UMAP or PCA tables "
+ "this reads (analysis/clustering/gene_expression_graphclust, ...)."
+ )
+ elif not analysis.is_dir():
+ warnings.append("No analysis/ or analysis.tar.gz: no clusters, UMAP or PCA to bring in.")
+
+ described = describe_import(matrix.n_barcodes, matrix.n_features, matrix.feature_types, targets)
+ return ConversionPlan(
+ source=source,
+ src_format=src_format,
+ dst_format=dst_format,
+ outputs=outputs,
+ summary=(
+ f"Xenium → {dst_format}, {described}, cells from {cell_file.name}"
+ + (f", analysis: {', '.join(found)}" if found else "")
+ ),
+ warnings=warnings,
+ detail=_Plan(
+ matrix=matrix,
+ matrix_path=matrix_path,
+ targets=targets,
+ cells=cells,
+ clusters=clusters,
+ projections=projections,
+ zarr_format=(options.zarr_format or 3) if dst_format == "zarr" else None,
+ unpacked=unpacked,
+ ),
+ )
+
+ def write(self, plan: ConversionPlan, staged: Sequence[Path], console: Console) -> None:
+ detail: _Plan = plan.detail
+ matrix = detail.matrix
+
+ with h5py.File(detail.matrix_path, "r") as handle, ExitStack() as stack:
+ group = handle[matrix.group]
+ roots = [
+ stack.enter_context(open_store(path, "w", zarr_format=detail.zarr_format)).root
+ for path in staged
+ ]
+ for root, target in zip(roots, detail.targets):
+ with console.status("[cyan]Writing obs, var and the cell table...[/]"):
+ obs = _write_obs_index(root, group)
+ ensure_anndata_skeleton(root)
+ columns = _write_cells(root, obs, detail.cells)
+ if detail.clusters is not None:
+ _write_clusters(root, obs, detail.clusters)
+ columns.append("graphclust")
+ set_column_order(obs, columns)
+ for key, table in detail.projections.items():
+ _write_projection(root, key, table)
+ write_var(root, target)
+ tenx._write_uns(root, matrix)
+
+ with console.status("[cyan]Copying X...[/]"):
+ fan_out_sparse_matrix(
+ tenx._AsCsr(group, (matrix.n_barcodes, matrix.n_features)),
+ "X",
+ [(root, None, target.keep) for root, target in zip(roots, detail.targets)],
+ )
+
+
+def _write_obs_index(root: Any, group: h5py.Group) -> Any:
+ barcodes = group["barcodes"]
+ n = barcodes.shape[0]
+ return write_obs_index(
+ root, n,
+ (
+ [str(v) for v in decode_str_array(np.asarray(barcodes[lo : lo + STRING_STEP]))]
+ for lo in range(0, n, STRING_STEP)
+ ),
+ )
+
+
+def _write_cells(root: Any, obs: Any, table: _Table) -> List[str]:
+ """The cell table's columns into obs, its centroids into `obsm["spatial"]`."""
+ scan, join = table.scan, table.join
+ n = join.n_out
+ columns, lookups = create_columns(obs, scan, n, exclude=_CENTROIDS)
+ spatial = create_dataset(root["obsm"], "spatial", shape=(n, 2), dtype=np.float64)
+ spec.set_encoding(spatial, spec.ARRAY)
+
+ by_name = {c.name: c for c in scan.columns}
+ sinks = [
+ Sink(
+ dest,
+ lambda block, column=by_name[name]: column_values(
+ scan, column, block[column.name], lookups.get(column.name, {})
+ ),
+ )
+ for name, dest in columns.items()
+ ]
+ sinks.append(Sink(spatial, lambda block: numbers(scan, _CENTROIDS, block)))
+ write_rows(scan, join, sinks, root)
+
+ if join.n_extra:
+ # Categories only the left-out cells used go, as anndata's
+ # subsetting drops them.
+ for name in lookups:
+ used = _used_codes(columns[name], len(lookups[name]), n)
+ if not used.all():
+ drop_unused(obs, name, used, n, block_rows=STRING_STEP)
+ return list(columns)
+
+
+def _used_codes(codes: Any, n_categories: int, n: int) -> np.ndarray:
+ used = np.zeros(n_categories, dtype=bool)
+ for lo in range(0, n, STRING_STEP):
+ block = np.asarray(codes[lo : lo + STRING_STEP], dtype=np.int64)
+ used[block[block >= 0]] = True
+ return used
+
+
+def _write_clusters(root: Any, obs: Any, table: _Table) -> None:
+ """Xenium's graph-based clusters, `obs["graphclust"]`, a categorical."""
+ scan, join = table.scan, table.join
+ column = next(c for c in scan.columns if c.name == "Cluster")
+ categories = natural_sorted(column.seen or set())
+ lookup = {c: i for i, c in enumerate(categories)}
+ codes = create_categorical(obs, "graphclust", categories, join.n_out)
+ write_rows(
+ scan, join,
+ [Sink(codes, lambda block: np.array([lookup.get(v, -1) for v in block["Cluster"]], dtype=np.int64), -1)],
+ root,
+ )
+ if join.n_extra:
+ used = _used_codes(codes, len(categories), join.n_out)
+ if not used.all():
+ drop_unused(obs, "graphclust", used, join.n_out, block_rows=STRING_STEP)
+
+
+def _write_projection(root: Any, key: str, table: _Table) -> None:
+ """An analysis projection (UMAP, PCA) into obsm, NaN for a cell it lacks."""
+ scan, join = table.scan, table.join
+ names = [c.name for c in scan.columns]
+ out = create_dataset(root["obsm"], key, shape=(join.n_out, len(names)), dtype=np.float64)
+ spec.set_encoding(out, spec.ARRAY)
+ write_rows(
+ scan, join,
+ [Sink(out, lambda block: numbers(scan, names, block), np.nan)],
+ root,
+ )
diff --git a/src/adata/core/scatter.py b/src/adata/core/scatter.py
new file mode 100644
index 0000000..5cf2b0b
--- /dev/null
+++ b/src/adata/core/scatter.py
@@ -0,0 +1,148 @@
+"""Rows read in one order, written in another, a block at a time.
+
+A table beside a matrix (Xenium's `cells.csv.gz`, an analysis projection)
+lists its rows in its own order, and AnnData wants them in the matrix's. The
+table can only be read front to back, and writing each row where it belongs
+is a fancy-indexed write per block, which HDF5 and Zarr do some 500 times
+slower than a slice.
+
+So the rows go through scratch: each block read is sorted by the output
+block its rows belong to and appended there, then each output block is read
+back from scratch in one slice, placed, and written in one slice. Memory is
+one block per column; the scratch costs one extra write and read of the
+values. Several columns share the layout, and are fed from one read of the
+table.
+"""
+
+from __future__ import annotations
+
+from contextlib import ExitStack
+from typing import Any, List
+
+import numpy as np
+
+from adata.core.csr_build import scratch_group
+from adata.storage import create_dataset
+
+#: Output rows per block.
+ROW_STEP = 1 << 16
+
+
+class RowScatter:
+ """Where each source row goes, shared by every column fed through it.
+
+ `target[i]` is the output row of source row `i`, or -1 to drop it; no
+ two source rows may share an output row. Output rows no source row
+ reaches keep the column's fill value.
+ """
+
+ def __init__(self, target: np.ndarray, n_out: int, scratch_parent: Any, row_step: int = ROW_STEP) -> None:
+ self.target = np.asarray(target, dtype=np.int64)
+ self.n_out = n_out
+ self.row_step = row_step
+ kept = np.flatnonzero(self.target >= 0)
+ block_of = self.target[kept] // row_step
+ self.n_blocks = -(-n_out // row_step)
+ order = np.argsort(block_of, kind="stable")
+ self.starts = np.concatenate(
+ ([0], np.cumsum(np.bincount(block_of, minlength=self.n_blocks)))
+ )
+ #: The output row of each scratch row.
+ self.placed = self.target[kept][order]
+ self.n_kept = int(kept.size)
+ self._stack = ExitStack()
+ self._scratch = self._stack.enter_context(scratch_group(scratch_parent, "__scatter_scratch__"))
+ self._columns: List["_Column"] = []
+
+ def column(self, dest: Any, fill: Any = None) -> "_Column":
+ """A column to fill: `dest` is created at `(n_out, ...)` already.
+
+ `fill` is what output rows no source row reaches hold; by default NaN
+ for floats, -1 for integers (a categorical's missing code), False
+ and the empty string."""
+ column = _Column(self, dest, fill, len(self._columns))
+ self._columns.append(column)
+ return column
+
+ def finish(self) -> None:
+ """Write every column out of scratch, and remove the scratch."""
+ try:
+ for column in self._columns:
+ column.finish()
+ finally:
+ self._stack.close()
+
+ def __enter__(self) -> "RowScatter":
+ return self
+
+ def __exit__(self, *exc: Any) -> None:
+ if exc[0] is None:
+ self.finish()
+ else:
+ self._stack.close()
+
+
+class _Column:
+ def __init__(self, layout: RowScatter, dest: Any, fill: Any, number: int) -> None:
+ self.layout = layout
+ self.dest = dest
+ self.fill = _default_fill(dest.dtype) if fill is None else fill
+ shape = (layout.n_kept,) + tuple(dest.shape[1:])
+ # HDF5 refuses a chunk larger than the data, so an empty scratch --
+ # no source row kept -- has none.
+ chunks = {"chunks": (min(layout.n_kept, layout.row_step),) + shape[1:]} if layout.n_kept else {}
+ self.scratch = create_dataset(
+ layout._scratch, f"c{number}", shape=shape, dtype=dest.dtype, **chunks
+ )
+ self.cursor = layout.starts[:-1].copy()
+ self.done = 0
+
+ def add(self, values: np.ndarray) -> None:
+ """The next block of source rows, in source order."""
+ layout = self.layout
+ lo, hi = self.done, self.done + len(values)
+ self.done = hi
+ here = layout.target[lo:hi]
+ have = np.flatnonzero(here >= 0)
+ if have.size == 0:
+ return
+ values = np.asarray(values)[have]
+ blocks = here[have] // layout.row_step
+ order = np.argsort(blocks, kind="stable")
+ values, blocks = values[order], blocks[order]
+ edges = np.flatnonzero(np.diff(blocks)) + 1
+ for run, b in zip(np.split(values, edges), blocks[np.r_[0, edges]]):
+ self.scratch[self.cursor[b] : self.cursor[b] + len(run)] = run
+ self.cursor[b] += len(run)
+
+ def finish(self) -> None:
+ layout = self.layout
+ if self.done != layout.target.size:
+ raise ValueError(
+ f"{self.done:,} source rows were given, where {layout.target.size:,} "
+ "were planned."
+ )
+ dtype = self.dest.dtype
+ object_dtype = dtype.kind not in "biufc"
+ for b in range(layout.n_blocks):
+ lo = b * layout.row_step
+ rows = min(layout.row_step, layout.n_out - lo)
+ block = np.full(
+ (rows,) + tuple(self.dest.shape[1:]), self.fill,
+ dtype=object if object_dtype else dtype,
+ )
+ a, z = layout.starts[b], layout.starts[b + 1]
+ if a < z:
+ block[layout.placed[a:z] - lo] = self.scratch[a:z]
+ self.dest[lo : lo + rows] = block
+
+
+def _default_fill(dtype: np.dtype) -> Any:
+ kind = np.dtype(dtype).kind
+ if kind in "fc":
+ return np.nan
+ if kind in "iu":
+ return -1 if kind == "i" else 0
+ if kind == "b":
+ return False
+ return ""
diff --git a/src/adata/formats/json_data.py b/src/adata/formats/json_data.py
index 1082c79..db4892b 100644
--- a/src/adata/formats/json_data.py
+++ b/src/adata/formats/json_data.py
@@ -136,17 +136,17 @@ def import_json(
if name in parent:
del parent[name]
- _write_json_to_group(parent, name, payload)
+ write_json_value(parent, name, payload)
console.print(f"[green]Imported[/] JSON data into '{obj}'")
-def _write_json_to_group(parent: Any, name: str, value: Any) -> None:
+def write_json_value(parent: Any, name: str, value: Any) -> None:
"""Write one JSON value as the AnnData element that best represents it."""
if isinstance(value, dict):
group = write_mapping(parent, name, replace=True)
for k, v in value.items():
- _write_json_to_group(group, k, v)
+ write_json_value(group, k, v)
return
if value is None:
@@ -194,3 +194,4 @@ def _write_json_list(parent: Any, name: str, value: list) -> None:
return
write_scalar(parent, name, json.dumps(value), replace=True)
+
diff --git a/src/adata/formats/table.py b/src/adata/formats/table.py
index 2a56d7d..3163068 100644
--- a/src/adata/formats/table.py
+++ b/src/adata/formats/table.py
@@ -83,6 +83,8 @@ class TableScan:
columns: List[Column]
#: Rows (0-based, after the header) where `select` matched, if asked.
selected: Optional[np.ndarray] = None
+ #: False for a file with no header line, its names given by the caller.
+ has_header: bool = True
def dtype_of(self, column: Column) -> str:
"""What the column is written as: int64, float64, bool, category or str."""
@@ -110,9 +112,13 @@ def _rows(
) -> Iterator[Tuple[int, List[List[str]]]]:
"""Blocks of parsed rows, each with the file line number of its first."""
line = 1
+ first = True
with open_binary(path) as fh:
for block in iter_line_blocks(fh, block_bytes):
- reader = csv.reader(io.StringIO(block.decode("utf-8")), delimiter=delimiter, strict=True)
+ # A byte order mark, as Excel saves CSV, is not part of the first field.
+ text = block.decode("utf-8-sig" if first else "utf-8")
+ first = False
+ reader = csv.reader(io.StringIO(text), delimiter=delimiter, strict=True)
rows: List[List[str]] = []
start = line
try:
@@ -158,40 +164,50 @@ def scan_table(
index: Optional[str] = None,
delimiter: str = ",",
select: Optional[Tuple[str, str]] = None,
+ names: Optional[Sequence[str]] = None,
+ on_index: Optional[Callable[[List[str]], None]] = None,
block_bytes: int = TABLE_BLOCK_BYTES,
) -> TableScan:
"""Pass 1: the header, the row count and each column's type.
`index` names the column that becomes the row names (the first, if not
given). `select=(column, value)` also records the rows where that column
- holds that value.
+ holds that value. `names` are the columns of a file with no header line,
+ such as Space Ranger 1's `tissue_positions_list.csv`. `on_index` is
+ given each block's index values, in file order, so a caller can match
+ them against something else in the same pass.
"""
path = Path(path)
header: Optional[List[str]] = None
+ if names is not None:
+ header = list(names)
columns: List[Column] = []
positions: List[int] = []
pick = -1
selected: List[np.ndarray] = []
n_rows = 0
+ at = -1
for start, rows in _rows(path, delimiter, block_bytes):
- if header is None:
- rows = [row for row in rows if row]
- if not rows:
- continue
- header, rows = rows[0], rows[1:]
+ if not columns:
+ if header is None:
+ rows = [row for row in rows if row]
+ if not rows:
+ continue
+ header, rows = rows[0], rows[1:]
+ start += 1
if len(set(header)) != len(header):
repeated = sorted({h for h in header if header.count(h) > 1})
raise ValueError(f"{path.name} repeats the column {', '.join(map(repr, repeated))}.")
index = header[0] if index is None else index
if index not in header:
raise ValueError(f"{path.name} has no {index!r} column; it has: {', '.join(header)}.")
+ at = header.index(index)
positions = [i for i, name in enumerate(header) if name != index]
columns = [Column(header[i]) for i in positions]
if select is not None:
if select[0] not in header:
raise ValueError(f"{path.name} has no {select[0]!r} column.")
pick = header.index(select[0])
- start += 1
_check_widths(path, rows, len(header), start)
rows = [row for row in rows if row]
if not rows:
@@ -199,12 +215,18 @@ def scan_table(
by_field = list(zip(*rows))
for column, i in zip(columns, positions):
_narrow(column, list(by_field[i]))
+ if on_index is not None:
+ on_index(list(by_field[at]))
if pick >= 0:
hits = np.flatnonzero(np.asarray(by_field[pick], dtype=object) == select[1])
selected.append(hits + n_rows)
n_rows += len(rows)
if header is None:
raise ValueError(f"{path.name} is empty: it has no header line.")
+ if not columns:
+ # A file of no rows: its header (or the names given) still stands.
+ index = header[0] if index is None else index
+ columns = [Column(name) for name in header if name != index]
return TableScan(
path=path,
delimiter=delimiter,
@@ -212,6 +234,7 @@ def scan_table(
index=index,
n_rows=n_rows,
columns=columns,
+ has_header=names is None,
selected=(
np.concatenate(selected) if selected else np.zeros(0, dtype=np.int64)
) if select is not None else None,
@@ -239,7 +262,7 @@ def iter_table_blocks(
"""
at = scan.header.index(scan.index)
done = 0
- first = True
+ first = scan.has_header
for _, block in _rows(scan.path, scan.delimiter, block_bytes):
block = [row for row in block if row]
if first and block:
@@ -261,7 +284,9 @@ def iter_table_blocks(
}
-def _values(scan: TableScan, column: Column, values: List[str], lookup: Dict[str, int]) -> np.ndarray:
+def column_values(scan: TableScan, column: Column, values: List[str], lookup: Dict[str, int]) -> np.ndarray:
+ """A block of one column's text, as the values `create_columns` stores:
+ numbers, booleans, category codes (-1 for a gap) or strings."""
kind = scan.dtype_of(column)
if kind == "int64":
return np.array(values, dtype=np.int64)
@@ -296,20 +321,7 @@ def write_table(
index = create_dataset(group, "_index", shape=(n,), dtype=string_dtype_for(backend))
spec.set_encoding(index, spec.STRING_ARRAY)
- targets: Dict[str, Any] = {}
- lookups: Dict[str, Dict[str, int]] = {}
- for column in scan.columns:
- kind = scan.dtype_of(column)
- if kind == "category":
- categories = column.categories or []
- targets[column.name] = create_categorical(group, column.name, categories, n)
- lookups[column.name] = {c: i for i, c in enumerate(categories)}
- continue
- dtype = string_dtype_for(backend) if kind == "str" else np.dtype(kind)
- targets[column.name] = create_dataset(group, column.name, shape=(n,), dtype=dtype)
- spec.set_encoding(
- targets[column.name], spec.STRING_ARRAY if kind == "str" else spec.ARRAY
- )
+ targets, lookups = create_columns(group, scan, n)
used = {name: np.zeros(len(lookup), dtype=bool) for name, lookup in lookups.items()}
written = 0
@@ -317,7 +329,7 @@ def write_table(
end = written + len(names)
index[written:end] = np.asarray(names, dtype=object)
for column in scan.columns:
- values = _values(scan, column, block[column.name], lookups.get(column.name, {}))
+ values = column_values(scan, column, block[column.name], lookups.get(column.name, {}))
targets[column.name][written:end] = values
if column.name in used:
used[column.name][values[values >= 0]] = True
@@ -327,11 +339,40 @@ def write_table(
if rows is not None:
for name, seen in used.items():
if not seen.all():
- _drop_unused(group, name, seen, n, block_rows=max(1, block_bytes // 8))
+ drop_unused(group, name, seen, n, block_rows=max(1, block_bytes // 8))
return group
-def _drop_unused(group: Any, name: str, used: np.ndarray, n: int, block_rows: int) -> None:
+def create_columns(
+ group: Any, scan: TableScan, n: int, *, exclude: Sequence[str] = ()
+) -> Tuple[Dict[str, Any], Dict[str, Dict[str, int]]]:
+ """Each column of `scan` but those in `exclude`, created at length `n`
+ in `group` as anndata would store it, to be filled a block at a time.
+
+ Returns the arrays by name, and each categorical's category lookup for
+ `column_values`.
+ """
+ backend = "zarr" if is_zarr_group(group) else "hdf5"
+ targets: Dict[str, Any] = {}
+ lookups: Dict[str, Dict[str, int]] = {}
+ for column in scan.columns:
+ if column.name in exclude:
+ continue
+ kind = scan.dtype_of(column)
+ if kind == "category":
+ categories = column.categories or []
+ targets[column.name] = create_categorical(group, column.name, categories, n)
+ lookups[column.name] = {c: i for i, c in enumerate(categories)}
+ continue
+ dtype = string_dtype_for(backend) if kind == "str" else np.dtype(kind)
+ targets[column.name] = create_dataset(group, column.name, shape=(n,), dtype=dtype)
+ spec.set_encoding(
+ targets[column.name], spec.STRING_ARRAY if kind == "str" else spec.ARRAY
+ )
+ return targets, lookups
+
+
+def drop_unused(group: Any, name: str, used: np.ndarray, n: int, block_rows: int) -> None:
"""Remove categories no written row uses, as anndata does on subsetting."""
remap = np.full(used.size + 1, -1, dtype=np.int64) # the last slot maps -1
remap[:-1][used] = np.arange(int(used.sum()))
@@ -345,10 +386,12 @@ def _drop_unused(group: Any, name: str, used: np.ndarray, n: int, block_rows: in
codes[lo : lo + block.size] = remap[block]
-def read_small_table(path: Path, delimiter: str = ",") -> Tuple[List[str], List[List[str]]]:
+def read_small_table(
+ path: Path, delimiter: str = ",", *, names: Optional[Sequence[str]] = None
+) -> Tuple[List[str], List[List[str]]]:
"""The header and every row of a table small enough to hold, such as a
- gene list."""
- header: Optional[List[str]] = None
+ gene list. `names` are the columns of a file with no header line."""
+ header: Optional[List[str]] = None if names is None else list(names)
body: List[List[str]] = []
for start, rows in _rows(Path(path), delimiter, TABLE_BLOCK_BYTES):
if header is None:
diff --git a/tests/conftest.py b/tests/conftest.py
index fc0ec0c..4f7019f 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -1,6 +1,7 @@
"""Pytest configuration and fixtures for adata tests."""
import base64
+import gzip
import json
import zlib
from contextlib import contextmanager
@@ -1327,3 +1328,218 @@ def write_h5mu(
with open(path, "r+b") as raw:
raw.write(header.ljust(512, b"\0"))
return MudataTruth(path, mods, obs, obsm, obsp, uns)
+
+
+# ---------------------------------------------------------------------------
+# Space Ranger output for standard Visium
+#
+# As Space Ranger writes it (checked against 10x's SR 1.1 and SR 2.1 sample
+# outputs): the matrix, and in spatial/ every spot's position -- in no
+# particular order, and including spots the filtered matrix leaves out --
+# the scale factors and two PNGs. SR 1 wrote the positions with no header.
+
+
+@dataclass
+class VisiumTruth:
+ outs: Path
+ barcodes: List[str] # the filtered matrix's, in its order
+ in_tissue: np.ndarray
+ array_row: np.ndarray
+ array_col: np.ndarray
+ spatial: np.ndarray # [pxl_col, pxl_row], per barcode
+ images: Dict[str, np.ndarray]
+ scalefactors: dict
+ counts: np.ndarray
+
+
+def write_visium_outs(
+ folder: Path,
+ *,
+ version: str = "sr2",
+ n_spots: int = 12,
+ n_extra: int = 5,
+ raw: bool = True,
+ images: bool = True,
+ seed: int = 0,
+) -> VisiumTruth:
+ """An `outs/` folder. `version` is sr1 (no header), sr2, or sr3_float
+ (pixel positions with decimals, as Space Ranger 3 and 4 may write)."""
+ from PIL import Image
+
+ rng = np.random.default_rng(seed)
+ outs = folder / "outs"
+ spatial = outs / "spatial"
+ spatial.mkdir(parents=True)
+ n_features = 6
+ barcodes = [f"AAAC{i:04d}-1" for i in range(n_spots)]
+ extra = [f"TTTG{i:04d}-1" for i in range(n_extra)]
+ counts = rng.integers(0, 5, size=(n_spots, n_features))
+ ids = [f"ENSG{i:05d}" for i in range(n_features)]
+ names = [f"G{i}" for i in range(n_features)]
+ types = ["Gene Expression"] * (n_features - 1) + ["Antibody Capture"]
+
+ def matrix(path: Path, rows: List[str], values: np.ndarray) -> None:
+ write_10x_h5(path, values, barcodes=rows, ids=ids, names=names, types=types)
+ with h5py.File(path, "a") as f:
+ f.attrs["library_ids"] = np.array([b"V1_Test_Section"])
+ f.attrs["chemistry_description"] = "Spatial 3' v1"
+ f.attrs["software_version"] = "spaceranger-2.1.0"
+
+ matrix(outs / "filtered_feature_bc_matrix.h5", barcodes, counts)
+ if raw:
+ matrix(
+ outs / "raw_feature_bc_matrix.h5", barcodes + extra,
+ np.vstack([counts, np.zeros((n_extra, n_features), dtype=counts.dtype)]),
+ )
+
+ every = barcodes + extra
+ n = len(every)
+ in_tissue = np.array([1] * n_spots + [0] * n_extra)
+ array_row = rng.integers(0, 78, n)
+ array_col = rng.integers(0, 128, n)
+ pxl_row = rng.integers(-50, 2000, n).astype(float)
+ pxl_col = rng.integers(0, 2000, n).astype(float)
+ if version == "sr3_float":
+ pxl_row += 0.25
+ pxl_col += 0.5
+ order = rng.permutation(n)
+ lines = []
+ if version != "sr1":
+ lines.append(",".join(
+ ("barcode", "in_tissue", "array_row", "array_col", "pxl_row_in_fullres", "pxl_col_in_fullres")
+ ))
+ fmt = (lambda v: f"{v:g}") if version == "sr3_float" else (lambda v: str(int(v)))
+ for i in order:
+ lines.append(
+ f"{every[i]},{in_tissue[i]},{array_row[i]},{array_col[i]},{fmt(pxl_row[i])},{fmt(pxl_col[i])}"
+ )
+ name = "tissue_positions_list.csv" if version == "sr1" else "tissue_positions.csv"
+ (spatial / name).write_text("\n".join(lines) + "\n")
+
+ scalefactors = {
+ "spot_diameter_fullres": 89.4,
+ "tissue_hires_scalef": 0.17,
+ "fiducial_diameter_fullres": 144.5,
+ "tissue_lowres_scalef": 0.05,
+ }
+ (spatial / "scalefactors_json.json").write_text(json.dumps(scalefactors))
+ pictures: Dict[str, np.ndarray] = {}
+ if images:
+ for key, size in (("hires", (20, 16)), ("lowres", (6, 5))):
+ pixels = rng.integers(0, 256, size=size + (3,), dtype=np.uint8)
+ Image.fromarray(pixels).save(spatial / f"tissue_{key}_image.png")
+ pictures[key] = pixels
+
+ take = slice(0, n_spots)
+ return VisiumTruth(
+ outs=outs,
+ barcodes=barcodes,
+ in_tissue=in_tissue[take],
+ array_row=array_row[take],
+ array_col=array_col[take],
+ spatial=np.stack([pxl_col, pxl_row], axis=1)[take],
+ images=pictures,
+ scalefactors=scalefactors,
+ counts=counts,
+ )
+
+
+# ---------------------------------------------------------------------------
+# Xenium output bundles
+#
+# As Xenium Onboard Analysis writes them (columns checked against 10x's XOA
+# 1.0.1, 3.0 and 4.0 sample bundles): the cell matrix, cells.csv.gz in the
+# matrix's order, experiment.xenium, and analysis/ CSVs keyed by Barcode.
+# XOA 1 numbered cells; later versions name them like "aaafiiei-1".
+
+
+@dataclass
+class XeniumTruth:
+ bundle: Path
+ cell_ids: List[str] # the matrix's, in its order
+ columns: Dict[str, np.ndarray] # cells.csv.gz's other columns, per cell
+ spatial: np.ndarray
+ clusters: Optional[np.ndarray] # per cell, "" where analysis lacks it
+ umap: Optional[np.ndarray] # NaN where analysis lacks it
+ counts: np.ndarray
+
+
+def write_xenium_bundle(
+ folder: Path,
+ *,
+ ids: str = "string",
+ n_cells: int = 15,
+ shuffled: bool = False,
+ extra_cells: int = 0,
+ analysis: bool = True,
+ missing_in_analysis: int = 0,
+ seed: int = 0,
+) -> XeniumTruth:
+ rng = np.random.default_rng(seed)
+ folder.mkdir(parents=True)
+ n_features = 5
+ if ids == "int":
+ cell_ids = [str(i + 1) for i in range(n_cells + extra_cells)]
+ else:
+ cell_ids = [
+ "".join("abcdefghijklmnop"[d] for d in rng.integers(0, 16, 8)) + "-1"
+ for _ in range(n_cells + extra_cells)
+ ]
+ counts = rng.integers(0, 6, size=(n_cells, n_features))
+ write_10x_h5(
+ folder / "cell_feature_matrix.h5", counts, barcodes=cell_ids[:n_cells],
+ ids=[f"ENSG{i:05d}" for i in range(n_features - 1)] + ["NegControlProbe_00001"],
+ names=[f"G{i}" for i in range(n_features - 1)] + ["NegControlProbe_00001"],
+ types=["Gene Expression"] * (n_features - 1) + ["Negative Control Probe"],
+ )
+ (folder / "experiment.xenium").write_text(json.dumps({"pixel_size": 0.2125, "major_version": 5}))
+
+ n = n_cells + extra_cells
+ table = {
+ "x_centroid": rng.uniform(0, 5000, n).round(4),
+ "y_centroid": rng.uniform(0, 5000, n).round(4),
+ "transcript_counts": rng.integers(0, 300, n),
+ "control_probe_counts": rng.integers(0, 3, n),
+ "cell_area": rng.uniform(10, 200, n).round(3),
+ "nucleus_area": rng.uniform(5, 50, n).round(3),
+ "segmentation_method": rng.choice(["Segmented by boundary stain", "Segmented by nucleus expansion"], n),
+ }
+ if extra_cells:
+ # A method only the left-out cells use: its category must go.
+ table["segmentation_method"][n_cells:] = "Imported cell segmentation"
+ order = rng.permutation(n) if shuffled else np.arange(n)
+ names = ["cell_id", *table]
+ lines = [",".join(names)]
+ for i in order:
+ values = [cell_ids[i]] + [str(table[k][i]) for k in table]
+ lines.append(",".join(f'"{v}"' if k in ("cell_id", "segmentation_method") else v for k, v in zip(names, values)))
+ (folder / "cells.csv.gz").write_bytes(gzip.compress(("\n".join(lines) + "\n").encode()))
+
+ clusters = umap = None
+ if analysis:
+ listed = rng.permutation(n_cells)[missing_in_analysis:]
+ clusters = np.full(n_cells, "", dtype=object)
+ clusters[listed] = [str(v) for v in rng.integers(1, 12, listed.size)]
+ umap = np.full((n_cells, 2), np.nan)
+ umap[listed] = rng.normal(size=(listed.size, 2)).round(5)
+ path = folder / "analysis" / "clustering" / "gene_expression_graphclust"
+ path.mkdir(parents=True)
+ (path / "clusters.csv").write_text(
+ "Barcode,Cluster\n" + "".join(f"{cell_ids[i]},{clusters[i]}\n" for i in listed)
+ )
+ path = folder / "analysis" / "umap" / "gene_expression_2_components"
+ path.mkdir(parents=True)
+ (path / "projection.csv").write_text(
+ "Barcode,UMAP-1,UMAP-2\n" + "".join(f"{cell_ids[i]},{umap[i, 0]},{umap[i, 1]}\n" for i in listed)
+ )
+
+ keep = slice(0, n_cells)
+ return XeniumTruth(
+ bundle=folder,
+ cell_ids=cell_ids[:n_cells],
+ columns={k: v[keep] for k, v in table.items() if k not in ("x_centroid", "y_centroid")},
+ spatial=np.stack([table["x_centroid"], table["y_centroid"]], axis=1)[keep],
+ clusters=clusters,
+ umap=umap,
+ counts=counts,
+ )
diff --git a/tests/test_external_readers.py b/tests/test_external_readers.py
index a4e367f..828678d 100644
--- a/tests/test_external_readers.py
+++ b/tests/test_external_readers.py
@@ -25,6 +25,8 @@
from adata.cli import app
from tests.conftest import (
+ write_visium_outs,
+ write_xenium_bundle,
write_10x_h5,
write_cellbender_h5,
write_hdf5_matrix,
@@ -1022,3 +1024,125 @@ def test_mudata_reads_what_we_export(tmp_path):
obj.write_h5ad(tmp_path / "multi.h5ad")
_convert(str(tmp_path / "multi.h5ad"), "-o", str(tmp_path / "multi.h5mu"))
assert _run(_MUDATA, _COMPARE_MUDATA_EXPORT, str(tmp_path / "multi.h5ad"), str(tmp_path / "multi.h5mu")) == []
+
+
+# ---------------------------------------------------------------------------
+# spatial: Visium against squidpy's reader, Xenium against 10x's recipe
+
+_COMPARE_VISIUM = r"""
+import json, sys, warnings
+warnings.filterwarnings("ignore")
+import numpy as np, anndata as ad, squidpy as sq
+
+outs, h5ad = sys.argv[1], sys.argv[2]
+ref = sq.read.visium(outs, counts_file="filtered_feature_bc_matrix.h5", load_images=True)
+got = ad.read_h5ad(h5ad)
+problems = []
+if ref.shape != got.shape:
+ problems.append(f"shape {ref.shape} != {got.shape}")
+else:
+ if not np.array_equal(ref.X.toarray(), got.X.toarray()):
+ problems.append("X differs")
+ if list(ref.obs_names) != list(got.obs_names):
+ problems.append("obs_names differ")
+ if list(ref.var_names) != list(got.var_names):
+ problems.append("var_names differ")
+ for column in ("in_tissue", "array_row", "array_col"):
+ if ref.obs[column].dtype != got.obs[column].dtype:
+ problems.append(f"obs[{column!r}] is {got.obs[column].dtype}, not {ref.obs[column].dtype}")
+ if not np.array_equal(ref.obs[column], got.obs[column]):
+ problems.append(f"obs[{column!r}] differs")
+ a, b = ref.obsm["spatial"], got.obsm["spatial"]
+ if a.dtype != b.dtype:
+ problems.append(f"obsm['spatial'] is {b.dtype}, not {a.dtype}")
+ if not np.array_equal(a, b):
+ problems.append("obsm['spatial'] differs")
+ (lib,) = ref.uns["spatial"]
+ if lib not in got.uns["spatial"]:
+ problems.append(f"no library {lib!r}")
+ else:
+ r, g = ref.uns["spatial"][lib], got.uns["spatial"][lib]
+ for name, image in r["images"].items():
+ mine = g["images"].get(name)
+ if mine is None or mine.dtype != image.dtype or not np.array_equal(mine, image):
+ problems.append(f"image {name} differs")
+ if {k: float(v) for k, v in r["scalefactors"].items()} != {k: float(v) for k, v in g["scalefactors"].items()}:
+ problems.append("scalefactors differ")
+ for key in ("chemistry_description", "software_version"):
+ if str(r["metadata"].get(key)) != str(g["metadata"].get(key)):
+ problems.append(f"metadata[{key!r}] differs")
+print("RESULT " + json.dumps(problems))
+"""
+
+
+@pytest.mark.parametrize("version", ["sr1", "sr2", "sr3_float"])
+def test_visium_imports_as_squidpy_reads_it(tmp_path, version):
+ truth = write_visium_outs(tmp_path / "run", version=version)
+ out = tmp_path / "visium.h5ad"
+ # squidpy keeps Gene Expression only and names var by symbol.
+ _convert(
+ str(truth.outs), "-o", str(out),
+ "--feature-type", "Gene Expression", "--var-names", "gene_symbols",
+ )
+ assert _run(["squidpy"], _COMPARE_VISIUM, str(truth.outs), str(out)) == []
+
+
+# 10x's documented recipe for Xenium in scanpy: the matrix, the cell table
+# reindexed to it, the centroids as obsm["spatial"], and the analysis CSVs.
+# spatialdata-io's reader needs the boundaries and morphology images too,
+# which a fixture of this size does not carry.
+_COMPARE_XENIUM = r"""
+import json, sys, warnings
+from pathlib import Path
+warnings.filterwarnings("ignore")
+import numpy as np, pandas as pd, anndata as ad, scanpy as sc
+
+bundle, h5ad = Path(sys.argv[1]), sys.argv[2]
+ref = sc.read_10x_h5(bundle / "cell_feature_matrix.h5", gex_only=False)
+cells = pd.read_csv(bundle / "cells.csv.gz", index_col="cell_id")
+cells.index = cells.index.astype(str)
+cells = cells.reindex(ref.obs_names)
+ref.obsm["spatial"] = cells[["x_centroid", "y_centroid"]].to_numpy()
+ref.obs = cells.drop(columns=["x_centroid", "y_centroid"])
+analysis = bundle / "analysis"
+clusters = pd.read_csv(analysis / "clustering/gene_expression_graphclust/clusters.csv", index_col="Barcode")
+clusters.index = clusters.index.astype(str)
+ref.obs["graphclust"] = clusters["Cluster"].astype(str).reindex(ref.obs_names).astype("category")
+umap = pd.read_csv(analysis / "umap/gene_expression_2_components/projection.csv", index_col="Barcode")
+umap.index = umap.index.astype(str)
+ref.obsm["X_umap"] = umap.reindex(ref.obs_names).to_numpy()
+ref.strings_to_categoricals()
+
+got = ad.read_h5ad(h5ad)
+problems = []
+if ref.shape != got.shape:
+ problems.append(f"shape {ref.shape} != {got.shape}")
+else:
+ if not np.array_equal(ref.X.toarray(), got.X.toarray()):
+ problems.append("X differs")
+ if list(ref.obs_names) != list(got.obs_names):
+ problems.append("obs_names differ")
+ if list(ref.obs.columns) != list(got.obs.columns):
+ problems.append(f"obs columns {list(ref.obs.columns)} != {list(got.obs.columns)}")
+ for column in ref.obs.columns:
+ if column not in got.obs:
+ continue
+ a, b = ref.obs[column], got.obs[column]
+ if str(a.dtype) != str(b.dtype):
+ problems.append(f"obs[{column!r}] is {b.dtype}, not {a.dtype}")
+ elif a.astype(str).tolist() != b.astype(str).tolist():
+ problems.append(f"obs[{column!r}] differs")
+ for key in ("spatial", "X_umap"):
+ if not np.allclose(ref.obsm[key], got.obsm[key], equal_nan=True):
+ problems.append(f"obsm[{key!r}] differs")
+print("RESULT " + json.dumps(problems))
+"""
+
+
+@pytest.mark.parametrize("ids,shuffled", [("int", False), ("string", True)])
+def test_xenium_imports_as_10x_s_scanpy_recipe_reads_it(tmp_path, ids, shuffled):
+ truth = write_xenium_bundle(tmp_path / "xenium", ids=ids, shuffled=shuffled, missing_in_analysis=2)
+ out = tmp_path / "xenium.h5ad"
+ _convert(str(truth.bundle), "-o", str(out), "--var-names", "gene_symbols")
+ assert _run(["scanpy"], _COMPARE_XENIUM, str(truth.bundle), str(out)) == []
+
diff --git a/tests/test_formats.py b/tests/test_formats.py
index 9dab72b..0dd61d6 100644
--- a/tests/test_formats.py
+++ b/tests/test_formats.py
@@ -476,9 +476,9 @@ def test_json_export_to_stdout(store, capsys):
def test_json_import_rejects_an_unrepresentable_value(store, temp_dir):
- from adata.formats.json_data import _write_json_to_group
+ from adata.formats.json_data import write_json_value
path, opener = store
with opener("a") as root:
with pytest.raises(ValueError, match="Cannot convert"):
- _write_json_to_group(root["uns"], "bad", {1, 2})
+ write_json_value(root["uns"], "bad", {1, 2})
diff --git a/tests/test_performance.py b/tests/test_performance.py
index 9951cad..d249a6d 100644
--- a/tests/test_performance.py
+++ b/tests/test_performance.py
@@ -980,6 +980,45 @@ def measure(n: int) -> int:
assert_independent_of(measure, what="convert CellBender write calls", axis="n_barcodes")
+def _spatial_bundle(path: Path, kind: str, n: int, shuffled: bool) -> Path:
+ from tests.conftest import write_visium_outs, write_xenium_bundle
+
+ if kind == "visium":
+ return write_visium_outs(path, n_spots=n, n_extra=n // 4).outs
+ return write_xenium_bundle(path, n_cells=n, shuffled=shuffled, extra_cells=n // 8 if shuffled else 0).bundle
+
+
+@pytest.mark.parametrize(
+ "kind,shuffled", [("visium", True), ("xenium", False), ("xenium", True)]
+)
+def test_spatial_import_work_grows_linearly(tmp_path, kind, shuffled):
+ """The positions or cell table are matched to the matrix while scanned,
+ and written straight through when in its order (Xenium's own) or through
+ the scatter buckets when not: linear either way, never a lookup per row."""
+ from adata.commands.convert import convert_store
+
+ def measure(n: int) -> int:
+ source = _spatial_bundle(tmp_path / f"sp{kind}{shuffled}{n}", kind, n, shuffled)
+ with count_io() as io:
+ convert_store(source, tmp_path / f"osp{kind}{shuffled}{n}.h5ad", QUIET)
+ return io.work
+
+ assert_grows_linearly(measure, what=f"convert {kind} (shuffled={shuffled})", axis="n_cells")
+
+
+@pytest.mark.parametrize("shuffled", [False, True])
+def test_xenium_cells_are_written_in_blocks_not_per_cell(tmp_path, shuffled):
+ from adata.commands.convert import convert_store
+
+ def measure(n: int) -> int:
+ source = _spatial_bundle(tmp_path / f"xw{shuffled}{n}", "xenium", n, shuffled)
+ with count_io() as io:
+ convert_store(source, tmp_path / f"oxw{shuffled}{n}.h5ad", QUIET)
+ return io.h5_write_calls
+
+ assert_independent_of(measure, what=f"convert Xenium write calls (shuffled={shuffled})", axis="n_cells")
+
+
def test_split_reads_grow_linearly_in_group_count(tmp_path, monkeypatch):
"""Splitting into k groups must stay linear in k at fixed store size.
diff --git a/tests/test_scatter.py b/tests/test_scatter.py
new file mode 100644
index 0000000..3f8dbd5
--- /dev/null
+++ b/tests/test_scatter.py
@@ -0,0 +1,74 @@
+"""Tests for writing rows in another order through scratch (`core/scatter.py`)."""
+
+from __future__ import annotations
+
+import numpy as np
+import pytest
+
+from adata.core.scatter import RowScatter
+from adata.storage import create_dataset
+
+
+def _feed(column, values, step):
+ for lo in range(0, len(values), step):
+ column.add(values[lo : lo + step])
+
+
+@pytest.mark.parametrize("row_step", [3, 1 << 16])
+def test_every_kept_row_lands_on_its_target_and_the_rest_keep_the_fill(new_store, row_step):
+ path, opener = new_store()
+ rng = np.random.default_rng(0)
+ n_source, n_out = 40, 30
+ # A permutation of 25 output rows, five source rows dropped, five output
+ # rows left unreached.
+ target = np.full(n_source, -1, dtype=np.int64)
+ sources = rng.choice(n_source, size=25, replace=False)
+ target[sources] = rng.choice(n_out, size=25, replace=False)
+ numbers = rng.random((n_source, 2))
+ codes = rng.integers(0, 4, n_source)
+
+ with opener("a") as root:
+ a = create_dataset(root["uns"], "a", shape=(n_out, 2), dtype=np.float64)
+ b = create_dataset(root["uns"], "b", shape=(n_out,), dtype=np.int64)
+ with RowScatter(target, n_out, root, row_step=row_step) as scatter:
+ first, second = scatter.column(a, np.nan), scatter.column(b, -1)
+ # Fed together, in one pass, in blocks that do not line up with
+ # the output blocks.
+ for lo in range(0, n_source, 7):
+ first.add(numbers[lo : lo + 7])
+ second.add(codes[lo : lo + 7])
+
+ with opener("r") as root:
+ got_a = np.asarray(root["uns"]["a"][...])
+ got_b = np.asarray(root["uns"]["b"][...])
+ assert not any(k.startswith("__scatter") for k in root.keys())
+
+ expected_a = np.full((n_out, 2), np.nan)
+ expected_b = np.full(n_out, -1)
+ expected_a[target[sources]] = numbers[sources]
+ expected_b[target[sources]] = codes[sources]
+ np.testing.assert_array_equal(got_a, expected_a)
+ np.testing.assert_array_equal(got_b, expected_b)
+
+
+def test_a_column_fed_fewer_rows_than_planned_is_refused(new_store):
+ path, opener = new_store()
+ with opener("a") as root:
+ dest = create_dataset(root["uns"], "a", shape=(3,), dtype=np.float64)
+ scatter = RowScatter(np.array([2, 0, 1]), 3, root)
+ _feed(scatter.column(dest, np.nan), np.arange(2.0), 1)
+ with pytest.raises(ValueError, match="2 source rows were given, where 3"):
+ scatter.finish()
+
+
+def test_a_layout_that_keeps_no_row_writes_only_the_fill(new_store):
+ """HDF5 refused the scratch's chunk shape, larger than its empty data,
+ when no source row was kept: an analysis table matching no cell, or an
+ empty matrix beside a table, crashed the write on h5ad alone."""
+ path, opener = new_store()
+ with opener("a") as root:
+ dest = create_dataset(root["uns"], "a", shape=(4,), dtype=np.float64)
+ with RowScatter(np.full(3, -1), 4, root) as scatter:
+ _feed(scatter.column(dest), np.arange(3.0), 2)
+ with opener("r") as root:
+ assert np.isnan(np.asarray(root["uns"]["a"][...])).all()
diff --git a/tests/test_table.py b/tests/test_table.py
index 491081f..60672a1 100644
--- a/tests/test_table.py
+++ b/tests/test_table.py
@@ -162,3 +162,26 @@ def test_natural_order_matches_natsort():
"", "A", "A01b", "A1b", "GRCh38", "GRCh37", "mm10", "s10_r2", "s2_r10",
]
assert natural_sorted(values) == natsort.natsorted(sorted(values))
+
+
+def test_a_file_with_no_header_line_reads_with_the_names_given(tmp_path):
+ """Space Ranger 1 writes `tissue_positions_list.csv` with no header."""
+ body = "AAAC-1,1,0,0,1251,1210\nAAAG-1,0,1,1,1260.5,-30\nAAAT-1,1,2,0,1270,1300\n"
+ names = ["barcode", "in_tissue", "array_row", "array_col", "pxl_row_in_fullres", "pxl_col_in_fullres"]
+ path = _write(tmp_path, body)
+
+ got = _read_back(tmp_path, scan_table(path, names=names, block_bytes=16))
+
+ expected = pd.read_csv(path, header=None, names=names).set_index("barcode")
+ expected.index.name = None
+ _same(got, expected)
+ assert table.read_small_table(path, names=names)[1][0][0] == "AAAC-1"
+
+
+def test_on_index_sees_every_index_value_once_and_in_order(tmp_path):
+ path = _write(tmp_path, CSV)
+ seen = []
+
+ scan_table(path, on_index=seen.extend, block_bytes=24)
+
+ assert seen == ["c1", "c2", "c3", "c4"]
diff --git a/tests/test_visium.py b/tests/test_visium.py
new file mode 100644
index 0000000..d375e86
--- /dev/null
+++ b/tests/test_visium.py
@@ -0,0 +1,266 @@
+"""Tests for `adata convert` of Space Ranger output for standard Visium.
+
+Fixtures are written as Space Ranger writes its `outs/` folder (see
+conftest): positions shuffled and including spots off the matrix, Space
+Ranger 1's headerless file and the later headed one. What anndata reads
+back is checked against what each fixture was built from, laid out as
+`squidpy.read.visium` lays it out.
+"""
+
+from __future__ import annotations
+
+import json
+import re
+from pathlib import Path
+
+import h5py
+import numpy as np
+import pytest
+from typer.testing import CliRunner
+
+from adata.cli import app
+from tests.conftest import write_visium_outs
+
+ad = pytest.importorskip("anndata")
+
+runner = CliRunner()
+_ANSI = re.compile(r"\x1b\[[0-9;]*[a-zA-Z]")
+
+
+def _out(result) -> str:
+ text = result.stdout + (result.stderr or "")
+ return " ".join(_ANSI.sub("", text).split())
+
+
+def _read(path: Path):
+ return ad.read_zarr(path) if path.suffix == ".zarr" else ad.read_h5ad(path)
+
+
+def _convert(*args):
+ return runner.invoke(app, ["convert", *map(str, args)])
+
+
+def _ok(*args):
+ result = _convert(*args)
+ assert result.exit_code == 0, _out(result)
+ return result
+
+
+def _refused(*args) -> str:
+ result = _convert(*args)
+ assert result.exit_code == 1, _out(result)
+ return _out(result)
+
+
+# ---------------------------------------------------------------------------
+# the mapping
+
+
+@pytest.mark.parametrize("version", ["sr1", "sr2", "sr3_float"])
+@pytest.mark.parametrize("suffix", [".h5ad", ".zarr"])
+def test_spots_images_and_scale_factors_land_where_squidpy_puts_them(tmp_path, version, suffix):
+ truth = write_visium_outs(tmp_path / "run", version=version)
+ out = tmp_path / f"out{suffix}"
+
+ result = _ok(truth.outs, "-o", out)
+
+ assert "Space Ranger (Visium)" in _out(result)
+ adata = _read(out)
+ assert list(adata.obs_names) == truth.barcodes
+ assert list(adata.obs.columns) == ["in_tissue", "array_row", "array_col"]
+ for column in ("in_tissue", "array_row", "array_col"):
+ assert adata.obs[column].dtype == np.int64, column
+ np.testing.assert_array_equal(adata.obs[column], getattr(truth, column))
+ spatial = adata.obsm["spatial"]
+ assert spatial.dtype == (np.float64 if version == "sr3_float" else np.int64)
+ np.testing.assert_array_equal(spatial, truth.spatial)
+
+ library = adata.uns["spatial"]["V1_Test_Section"]
+ for name, pixels in truth.images.items():
+ assert library["images"][name].dtype == np.uint8
+ np.testing.assert_array_equal(library["images"][name], pixels)
+ assert library["scalefactors"] == truth.scalefactors
+ assert library["metadata"] == {
+ "chemistry_description": "Spatial 3' v1",
+ "software_version": "spaceranger-2.1.0",
+ }
+ np.testing.assert_array_equal(adata.X.toarray(), truth.counts)
+ assert "cellranger" in adata.uns
+
+
+def test_the_run_folder_holding_outs_is_read_too(tmp_path):
+ truth = write_visium_outs(tmp_path / "run")
+ _ok(tmp_path / "run", "-o", tmp_path / "out.h5ad")
+ assert list(_read(tmp_path / "out.h5ad").obs_names) == truth.barcodes
+
+
+def test_unfiltered_reads_every_spot_with_its_position(tmp_path):
+ truth = write_visium_outs(tmp_path / "run", n_spots=8, n_extra=4)
+ _ok(truth.outs, "-o", tmp_path / "out.h5ad", "--unfiltered")
+ adata = _read(tmp_path / "out.h5ad")
+ assert adata.n_obs == 12
+ assert adata.obs["in_tissue"].tolist() == [1] * 8 + [0] * 4
+
+
+def test_a_matrix_given_directly_takes_its_folder_as_the_output(tmp_path):
+ truth = write_visium_outs(tmp_path / "run")
+ _ok(
+ truth.outs / "filtered_feature_bc_matrix.h5", "-o", tmp_path / "out.h5ad",
+ "--from", "visium",
+ )
+ assert "spatial" in _read(tmp_path / "out.h5ad").obsm
+
+
+def test_a_10x_download_with_its_sample_prefix_is_found(tmp_path):
+ truth = write_visium_outs(tmp_path / "run", raw=False)
+ (truth.outs / "filtered_feature_bc_matrix.h5").rename(
+ truth.outs / "V1_Mouse_Brain_filtered_feature_bc_matrix.h5"
+ )
+ _ok(truth.outs, "-o", tmp_path / "out.h5ad")
+ assert list(_read(tmp_path / "out.h5ad").obs_names) == truth.barcodes
+
+
+def test_library_id_and_no_images(tmp_path):
+ truth = write_visium_outs(tmp_path / "run")
+ _ok(truth.outs, "-o", tmp_path / "out.h5ad", "--library-id", "brain", "--no-images")
+ library = _read(tmp_path / "out.h5ad").uns["spatial"]["brain"]
+ assert "images" not in library
+ assert library["scalefactors"] == truth.scalefactors
+
+
+def test_a_missing_image_is_a_warning_not_a_refusal(tmp_path):
+ truth = write_visium_outs(tmp_path / "run")
+ (truth.outs / "spatial" / "tissue_hires_image.png").unlink()
+ result = _ok(truth.outs, "-o", tmp_path / "out.h5ad")
+ assert "no hires image" in _out(result)
+ assert list(_read(tmp_path / "out.h5ad").uns["spatial"]["V1_Test_Section"]["images"]) == ["lowres"]
+
+
+def test_split_feature_types_puts_the_positions_in_every_output(tmp_path):
+ truth = write_visium_outs(tmp_path / "run")
+ _ok(truth.outs, "-o", tmp_path / "out.h5ad", "--split-feature-types")
+ for name in ("out.Gene_Expression.h5ad", "out.Antibody_Capture.h5ad"):
+ adata = _read(tmp_path / name)
+ np.testing.assert_array_equal(adata.obsm["spatial"], truth.spatial)
+ assert "V1_Test_Section" in adata.uns["spatial"]
+
+
+# ---------------------------------------------------------------------------
+# detection
+
+
+def test_an_outs_folder_is_visium_not_a_cell_ranger_2_mtx_layout(tmp_path):
+ """A Space Ranger outs/ holds filtered_ and raw_feature_bc_matrix/ MEX
+ folders, which 10x-mtx took for a Cell Ranger 2 file's genome folders
+ and refused with "pick one with --genome"."""
+ truth = write_visium_outs(tmp_path / "run")
+ for name in ("filtered_feature_bc_matrix", "raw_feature_bc_matrix"):
+ mex = truth.outs / name
+ mex.mkdir()
+ for file in ("matrix.mtx.gz", "barcodes.tsv.gz", "features.tsv.gz"):
+ (mex / file).write_bytes(b"")
+
+ result = _ok(truth.outs, "-o", tmp_path / "out.h5ad")
+
+ assert "Space Ranger (Visium)" in _out(result)
+
+
+def test_a_bare_matrix_is_still_a_plain_10x_matrix(tmp_path):
+ truth = write_visium_outs(tmp_path / "run")
+ result = _ok(truth.outs / "filtered_feature_bc_matrix.h5", "-o", tmp_path / "out.h5ad")
+ assert "10x h5" in _out(result)
+ assert "spatial" not in _read(tmp_path / "out.h5ad").obsm
+
+
+# ---------------------------------------------------------------------------
+# refusals, each before anything is written
+
+
+def test_a_barcode_with_no_position_is_refused(tmp_path):
+ truth = write_visium_outs(tmp_path / "run")
+ positions = truth.outs / "spatial" / "tissue_positions.csv"
+ lines = positions.read_text().splitlines()
+ positions.write_text("\n".join(l for l in lines if not l.startswith(truth.barcodes[3])) + "\n")
+
+ message = _refused(truth.outs, "-o", tmp_path / "out.h5ad")
+
+ assert f"1 of the matrix's 12 barcodes have no row in spatial/tissue_positions.csv (the first is '{truth.barcodes[3]}')" in message
+ assert not list(tmp_path.glob("out*"))
+
+
+def test_a_repeated_position_is_refused(tmp_path):
+ truth = write_visium_outs(tmp_path / "run")
+ positions = truth.outs / "spatial" / "tissue_positions.csv"
+ lines = positions.read_text().splitlines()
+ positions.write_text("\n".join(lines + [lines[1]]) + "\n")
+ assert "more than one row for 1 of the matrix's barcodes" in _refused(truth.outs, "-o", tmp_path / "out.h5ad")
+
+
+def test_visium_hd_is_refused_by_name(tmp_path):
+ outs = tmp_path / "hd" / "outs"
+ (outs / "binned_outputs" / "square_008um").mkdir(parents=True)
+ message = _refused(outs, "-o", tmp_path / "out.h5ad")
+ assert "Visium HD output, which is not read yet" in message
+
+
+def test_aggr_output_is_refused(tmp_path):
+ truth = write_visium_outs(tmp_path / "run")
+ with h5py.File(truth.outs / "filtered_feature_bc_matrix.h5", "a") as f:
+ f.attrs["library_ids"] = np.array([b"a", b"b"])
+ assert "holds 2 libraries (a, b)" in _refused(truth.outs, "-o", tmp_path / "out.h5ad")
+
+
+def test_a_bundle_with_only_mex_folders_points_at_10x_mtx(tmp_path):
+ truth = write_visium_outs(tmp_path / "run", raw=False)
+ (truth.outs / "filtered_feature_bc_matrix.h5").unlink()
+ (truth.outs / "filtered_feature_bc_matrix").mkdir()
+ message = _refused(truth.outs, "-o", tmp_path / "out.h5ad")
+ assert "Convert the folder with --from 10x-mtx" in message
+
+
+@pytest.mark.parametrize("flag", [["--unfiltered"], ["--library-id", "x"], ["--no-images"]])
+def test_the_visium_flags_are_refused_for_other_formats(tmp_path, flag):
+ truth = write_visium_outs(tmp_path / "run")
+ message = _refused(truth.outs / "filtered_feature_bc_matrix.h5", "-o", tmp_path / "out.h5ad", *flag)
+ assert "does not apply when converting 10x-h5 to h5ad" in message
+
+
+def test_unfiltered_with_a_matrix_given_is_refused(tmp_path):
+ truth = write_visium_outs(tmp_path / "run")
+ message = _refused(
+ truth.outs / "filtered_feature_bc_matrix.h5", "-o", tmp_path / "out.h5ad",
+ "--from", "visium", "--unfiltered",
+ )
+ assert "is a matrix already" in message
+
+
+def test_a_positions_file_saved_with_a_byte_order_mark_keeps_its_header(tmp_path):
+ """Excel's CSV starts with a BOM, which made the first field
+ 'barcode': the header was taken for a data row, and the write
+ failed on 'could not convert string to float'."""
+ truth = write_visium_outs(tmp_path / "run")
+ positions = truth.outs / "spatial" / "tissue_positions.csv"
+ positions.write_bytes(b"\xef\xbb\xbf" + positions.read_bytes())
+ _ok(truth.outs, "-o", tmp_path / "out.h5ad")
+ np.testing.assert_array_equal(_read(tmp_path / "out.h5ad").obsm["spatial"], truth.spatial)
+
+
+def test_a_position_that_is_not_a_number_is_refused_before_writing(tmp_path):
+ """It used to pass the plan and fail in the write, after staging began."""
+ truth = write_visium_outs(tmp_path / "run")
+ positions = truth.outs / "spatial" / "tissue_positions.csv"
+ lines = positions.read_text().splitlines()
+ fields = lines[1].split(",")
+ fields[5] = "left"
+ positions.write_text("\n".join([lines[0], ",".join(fields), *lines[2:]]) + "\n")
+ message = _refused(truth.outs, "-o", tmp_path / "out.h5ad")
+ assert "pxl_col_in_fullres holds values that are not numbers" in message
+ assert not list(tmp_path.glob("out*"))
+
+
+@pytest.mark.parametrize("name", ["", "a/b"])
+def test_a_library_id_that_cannot_be_a_key_is_refused(tmp_path, name):
+ """An empty one stopped with a raw h5py error; one with '/' silently
+ nested as uns['spatial']['a']['b']."""
+ truth = write_visium_outs(tmp_path / "run")
+ assert "cannot be a key of uns" in _refused(truth.outs, "-o", tmp_path / "out.h5ad", "--library-id", name)
diff --git a/tests/test_xenium.py b/tests/test_xenium.py
new file mode 100644
index 0000000..0f3b1b8
--- /dev/null
+++ b/tests/test_xenium.py
@@ -0,0 +1,286 @@
+"""Tests for `adata convert` of Xenium output bundles.
+
+Fixtures are written as Xenium Onboard Analysis writes its bundle (see
+conftest): XOA 1's numbered cells and later versions' named ones, the cell
+table in the matrix's order or shuffled, with cells the matrix lacks, and
+analysis/ with gaps. What anndata reads back is checked against what each
+fixture was built from.
+"""
+
+from __future__ import annotations
+
+import re
+import shutil
+from pathlib import Path
+
+import numpy as np
+import pytest
+from typer.testing import CliRunner
+
+from adata.cli import app
+from tests.conftest import write_xenium_bundle
+
+ad = pytest.importorskip("anndata")
+pd = pytest.importorskip("pandas")
+
+runner = CliRunner()
+_ANSI = re.compile(r"\x1b\[[0-9;]*[a-zA-Z]")
+
+
+def _out(result) -> str:
+ text = result.stdout + (result.stderr or "")
+ return " ".join(_ANSI.sub("", text).split())
+
+
+def _read(path: Path):
+ return ad.read_zarr(path) if path.suffix == ".zarr" else ad.read_h5ad(path)
+
+
+def _convert(*args):
+ return runner.invoke(app, ["convert", *map(str, args)])
+
+
+def _ok(*args):
+ result = _convert(*args)
+ assert result.exit_code == 0, _out(result)
+ return result
+
+
+def _refused(*args) -> str:
+ result = _convert(*args)
+ assert result.exit_code == 1, _out(result)
+ return _out(result)
+
+
+def _check_cells(adata, truth) -> None:
+ assert list(adata.obs_names) == truth.cell_ids
+ for name, values in truth.columns.items():
+ got = adata.obs[name]
+ if name == "segmentation_method":
+ assert isinstance(got.dtype, pd.CategoricalDtype)
+ assert got.astype(str).tolist() == list(values)
+ else:
+ np.testing.assert_allclose(got.to_numpy(dtype=float), values.astype(float), err_msg=name)
+ np.testing.assert_allclose(adata.obsm["spatial"], truth.spatial)
+ np.testing.assert_array_equal(adata.X.toarray(), truth.counts)
+
+
+# ---------------------------------------------------------------------------
+# the mapping
+
+
+@pytest.mark.parametrize("ids", ["int", "string"])
+@pytest.mark.parametrize("suffix", [".h5ad", ".zarr"])
+def test_cells_centroids_and_analysis_land_where_spatialdata_io_puts_them(tmp_path, ids, suffix):
+ truth = write_xenium_bundle(tmp_path / "xenium", ids=ids)
+ out = tmp_path / f"out{suffix}"
+
+ result = _ok(truth.bundle, "-o", out)
+
+ assert "Xenium →" in _out(result)
+ assert "analysis: clusters, X_umap" in _out(result)
+ adata = _read(out)
+ _check_cells(adata, truth)
+ assert adata.obs["transcript_counts"].dtype == np.int64
+ assert list(adata.obs.columns) == [*truth.columns, "graphclust"]
+ assert isinstance(adata.obs["graphclust"].dtype, pd.CategoricalDtype)
+ assert adata.obs["graphclust"].astype(str).tolist() == list(truth.clusters)
+ np.testing.assert_allclose(adata.obsm["X_umap"], truth.umap)
+ assert adata.var["feature_types"].tolist()[-1] == "Negative Control Probe"
+
+
+def test_a_cell_table_in_another_order_is_written_in_the_matrix_s(tmp_path):
+ truth = write_xenium_bundle(tmp_path / "xenium", shuffled=True)
+ _ok(truth.bundle, "-o", tmp_path / "out.h5ad")
+ _check_cells(_read(tmp_path / "out.h5ad"), truth)
+
+
+def test_cells_the_matrix_lacks_are_left_out_with_their_categories(tmp_path):
+ truth = write_xenium_bundle(tmp_path / "xenium", extra_cells=4, shuffled=True)
+
+ result = _ok(truth.bundle, "-o", tmp_path / "out.h5ad")
+
+ assert "4 cells of cells.csv.gz are not in the matrix; they are left out" in _out(result)
+ adata = _read(tmp_path / "out.h5ad")
+ _check_cells(adata, truth)
+ assert "Imported cell segmentation" not in adata.obs["segmentation_method"].cat.categories
+
+
+def test_cells_missing_from_the_analysis_are_left_empty(tmp_path):
+ truth = write_xenium_bundle(tmp_path / "xenium", missing_in_analysis=3)
+
+ result = _ok(truth.bundle, "-o", tmp_path / "out.h5ad")
+
+ assert "3 cells have no clusters in analysis/" in _out(result)
+ adata = _read(tmp_path / "out.h5ad")
+ assert adata.obs["graphclust"].isna().sum() == 3
+ assert adata.obs["graphclust"].astype(object).fillna("").tolist() == list(truth.clusters)
+ np.testing.assert_allclose(adata.obsm["X_umap"], truth.umap)
+
+
+def test_no_analysis_is_a_warning(tmp_path):
+ truth = write_xenium_bundle(tmp_path / "xenium", analysis=False)
+ result = _ok(truth.bundle, "-o", tmp_path / "out.h5ad")
+ assert "No analysis/ or analysis.tar.gz" in _out(result)
+ adata = _read(tmp_path / "out.h5ad")
+ assert "graphclust" not in adata.obs and "X_umap" not in adata.obsm
+
+
+def test_the_analysis_is_read_from_the_archive_xenium_ships_it_in(tmp_path):
+ """Xenium packs analysis/ as analysis.tar.gz, so a bundle as downloaded
+ has no analysis/ folder at all."""
+ truth = write_xenium_bundle(tmp_path / "xenium", missing_in_analysis=2)
+ shutil.make_archive(str(truth.bundle / "analysis"), "gztar", truth.bundle, "analysis")
+ shutil.rmtree(truth.bundle / "analysis")
+
+ result = _ok(truth.bundle, "-o", tmp_path / "out.h5ad")
+
+ assert "analysis: clusters, X_umap" in _out(result)
+ adata = _read(tmp_path / "out.h5ad")
+ assert adata.obs["graphclust"].astype(object).fillna("").tolist() == list(truth.clusters)
+ np.testing.assert_allclose(adata.obsm["X_umap"], truth.umap)
+
+
+def test_an_archive_cannot_write_outside_its_temporary_folder(tmp_path):
+ """Members are copied to paths built from the names this wants, never
+ extracted where the archive says."""
+ import io
+ import tarfile
+
+ truth = write_xenium_bundle(tmp_path / "xenium", analysis=False)
+ with tarfile.open(truth.bundle / "analysis.tar.gz", "w:gz") as tar:
+ data = b"Barcode,Cluster\n"
+ info = tarfile.TarInfo("../../escaped.csv")
+ info.size = len(data)
+ tar.addfile(info, io.BytesIO(data))
+ _ok(truth.bundle, "-o", tmp_path / "out.h5ad")
+ assert not list(tmp_path.rglob("escaped.csv"))
+
+
+def test_feature_type_keeps_the_cells_and_their_positions(tmp_path):
+ truth = write_xenium_bundle(tmp_path / "xenium")
+ _ok(truth.bundle, "-o", tmp_path / "out.h5ad", "--feature-type", "Gene Expression")
+ adata = _read(tmp_path / "out.h5ad")
+ assert adata.n_vars == 4
+ np.testing.assert_allclose(adata.obsm["spatial"], truth.spatial)
+
+
+def test_a_bundle_is_xenium_not_a_mtx_layout(tmp_path):
+ """A Xenium bundle holds cell_feature_matrix/ as a MEX folder too, which
+ the 10x-mtx sniffer could take for the folder to read."""
+ truth = write_xenium_bundle(tmp_path / "xenium")
+ mex = truth.bundle / "cell_feature_matrix"
+ mex.mkdir()
+ for file in ("matrix.mtx.gz", "barcodes.tsv.gz", "features.tsv.gz"):
+ (mex / file).write_bytes(b"")
+ assert "Xenium →" in _out(_ok(truth.bundle, "-o", tmp_path / "out.h5ad"))
+
+
+# ---------------------------------------------------------------------------
+# refusals, each before anything is written
+
+
+def test_a_cell_missing_from_the_table_is_refused(tmp_path):
+ truth = write_xenium_bundle(tmp_path / "xenium")
+ path = truth.bundle / "cells.csv.gz"
+ import gzip
+
+ lines = gzip.decompress(path.read_bytes()).decode().splitlines()
+ path.write_bytes(gzip.compress(("\n".join(lines[:-1]) + "\n").encode()))
+
+ message = _refused(truth.bundle, "-o", tmp_path / "out.h5ad")
+
+ assert "1 of the matrix's 15 barcodes have no row in cells.csv.gz" in message
+ assert not list(tmp_path.glob("out*"))
+
+
+def test_a_bundle_with_only_parquet_cells_says_so(tmp_path):
+ truth = write_xenium_bundle(tmp_path / "xenium")
+ (truth.bundle / "cells.csv.gz").rename(truth.bundle / "cells.parquet")
+ assert "has cells.parquet only, which is not read yet" in _refused(truth.bundle, "-o", tmp_path / "out.h5ad")
+
+
+@pytest.mark.parametrize("flag", [["--unfiltered"], ["--library-id", "x"], ["--no-images"]])
+def test_visium_s_flags_do_not_apply(tmp_path, flag):
+ truth = write_xenium_bundle(tmp_path / "xenium")
+ message = _refused(truth.bundle, "-o", tmp_path / "out.h5ad", *flag)
+ assert "does not apply when converting xenium to h5ad" in message
+
+
+@pytest.mark.parametrize("shuffled", [False, True])
+def test_a_table_in_the_matrix_s_order_is_written_straight_through(tmp_path, shuffled):
+ """Xenium writes cells.csv.gz in the matrix's order; only a table in
+ another order should pay for the scatter's scratch copy."""
+ from adata.conversion import Options
+ from adata.conversion.xenium import XeniumImport
+
+ truth = write_xenium_bundle(tmp_path / "xenium", shuffled=shuffled)
+ plan = XeniumImport().plan(truth.bundle, tmp_path / "out.h5ad", "xenium", "h5ad", Options())
+ assert (plan.detail.cells.join.target is None) is not shuffled
+
+
+def _rewrite(path: Path, edit) -> None:
+ path.write_text(edit(path.read_text()))
+
+
+def test_clusters_matching_no_cell_are_refused_as_another_run_s(tmp_path):
+ """Planned as 'every cell left empty', then crashed the h5ad write."""
+ truth = write_xenium_bundle(tmp_path / "xenium")
+ clusters = truth.bundle / "analysis/clustering/gene_expression_graphclust/clusters.csv"
+ _rewrite(clusters, lambda text: text.replace("-1,", "-2,"))
+ message = _refused(truth.bundle, "-o", tmp_path / "out.h5ad")
+ assert "has a row for none of the matrix's 15 barcodes" in message
+
+
+def test_a_clusters_table_without_its_column_is_refused_before_writing(tmp_path):
+ """It used to fail inside the write, after staging began."""
+ truth = write_xenium_bundle(tmp_path / "xenium")
+ clusters = truth.bundle / "analysis/clustering/gene_expression_graphclust/clusters.csv"
+ _rewrite(clusters, lambda text: text.replace("Barcode,Cluster", "Barcode,cluster", 1))
+ assert "clusters.csv in analysis has no Cluster column" in _refused(truth.bundle, "-o", tmp_path / "out.h5ad")
+ assert not list(tmp_path.glob("out*"))
+
+
+def test_a_missing_centroid_is_nan_and_a_word_is_refused(tmp_path):
+ import gzip
+
+ truth = write_xenium_bundle(tmp_path / "xenium")
+ path = truth.bundle / "cells.csv.gz"
+ lines = gzip.decompress(path.read_bytes()).decode().splitlines()
+ fields = lines[1].split(",")
+ fields[1] = "NA"
+ path.write_bytes(gzip.compress(("\n".join([lines[0], ",".join(fields), *lines[2:]]) + "\n").encode()))
+ _ok(truth.bundle, "-o", tmp_path / "out.h5ad")
+ assert np.isnan(_read(tmp_path / "out.h5ad").obsm["spatial"][0, 0])
+
+ fields[1] = "left"
+ path.write_bytes(gzip.compress(("\n".join([lines[0], ",".join(fields), *lines[2:]]) + "\n").encode()))
+ assert "x_centroid holds values that are not numbers" in _refused(truth.bundle, "-o", tmp_path / "out2.h5ad")
+
+
+def test_clusters_of_cells_the_matrix_lacks_leave_no_category(tmp_path):
+ truth = write_xenium_bundle(tmp_path / "xenium")
+ clusters = truth.bundle / "analysis/clustering/gene_expression_graphclust/clusters.csv"
+ _rewrite(clusters, lambda text: text + "zzzzzzzz-1,99\n")
+ _ok(truth.bundle, "-o", tmp_path / "out.h5ad")
+ assert "99" not in _read(tmp_path / "out.h5ad").obs["graphclust"].cat.categories
+
+
+def test_an_archive_without_the_tables_says_so(tmp_path):
+ """Its empty unpacked folder passed for an analysis folder, and the
+ import went ahead without a word."""
+ import tarfile
+
+ truth = write_xenium_bundle(tmp_path / "xenium")
+ shutil.rmtree(truth.bundle / "analysis")
+ with tarfile.open(truth.bundle / "analysis.tar.gz", "w:gz") as tar:
+ tar.add(truth.bundle / "experiment.xenium", arcname="outs/other.txt")
+ assert "holds none of the clustering, UMAP or PCA tables" in _out(_ok(truth.bundle, "-o", tmp_path / "out.h5ad"))
+ assert not list(tmp_path.glob(".adata-xenium-*"))
+
+
+def test_from_xenium_names_the_matrix_not_another_file(tmp_path):
+ """Any file was taken for its folder, and the matrix beside it read."""
+ truth = write_xenium_bundle(tmp_path / "xenium")
+ message = _refused(truth.bundle / "cells.csv.gz", "-o", tmp_path / "out.h5ad", "--from", "xenium")
+ assert "give its folder, or its cell_feature_matrix.h5" in message