convert: loom, import and export

AGENTS.md

@@ -102,7 +102,9 @@ Each layer calls only the layers below it:
   passes in, block writes out) shared by `import`/`export sparse` and the
   `10x-mtx` and `parse` converters; `textio.py` holds its gzip-aware block
   I/O, and `table.py` a streaming CSV reader that types columns as pandas
-  does (Parse's cell metadata; meant for #25's delimited tables too).
+  does (Parse's cell metadata; meant for #25's delimited tables too), plus
+  `write_text_column`, which writes any text column as anndata would
+  (categorical when values repeat); the loom import uses it.
 - **`util/`**: `staging.py` and `path.py` (below), plus `names.py`.
 
 ## Rules the code depends on

README.md

@@ -102,6 +102,7 @@ adata convert filtered_feature_bc_matrix.h5 -o counts.h5ad --feature-type "Gene
 adata convert counts.h5ad -o for_cellbender.h5 --layer counts
 adata convert Solo.out/Velocyto/raw -o velocyto.h5ad
 adata convert splitpipe_out -o parse.h5ad --unfiltered
+adata convert velocyto.loom -o velocyto.h5ad
 ```
 
 ## Documentation

docs/COMMANDS.md

@@ -106,6 +106,8 @@ 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 splitpipe_out -o parse.h5ad
+adata convert velocyto.loom -o velocyto.h5ad
+adata convert counts.h5ad -o counts.loom
 ```
 
 | Flag | Meaning |
@@ -126,6 +128,7 @@ adata convert splitpipe_out -o parse.h5ad
 | `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 |
 | `parse` | a folder | Parse Biosciences split-pipe: a run, a sample or a `DGE_*` folder | no |
+| `loom` | `.loom` | loompy, velocyto, kallisto\|bustools; spec 2 and 3 | spec 3.0.0, as loompy 3 writes it |
 
 Formats other than h5ad and Zarr convert to and from AnnData only: to turn a
 10x file into another format, convert it to `.h5ad` first.
@@ -343,6 +346,58 @@ Both files are read in two streaming passes, as for 10x MTX folders, so a
 run's size is not a limit: the first checks the matrix against both tables
 before any output exists, and the second writes obs and X.
 
+### loom
+
+A loom holds a dense genes x cells matrix, layers of the same shape, an
+attribute per cell (`col_attrs`) and per gene (`row_attrs`), graphs, and
+global attributes. velocyto writes its counts this way, as do loompy and
+kallisto|bustools.
+
+```bash
+adata convert velocyto.loom -o velocyto.h5ad
+adata convert velocyto.loom -o unspliced.h5ad --layer unspliced
+adata convert scope.loom -o dense.h5ad --dense --var-names gene_symbols
+adata convert counts.h5ad -o counts.loom
+```
+
+| loom | AnnData |
+|---|---|
+| `/matrix`, or `/layers/spliced` when there is one | `X`, CSR (`--dense` keeps it dense) |
+| `/layers/<name>`, and `/matrix` when X is `spliced` | `layers[<name>]`, `layers["matrix"]` |
+| `col_attrs` with one value per cell / several | obs columns / obsm |
+| `row_attrs`, the same way | var columns / varm |
+| `col_graphs`, `row_graphs` | obsp, varp |
+| global attributes | `uns["loom"]["attrs"]` |
+
+| Flag | Meaning |
+|---|---|
+| `--layer` | The loom layer that becomes X: `matrix` (the main one) or a layer name |
+| `--dense` | Keep X and the layers dense |
+| `--obs-names-attr` | The attribute naming cells (default `CellID`, else `cell_id`, `obs_names` or `barcode`) |
+| `--var-names-attr` | The attribute naming genes (default the ids: `Accession` or `gene_id`) |
+
+Reading gives what `anndata.read_loom` gives, checked against it in CI, with
+three differences: values keep their dtype (velocyto's integer counts stay
+integers, where anndata makes X float32); graphs and global attributes are
+kept, where anndata drops them; and genes are named by id (`Accession`),
+the symbols (`Gene`) staying a column, unless `--var-names gene_symbols` is
+given. Strings are read as loompy reads them, with XML references such as
+`&#233;` turned back into characters. SCope's compound attributes
+(`RegulonsAUC`) are left out with a warning; #27 will read them.
+
+The matrix is read a slab of cells at a time, a whole number of the file's
+64-cell chunks, so each chunk is decompressed once and memory is set by the
+slab, not the file.
+
+Writing produces a spec 3.0.0 loom that loompy reads and validates: X (or
+the matrix `--layer` names) as `/matrix`, every layer, obs and var as
+attributes (`CellID`; `Gene` for symbols and `Accession` for ids), 2-D obsm
+and varm, obsp and varp as graphs, and uns scalars as global attributes.
+Strings are written as loompy writes them, ASCII with XML references. A
+loom's matrices are dense by definition: the summary gives their size before
+compression, and above 10 GB a warning says it will take a while. Chunks
+holding only zeros compress to almost nothing.
+
 ## `recast`
 
 Change a matrix's dtype, layout or density. Counts held as float64 cost twice

docs/TESTING.md

@@ -30,9 +30,10 @@ separate jobs.
 | `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_parse.py` | Parse split-pipe runs, as v1.6.1 writes them: run, sample and DGE folders, `--sample`, `--all-samples`, `--unfiltered`, `--genome`, refusals |
+| `test_loom.py` | loom import and export against loompy- and velocyto-shaped fixtures (spec 2 and 3): X choice, names, attributes, graphs, `--dense`, round trips |
 | `test_table.py` | The streaming CSV reader against pandas: column types, quoting, row subsets, blocks, 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, and Parse runs, checked by scanpy (Parse's own recipe for Parse) and PyTables in their own uv environments (`integration`) |
+| `test_external_readers.py` | 10x conversions, HDF5 and MTX, Parse runs and looms, checked by scanpy (Parse's own recipe for Parse), anndata with loompy, and PyTables in their own uv environments (`integration`) |
 | `test_performance.py` | Complexity guards -- what an operation costs, not how long it takes (see below) |
 | `test_benchmark_harness.py` | That the benchmark's measurement is trustworthy |
 | `test_subset.py`, `test_export.py`, `test_import.py`, `test_info_read.py`, `test_zarr.py`, `test_query.py` | Per-feature unit tests |

docs/index.md

@@ -71,8 +71,8 @@ in.h5ad -o out.zarr` converts a whole store, and `subset`, `split` and
 `*_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 (features named by gene id unless `--var-names gene_symbols`),
-CellBender's `remove-background` output with its latent variables, and
-Parse Biosciences split-pipe runs.
+CellBender's `remove-background` output with its latent variables, Parse
+Biosciences split-pipe runs, and looms (velocyto's included), both ways.
 
 | Element | Read | Write |
 |---|---|---|

src/adata/cli.py

@@ -527,6 +527,9 @@ def concat(
     "sample": "--sample",
     "all_samples": "--all-samples",
     "unfiltered": "--unfiltered",
+    "dense": "--dense",
+    "obs_names_attr": "--obs-names-attr",
+    "var_names_attr": "--var-names-attr",
 }
 
 #: Flags that only make sense for format conversion; with recast's flags,
@@ -634,7 +637,9 @@ def convert(
     layer: Optional[str] = typer.Option(
         None,
         "--layer",
-        help="10x output: the matrix to write: X (default), raw/X, or a layer",
+        help="10x or loom output: the matrix to write: X (default), raw/X, or "
+        "a layer; loom input: the loom layer that becomes X (default spliced, "
+        "else the main matrix)",
     ),
     gene_ids_column: Optional[str] = typer.Option(
         None,
@@ -681,6 +686,22 @@ def convert(
         "--unfiltered",
         help="Parse input: read DGE_unfiltered (every barcode) not DGE_filtered",
     ),
+    dense: bool = typer.Option(
+        False,
+        "--dense",
+        help="loom input: keep X and the layers dense (default: sparse CSR)",
+    ),
+    obs_names_attr: Optional[str] = typer.Option(
+        None,
+        "--obs-names-attr",
+        help="loom input: the column attribute naming cells (default CellID)",
+    ),
+    var_names_attr: Optional[str] = typer.Option(
+        None,
+        "--var-names-attr",
+        help="loom input: the row attribute naming genes (default Accession, "
+        "or Gene with --var-names gene_symbols)",
+    ),
     # 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),
@@ -694,9 +715,9 @@ def convert(
     """
     Convert a store to another format.
 
-    The output's extension picks the format: .h5ad, .zarr, or .h5 for a 10x
-    Genomics matrix (Cell Ranger's *_feature_bc_matrix.h5); --to 10x-mtx
-    writes a Matrix Market folder. Between h5ad and zarr every element is
+    The output's extension picks the format: .h5ad, .zarr, .loom, or .h5 for
+    a 10x 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. CellBender output keeps its latent variables, as
     CellBender's own loader names them. A 10x matrix is read as scanpy's
@@ -728,6 +749,8 @@ def convert(
         adata convert cellbender_output.h5 -o cells.h5ad --min-cell-probability 0.5
         adata convert splitpipe_out -o parse.h5ad
         adata convert splitpipe_out -o parse.h5ad --all-samples --unfiltered
+        adata convert velocyto.loom -o velocyto.h5ad
+        adata convert counts.h5ad -o counts.loom
     """
     entries = list(ctx.args)
     recast_flags = _given(ctx, _RECAST_FLAGS)
@@ -816,6 +839,9 @@ def convert(
                 sample=sample,
                 all_samples=all_samples,
                 unfiltered=unfiltered,
+                dense=dense,
+                obs_names_attr=obs_names_attr,
+                var_names_attr=var_names_attr,
             ),
         )
     except Exception as e:

src/adata/conversion/__init__.py

@@ -17,7 +17,7 @@
 
 import h5py
 
-from adata.conversion import cellbender, mtx, parse, tenx
+from adata.conversion import cellbender, loom, mtx, parse, 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
@@ -83,6 +83,14 @@ def _sniff_h5ad(path: Path) -> bool:
         importer=tenx.TenxImport(),
         exporter=tenx.TenxExport(),
     ),
+    "loom": Format(
+        name="loom",
+        description="Loom (loompy, velocyto, kallisto|bustools; spec 2 and 3)",
+        extensions=(".loom",),
+        sniff=loom.sniff,
+        importer=loom.LoomImport(),
+        exporter=loom.LoomExport(),
+    ),
     # A folder, so no extension: chosen as output with `--to 10x-mtx`.
     "10x-mtx": Format(
         name="10x-mtx",

src/adata/conversion/base.py

@@ -44,6 +44,11 @@ class Options:
     all_samples: bool = False
     unfiltered: bool = False
 
+    # Reading a loom.
+    dense: bool = False
+    obs_names_attr: Optional[str] = None
+    var_names_attr: Optional[str] = None
+
 
 @dataclass
 class ConversionPlan: