diff --git a/docs/COMMANDS.md b/docs/COMMANDS.md index 66b4471..2a40836 100644 --- a/docs/COMMANDS.md +++ b/docs/COMMANDS.md @@ -383,7 +383,10 @@ 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. +(`RegulonsAUC`) are left out with a warning; #27 will read them. Looms from +other writers read too: R's LoomExperiment (graph vertices stored as +floats), Baysor (genes in a `Name` attribute, and a `.loom` holding the +matrix alone, which is recognised by its name). 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 diff --git a/src/adata/conversion/loom.py b/src/adata/conversion/loom.py index 2fcdbf7..9537f54 100644 --- a/src/adata/conversion/loom.py +++ b/src/adata/conversion/loom.py @@ -86,7 +86,7 @@ #: 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") +GENE_SYMBOLS = ("Gene", "gene_name", "gene_names", "gene_symbols", "Name") #: Bytes of matrix held at once: one slab of cells across every gene. SLAB_BYTES = 64 << 20 @@ -108,15 +108,17 @@ def sniff(path: Path) -> bool: + """A 2-D `/matrix` dataset with `row_attrs` and `col_attrs` -- or, in a + file named `.loom`, the matrix alone, as some tools (Baysor) write it.""" 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) + if not isinstance(matrix, h5py.Dataset) or matrix.ndim != 2: + return False + return path.suffix.lower() == ".loom" or ( + isinstance(handle.get("row_attrs"), h5py.Group) and isinstance(handle.get("col_attrs"), h5py.Group) ) @@ -138,11 +140,15 @@ def decode(values: Any) -> Any: def _version(handle: h5py.File) -> str: + value = None 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" + value = handle["attrs"]["LOOM_SPEC_VERSION"][()] + elif "LOOM_SPEC_VERSION" in handle.attrs: + value = handle.attrs["LOOM_SPEC_VERSION"] + if value is None: + return "pre-2" + # R's LoomExperiment writes it as a one-element array. + return str(np.asarray(decode(value), dtype=object).reshape(-1)[0]) @dataclass @@ -156,7 +162,7 @@ class _Axis: skipped: List[str] = field(default_factory=list) -def _axis(group: h5py.Group, n: int, what: str) -> _Axis: +def _axis(group: Any, n: int, what: str) -> _Axis: axis = _Axis(n) for name, ds in group.items(): if not isinstance(ds, h5py.Dataset): @@ -262,8 +268,9 @@ def plan( # 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") + empty: Dict[str, Any] = {} + obs = _axis(handle.get("col_attrs", empty), n_cells, "Column") + var = _axis(handle.get("row_attrs", empty), n_genes, "Row") obs.names = _pick(obs, options.obs_names_attr, OBS_NAMES, "--obs-names-attr") if obs.names is None: warnings.append( @@ -298,7 +305,7 @@ def plan( ) var_names = ( - [str(v) for v in decode(handle["row_attrs"][var.names][...]).reshape(-1)] + [str(v) for v in decode(handle.get("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: @@ -364,12 +371,12 @@ def write( 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) + _write_obs(root, handle.get("col_attrs", {}), detail.obs) + _write_var(root, handle.get("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"), + (detail.obs, handle.get("col_attrs", {}), "obsm"), + (detail.var, handle.get("row_attrs", {}), "varm"), ): for name in axis.matrices: _copy_rows(group[name], root[key], name) @@ -458,10 +465,16 @@ def _write_graph(graph: h5py.Group, root: Any, key: str, name: str, n: int) -> N raise ValueError(f"Graph {name!r} has a, b and w of different lengths.") step = 1 << 20 + def coordinates(ds: h5py.Dataset, lo: int) -> np.ndarray: + values = np.asarray(ds[lo : lo + step]) + if values.dtype.kind == "f" and not np.array_equal(values, np.trunc(values)): + raise ValueError(f"Graph {name!r} has a vertex that is not a whole number.") + return values.astype(np.int64) + 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) + # R's LoomExperiment stores the vertices as floats. + rows, cols = coordinates(a, lo), coordinates(b, lo) 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]) diff --git a/tests/test_loom.py b/tests/test_loom.py index b4aba36..42da281 100644 --- a/tests/test_loom.py +++ b/tests/test_loom.py @@ -332,3 +332,63 @@ def test_explicit_gene_ids_are_the_accession_whatever_the_var_names(tmp_path): 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"] + + +# --------------------------------------------------------------------------- +# looms other tools write, as found on real files + + +def test_a_matrix_only_loom_is_read_by_its_name(tmp_path): + """Baysor writes some looms as /matrix alone, with no attribute groups; + these were not recognised as looms at all.""" + path = tmp_path / "ncvs.loom" + values = np.arange(12, dtype=np.float32).reshape(3, 4) # genes x cells + with h5py.File(path, "w") as f: + f.create_dataset("matrix", data=values) + out = tmp_path / "out.h5ad" + result = _ok(path, "-o", out) + assert "loom (spec pre-2) → h5ad, 4 cells x 3 genes" in _out(result) + got = _read(out) + assert np.array_equal(got.X.toarray(), values.T) + assert list(got.obs_names) == ["0", "1", "2", "3"] + + renamed = tmp_path / "ncvs.h5" + path.rename(renamed) + result = _convert(renamed, "-o", tmp_path / "x.h5ad") + assert result.exit_code == 1, "without the name, a bare matrix is no evidence" + + +def test_an_r_loomexperiment_file_reads_its_version_and_float_graphs(tmp_path): + """R's LoomExperiment writes LOOM_SPEC_VERSION as a one-element array + (shown as "['2.0.1']") and graph vertices as floats.""" + truth = write_loom(tmp_path / "r.loom", spec="2.0.1") + with h5py.File(truth.path, "a") as f: + f.attrs["LOOM_SPEC_VERSION"] = np.array([b"2.0.1"]) + for key in ("a", "b"): + values = f["col_graphs/knn"][key][...].astype(np.float64) + del f["col_graphs/knn"][key] + f["col_graphs/knn"].create_dataset(key, data=values) + out = tmp_path / "out.h5ad" + result = _ok(truth.path, "-o", out) + assert "loom (spec 2.0.1)" in _out(result) + assert _read(out).obsp["knn"].nnz == 60 + + with h5py.File(truth.path, "a") as f: + values = f["col_graphs/knn/a"][...] + values[0] += 0.5 + f["col_graphs/knn/a"][...] = values + result = _convert(truth.path, "-o", tmp_path / "bad.h5ad") + assert result.exit_code == 1 + assert "not a whole number" in _out(result) + + +def test_baysor_names_genes_by_name(tmp_path): + truth = write_loom(tmp_path / "b.loom") + with h5py.File(truth.path, "a") as f: + names = f["row_attrs/Gene"][...] + del f["row_attrs/Gene"], f["row_attrs/Accession"] + f["row_attrs"].create_dataset("Name", data=names) + out = tmp_path / "out.h5ad" + result = _ok(truth.path, "-o", out) + assert "genes are named by 'Name'" in _out(result) + assert list(_read(out).var_names) == truth.gene_names