Merge pull request #51 from cellgeni/feat/convert-cellbender

convert: CellBender remove-background output, with its latents

Author: Claptar

AGENTS.md

@@ -93,7 +93,10 @@ Each layer calls only the layers below it:
 
   Content sniffing runs before extensions, because extensions lie (`.h5`).
   A format with no extension, such as an MTX folder (`10x-mtx`), is chosen as
-  an output only with `--to`.
+  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.
 - **`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

README.md

@@ -8,7 +8,7 @@ A command-line tool for exploring huge AnnData stores (`.h5ad` and `.zarr`) with
 - Auto-detects `.h5ad` files vs `.zarr` directories
 - Chunked processing for dense and sparse matrices (CSR/CSC)
 - Reads every AnnData on-disk layout, from 0.7.x through the current spec, and always writes the current one
-- Converts between HDF5 and Zarr (v2 and v3) in either direction with `adata convert`, and to and from 10x Genomics HDF5 matrices and 10x / STARsolo MTX folders
+- Converts between HDF5 and Zarr (v2 and v3) in either direction with `adata convert`, and to and from 10x Genomics HDF5 matrices and 10x / STARsolo MTX folders, and from CellBender output
 - Rich terminal output with progress indicators, kept on stderr so results pipe cleanly
 
 **Documentation: [cellgeni.github.io/adata-cli](https://cellgeni.github.io/adata-cli/)**
@@ -56,7 +56,7 @@ Run help at any level (e.g. `adata --help`, `adata export --help`).
 - `subset` – stream and write a filtered copy, selected by obs/var name lists (`--obs`/`--var`) or by expression (`--obs-query`/`--var-query`).
 - `split` – write one store per distinct value of an annotation column, with a CSV manifest.
 - `concat` – concatenate stores along the obs axis, with `--join inner|outer` and merge strategies for var and uns.
-- `convert` – change a store's format: h5ad ↔ Zarr v2/v3, 10x Genomics HDF5 (Cell Ranger's `*_feature_bc_matrix.h5`) and 10x / STARsolo MTX folders in and out, streaming, with the output written beside its destination and moved into place only once complete.
+- `convert` – change a store's format: h5ad ↔ Zarr v2/v3, 10x Genomics HDF5 (Cell Ranger's `*_feature_bc_matrix.h5`) and 10x / STARsolo MTX folders in and out, CellBender output in, streaming, with the output written beside its destination and moved into place only once complete.
 - `recast` – change a matrix's dtype, layout (CSR/CSC/dense) or density, streaming; refuses a lossy cast or a large size increase unless forced. Called `convert` until 0.7.0.
 - `export` – extract data from a store; subcommands: `dataframe` (any dataframe group to CSV), `array` (dense to `.npy`), `sparse` (CSR/CSC to `.mtx`), `dict` (JSON), `image` (PNG). Results go to stdout when no `--output` is given.
 - `import` – write new data into a store at any path; subcommands: `dataframe` (CSV), `array` (`.npy`), `sparse` (`.mtx`), `dict` (JSON), `image` (PNG/JPEG/TIFF).

docs/COMMANDS.md

@@ -122,6 +122,7 @@ adata convert counts.h5ad -o counts_mtx --to 10x-mtx
 | `h5ad` | `.h5ad` | yes | yes |
 | `zarr` | `.zarr` (v2 or v3) | yes | yes |
 | `10x-h5` | `.h5` | Cell Ranger 2 and 3+, Space Ranger, Xenium | Cell Ranger 3+ layout |
+| `cellbender` | `.h5` | CellBender `remove-background` 0.2 and later | no |
 | `10x-mtx` | a folder (`--to 10x-mtx`) | Cell Ranger 2 and 3+, STARsolo (velocyto too) | Cell Ranger 3+ layout, gzipped |
 
 Formats other than h5ad and Zarr convert to and from AnnData only: to turn a
@@ -196,6 +197,42 @@ is refused, and `--force` writes the names as ids and `Gene Expression` as
 the type. A matrix of floats is written as it is, with a warning. A file
 imported by `convert` exports with its run metadata and tag fields intact.
 
+### CellBender
+
+The output of CellBender's `remove-background` (`output.h5` and
+`output_filtered.h5`, 0.2 and later) is recognised by its latent variables,
+whatever its name. Its matrix is a 10x one and is read as above; the latents
+are kept, under the names CellBender's own loaders give them, so code written
+against those keeps working:
+
+| CellBender | AnnData |
+|---|---|
+| denoised counts | `X` |
+| per-droplet latents (`cell_probability`, `cell_size`, `background_fraction`, `droplet_efficiency`; `latent_*` in 0.2) | obs columns, NaN for droplets CellBender did not analyse |
+| which droplets were analysed | `obs["cellbender_analyzed"]` |
+| `gene_expression_encoding` (`latent_gene_encoding` in 0.2) | obsm, same name |
+| `ambient_expression` | `var["ambient_expression"]` |
+| which features were analysed (`features_analyzed_inds`, 0.3+) | `var["cellbender_analyzed"]` |
+| everything else (fitted constants, learning curve, settings) | `uns["cellbender"]` |
+
+```bash
+adata convert output.h5 -o cellbender.h5ad
+adata convert output.h5 -o cells.h5ad --min-cell-probability 0.5
+adata convert output_filtered.h5 -o cb.h5ad --raw raw_feature_bc_matrix.h5
+```
+
+| Flag | Meaning |
+|---|---|
+| `--min-cell-probability` | Keep droplets whose `cell_probability` is at least this; CellBender's own cutoff is 0.5 |
+| `--raw` | Add the raw Cell Ranger matrix CellBender ran on as `layers["cellranger"]`, and the denoised counts as `layers["cellbender"]`, as `load_anndata_from_input_and_output` does |
+
+The full `output.h5` keeps every droplet, as `read_10x_h5` would. `--raw`
+matches barcodes, so it works for `output_filtered.h5` too, whose cells
+CellBender orders by count; the raw file must list the same features in the
+same order, and every output barcode, or it is refused. The 10x flags above
+apply as well. `--from 10x-h5` reads a CellBender file as a plain 10x matrix,
+without the latents.
+
 ### 10x Genomics and STARsolo MTX folders
 
 A folder holding `matrix.mtx`, `barcodes.tsv` and `features.tsv` (or

docs/TESTING.md

@@ -28,6 +28,7 @@ separate jobs.
 | `test_recast.py` | `recast` against scipy: casts, transposes, density changes and every refusal |
 | `test_tenx.py` | 10x HDF5 import and export against Cell Ranger-shaped fixtures (v2, v3, multimodal, probe-based): mapping, filters, round trips, refusals |
 | `test_mtx.py` | 10x / STARsolo MTX folders against hand-written fixtures (Cell Ranger 2 and 3, STARsolo, velocyto, GEO prefixes): mapping, entry orders, refusals that leave nothing, round trips |
+| `test_cellbender.py` | CellBender 0.2 and 0.3+ output, full and filtered: each latent at its barcode, `--min-cell-probability`, `--raw` alignment, refusals |
 | `test_mtx_codec.py`, `test_csr_build.py` | The streaming Matrix Market reader and writer against scipy, and the bucket sort they and `recast` share |
 | `test_external_readers.py` | 10x conversions, HDF5 and MTX, checked by scanpy and PyTables in their own uv environments (`integration`) |
 | `test_performance.py` | Complexity guards -- what an operation costs, not how long it takes (see below) |
@@ -245,6 +246,7 @@ should scale with, scale only that, and hold everything else fixed.
 | `convert` 10x import (whole, `--split-feature-types`), export | n_barcodes | linear |
 | `convert` 10x MTX import (in order, shuffled, `--split-feature-types`) | n_barcodes | linear; **flat** write calls |
 | `convert` to 10x MTX | n_obs | **flat** read calls |
+| `convert` CellBender (plain, `--min-cell-probability`, `--raw` full and filtered) | n_barcodes | linear; **flat** write calls |
 | `recast --dtype`, `--layout csc`, `--layout dense` | nnz, n_obs | linear; the streaming transpose's **peak memory** is bounded |
 
 The `slow` marker is on the three guards that build a 65,536-row store; they still gate

docs/index.md

@@ -70,7 +70,7 @@ in.h5ad -o out.zarr` converts a whole store, and `subset`, `split` and
 `convert` also reads and writes 10x Genomics HDF5 matrices (Cell Ranger's
 `*_feature_bc_matrix.h5`) and 10x / STARsolo Matrix Market folders,
 velocyto's included, matching what scanpy's `read_10x_h5` and `read_10x_mtx`
-read.
+read, and CellBender's `remove-background` output with its latent variables.
 
 | Element | Read | Write |
 |---|---|---|

src/adata/cli.py

@@ -522,6 +522,8 @@ def concat(
     "gene_ids_column": "--gene-ids-column",
     "gene_names_column": "--gene-names-column",
     "feature_types_column": "--feature-types-column",
+    "raw": "--raw",
+    "min_cell_probability": "--min-cell-probability",
 }
 
 #: Flags that only make sense for format conversion; with recast's flags,
@@ -645,6 +647,21 @@ def convert(
         "--feature-types-column",
         help="10x output: var column holding feature types (default feature_types)",
     ),
+    raw: Optional[Path] = typer.Option(
+        None,
+        "--raw",
+        help="CellBender input: the raw Cell Ranger .h5 it ran on, added as "
+        "layers['cellranger']",
+        exists=True,
+        dir_okay=False,
+        readable=True,
+    ),
+    min_cell_probability: Optional[float] = typer.Option(
+        None,
+        "--min-cell-probability",
+        help="CellBender input: keep droplets with cell_probability at least "
+        "this (CellBender's own cutoff is 0.5)",
+    ),
     # Deprecated: recast's flags, so `convert` keeps working as it did until
     # 1.0.0. Hidden from --help; the old matrix paths arrive in ctx.args.
     inplace: bool = typer.Option(False, "--inplace", hidden=True),
@@ -662,7 +679,8 @@ def convert(
     Genomics matrix (Cell Ranger's *_feature_bc_matrix.h5); --to 10x-mtx
     writes a Matrix Market folder. Between h5ad and zarr every element is
     copied, streamed, keeping the source's chunking and compression where the
-    target can express them. A 10x matrix is read as scanpy's read_10x_h5 or
+    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 Ranger and
     STARsolo, velocyto's included, are read in two streaming passes.
@@ -684,6 +702,8 @@ def convert(
         adata convert Solo.out/Gene/filtered -o counts.h5ad
         adata convert Solo.out/Velocyto/raw -o velocyto.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
     """
     entries = list(ctx.args)
     recast_flags = _given(ctx, _RECAST_FLAGS)
@@ -767,6 +787,8 @@ def convert(
                 gene_ids_column=gene_ids_column,
                 gene_names_column=gene_names_column,
                 feature_types_column=feature_types_column,
+                raw=raw,
+                min_cell_probability=min_cell_probability,
             ),
         )
     except Exception as e:

src/adata/conversion/__init__.py

@@ -17,7 +17,7 @@
 
 import h5py
 
-from adata.conversion import mtx, tenx
+from adata.conversion import cellbender, mtx, tenx
 from adata.conversion.anndata_copy import AnnDataCopy
 from adata.conversion.base import ConversionPlan, Converter, Format, Options
 from adata.storage import has_valid_anndata_root_attrs, is_zarr_path
@@ -66,6 +66,15 @@ def _sniff_h5ad(path: Path) -> bool:
         anndata=True,
         sniff=is_zarr_path,
     ),
+    # Before 10x-h5: a CellBender file is a 10x matrix too, and the first
+    # sniffer to match wins. `--from 10x-h5` still reads one as plain 10x.
+    "cellbender": Format(
+        name="cellbender",
+        description="CellBender remove-background output (0.2 and later)",
+        extensions=(".h5",),
+        sniff=cellbender.sniff,
+        importer=cellbender.CellbenderImport(),
+    ),
     "10x-h5": Format(
         name="10x-h5",
         description="10x Genomics feature-barcode matrix (Cell Ranger v2 or v3)",
@@ -136,10 +145,21 @@ def output_format(path: Path, requested: Optional[str] = None) -> str:
     path = Path(path)
     if requested is not None and requested not in FORMATS:
         raise ValueError(f"Unknown format {requested!r}; use one of: {_names()}.")
+    if requested is not None and not (
+        FORMATS[requested].anndata or FORMATS[requested].exporter is not None
+    ):
+        raise ValueError(f"{requested} can be read but not written.")
 
     suffix = path.suffix.lower()
+    # Only formats that can be written: `.h5` names CellBender output too,
+    # but an `.h5` output is a 10x matrix.
     by_extension = next(
-        (fmt.name for fmt in FORMATS.values() if suffix in fmt.extensions), None
+        (
+            fmt.name
+            for fmt in FORMATS.values()
+            if suffix in fmt.extensions and (fmt.anndata or fmt.exporter is not None)
+        ),
+        None,
     )
     if requested is None:
         if by_extension is None:

src/adata/conversion/base.py

@@ -35,6 +35,10 @@ class Options:
     gene_names_column: Optional[str] = None
     feature_types_column: Optional[str] = None
 
+    # Reading CellBender output.
+    raw: Optional[Path] = None
+    min_cell_probability: Optional[float] = None
+
 
 @dataclass
 class ConversionPlan: