diff --git a/AGENTS.md b/AGENTS.md index 97049d6..f272314 100644 --- a/AGENTS.md +++ b/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 diff --git a/README.md b/README.md index bc1ac80..6b47697 100644 --- a/README.md +++ b/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 diff --git a/docs/COMMANDS.md b/docs/COMMANDS.md index a30743c..66b4471 100644 --- a/docs/COMMANDS.md +++ b/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/`, and `/matrix` when X is `spliced` | `layers[]`, `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 +`é` 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 diff --git a/docs/TESTING.md b/docs/TESTING.md index 037c160..5130170 100644 --- a/docs/TESTING.md +++ b/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 | diff --git a/docs/index.md b/docs/index.md index d27972b..604151a 100644 --- a/docs/index.md +++ b/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 | |---|---|---| diff --git a/src/adata/cli.py b/src/adata/cli.py index 0115286..5775de8 100644 --- a/src/adata/cli.py +++ b/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: diff --git a/src/adata/conversion/__init__.py b/src/adata/conversion/__init__.py index 2f01d7c..4b4f287 100644 --- a/src/adata/conversion/__init__.py +++ b/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", diff --git a/src/adata/conversion/base.py b/src/adata/conversion/base.py index 66fc806..915b95f 100644 --- a/src/adata/conversion/base.py +++ b/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: diff --git a/src/adata/conversion/loom.py b/src/adata/conversion/loom.py new file mode 100644 index 0000000..2fcdbf7 --- /dev/null +++ b/src/adata/conversion/loom.py @@ -0,0 +1,956 @@ +"""Loom files: loompy's, velocyto's, kallisto|bustools'. + +A loom is one HDF5 file: + +* `/matrix`, genes x cells, dense, usually chunked 64 x 64 and gzipped, and + `/layers/` of the same shape; +* `/col_attrs/`, one value (or row) per cell, and `/row_attrs/` + per gene; +* `/col_graphs/` and `/row_graphs/`, COO graphs (`a`, `b`, `w`); +* global attributes: datasets under `/attrs` in spec 3.0.0, root HDF5 + attributes before it. + +Import follows `anndata.read_loom`: X is the `spliced` layer when there is +one (with the main matrix kept as `layers["matrix"]`), else the main matrix; +attributes with one value per cell are obs columns and those with several +are obsm. Three things differ, on purpose: + +* values keep their dtype, where anndata casts X to float32; +* graphs become obsp/varp, and global attributes `uns["loom"]`, where + anndata drops them; +* var is named by gene id when the file records one, as `convert` names + every other format. + +Strings are read as loompy reads them -- XML character references such as +`é` unescaped -- and written as loompy writes them, ASCII with those +references, so loompy reads our files back unchanged. + +The matrix is read and written a slab of cells at a time, aligned to the +file's chunks, so each chunk is decompressed or compressed once and memory +is set by the slab, not the file. +""" + +from __future__ import annotations + +import html +from contextlib import ExitStack +from dataclasses import dataclass, field +from datetime import datetime, timezone +from pathlib import Path +from typing import Any, Dict, Iterator, List, Optional, Sequence, Tuple + +import h5py +import numpy as np +from rich.console import Console + +from adata.conversion.base import ConversionPlan, Converter, Options +from adata.conversion.features import ( + STRING_STEP, + Target, + check_import_options, + resolve, + write_obs_index, + write_var, +) +from adata.core.csr_build import CsrAppender, count_keys, scratch_group, write_runs +from adata.core.recast import Matrix, describe, transpose_sparse_streaming +from adata.elements import spec +from adata.elements.read import ( + dataframe_columns, + element_len, + read_str_chunk, + resolve_index, +) +from adata.elements.write import ( + ensure_anndata_skeleton, + set_column_order, + write_dense, + write_mapping, + write_scalar, + write_string_array, +) +from adata.formats.table import write_text_column +from adata.storage import ( + COMPRESSIONS, + create_dataset, + is_dataset, + is_group, + is_zarr_group, + open_store, + with_compression, +) +from adata.util.names import count_duplicates, make_index_unique + +#: Attributes that name the cells, the gene ids and the gene symbols, in the +#: order they are looked for: loompy and velocyto's names first, then +#: kallisto|bustools'. +OBS_NAMES = ("CellID", "cell_id", "obs_names", "barcode", "barcodes") +GENE_IDS = ("Accession", "gene_id", "gene_ids") +GENE_SYMBOLS = ("Gene", "gene_name", "gene_names", "gene_symbols") + +#: Bytes of matrix held at once: one slab of cells across every gene. +SLAB_BYTES = 64 << 20 + +#: Above this many bytes, dense, an export says how large the matrix is. +LARGE_DENSE = 10 << 30 + +_MATRIX_TYPES = frozenset( + "float16 float32 float64 int8 int16 int32 int64 uint8 uint16 uint32 uint64".split() +) + +#: How loompy writes the matrix and layers. +_CHUNK = 64 +_GZIP_LEVEL = 2 + + +# --------------------------------------------------------------------------- +# reading + + +def sniff(path: Path) -> bool: + path = Path(path) + if path.is_dir() or not h5py.is_hdf5(path): + return False + with h5py.File(path, "r") as handle: + matrix = handle.get("matrix") + return ( + isinstance(matrix, h5py.Dataset) + and matrix.ndim == 2 + and isinstance(handle.get("row_attrs"), h5py.Group) + and isinstance(handle.get("col_attrs"), h5py.Group) + ) + + +def decode(values: Any) -> Any: + """Text as loompy reads it: bytes decoded, XML references unescaped.""" + if isinstance(values, (str, bytes, np.str_, np.bytes_)): + text = values.decode("utf-8", "replace") if isinstance(values, bytes) else str(values) + return html.unescape(text) + arr = np.asarray(values) + if arr.dtype.kind not in "SOU": + return values + flat = [ + html.unescape(v.decode("utf-8", "replace") if isinstance(v, bytes) else str(v)) + for v in arr.reshape(-1) + ] + out = np.array(flat, dtype=object).reshape(arr.shape) + return out if arr.ndim else out[()] + + +def _version(handle: h5py.File) -> str: + if isinstance(handle.get("attrs"), h5py.Group) and "LOOM_SPEC_VERSION" in handle["attrs"]: + return str(decode(handle["attrs"]["LOOM_SPEC_VERSION"][()])) + if "LOOM_SPEC_VERSION" in handle.attrs: + return str(decode(handle.attrs["LOOM_SPEC_VERSION"])) + return "pre-2" + + +@dataclass +class _Axis: + """One axis's attributes, sorted by where they go.""" + + n: int + names: Optional[str] = None + columns: List[str] = field(default_factory=list) + matrices: List[str] = field(default_factory=list) + skipped: List[str] = field(default_factory=list) + + +def _axis(group: h5py.Group, n: int, what: str) -> _Axis: + axis = _Axis(n) + for name, ds in group.items(): + if not isinstance(ds, h5py.Dataset): + axis.skipped.append(name) + continue + if ds.dtype.names is not None or ds.dtype.kind == "V" or ds.ndim == 0: + axis.skipped.append(name) + continue + if ds.shape[0] != n: + raise ValueError( + f"{what} attribute {name!r} has {ds.shape[0]:,} entries, but the " + f"matrix has {n:,} {'cells' if what == 'Column' else 'genes'}." + ) + if ds.ndim == 1 or (ds.ndim == 2 and ds.shape[1] == 1): + axis.columns.append(name) + elif ds.ndim == 2 and ds.dtype.kind in "biuf": + axis.matrices.append(name) + else: + axis.skipped.append(name) + return axis + + +def _pick(axis: _Axis, explicit: Optional[str], candidates: Sequence[str], flag: str) -> Optional[str]: + if explicit is not None: + if explicit not in axis.columns: + raise ValueError( + f"{flag} {explicit!r}: no such attribute; there are " + f"{', '.join(axis.columns) or 'none'}." + ) + return explicit + return next((name for name in candidates if name in axis.columns), None) + + +@dataclass +class _Import: + version: str + velocyto: bool + n_cells: int + n_genes: int + #: Loom matrix ("" for the main one) -> the AnnData places it fills. + sources: Dict[str, List[str]] + obs: _Axis + var: _Axis + var_names: List[str] + graphs: Dict[str, List[str]] + dense: bool + compression: Optional[str] + zarr_format: Optional[int] + + +class LoomImport(Converter): + options = frozenset({ + "--layer", "--var-names", "--make-unique", "--compression", "--dense", + "--obs-names-attr", "--var-names-attr", + }) + + def plan( + self, + source: Path, + output: Path, + src_format: str, + dst_format: str, + options: Options, + ) -> ConversionPlan: + check_import_options(options) + backend = "zarr" if dst_format == "zarr" else "hdf5" + if options.compression is not None and options.compression not in COMPRESSIONS[backend]: + raise ValueError( + f"--compression {options.compression} is not available for " + f"{dst_format}; use one of: {', '.join(COMPRESSIONS[backend])}." + ) + warnings: List[str] = [] + with h5py.File(source, "r") as handle: + matrix = handle["matrix"] + n_genes, n_cells = matrix.shape + layers = handle["layers"] if isinstance(handle.get("layers"), h5py.Group) else {} + for name, layer in layers.items(): + if not isinstance(layer, h5py.Dataset) or layer.shape != matrix.shape: + raise ValueError( + f"Layer {name!r} is not a {n_genes:,} x {n_cells:,} matrix " + "like the main one." + ) + for name, ds in (("matrix", matrix), *layers.items()): + if ds.dtype.kind not in "biuf": + raise ValueError(f"{name} holds {ds.dtype} values, which a loom may not.") + + x = options.layer + if x is None: + x = "spliced" if "spliced" in layers else "" + elif x == "matrix": + x = "" + elif x not in layers: + raise ValueError( + f"--layer {x!r}: no such layer; there are matrix (the main one)" + + "".join(f", {name}" for name in layers) + "." + ) + sources: Dict[str, List[str]] = {"": [], **{name: [] for name in layers}} + sources[x].append("X") + for name in layers: + sources[name].append(f"layers/{name}") + if x != "" and "matrix" not in layers: + # anndata's place for the main matrix, unless a layer of + # that name -- from an exported AnnData -- takes it. + sources[""].append("layers/matrix") + + obs = _axis(handle["col_attrs"], n_cells, "Column") + var = _axis(handle["row_attrs"], n_genes, "Row") + obs.names = _pick(obs, options.obs_names_attr, OBS_NAMES, "--obs-names-attr") + if obs.names is None: + warnings.append( + "No attribute names the cells (looked for " + f"{', '.join(OBS_NAMES)}); they are numbered 0..{n_cells - 1}, " + "as anndata numbers them. Name one with --obs-names-attr." + ) + first, second = ( + (GENE_IDS, GENE_SYMBOLS) if options.var_names == "gene_ids" + else (GENE_SYMBOLS, GENE_IDS) + ) + var.names = _pick(var, options.var_names_attr, first, "--var-names-attr") + if var.names is None and options.var_names_attr is None: + var.names = _pick(var, None, second, "--var-names-attr") + if var.names is not None: + warnings.append( + f"No gene {'ids' if first is GENE_IDS else 'symbols'} recorded; " + f"genes are named by {var.names!r}." + ) + if var.names is None: + warnings.append( + f"No attribute names the genes; they are numbered 0..{n_genes - 1}. " + "Name one with --var-names-attr." + ) + for axis, what in ((obs, "col_attrs"), (var, "row_attrs")): + if axis.names in axis.columns: + axis.columns.remove(axis.names) + if axis.skipped: + warnings.append( + f"Not read: {what} {', '.join(axis.skipped)} -- compound or " + "higher-dimensional values, such as SCope's regulons (#27)." + ) + + var_names = ( + [str(v) for v in decode(handle["row_attrs"][var.names][...]).reshape(-1)] + if var.names is not None else [str(i) for i in range(n_genes)] + ) + if options.make_unique: + var_names = make_index_unique(var_names) + elif count_duplicates(var_names): + warnings.append( + f"{count_duplicates(var_names):,} var names repeat an earlier one. " + "Pass --make-unique to suffix them as anndata's " + "var_names_make_unique does." + ) + + graphs = { + kind: [ + name for name, graph in handle[kind].items() + if isinstance(graph, h5py.Group) and {"a", "b", "w"} <= set(graph) + ] + for kind in ("col_graphs", "row_graphs") + if isinstance(handle.get(kind), h5py.Group) + } + version = _version(handle) + velocyto = any( + key.startswith("velocyto") for key in _global_attrs(handle) + ) or {"spliced", "unspliced"} <= set(layers) + + detail = _Import( + version=version, + velocyto=velocyto, + n_cells=n_cells, + n_genes=n_genes, + sources=sources, + obs=obs, + var=var, + var_names=var_names, + graphs=graphs, + dense=options.dense, + compression=options.compression or ("gzip" if backend == "hdf5" else None), + zarr_format=(options.zarr_format or 3) if dst_format == "zarr" else None, + ) + kind = "velocyto" if velocyto else f"spec {version}" + named = f"X = {x or 'matrix'}" + others = [p.split("/", 1)[1] for places in sources.values() for p in places if p != "X"] + if others: + named += f", layers: {', '.join(others)}" + return ConversionPlan( + source=source, + src_format=src_format, + dst_format=dst_format, + outputs=[output], + summary=( + f"loom ({kind}) → {dst_format}, {n_cells:,} cells x {n_genes:,} " + f"genes, {named}{', dense' if options.dense else ''}" + ), + warnings=warnings, + detail=detail, + ) + + def write( + self, plan: ConversionPlan, staged: Sequence[Path], console: Console + ) -> None: + detail: _Import = plan.detail + with h5py.File(plan.source, "r") as handle, ExitStack() as stack: + root = stack.enter_context( + open_store(staged[0], "w", zarr_format=detail.zarr_format) + ).root + with console.status("[cyan]Writing obs, var and the attributes...[/]"): + _write_obs(root, handle["col_attrs"], detail.obs) + _write_var(root, handle["row_attrs"], detail.var, detail.var_names) + ensure_anndata_skeleton(root) + for axis, group, key in ( + (detail.obs, handle["col_attrs"], "obsm"), + (detail.var, handle["row_attrs"], "varm"), + ): + for name in axis.matrices: + _copy_rows(group[name], root[key], name) + for kind, key in (("col_graphs", "obsp"), ("row_graphs", "varp")): + n = detail.n_cells if kind == "col_graphs" else detail.n_genes + for name in detail.graphs.get(kind, []): + _write_graph(handle[kind][name], root, key, name, n) + _write_uns(root, handle, detail) + + for loom_name, places in detail.sources.items(): + if not places: + continue + ds = handle["matrix"] if loom_name == "" else handle["layers"][loom_name] + with console.status(f"[cyan]Reading {loom_name or 'matrix'}...[/]"): + _write_matrix(ds, root, places, detail) + + +def _global_attrs(handle: h5py.File) -> Dict[str, Any]: + """Global attributes, spec 3's datasets and older root attributes alike.""" + found: Dict[str, Any] = {key: value for key, value in handle.attrs.items()} + attrs = handle.get("attrs") + if isinstance(attrs, h5py.Group): + for key, ds in attrs.items(): + if isinstance(ds, h5py.Dataset): + found[key] = ds[()] + return found + + +def _write_obs(root: Any, group: h5py.Group, axis: _Axis) -> None: + n = axis.n + if axis.names is not None: + names = group[axis.names] + blocks = ( + [str(v) for v in decode(names[lo : lo + STRING_STEP]).reshape(-1)] + for lo in range(0, n, STRING_STEP) + ) + else: + blocks = ( + [str(i) for i in range(lo, min(lo + STRING_STEP, n))] + for lo in range(0, n, STRING_STEP) + ) + obs = write_obs_index(root, n, blocks) + for name in axis.columns: + _write_column(obs, name, group[name], n) + set_column_order(obs, axis.columns) + + +def _write_column(parent: Any, name: str, ds: h5py.Dataset, n: int) -> None: + """One attribute as a dataframe column, a block of rows at a time.""" + if ds.dtype.kind in "SOU": + def blocks() -> Iterator[List[str]]: + for lo in range(0, n, STRING_STEP): + yield [str(v) for v in decode(ds[lo : lo + STRING_STEP]).reshape(-1)] + + write_text_column(parent, name, n, blocks) + return + target = create_dataset(parent, name, shape=(n,), dtype=ds.dtype) + spec.set_encoding(target, spec.ARRAY) + for lo in range(0, n, STRING_STEP): + target[lo : lo + STRING_STEP] = np.asarray(ds[lo : lo + STRING_STEP]).reshape(-1) + + +def _write_var(root: Any, group: h5py.Group, axis: _Axis, var_names: List[str]) -> None: + columns = {} + for name in axis.columns: + values = group[name][...] + values = values.reshape(-1) if values.ndim == 2 else values + columns[name] = decode(values) if values.dtype.kind in "SOU" else values + write_var(root, Target(label=None, keep=None, var_names=var_names, columns=columns)) + + +def _copy_rows(ds: h5py.Dataset, parent: Any, name: str) -> None: + """A 2-D attribute as obsm/varm, a block of rows at a time.""" + n = ds.shape[0] + target = create_dataset(parent, name, shape=ds.shape, dtype=ds.dtype) + spec.set_encoding(target, spec.ARRAY) + for lo in range(0, n, STRING_STEP): + target[lo : lo + STRING_STEP] = ds[lo : lo + STRING_STEP] + + +def _write_graph(graph: h5py.Group, root: Any, key: str, name: str, n: int) -> None: + """A COO graph as an n x n CSR matrix, its edges sorted through buckets.""" + a, b, w = graph["a"], graph["b"], graph["w"] + m = a.shape[0] + if b.shape[0] != m or w.shape[0] != m: + raise ValueError(f"Graph {name!r} has a, b and w of different lengths.") + step = 1 << 20 + + def blocks() -> Iterator[Tuple[np.ndarray, np.ndarray, np.ndarray]]: + for lo in range(0, m, step): + rows = np.asarray(a[lo : lo + step], dtype=np.int64) + cols = np.asarray(b[lo : lo + step], dtype=np.int64) + if rows.size and (rows.min() < 0 or rows.max() >= n or cols.min() < 0 or cols.max() >= n): + raise ValueError(f"Graph {name!r} has an edge outside 0..{n - 1}.") + yield rows, cols, np.asarray(w[lo : lo + step]) + + counts = np.zeros(n, dtype=np.int64) + grouped, last = True, -1 + for rows, _, _ in blocks(): + count_keys(counts, rows) + if rows.size: + grouped = grouped and rows[0] >= last and bool((np.diff(rows) >= 0).all()) + last = int(rows[-1]) + appender = CsrAppender(root[key], name, shape=(n, n), data_dtype=w.dtype) + write_runs( + blocks(), n_rows=n, n_targets=1, destinations=[[appender]], scratch_parent=root, + grouped=grouped, canonical=False, counts=counts, value_dtype=w.dtype, chunk=step, + ) + appender.finish() + + +def _write_uns(root: Any, handle: h5py.File, detail: _Import) -> None: + meta = write_mapping(root["uns"], "loom") + write_scalar(meta, "spec", detail.version) + for key, name in (("obs_names", detail.obs.names), ("var_names", detail.var.names)): + if name is not None: + write_scalar(meta, key, name) + attrs = write_mapping(meta, "attrs") + for key, value in _global_attrs(handle).items(): + value = decode(value) + arr = np.asarray(value) + if arr.ndim == 0: + item = arr.item() if arr.dtype.kind != "O" else value + if isinstance(item, (str, int, float, np.integer, np.floating)): + write_scalar(attrs, key, item) + elif arr.dtype.kind == "O": + write_string_array(attrs, key, [str(v) for v in arr.reshape(-1)]) + elif arr.dtype.kind in "biuf": + write_dense(attrs, key, arr) + + +def slab_cells(n_genes: int, itemsize: int, chunk_width: int) -> int: + """Cells per slab: as many whole chunk columns as `SLAB_BYTES` holds, one at least.""" + per_cell = max(1, n_genes * itemsize) + fit = max(1, SLAB_BYTES // per_cell) + return max(chunk_width, fit // chunk_width * chunk_width) + + +def _write_matrix(ds: h5py.Dataset, root: Any, places: List[str], detail: _Import) -> None: + """One loom matrix, read once, written to every place that holds it.""" + n_genes, n_cells = ds.shape + width = ds.chunks[1] if ds.chunks else _CHUNK + slab = slab_cells(n_genes, ds.dtype.itemsize, width) + parents = [ + (root, "X") if place == "X" else (root["layers"], place.split("/", 1)[1]) + for place in places + ] + backend = "zarr" if is_zarr_group(root) else "hdf5" + if detail.dense: + targets = [] + for parent, name in parents: + kw: Dict[str, Any] = {"chunks": (min(max(n_cells, 1), slab), max(n_genes, 1))} + if detail.compression is not None: + kw = with_compression( + kw, detail.compression, shape=(n_cells, n_genes), dtype=ds.dtype, + backend=backend, zarr_format=detail.zarr_format, + ) + target = create_dataset(parent, name, shape=(n_cells, n_genes), dtype=ds.dtype, **kw) + spec.set_encoding(target, spec.ARRAY) + targets.append(target) + for lo in range(0, n_cells, slab): + block = np.asarray(ds[:, lo : lo + slab]).T + for target in targets: + target[lo : lo + block.shape[0]] = block + return + + appenders = [ + CsrAppender( + parent, name, shape=(n_cells, n_genes), data_dtype=ds.dtype, + compression=detail.compression, + ) + for parent, name in parents + ] + for lo in range(0, n_cells, slab): + block = np.asarray(ds[:, lo : lo + slab]).T + rows, cols = np.nonzero(block) + values = block[rows, cols] + for appender in appenders: + appender.append(rows + lo, cols, values) + for appender in appenders: + appender.finish() + + +# --------------------------------------------------------------------------- +# export + + +def _ascii(values: Sequence[str]) -> np.ndarray: + """Text as loompy writes it: ASCII, anything else as XML references.""" + return np.array([str(v).encode("ascii", "xmlcharrefreplace") for v in values], dtype=object) + + +_ASCII = h5py.string_dtype(encoding="ascii") + + +@dataclass +class _Export: + matrix: str + layers: List[str] + obs_columns: List[str] + var_columns: List[str] + obsm: List[str] + varm: List[str] + obsp: List[str] + varp: List[str] + #: Where Gene and Accession come from: a var column, "" for var_names, + #: None (Accession only) for nowhere. + gene: str + accession: Optional[str] + uns: List[str] + #: var, or raw/var when the matrix is raw/X. + var: str = "var" + varm_key: str = "varm" + + +class LoomExport(Converter): + options = frozenset({"--layer"}) + + def plan( + self, + source: Path, + output: Path, + src_format: str, + dst_format: str, + options: Options, + ) -> ConversionPlan: + warnings: List[str] = [] + path = options.layer or "X" + if path not in ("X", "raw/X") and not path.startswith("layers/"): + path = f"layers/{path}" + with open_store(source, "r", require_anndata=False) as store: + root = store.root + raw = path.startswith("raw/") + var_path = "raw/var" if raw else "var" + var_group = resolve(root, var_path) + obs_index, obs_name = resolve_index(root["obs"], "obs") + var_index, var_name = resolve_index(var_group, "var") + n_obs, n_var = element_len(obs_index), element_len(var_index) + main = _matrix(root, path, (n_obs, n_var)) + layers = [] + if raw: + # raw/X has raw/var's genes; X and the layers have var's. + if "layers" in root and len(root["layers"]) or path != "X": + warnings.append( + "Not written: X and the layers -- they follow var, and " + "raw/X follows raw/var." + ) + elif "layers" in root: + for name in root["layers"]: + if f"layers/{name}" == path: + continue + try: + _matrix(root, f"layers/{name}", (n_obs, n_var)) + layers.append(name) + except ValueError as exc: + warnings.append(f"Not written: layers/{name}: {exc}") + if path != "X" and not raw: + layers.insert(0, "X") + obs_columns = _columns(root["obs"], obs_name, warnings, "obs") + var_columns = _columns(var_group, var_name, warnings, var_path) + gene, accession = _gene_attrs(var_columns) + varm_key = "raw/varm" if raw else "varm" + obsm = _arrays(root, "obsm", n_obs, warnings) + varm = _arrays(root, varm_key, n_var, warnings) + obsp = _graphs(root, "obsp", n_obs, warnings) + varp = [] if raw else _graphs(root, "varp", n_var, warnings) + uns = _uns_scalars(root, warnings) + + dense = n_obs * n_var * max(1, main.dtype.itemsize) * (1 + len(layers)) + summary = ( + f"{src_format} → loom 3.0.0, {n_var:,} genes x {n_obs:,} cells from " + f"{path}{', layers: ' + ', '.join(layers) if layers else ''}; " + f"dense {_size(dense)} before compression" + ) + if dense > LARGE_DENSE: + warnings.append( + f"A loom holds its matrices dense: {_size(dense)} to compress and " + "write. It is written, but may take a while; h5ad or zarr keep " + "them sparse." + ) + return ConversionPlan( + source=source, + src_format=src_format, + dst_format=dst_format, + outputs=[output], + summary=summary, + warnings=warnings, + detail=_Export( + matrix=path, layers=layers, obs_columns=obs_columns, + var_columns=var_columns, obsm=obsm, varm=varm, obsp=obsp, varp=varp, + gene=gene, accession=accession, uns=uns, var=var_path, + varm_key=varm_key, + ), + ) + + def write( + self, plan: ConversionPlan, staged: Sequence[Path], console: Console + ) -> None: + detail: _Export = plan.detail + with open_store(plan.source, "r", require_anndata=False) as store, h5py.File( + staged[0], "w" + ) as out: + root = store.root + with console.status(f"[cyan]Writing {detail.matrix} dense...[/]"): + _export_matrix(resolve(root, detail.matrix), out, "matrix") + layers = out.create_group("layers") + for name in detail.layers: + with console.status(f"[cyan]Writing layer {name} dense...[/]"): + obj = root["X"] if name == "X" else root["layers"][name] + _export_matrix(obj, layers, name) + + with console.status("[cyan]Writing the attributes...[/]"): + cols = out.create_group("col_attrs") + rows = out.create_group("row_attrs") + obs, var = root["obs"], resolve(root, detail.var) + obs_index, _ = resolve_index(obs, "obs") + var_index, _ = resolve_index(var, "var") + _export_strings(cols, "CellID", obs_index, obs) + for name in detail.obs_columns: + if name != "CellID": + _export_column(cols, name, obs[name], obs) + for attr, where in (("Gene", detail.gene), ("Accession", detail.accession)): + if where == "": + _export_strings(rows, attr, var_index, var) + elif where is not None: + _export_column(rows, attr, var[where], var) + skip = {detail.gene, detail.accession, "Gene", "Accession"} + for name in detail.var_columns: + if name not in skip: + _export_column(rows, name, var[name], var) + for key, group, where in ( + ("obsm", cols, "obsm"), ("varm", rows, detail.varm_key) + ): + for name in getattr(detail, key): + _export_array(resolve(root, where)[name], group, name) + for key, kind in (("obsp", "col_graphs"), ("varp", "row_graphs")): + graphs = out.create_group(kind) + for name in getattr(detail, key): + _export_graph(describe(root[key][name]), graphs.create_group(name)) + attrs = out.create_group("attrs") + _attr(attrs, "LOOM_SPEC_VERSION", "3.0.0") + _attr(attrs, "CreationDate", datetime.now(timezone.utc).strftime("%Y%m%dT%H%M%S.%fZ")) + for name in detail.uns: + value = _uns_value(root, name) + if name not in ("LOOM_SPEC_VERSION", "CreationDate"): + _attr(attrs, name, value) + + +def _size(n: int) -> str: + for unit, scale in (("TB", 1 << 40), ("GB", 1 << 30), ("MB", 1 << 20)): + if n >= scale: + return f"{n / scale:,.1f} {unit}" + return f"{n:,} bytes" + + +def _matrix(root: Any, path: str, shape: Tuple[int, int]) -> Matrix: + obj = resolve(root, path) + matrix = describe(obj) + if matrix.shape != shape: + raise ValueError( + f"{path} is {matrix.shape[0]:,} x {matrix.shape[1]:,}, not " + f"{shape[0]:,} x {shape[1]:,}." + ) + if matrix.dtype.kind not in "biuf": + raise ValueError(f"{path} holds {matrix.dtype} values, which a loom may not.") + return matrix + + +def _column_kind(col: Any) -> Optional[str]: + """'number', 'text', or None for what a loom attribute cannot hold.""" + if is_dataset(col): + if col.ndim != 1: + return None + kind = col.dtype.kind + if kind in "biuf": + return "number" + if kind in "SOUT": + return "text" + return None + if is_group(col): + enc = spec.encoding_type(col) + if enc in (spec.CATEGORICAL, spec.NULLABLE_STRING_ARRAY): + return "text" + if enc in (spec.NULLABLE_INTEGER, spec.NULLABLE_BOOLEAN): + return "masked" + return None + + +def _columns(group: Any, index_name: str, warnings: List[str], axis: str) -> List[str]: + kept, dropped = [], [] + for name in dataframe_columns(group, index_name): + (kept if _column_kind(group[name]) else dropped).append(name) + if dropped: + warnings.append( + f"Not written: {axis} {', '.join(dropped)} -- a loom attribute holds " + "numbers or text only." + ) + return kept + + +def _gene_attrs(columns: List[str]) -> Tuple[str, Optional[str]]: + """Where loom's Gene (symbol) and Accession (id) come from. + + Each is chosen on its own: an explicit column first (`gene_symbols`, + `Gene`, ... and `gene_ids`, `Accession`), and var_names ("") only for the + one no column holds -- ids when there is a symbols column (as `convert` + imports by default), symbols otherwise. Accession is None when nothing + can be taken for ids. + """ + ids = next((n for n in ("gene_ids", "Accession") if n in columns), None) + symbols = next( + (n for n in ("gene_symbols", "Gene", "gene_name", "gene_names") if n in columns), None + ) + if symbols is None: + return "", ids + return symbols, "" if ids is None else ids + + +def _arrays(root: Any, key: str, n: int, warnings: List[str]) -> List[str]: + kept, dropped = [], [] + for name in (root[key] if key in root else []): + obj = root[key][name] + if is_dataset(obj) and obj.ndim == 2 and obj.dtype.kind in "biuf" and obj.shape[0] == n: + kept.append(name) + else: + dropped.append(name) + if dropped: + warnings.append(f"Not written: {key} {', '.join(dropped)} -- only numeric 2-D arrays fit a loom.") + return kept + + +def _graphs(root: Any, key: str, n: int, warnings: List[str]) -> List[str]: + kept, dropped = [], [] + for name in (root[key] if key in root else []): + try: + matrix = describe(root[key][name]) + if not matrix.sparse or matrix.shape != (n, n): + raise ValueError + kept.append(name) + except (ValueError, KeyError): + dropped.append(name) + if dropped: + warnings.append(f"Not written: {key} {', '.join(dropped)} -- only sparse graphs fit a loom.") + return kept + + +def _uns_scalars(root: Any, warnings: List[str]) -> List[str]: + """uns entries a global attribute can hold: scalars and 1-D arrays, top + level or, from a loom import, under uns["loom"]["attrs"].""" + found = [] + uns = root["uns"] if "uns" in root else {} + for name in uns: + if name == "loom": + attrs = uns["loom"].get("attrs") if is_group(uns["loom"]) else None + found += [f"loom/attrs/{k}" for k in (attrs or []) if is_dataset(attrs[k]) and attrs[k].ndim <= 1] + elif is_dataset(uns[name]) and uns[name].ndim <= 1: + found.append(name) + return found + + +def _uns_value(root: Any, name: str) -> Any: + value = root["uns"][name][()] + return decode(value) if np.asarray(value).dtype.kind in "SOUT" else value + + +def _attr(group: h5py.Group, name: str, value: Any) -> None: + key = name.rsplit("/", 1)[-1] + if key in group: + return + arr = np.asarray(value) + if arr.dtype.kind in "OSUT": + if arr.ndim == 0: + group.create_dataset(key, data=str(value).encode("ascii", "xmlcharrefreplace"), dtype=_ASCII) + else: + group.create_dataset(key, data=_ascii(arr.reshape(-1)), dtype=_ASCII) + elif arr.dtype.kind == "b": + group.create_dataset(key, data=arr.astype(np.uint8)) + else: + group.create_dataset(key, data=arr) + + +def _export_strings(group: h5py.Group, name: str, column: Any, parent: Any) -> None: + n = element_len(column) + target = group.create_dataset(name, shape=(n,), dtype=_ASCII) + cache: Dict[str, np.ndarray] = {} + for lo in range(0, n, STRING_STEP): + hi = min(lo + STRING_STEP, n) + target[lo:hi] = _ascii(read_str_chunk(column, lo, hi, cache, parent_group=parent)) + + +def _export_column(group: h5py.Group, name: str, column: Any, parent: Any) -> None: + kind = _column_kind(column) + if kind == "text": + _export_strings(group, name, column, parent) + return + n = element_len(column) + if kind == "masked": + # A loom number cannot be missing but as NaN, so these become floats. + target = group.create_dataset(name, shape=(n,), dtype=np.float64) + for lo in range(0, n, STRING_STEP): + values = np.asarray(column["values"][lo : lo + STRING_STEP], dtype=np.float64) + values[np.asarray(column["mask"][lo : lo + STRING_STEP], dtype=bool)] = np.nan + target[lo : lo + values.size] = values + return + dtype = np.dtype(np.uint8) if column.dtype.kind == "b" else column.dtype + target = group.create_dataset(name, shape=(n,), dtype=dtype) + for lo in range(0, n, STRING_STEP): + target[lo : lo + STRING_STEP] = np.asarray(column[lo : lo + STRING_STEP]).astype(dtype, copy=False) + + +def _export_array(source: Any, group: h5py.Group, name: str) -> None: + target = group.create_dataset(name, shape=source.shape, dtype=source.dtype) + for lo in range(0, source.shape[0], STRING_STEP): + target[lo : lo + STRING_STEP] = source[lo : lo + STRING_STEP] + + +def _export_graph(matrix: Matrix, group: h5py.Group) -> None: + """A sparse graph as loom's COO vectors: a (from), b (to), w (weight).""" + weight = matrix.dtype if matrix.dtype.kind == "f" else np.dtype(np.float32) + nnz = matrix.nnz + a = group.create_dataset("a", shape=(nnz,), dtype=np.int64) + b = group.create_dataset("b", shape=(nnz,), dtype=np.int64) + w = group.create_dataset("w", shape=(nnz,), dtype=weight) + indptr = np.asarray(matrix.obj["indptr"][...], dtype=np.int64) + step = 1 << 20 + for lo in range(0, nnz, step): + hi = min(lo + step, nnz) + minor = np.asarray(matrix.obj["indices"][lo:hi], dtype=np.int64) + # The major index of each entry, for this block only. + major = np.searchsorted(indptr, np.arange(lo, hi, dtype=np.int64), side="right") - 1 + rows, cols = (major, minor) if matrix.kind == spec.CSR_MATRIX else (minor, major) + a[lo:hi], b[lo:hi] = rows, cols + w[lo:hi] = np.asarray(matrix.obj["data"][lo:hi]).astype(weight, copy=False) + + +def _loom_dtype(dtype: np.dtype) -> np.dtype: + return np.dtype(np.uint8) if dtype.kind == "b" else dtype + + +def _export_matrix(obj: Any, parent: h5py.Group, name: str) -> None: + """An AnnData matrix (cells x genes, any layout) as a dense genes x cells + dataset, a slab of cells at a time, each slab a whole number of chunks.""" + matrix = describe(obj) + n_cells, n_genes = matrix.shape + dtype = _loom_dtype(matrix.dtype) + kw: Dict[str, Any] = {} + if n_cells and n_genes: + kw = { + "chunks": (min(_CHUNK, n_genes), min(_CHUNK, n_cells)), + "compression": "gzip", "compression_opts": _GZIP_LEVEL, + } + target = parent.create_dataset(name, shape=(n_genes, n_cells), dtype=dtype, fillvalue=0, **kw) + slab = slab_cells(n_genes, dtype.itemsize, _CHUNK) + + if not matrix.sparse: + for lo in range(0, n_cells, slab): + target[:, lo : lo + slab] = np.asarray(matrix.obj[lo : lo + slab]).astype(dtype, copy=False).T + return + if matrix.kind == spec.CSC_MATRIX: + # Cells are the minor axis: sort them to the major one on disk first, + # rather than hold every cell of a slab of genes. + with scratch_group(parent, "__loom_scratch__") as scratch: + transpose_sparse_streaming( + matrix, scratch, "csr", data_dtype=matrix.dtype, index_dtype=np.dtype(np.int64) + ) + _dense_from_csr(describe(scratch["csr"]), target, slab, dtype) + return + _dense_from_csr(matrix, target, slab, dtype) + + +def _dense_from_csr(matrix: Matrix, target: h5py.Dataset, slab: int, dtype: np.dtype) -> None: + n_cells, n_genes = matrix.shape + indptr = np.asarray(matrix.obj["indptr"][...], dtype=np.int64) + indices, data = matrix.obj["indices"], matrix.obj["data"] + for lo in range(0, n_cells, slab): + hi = min(lo + slab, n_cells) + block = np.zeros((n_genes, hi - lo), dtype=dtype) + start, end = int(indptr[lo]), int(indptr[hi]) + if end > start: + genes = np.asarray(indices[start:end], dtype=np.int64) + cells = np.repeat(np.arange(hi - lo, dtype=np.int64), np.diff(indptr[lo : hi + 1])) + # Repeated coordinates mean their sum, as scipy's toarray gives. + np.add.at(block, (genes, cells), np.asarray(data[start:end]).astype(dtype, copy=False)) + target[:, lo:hi] = block diff --git a/src/adata/formats/table.py b/src/adata/formats/table.py index 1b4c3bb..2a56d7d 100644 --- a/src/adata/formats/table.py +++ b/src/adata/formats/table.py @@ -27,7 +27,7 @@ import io from dataclasses import dataclass, field from pathlib import Path -from typing import Any, Dict, Iterator, List, Optional, Set, Tuple +from typing import Any, Callable, Dict, Iterable, Iterator, List, Optional, Sequence, Set, Tuple import numpy as np @@ -362,3 +362,45 @@ def read_small_table(path: Path, delimiter: str = ",") -> Tuple[List[str], List[ if header is None: raise ValueError(f"{Path(path).name} is empty: it has no header line.") return header, body + + +def write_text_column( + parent: Any, + name: str, + n: int, + blocks: Callable[[], Iterable[Sequence[str]]], + *, + missing: Sequence[str] = (), +) -> Any: + """A text column of `n` values, from `blocks()` called twice, as anndata + would write it: a categorical when it repeats a value, else strings. + + Pass 1 collects the distinct values, up to `CATEGORY_CAP`; pass 2 writes + codes in anndata's (natural) category order, or the strings themselves. + Values in `missing` are gaps: no category, code -1. + """ + gaps = frozenset(missing) + seen: Optional[Set[str]] = set() + for block in blocks(): + seen.update(v for v in block if v not in gaps) + if len(seen) > CATEGORY_CAP: + seen = None + break + backend = "zarr" if is_zarr_group(parent) else "hdf5" + if seen is not None and len(seen) < n: + categories = natural_sorted(seen) + lookup = {c: i for i, c in enumerate(categories)} + target = create_categorical(parent, name, categories, n) + convert = lambda block: np.array([lookup.get(v, -1) for v in block], dtype=np.int64) # noqa: E731 + else: + target = create_dataset(parent, name, shape=(n,), dtype=string_dtype_for(backend)) + spec.set_encoding(target, spec.STRING_ARRAY) + convert = lambda block: np.asarray(list(block), dtype=object) # noqa: E731 + written = 0 + for block in blocks(): + values = convert(block) + target[written : written + values.size] = values + written += values.size + if written != n: + raise ValueError(f"Column {name!r} has {written:,} values, not {n:,}.") + return parent[name] diff --git a/tests/conftest.py b/tests/conftest.py index d0f7141..016d4a1 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -739,3 +739,102 @@ def write_parse_run( ) dges[(sample, which)] = dge return ParseRun(root, ids, names, gene_genomes, dges) + + +# --------------------------------------------------------------------------- +# Loom, as loompy 3 and velocyto write it (checked on velocyto 0.17 and +# kallisto|bustools looms): genes x cells, chunked 64 x 64 and gzipped; +# strings as variable-length ASCII with XML character references; booleans +# as uint8; spec 3.0.0 keeps global attributes as datasets under /attrs, +# spec 2.0.1 as root attributes. + + +@dataclass +class LoomTruth: + path: Path + counts: np.ndarray # cells x genes, the main matrix + layers: Dict[str, np.ndarray] + cell_ids: List[str] + gene_ids: List[str] + gene_names: List[str] + + +def _loom_strings(values: Sequence[str]) -> np.ndarray: + return np.array([str(v).encode("ascii", "xmlcharrefreplace") for v in values], dtype=object) + + +def write_loom( + path: Path, + *, + n_cells: int = 30, + n_genes: int = 10, + spec: str = "3.0.0", + velocyto: bool = False, + graph: bool = True, + compound: bool = False, + seed: int = 0, +) -> LoomTruth: + rng = np.random.default_rng(seed) + ascii_str = h5py.string_dtype(encoding="ascii") + counts = rng.poisson(0.7, (n_cells, n_genes)).astype(np.float32) + layers: Dict[str, np.ndarray] = {} + if velocyto: + for name in ("spliced", "unspliced", "ambiguous"): + layers[name] = rng.poisson(0.5, (n_cells, n_genes)).astype(np.uint32) + else: + layers["counts"] = counts.astype(np.int32) * 2 + cell_ids = [f"sample:CELL{i:05d}x" for i in range(n_cells)] + gene_ids = [f"ENSMUSG{i:011d}" for i in range(n_genes)] + gene_names = [f"Gène{i % (n_genes - 1)}" for i in range(n_genes)] # one repeat, one accent + + def chunks(shape): + return (min(64, shape[0]), min(64, shape[1])) + + with h5py.File(path, "w") as f: + f.create_dataset( + "matrix", data=counts.T, chunks=chunks(counts.T.shape), + compression="gzip", compression_opts=2, + ) + group = f.create_group("layers") + for name, values in layers.items(): + group.create_dataset( + name, data=values.T, chunks=chunks(values.T.shape), + compression="gzip", compression_opts=2, + ) + col = f.create_group("col_attrs") + col.create_dataset("CellID", data=_loom_strings(cell_ids), dtype=ascii_str) + col.create_dataset("Clusters", data=np.arange(n_cells) % 3) + col.create_dataset("Tissue", data=_loom_strings(["cortex", "hippocampus"] * (n_cells // 2) + ["cortex"] * (n_cells % 2)), dtype=ascii_str) + col.create_dataset("Passed", data=(np.arange(n_cells) % 2).astype(np.uint8)) + col.create_dataset("_X", data=rng.normal(size=n_cells).astype(np.float32)) + col.create_dataset("embedding", data=rng.normal(size=(n_cells, 2))) + if compound: + dt = np.dtype([("Regulon_A", "f8"), ("Regulon_B", "f8")]) + col.create_dataset("RegulonsAUC", data=np.zeros(n_cells, dtype=dt)) + row = f.create_group("row_attrs") + row.create_dataset("Gene", data=_loom_strings(gene_names), dtype=ascii_str) + row.create_dataset("Accession", data=_loom_strings(gene_ids), dtype=ascii_str) + row.create_dataset("Start", data=np.arange(n_genes) * 100) + if velocyto: + row.create_dataset("Chromosome", data=_loom_strings(["chr1"] * n_genes), dtype=ascii_str) + cg, rgraphs = f.create_group("col_graphs"), f.create_group("row_graphs") + if graph: + g = cg.create_group("knn") + a = np.repeat(np.arange(n_cells), 2) + b = (a + np.tile([1, 2], n_cells)) % n_cells + order = rng.permutation(a.size) # COO need not be sorted + g.create_dataset("a", data=a[order]) + g.create_dataset("b", data=b[order]) + g.create_dataset("w", data=np.linspace(0.1, 1, a.size)[order].astype(np.float32)) + del rgraphs + global_attrs = {"LOOM_SPEC_VERSION": spec, "title": "Démo"} + if velocyto: + global_attrs["velocyto.__version__"] = "0.17.17" + if spec.startswith("3"): + attrs = f.create_group("attrs") + for key, value in global_attrs.items(): + attrs.create_dataset(key, data=str(value).encode("ascii", "xmlcharrefreplace"), dtype=ascii_str) + else: + for key, value in global_attrs.items(): + f.attrs[key] = str(value).encode("ascii", "xmlcharrefreplace") + return LoomTruth(path, counts, layers, cell_ids, gene_ids, gene_names) diff --git a/tests/test_convert.py b/tests/test_convert.py index 73a87db..199a570 100644 --- a/tests/test_convert.py +++ b/tests/test_convert.py @@ -235,8 +235,8 @@ def test_to_names_the_format_of_an_unfamiliar_extension(h5ad, tmp_path): (["-o", "{tmp}/o.h5ad", "--to", "zarr"], "contradicts the output's extension"), (["-o", "{tmp}/o.data", "--to", "zarr"], "needs a name ending in .zarr"), (["-o", "{tmp}/o.data"], "Cannot tell the output format"), - (["-o", "{tmp}/o.h5ad", "--to", "loom"], "Unknown format 'loom'"), - (["-o", "{tmp}/o.h5ad", "--from", "loom"], "Unknown format 'loom'"), + (["-o", "{tmp}/o.h5ad", "--to", "h5seurat"], "Unknown format 'h5seurat'"), + (["-o", "{tmp}/o.h5ad", "--from", "h5seurat"], "Unknown format 'h5seurat'"), (["-o", "-"], "cannot go to stdout"), ([], "Output path is required"), (["-o", "{tmp}/o.h5ad", "--in-memory"], "--in-memory applies to `adata recast`"), diff --git a/tests/test_external_readers.py b/tests/test_external_readers.py index 028baf2..b50fca9 100644 --- a/tests/test_external_readers.py +++ b/tests/test_external_readers.py @@ -29,6 +29,7 @@ write_cellbender_h5, write_10x_mtx, write_10x_v2_h5, + write_loom, write_parse_run, write_starsolo_velocyto, ) @@ -510,3 +511,101 @@ def test_a_parse_run_imports_as_parses_scanpy_recipe_reads_it(tmp_path, genome): problems.append(f"{column!r} differs") print("RESULT " + json.dumps(problems)) """ + + +@pytest.mark.parametrize("velocyto", [False, True], ids=["loompy", "velocyto"]) +def test_a_loom_imports_as_anndata_reads_it(tmp_path, velocyto): + """anndata.read_loom (with loompy) is the reference: X, layers, names and + attributes. Ours keeps dtypes where anndata casts X to float32, so values + are compared, not types.""" + truth = write_loom(tmp_path / "in.loom", velocyto=velocyto, compound=False) + out = tmp_path / "out.h5ad" + _convert(str(truth.path), "-o", str(out), "--var-names", "gene_symbols") + problems = _run(["anndata", "loompy", "scipy"], _COMPARE_LOOM, str(truth.path), str(out)) + assert problems == [], problems + + +_COMPARE_LOOM = r""" +import json, sys, warnings +warnings.filterwarnings("ignore") +import numpy as np, anndata as ad + +ref, got = ad.io.read_loom(sys.argv[1]), ad.read_h5ad(sys.argv[2]) +problems = [] +dense = lambda m: m.toarray() if hasattr(m, "toarray") else np.asarray(m) +names = lambda a: sorted(k for k in a.layers.keys() if k) +if ref.shape != got.shape: + problems.append(f"shape {ref.shape} != {got.shape}") +else: + if not np.array_equal(dense(ref.X), dense(got.X).astype(np.float32)): + 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") + if names(ref) != names(got): + problems.append(f"layers {names(ref)} != {names(got)}") + for key in names(ref): + if key in got.layers and not np.array_equal(dense(ref.layers[key]), dense(got.layers[key])): + problems.append(f"layers[{key!r}] differs") + for frame, mine, axis in ((ref.obs, got.obs, "obs"), (ref.var, got.var, "var")): + if sorted(frame.columns) != sorted(mine.columns): + problems.append(f"{axis} columns {sorted(frame.columns)} != {sorted(mine.columns)}") + for column in frame.columns: + if column in mine and list(frame[column].astype(str)) != list(mine[column].astype(str)): + problems.append(f"{axis}[{column!r}] differs") + for key in ref.obsm: + if key not in got.obsm or not np.allclose(ref.obsm[key], got.obsm[key]): + problems.append(f"obsm[{key!r}] differs") +print("RESULT " + json.dumps(problems)) +""" + + +@pytest.mark.parametrize("layout", ["csr", "dense"]) +def test_loompy_accepts_and_anndata_reads_what_we_export(tmp_path, layout): + """loompy's own validator passes our loom, and anndata.read_loom reads + back the matrix, layers, names and columns we wrote.""" + import anndata as ad + import pandas as pd + from scipy import sparse + + counts = _counts(40, 7).astype(np.float32) + X = sparse.csr_matrix(counts) if layout == "csr" else counts + source = tmp_path / "in.h5ad" + ad.AnnData( + X=X, + obs=pd.DataFrame({"cluster": pd.Categorical(["a", "b"] * 20)}, index=_barcodes(40)), + var=pd.DataFrame({"gene_symbols": [f"G{i}" for i in range(7)]}, index=[f"ENSG{i}" for i in range(7)]), + layers={"spliced": sparse.csr_matrix(counts * 2)}, + ).write_h5ad(source) + out = tmp_path / "out.loom" + _convert(str(source), "-o", str(out)) + problems = _run(["anndata", "loompy", "scipy"], _CHECK_LOOM_EXPORT, str(source), str(out)) + assert problems == [], problems + + +_CHECK_LOOM_EXPORT = r""" +import json, sys, warnings +warnings.filterwarnings("ignore") +import numpy as np, anndata as ad +from loompy.loom_validator import LoomValidator + +src = ad.read_h5ad(sys.argv[1]) +validator = LoomValidator() +problems = [] if validator.validate(sys.argv[2]) else [f"validator: {validator.errors}"] +got = ad.io.read_loom(sys.argv[2], var_names="Accession") +dense = lambda m: m.toarray() if hasattr(m, "toarray") else np.asarray(m) +if not np.array_equal(dense(got.X), dense(src.layers["spliced"]).astype(np.float32)): + problems.append("X (spliced, as anndata picks it) differs") +if not np.array_equal(dense(got.layers["matrix"]), dense(src.X)): + problems.append("the main matrix differs from X") +if list(got.obs_names) != list(src.obs_names): + problems.append("obs_names differ") +if list(got.var_names) != list(src.var_names): + problems.append("var_names differ") +if list(got.var["Gene"]) != list(src.var["gene_symbols"]): + problems.append("Gene differs from gene_symbols") +if list(got.obs["cluster"].astype(str)) != list(src.obs["cluster"].astype(str)): + problems.append("cluster differs") +print("RESULT " + json.dumps(problems)) +""" diff --git a/tests/test_loom.py b/tests/test_loom.py new file mode 100644 index 0000000..b4aba36 --- /dev/null +++ b/tests/test_loom.py @@ -0,0 +1,334 @@ +"""Tests for `adata convert` to and from loom. + +Fixtures are written with h5py as loompy 3 and velocyto write looms (see +conftest). Import is checked against what each fixture was built from, and +against `anndata.read_loom` in test_external_readers.py; export by reading +the file back, here and with loompy there. +""" + +from __future__ import annotations + +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_loom + +ad = pytest.importorskip("anndata") +pd = pytest.importorskip("pandas") +sparse = pytest.importorskip("scipy.sparse") + +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 _dense(m) -> np.ndarray: + return m.toarray() if sparse.issparse(m) else np.asarray(m) + + +# --------------------------------------------------------------------------- +# import + + +@pytest.mark.parametrize("suffix", [".h5ad", ".zarr"]) +@pytest.mark.parametrize("spec", ["3.0.0", "2.0.1"]) +def test_a_loom_imports_with_its_attributes_graphs_and_globals(tmp_path, spec, suffix): + truth = write_loom(tmp_path / "in.loom", spec=spec) + out = tmp_path / f"out{suffix}" + result = _ok(truth.path, "-o", out) + assert f"loom (spec {spec})" in _out(result) + assert "X = matrix, layers: counts" in _out(result) + + got = _read(out) + assert sparse.issparse(got.X) and got.X.format == "csr" + assert got.X.dtype == np.float32 + assert np.array_equal(got.X.toarray(), truth.counts) + assert got.layers["counts"].dtype == np.int32, "dtypes are kept" + assert np.array_equal(got.layers["counts"].toarray(), truth.layers["counts"]) + + assert list(got.obs_names) == truth.cell_ids + assert list(got.var_names) == truth.gene_ids, "named by gene id" + assert list(got.var["Gene"]) == truth.gene_names, "XML references unescaped" + assert got.obs["Tissue"].dtype == "category" + assert got.obs["Passed"].dtype == np.uint8, "loompy's booleans stay uint8" + assert got.obs["Clusters"].tolist() == (np.arange(30) % 3).tolist() + assert got.obsm["embedding"].shape == (30, 2) + + knn = got.obsp["knn"] + assert knn.shape == (30, 30) and knn.nnz == 60 and knn.has_sorted_indices + assert knn[0, 1] == pytest.approx(0.1) + + loom = got.uns["loom"] + assert loom["spec"] == spec + assert loom["obs_names"] == "CellID" and loom["var_names"] == "Accession" + assert loom["attrs"]["title"] == "Démo" + + +def test_a_velocyto_loom_makes_x_spliced_and_keeps_the_main_matrix(tmp_path): + truth = write_loom(tmp_path / "velo.loom", velocyto=True) + out = tmp_path / "velo.h5ad" + result = _ok(truth.path, "-o", out) + assert "loom (velocyto)" in _out(result) + got = _read(out) + assert got.X.dtype == np.uint32 + assert np.array_equal(got.X.toarray(), truth.layers["spliced"]) + assert np.array_equal(got.layers["matrix"].toarray(), truth.counts) + for name in ("spliced", "unspliced", "ambiguous"): + assert np.array_equal(got.layers[name].toarray(), truth.layers[name]) + + +def test_layer_picks_x(tmp_path): + truth = write_loom(tmp_path / "velo.loom", velocyto=True) + out = tmp_path / "u.h5ad" + _ok(truth.path, "-o", out, "--layer", "unspliced") + got = _read(out) + assert np.array_equal(got.X.toarray(), truth.layers["unspliced"]) + assert "matrix" in got.layers + + out = tmp_path / "m.h5ad" + _ok(truth.path, "-o", out, "--layer", "matrix") + assert np.array_equal(_read(out).X.toarray(), truth.counts) + + +def test_dense_keeps_the_matrices_dense(tmp_path): + truth = write_loom(tmp_path / "velo.loom", velocyto=True, n_cells=150) + out = tmp_path / "dense.h5ad" + _ok(truth.path, "-o", out, "--dense") + got = _read(out) + assert isinstance(got.X, np.ndarray) + assert np.array_equal(got.X, truth.layers["spliced"]) + assert np.array_equal(got.layers["matrix"], truth.counts) + + +def test_slabs_smaller_than_the_file_read_the_same(tmp_path, monkeypatch): + """Cells are read a slab at a time; with tiny slabs every slab boundary + falls inside the fixture.""" + from adata.conversion import loom + + monkeypatch.setattr(loom, "SLAB_BYTES", 1) + truth = write_loom(tmp_path / "in.loom", n_cells=200) + out = tmp_path / "out.h5ad" + _ok(truth.path, "-o", out) + assert np.array_equal(_read(out).X.toarray(), truth.counts) + + +def test_var_names_can_be_the_symbols_or_any_attribute(tmp_path): + truth = write_loom(tmp_path / "in.loom") + out = tmp_path / "sym.h5ad" + result = _ok(truth.path, "-o", out, "--var-names", "gene_symbols") + assert "var names repeat an earlier one" in _out(result) + got = _read(out) + assert list(got.var_names) == truth.gene_names + assert list(got.var["Accession"]) == truth.gene_ids + + out = tmp_path / "unique.h5ad" + _ok(truth.path, "-o", out, "--var-names", "gene_symbols", "--make-unique") + names = list(_read(out).var_names) + assert len(set(names)) == len(names) + + out = tmp_path / "attr.h5ad" + _ok(truth.path, "-o", out, "--obs-names-attr", "Tissue", "--var-names-attr", "Start") + got = _read(out) + assert got.obs_names[1] == "hippocampus" + assert list(got.var_names) == [str(i * 100) for i in range(10)] + + +def test_a_loom_without_names_is_numbered_as_anndata_numbers_it(tmp_path): + truth = write_loom(tmp_path / "in.loom") + with h5py.File(truth.path, "a") as f: + del f["col_attrs/CellID"], f["row_attrs/Accession"], f["row_attrs/Gene"] + out = tmp_path / "out.h5ad" + result = _ok(truth.path, "-o", out) + assert "they are numbered 0..29" in _out(result) + got = _read(out) + assert list(got.obs_names) == [str(i) for i in range(30)] + assert list(got.var_names) == [str(i) for i in range(10)] + + +def test_scope_compound_attributes_are_skipped_with_a_pointer(tmp_path): + truth = write_loom(tmp_path / "scope.loom", compound=True) + out = tmp_path / "out.h5ad" + result = _ok(truth.path, "-o", out) + assert "Not read: col_attrs RegulonsAUC" in _out(result) + assert "#27" in _out(result) + assert "RegulonsAUC" not in _read(out).obs + + +@pytest.mark.parametrize( + "spoil,extra,message", + [ + ("layer", [], "Layer 'bad' is not a 10 x 30 matrix"), + ("attr", [], "Column attribute 'short' has 3 entries, but the matrix has 30 cells"), + (None, ["--layer", "nope"], "--layer 'nope': no such layer; there are matrix (the main one), counts"), + (None, ["--var-names-attr", "nope"], "--var-names-attr 'nope': no such attribute"), + ], +) +def test_a_loom_that_cannot_be_read_is_refused_before_writing(tmp_path, spoil, extra, message): + truth = write_loom(tmp_path / "in.loom") + with h5py.File(truth.path, "a") as f: + if spoil == "layer": + f["layers"].create_dataset("bad", data=np.zeros((3, 3))) + elif spoil == "attr": + f["col_attrs"].create_dataset("short", data=np.zeros(3)) + out = tmp_path / "out.h5ad" + result = _convert(truth.path, "-o", out, *extra) + assert result.exit_code == 1 + assert message in _out(result) + assert not out.exists() + + +def test_loom_flags_are_refused_for_other_formats(tmp_path): + truth = write_loom(tmp_path / "in.loom") + store = tmp_path / "s.h5ad" + _ok(truth.path, "-o", store) + result = _convert(store, "-o", tmp_path / "o.zarr", "--dense") + assert result.exit_code == 1 + assert "does not apply" in _out(result) + + +# --------------------------------------------------------------------------- +# export + + +def _anndata_source(tmp_path: Path, layout: str): + rng = np.random.default_rng(3) + counts = rng.poisson(0.6, (130, 9)).astype(np.int32) + X = {"csr": sparse.csr_matrix, "csc": sparse.csc_matrix, "dense": np.asarray}[layout](counts) + obj = ad.AnnData( + X=X, + obs=pd.DataFrame( + {"cluster": pd.Categorical(["a", "b", "c"] * 43 + ["a"]), "score": rng.normal(size=130), + "flag": np.arange(130) % 2 == 0}, + index=[f"c{i}" for i in range(130)], + ), + var=pd.DataFrame({"gene_symbols": [f"Gène{i}" for i in range(9)]}, index=[f"ENSG{i}" for i in range(9)]), + layers={"spliced": sparse.csr_matrix(counts * 2)}, + obsm={"X_umap": rng.normal(size=(130, 2))}, + obsp={"connectivities": sparse.csr_matrix(np.eye(130, dtype=np.float32))}, + uns={"title": "tést", "n": 3}, + ) + path = tmp_path / f"{layout}.h5ad" + obj.write_h5ad(path) + return path, counts + + +@pytest.mark.parametrize("layout", ["csr", "csc", "dense"]) +def test_export_writes_a_loom_as_loompy_does(tmp_path, layout, monkeypatch): + from adata.conversion import loom + + monkeypatch.setattr(loom, "SLAB_BYTES", 1) # a slab per 64 cells: three + source, counts = _anndata_source(tmp_path, layout) + out = tmp_path / "out.loom" + result = _ok(source, "-o", out) + assert "dense" in _out(result) and "before compression" in _out(result) + with h5py.File(out, "r") as f: + assert f["matrix"].shape == (9, 130) + assert f["matrix"].chunks == (9, 64) and f["matrix"].compression == "gzip" + assert np.array_equal(f["matrix"][...], counts.T) + assert np.array_equal(f["layers/spliced"][...], counts.T * 2) + assert f["attrs/LOOM_SPEC_VERSION"][()] == b"3.0.0" + assert f["col_attrs/CellID"][0] == b"c0" + assert f["row_attrs/Accession"][0] == b"ENSG0" + assert f["row_attrs/Gene"][0] == b"Gène0", "non-ASCII as loompy escapes it" + assert f["col_attrs/flag"].dtype == np.uint8 + assert f["col_attrs/cluster"][1] == b"b" + assert f["col_attrs/X_umap"].shape == (130, 2) + assert f["col_graphs/connectivities/a"].shape == (130,) + assert f["attrs/title"][()] == b"tést" + + +def test_a_round_trip_keeps_what_a_loom_holds(tmp_path): + truth = write_loom(tmp_path / "velo.loom", velocyto=True) + h5ad, back, again = tmp_path / "a.h5ad", tmp_path / "back.loom", tmp_path / "b.h5ad" + _ok(truth.path, "-o", h5ad) + _ok(h5ad, "-o", back) + _ok(back, "-o", again) + first, second = _read(h5ad), _read(again) + assert np.array_equal(first.X.toarray(), second.X.toarray()) + assert list(first.obs_names) == list(second.obs_names) + assert list(first.var_names) == list(second.var_names) + for name in ("spliced", "unspliced", "ambiguous"): + assert np.array_equal(first.layers[name].toarray(), second.layers[name].toarray()) + assert list(first.var["Gene"]) == list(second.var["Gene"]) + assert first.obsp["knn"].nnz == second.obsp["knn"].nnz + assert second.uns["loom"]["attrs"]["velocyto.__version__"] == "0.17.17" + + +def test_a_large_export_warns_but_writes(tmp_path, monkeypatch): + from adata.conversion import loom + + monkeypatch.setattr(loom, "LARGE_DENSE", 10) + source, _ = _anndata_source(tmp_path, "csr") + result = _ok(source, "-o", tmp_path / "out.loom") + assert "A loom holds its matrices dense" in _out(result) + assert (tmp_path / "out.loom").exists() + + +def test_raw_x_exports_with_raw_var(tmp_path): + """raw/X has raw/var's genes, which filtering var leaves behind; it was + checked against var and refused, or written with var's names.""" + counts = np.arange(40, dtype=np.float32).reshape(8, 5) + full = ad.AnnData( + X=sparse.csr_matrix(counts), + obs=pd.DataFrame(index=[f"c{i}" for i in range(8)]), + var=pd.DataFrame({"gene_symbols": list("ABCDE")}, index=[f"E{i}" for i in range(5)]), + ) + full.raw = full + filtered = full[:, :3].copy() + filtered.layers["counts"] = filtered.X.copy() + source = tmp_path / "f.h5ad" + filtered.write_h5ad(source) + + out = tmp_path / "raw.loom" + result = _ok(source, "-o", out, "--layer", "raw/X") + assert "Not written: X and the layers" in _out(result) + with h5py.File(out, "r") as f: + assert f["matrix"].shape == (5, 8) + assert np.array_equal(f["matrix"][...], counts.T) + assert [v.decode() for v in f["row_attrs/Accession"][...]] == [f"E{i}" for i in range(5)] + assert [v.decode() for v in f["row_attrs/Gene"][...]] == list("ABCDE") + assert list(f["layers"]) == [] + + +def test_explicit_gene_ids_are_the_accession_whatever_the_var_names(tmp_path): + """With gene_symbols and gene_ids both columns, var_names were taken for + the ids and the gene_ids column dropped.""" + obj = ad.AnnData( + X=np.ones((2, 3), dtype=np.float32), + obs=pd.DataFrame(index=["a", "b"]), + var=pd.DataFrame( + {"gene_symbols": ["S1", "S2", "S3"], "gene_ids": ["E1", "E2", "E3"]}, + index=["key1", "key2", "key3"], + ), + ) + source = tmp_path / "s.h5ad" + obj.write_h5ad(source) + out = tmp_path / "s.loom" + _ok(source, "-o", out) + with h5py.File(out, "r") as f: + assert [v.decode() for v in f["row_attrs/Accession"][...]] == ["E1", "E2", "E3"] + assert [v.decode() for v in f["row_attrs/Gene"][...]] == ["S1", "S2", "S3"] diff --git a/tests/test_performance.py b/tests/test_performance.py index 379c2f9..21a6880 100644 --- a/tests/test_performance.py +++ b/tests/test_performance.py @@ -592,6 +592,59 @@ def measure(n: int) -> int: assert_independent_of(measure, what="convert h5ad to 10x mtx read calls", axis="n_obs") +@pytest.mark.parametrize("extra", [{}, {"dense": True}]) +def test_loom_import_work_grows_linearly(tmp_path, monkeypatch, extra): + """Slabs of cells: each chunk read once, whatever the cell count.""" + from adata.commands.convert import convert_store + from adata.conversion import Options, loom + from tests.conftest import write_loom + + monkeypatch.setattr(loom, "SLAB_BYTES", 1) + tag = "".join(extra) + + def measure(n: int) -> int: + source = write_loom(tmp_path / f"lm{tag}{n}.loom", n_cells=n, velocyto=True).path + with count_io() as io: + convert_store(source, tmp_path / f"olm{tag}{n}.h5ad", QUIET, options=Options(**extra)) + return io.work + + assert_grows_linearly(measure, what=f"convert loom ({extra}) to h5ad", axis="n_cells") + + +def test_loom_import_writes_in_blocks(tmp_path): + from adata.commands.convert import convert_store + from tests.conftest import write_loom + + def measure(n: int) -> int: + source = write_loom(tmp_path / f"lw{n}.loom", n_cells=n, velocyto=True).path + with count_io() as io: + convert_store(source, tmp_path / f"olw{n}.h5ad", QUIET) + return io.h5_write_calls + + assert_independent_of(measure, what="convert loom write calls", axis="n_cells") + + +@pytest.mark.parametrize("layout", ["csr", "csc", "dense"]) +def test_loom_export_work_grows_linearly(tmp_path, monkeypatch, layout): + from adata.commands.convert import convert_store + from adata.conversion import loom + + monkeypatch.setattr(loom, "SLAB_BYTES", 1) + + def measure(n: int) -> int: + from scipy import sparse + + counts = np.random.default_rng(0).poisson(0.5, (n, 8)).astype(np.float32) + X = {"csr": sparse.csr_matrix, "csc": sparse.csc_matrix, "dense": np.asarray}[layout](counts) + source = tmp_path / f"le{layout}{n}.h5ad" + ad.AnnData(X=X, obs=pd.DataFrame(index=[f"c{i}" for i in range(n)])).write_h5ad(source) + with count_io() as io: + convert_store(source, tmp_path / f"ole{layout}{n}.loom", QUIET) + return io.work + + assert_grows_linearly(measure, what=f"convert {layout} h5ad to loom", axis="n_obs") + + def _parse_run(path: Path, n: int) -> Path: from tests.conftest import write_parse_run