diff --git a/docs/COMMANDS.md b/docs/COMMANDS.md index 0b3310e..95137b1 100644 --- a/docs/COMMANDS.md +++ b/docs/COMMANDS.md @@ -131,7 +131,7 @@ adata convert cistopic_obj.pkl -o atac.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, SCope (pySCENIC, SCENIC+, pycistopic); spec 2 and 3 | spec 3.0.0, as loompy 3 writes it | +| `loom` | `.loom` | loompy, velocyto, kallisto\|bustools, SCope (pySCENIC, SCENIC+, pycistopic); spec 2 and 3 | spec 3.0.0, as loompy 3 writes it; a SCope loom with `--scope` | | `text` | `.csv`, `.tsv`, `.tab`, `.txt`, `.data`, each also `.gz` | a dense matrix, as anndata's `read_csv` and `read_text` read it | no | | `hdf5` | `.hdf5`, `.he5`, `.hdf`; a plain `.h5` with `--from hdf5` | one 2-D dataset, as anndata's `read_hdf` reads it | no | | `cistopic` | `.pkl`, `.pickle`, gzipped too | a pycisTopic `CistopicObject` or `CistopicImputedFeatures`, without pycisTopic (needs `pyadata-cli[pickle]`) | no | @@ -449,8 +449,56 @@ where genes go, so var is named by region, and SCENIC+'s `linkedGene` is a var column. The compound attributes are read a block of cells at a time, so memory is one block across all regulons, never the table. -Writing a loom writes a plain one: an AnnData's regulons, embeddings and -clusterings are not yet written in SCope's form (#61). +#### Writing SCope looms + +`--scope` writes a loom the SCope viewer reads: an AnnData's regulons, +embeddings and clusterings in SCope's form, with a `MetaData` naming them. +Without it a loom export is plain, and one of an AnnData that came from a +SCope loom says that `--scope` would keep its regulons. + +```bash +adata convert scenic.h5ad -o scenic.loom --scope +adata convert scenic.h5ad -o scenic.loom --scope --clusterings leiden,cell_type +adata convert scenic.h5ad -o scenic.loom --scope --default-embedding X_umap +``` + +| AnnData | SCope loom | +|---|---| +| a dataframe in obsm or varm (`RegulonsAUC`, `Regulons`, `ClusterMarkers_`) | a compound `col_attrs` or `row_attrs` of its numeric columns | +| `obsm` arrays of two columns | `Embeddings_X` and `Embeddings_Y`, a field per id; the default (id -1) as `Embedding` too | +| the categorical obs columns a SCope import recorded, and those `--clusterings` names | `Clusterings`, int codes per clustering (and not a plain attribute too), with `ClusterID` for the first; their clusters described in `MetaData` | +| other obs columns | plain attributes, as ever, listed in `MetaData` as metrics (numbers) and annotations (up to 245 values) | +| `uns["scope"]` | `MetaData`'s regulon thresholds (from the `regulon_thresholds` table, as edited), and the recorded ids and names | + +| Flag | Meaning | +|---|---| +| `--scope` | Write a SCope loom | +| `--clusterings` | Categorical obs columns to write as clusterings, comma separated | +| `--default-embedding` | The obsm key SCope shows first (default: the recorded default, else `X_umap`, else `X_tsne`, else the first) | + +An AnnData that came from a SCope loom keeps its embeddings' and +clusterings' ids and names, recorded in `uns["scope"]`, so a SCope loom +converts to h5ad and back with the same attributes. As SCope expects, +clusterings and their clusters are numbered 0..n-1 in order, regulon names +lose their spaces (`TP53 (5g)` becomes `TP53_(5g)`, alike in the AUC and +the Regulons, with a warning), and cluster markers are written only as +`ClusterMarkers_` with their metrics (`ClusterMarkers__avg_logFC` and the like, which +`clusterMarkerMetrics` then lists, under the names a SCope import recorded), their fields renumbered with the clusters; a +clustering `--clusterings` names owns `ClusterMarkers_`. `ClusterID` is always written afresh from the first +clustering, in place of an obs column of that name, so the two cannot +disagree. `MetaData` is plain JSON, which SCope reads as it +reads compressed. A clustering that is not categorical is refused, as is a +table with more columns than HDF5 can hold in one compound attribute (its +type must fit 64 KiB, a limit pySCENIC's own export shares), before any +output exists. Attributes are written a block of cells at a time. A loom's +attributes share one flat namespace per axis, with or without `--scope`: +where an obs column, an obsm entry or one of SCope's parts ask for one name, +the first keeps it (loom's and SCope's own attributes, then columns, SCope's +tables, arrays) and the others are left out with a warning. + +SCope opens looms without loompy's validation, which rejects compound +attributes; so do `loompy.connect(..., validate=False)` and the pySCENIC +tutorial's recipe, which CI uses to read our SCope looms back. ### Text tables diff --git a/docs/TESTING.md b/docs/TESTING.md index ed01ef4..c763039 100644 --- a/docs/TESTING.md +++ b/docs/TESTING.md @@ -36,7 +36,7 @@ separate jobs. | `test_cistopic.py` | pycisTopic pickles read through the allowlist: a CistopicObject's matrices, frames, model and projections, an imputed-features pickle, a pickle naming `os.system` refused unrun, the missing-extra hint | | `test_pandas_frame.py` | Each pandas dtype written as anndata writes it, on every backend | | `test_concat_axis_elements.py` | `concat` on obsm dataframes and varm, against `anndata.concat` for each join and merge strategy | -| `test_scope.py` | SCope looms as pySCENIC's `export2loom` writes them: regulons as dataframes, embeddings and clusterings named from `MetaData`, plain or compressed `MetaData`, region looms, R's logicals and factors, renamed fields, small blocks | +| `test_scope.py` | SCope looms, read as pySCENIC's `export2loom` writes them and written with `--scope` (round trips, plain AnnData, refusals, markers): regulons as dataframes, embeddings and clusterings named from `MetaData`, plain or compressed `MetaData`, region looms, R's logicals and factors, renamed fields, small blocks | | `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, Parse runs and looms, checked by scanpy (Parse's own recipe for Parse), anndata with loompy (the pySCENIC tutorial's recipe for SCope looms), anndata's `read_csv`, `read_text` and `read_hdf`, a pandas 1.5 pycisTopic pickle read back there, and PyTables in their own uv environments (`integration`) | diff --git a/src/adata/cli.py b/src/adata/cli.py index eb4dadd..1a1cd83 100644 --- a/src/adata/cli.py +++ b/src/adata/cli.py @@ -537,6 +537,9 @@ def concat( "key": "--key", "obs_names_key": "--obs-names-key", "var_names_key": "--var-names-key", + "scope": "--scope", + "clusterings": "--clusterings", + "default_embedding": "--default-embedding", } #: Flags that only make sense for format conversion; with recast's flags, @@ -738,6 +741,24 @@ def convert( help="HDF5 input: the dataset naming the cells (default row_names, " "col_names with --transpose)", ), + scope: bool = typer.Option( + False, + "--scope", + help="loom output: write a SCope loom: regulons, embeddings, clusterings " + "and MetaData, as the SCope viewer reads them", + ), + clusterings: Optional[str] = typer.Option( + None, + "--clusterings", + help="SCope output: categorical obs columns written as clusterings, " + "comma separated (besides those a SCope import recorded)", + ), + default_embedding: Optional[str] = typer.Option( + None, + "--default-embedding", + help="SCope output: the obsm key SCope shows first (default: the " + "recorded default, else X_umap, else X_tsne, else the first)", + ), var_names_key: Optional[str] = typer.Option( None, "--var-names-key", @@ -773,7 +794,7 @@ def convert( is read as anndata.read_loom reads it; a SCope loom (pySCENIC, SCENIC+, pycistopic) also gives its regulons as dataframes in obsm and varm, its embeddings and clusterings named from its MetaData, and the MetaData in - uns["scope"]. Features are named by gene id, unique but in probe-level + uns["scope"]; --scope writes them back as SCope reads them. Features are named by gene id, unique but in probe-level (Flex) matrices; --var-names gene_symbols names them as scanpy does. Nothing is written until every check has passed, and the output is moved @@ -800,6 +821,7 @@ def convert( adata convert velocyto.loom -o velocyto.h5ad adata convert pyscenic_output.loom -o scenic.h5ad adata convert counts.h5ad -o counts.loom + adata convert scenic.h5ad -o scenic.loom --scope --clusterings leiden adata convert GSE123_counts.txt.gz -o counts.h5ad --transpose adata convert matrix.hdf5 -o matrix.h5ad --key data/counts adata convert cistopic_obj.pkl -o atac.h5ad @@ -900,6 +922,9 @@ def convert( key=key, obs_names_key=obs_names_key, var_names_key=var_names_key, + scope=scope, + clusterings=clusterings, + default_embedding=default_embedding, ), ) except Exception as e: diff --git a/src/adata/conversion/base.py b/src/adata/conversion/base.py index bf6539d..187043a 100644 --- a/src/adata/conversion/base.py +++ b/src/adata/conversion/base.py @@ -59,6 +59,11 @@ class Options: obs_names_key: Optional[str] = None var_names_key: Optional[str] = None + # Writing a SCope loom. + scope: bool = False + clusterings: Optional[str] = None + default_embedding: Optional[str] = None + @dataclass class ConversionPlan: diff --git a/src/adata/conversion/loom.py b/src/adata/conversion/loom.py index afb41c3..90b724c 100644 --- a/src/adata/conversion/loom.py +++ b/src/adata/conversion/loom.py @@ -38,6 +38,7 @@ from __future__ import annotations import html +import json from contextlib import ExitStack from dataclasses import dataclass, field from datetime import datetime, timezone @@ -741,10 +742,12 @@ class _Export: #: var, or raw/var when the matrix is raw/X. var: str = "var" varm_key: str = "varm" + #: SCope's parts, with --scope. + scope: Optional[scope.ScopeExport] = None class LoomExport(Converter): - options = frozenset({"--layer"}) + options = frozenset({"--layer", "--scope", "--clusterings", "--default-embedding"}) def plan( self, @@ -755,6 +758,9 @@ def plan( options: Options, ) -> ConversionPlan: warnings: List[str] = [] + for flag in ("--clusterings", "--default-embedding"): + if flag in options.given and not options.scope: + raise ValueError(f"{flag} shapes a SCope loom: pass --scope too.") path = options.layer or "X" if path not in ("X", "raw/X") and not path.startswith("layers/"): path = f"layers/{path}" @@ -790,16 +796,54 @@ def plan( 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) + obsm = _arrays(root, "obsm", n_obs, warnings, frames=options.scope) + varm = _arrays(root, varm_key, n_var, warnings, frames=options.scope) obsp = _graphs(root, "obsp", n_obs, warnings) varp = [] if raw else _graphs(root, "varp", n_var, warnings) uns = _uns_scalars(root, warnings) + # A clustering goes out inside Clusterings alone (a plain copy too + # would come back from the loom as a second column), so it claims + # no attribute name of its own. + clustered = set(scope.clustering_columns(root, options.clusterings)) if options.scope else set() + obs_columns, var_columns, owners = _column_owners( + [c for c in obs_columns if c not in clustered], var_columns, gene, accession, + options.scope, warnings, + ) + found = None + if options.scope: + found = scope.plan_export( + root, obs_columns=[c for c in obs_columns if c != "CellID"], + varm_key=varm_key, clusterings=options.clusterings, + default_embedding=options.default_embedding, warnings=warnings, + owners=owners, + ) + if found.cluster_id and "ClusterID" in obs_columns: + # A stale copy would disagree with Clusterings once the + # first clustering is edited. + obs_columns = [c for c in obs_columns if c != "ClusterID"] + warnings.append("obs column 'ClusterID' is written afresh from the " + f"first clustering, {found.clusterings[0].column!r}.") + # Embeddings go out as SCope's, not again as plain attributes. + embedded = {e.key for e in found.embeddings} + obsm = [name for name in obsm if name not in embedded] + elif "uns" in root and "scope" in root["uns"]: + warnings.append( + "This came from a SCope loom: pass --scope to write its regulons, " + "embeddings and clusterings as SCope reads them; a plain loom " + "leaves them out." + ) + + # Arrays last: every name they could want is settled by now. + obsm = [k for k in obsm if scope.claim_attr(owners, "col_attrs", k, f"obsm/{k}", warnings)] + varm = [k for k in varm + if scope.claim_attr(owners, "row_attrs", k, f"{varm_key}/{k}", 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"{src_format} → loom 3.0.0{' (SCope)' if found is not None else ''}, " + f"{n_var:,} genes x {n_obs:,} cells from " + f"{path}{', layers: ' + ', '.join(layers) if layers else ''}" + f"{', ' + found.describe() if found is not None and found.describe() else ''}; " f"dense {_size(dense)} before compression" ) if dense > LARGE_DENSE: @@ -819,7 +863,7 @@ def plan( 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, + varm_key=varm_key, scope=found, ), ) @@ -870,12 +914,76 @@ def write( 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")) + if detail.scope is not None: + _export_scope(root, cols, rows, detail.scope, detail.varm_key) + # Before uns's own attributes, which a stale MetaData among + # them must not replace. + _attr(attrs, scope.METADATA, json.dumps(scope.json_safe(detail.scope.metadata))) for name in detail.uns: value = _uns_value(root, name) if name not in ("LOOM_SPEC_VERSION", "CreationDate"): _attr(attrs, name, value) +def _export_scope( + root: Any, cols: h5py.Group, rows: h5py.Group, found: scope.ScopeExport, varm_key: str +) -> None: + """SCope's compound attributes, a block of rows at a time across their + fields: the dataframes, the embeddings, the clusterings.""" + n_obs = element_len(resolve_index(root["obs"], "obs")[0]) + for table in found.tables: + source = resolve(root, table.source) + target = cols if table.axis == "col" else rows + n = element_len(resolve_index(source)[0]) + _write_records(target, table.key, n, table.dtype, [ + (field, (lambda lo, hi, c=column: np.asarray(source[c][lo:hi]))) + for field, column in zip(table.fields, table.columns) + ]) + obsm = root["obsm"] if "obsm" in root else {} + if found.embeddings: + default = found.embeddings[0] + for axis, attr in ((0, "Embeddings_X"), (1, "Embeddings_Y")): + dtype = np.dtype([(str(e.id), np.dtype(obsm[e.key].dtype.str)) for e in found.embeddings]) + _write_records(cols, attr, n_obs, dtype, [ + (str(e.id), (lambda lo, hi, k=e.key, a=axis: np.asarray(obsm[k][lo:hi, a]))) + for e in found.embeddings + ]) + if default.id == -1: + dtype = np.dtype([(f, np.dtype(obsm[default.key].dtype.str)) for f in ("_X", "_Y")]) + _write_records(cols, "Embedding", n_obs, dtype, [ + (f, (lambda lo, hi, a=axis: np.asarray(obsm[default.key][lo:hi, a]))) + for axis, f in enumerate(("_X", "_Y")) + ]) + if found.clusterings: + obs = root["obs"] + dtype = np.dtype([(str(c.id), np.int64) for c in found.clusterings]) + codes = [ + (str(c.id), (lambda lo, hi, col=c.column: np.asarray(scope.categorical_codes(obs[col])[lo:hi], dtype=np.int64))) + for c in found.clusterings + ] + _write_records(cols, "Clusterings", n_obs, dtype, codes) + if found.cluster_id: + first = codes[0][1] + target = cols.create_dataset("ClusterID", shape=(n_obs,), dtype=np.int64) + for lo in range(0, n_obs, STRING_STEP): + target[lo : lo + STRING_STEP] = first(lo, lo + STRING_STEP) + + +def _write_records( + group: h5py.Group, name: str, n: int, dtype: np.dtype, + fields: Sequence[Tuple[str, Callable[[int, int], np.ndarray]]], +) -> None: + """A compound attribute of `n` records, each field read a block at a time.""" + target = group.create_dataset(name, shape=(n,), dtype=dtype) + step = max(1, min(STRING_STEP, SLAB_BYTES // max(1, dtype.itemsize))) + for lo in range(0, n, step): + hi = min(lo + step, n) + block = np.empty(hi - lo, dtype=dtype) + for field, read in fields: + block[field] = read(lo, hi) + target[lo:hi] = block + + def _size(n: int) -> str: for unit, scale in (("TB", 1 << 40), ("GB", 1 << 30), ("MB", 1 << 20)): if n >= scale: @@ -928,6 +1036,36 @@ def _columns(group: Any, index_name: str, warnings: List[str], axis: str) -> Lis return kept +def _column_owners( + obs_columns: List[str], var_columns: List[str], gene: str, accession: Optional[str], + scope_export: bool, warnings: List[str], +) -> Tuple[List[str], List[str], Dict[str, Dict[str, str]]]: + """A loom's col_attrs and row_attrs are flat, so an obs column, an obsm + entry and one of SCope's parts can ask for one name, and the second + write failed with the output half written. Names go to the first that + asks: loom's and SCope's own attributes, then columns here; SCope's + tables (`scope.plan_export`) and arrays claim theirs after. The others + are left out, each with a warning.""" + owners: Dict[str, Dict[str, str]] = {"col_attrs": {"CellID": "the cell names"}, + "row_attrs": {"Gene": "the gene names"}} + if accession is not None: + owners["row_attrs"]["Accession"] = "the gene ids" + if scope_export: + owners["col_attrs"].update({ + scope.EMBEDDING: "SCope's default embedding", scope.EMBEDDINGS_X: "SCope's embeddings", + scope.EMBEDDINGS_Y: "SCope's embeddings", scope.CLUSTERINGS: "SCope's clusterings", + }) + # ClusterID is SCope's only once there are clusterings; the plan decides. + obs_columns = [c for c in obs_columns if c in ("CellID", "ClusterID") + or scope.claim_attr(owners, "col_attrs", c, f"obs column {c!r}", warnings)] + skip = {gene, accession, "Gene", "Accession"} + var_columns = [c for c in var_columns if c in skip + or scope.claim_attr(owners, "row_attrs", c, f"var column {c!r}", warnings)] + if scope_export: + owners["col_attrs"].setdefault("ClusterID", "SCope's first clustering") + return obs_columns, var_columns, owners + + def _gene_attrs(columns: List[str]) -> Tuple[str, Optional[str]]: """Where loom's Gene (symbol) and Accession (id) come from. @@ -946,13 +1084,15 @@ def _gene_attrs(columns: List[str]) -> Tuple[str, Optional[str]]: return symbols, "" if ids is None else ids -def _arrays(root: Any, key: str, n: int, warnings: List[str]) -> List[str]: +def _arrays(root: Any, key: str, n: int, warnings: List[str], *, frames: bool = False) -> List[str]: + """The numeric 2-D arrays a loom attribute can hold; with `frames`, the + dataframes are not reported, since a SCope loom writes them.""" 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: + elif not (frames and spec.encoding_type(obj) == spec.DATAFRAME): dropped.append(name) if dropped: warnings.append(f"Not written: {key} {', '.join(dropped)} -- only numeric 2-D arrays fit a loom.") diff --git a/src/adata/conversion/scope.py b/src/adata/conversion/scope.py index 911ec96..596bff2 100644 --- a/src/adata/conversion/scope.py +++ b/src/adata/conversion/scope.py @@ -26,6 +26,7 @@ import base64 import json +import re import zlib from dataclasses import dataclass, field from typing import Any, Dict, Iterable, Iterator, List, Optional, Sequence, Set, Tuple @@ -42,7 +43,7 @@ write_scalar, write_string_array, ) -from adata.storage import create_dataset +from adata.storage import create_dataset, is_dataset, is_group from adata.util.names import claim, make_index_unique, sanitize_label, storable METADATA = "MetaData" @@ -395,3 +396,500 @@ def _write_thresholds(parent: Any, entries: Any) -> None: else: columns[key] = ["" if row.get(key) is None else str(row[key]) for row in rows] _write_frame(parent, "regulon_thresholds", [str(row["regulon"]) for row in rows], columns) + + +# --------------------------------------------------------------------------- +# export: an AnnData as a SCope loom (`adata convert ... -o x.loom --scope`) +# +# What SCope's server reads, and must find consistent, decides the rules: +# compound attributes with numeric fields; embedding ids as field names, -1 +# the default; clusterings as int codes with every cluster described in +# MetaData, clustering and cluster ids 0..n-1 in list order; regulon names +# without spaces, alike in the AUC and the Regulons attributes; markers only +# with their metrics described. An AnnData that came from a SCope loom keeps +# its ids and names, recorded in uns["scope"] on import. + +#: Categories an obs column may have to be offered as a SCope annotation, the +#: length of SCope's colour list. +MAX_ANNOTATION_VALUES = 245 + +_MARKERS = re.compile(r"ClusterMarkers_(\d+)(?:_(.+))?$") + + +@dataclass +class ExportTable: + #: "col" or "row": which attribute group it goes to. + axis: str + key: str + #: The obsm or varm entry it comes from, and its numeric columns. + source: str + columns: List[str] + fields: List[str] + dtype: np.dtype + + +@dataclass +class ExportEmbedding: + id: int + name: str + key: str + + +@dataclass +class ExportClustering: + id: int + column: str + entry: Dict[str, Any] + + +@dataclass +class ScopeExport: + tables: List[ExportTable] + embeddings: List[ExportEmbedding] + clusterings: List[ExportClustering] + metadata: Dict[str, Any] + cluster_id: bool + regulons: int + + def describe(self) -> str: + counts = ((self.regulons, "regulon"), (len(self.embeddings), "embedding"), + (len(self.clusterings), "clustering")) + return ", ".join(f"{n:,} {noun}{'' if n == 1 else 's'}" for n, noun in counts if n) + + +def compound_fits(dtype: np.dtype) -> bool: + """Whether HDF5 can hold `dtype`: a compound type must fit one object + header message, 64 KiB, so a few hundred to a thousand fields.""" + import io + + with h5py.File(io.BytesIO(), "w") as handle: + try: + handle.create_dataset("probe", shape=(1,), dtype=dtype) + except ValueError: + return False + return True + + +def _frame_rows(group: Any) -> Dict[str, Dict[str, str]]: + """A small text dataframe (uns["scope"]["embeddings"]) as {index: {column: value}}.""" + from adata.elements.read import dataframe_columns, read_str_all, resolve_index + + index, name = resolve_index(group) + labels = read_str_all(index) + rows: Dict[str, Dict[str, str]] = {label: {} for label in labels} + for column in dataframe_columns(group, name): + for label, value in zip(labels, read_str_all(group[column], parent_group=group)): + rows[label][column] = value + return rows + + +def _recorded(root: Any) -> Tuple[Dict[str, Any], Dict[str, Tuple[int, str]], Dict[str, Tuple[int, str]]]: + """What a SCope import recorded: MetaData's JSON, and the obsm keys and obs + columns of its embeddings and clusterings, each with its id and name.""" + uns = root["uns"] if "uns" in root else {} + group = uns["scope"] if "scope" in uns else None + parsed: Dict[str, Any] = {} + embeddings: Dict[str, Tuple[int, str]] = {} + clusterings: Dict[str, Tuple[int, str]] = {} + if group is None: + return parsed, embeddings, clusterings + if "metadata" in group: + decoded, _ = read_metadata(_scalar_text(group["metadata"])) + parsed = decoded or {} + for key, target, slot in (("embeddings", embeddings, "obsm"), ("clusterings", clusterings, "obs")): + if key in group: + for id, row in _frame_rows(group[key]).items(): + number = _number(id) + if number is not None and row.get(slot): + target[row[slot]] = (int(number), row.get("name", "")) + return parsed, embeddings, clusterings + + +def _scalar_text(ds: Any) -> str: + value = ds[()] + return value.decode("utf-8") if isinstance(value, bytes) else str(value) + + +def _check_fits(fields: Sequence[Tuple[str, Any]], what: str, advice: str) -> None: + """Refuse, in the plan, a compound attribute HDF5 could not create.""" + if fields and not compound_fits(np.dtype(list(fields))): + raise ValueError( + f"{what} are more than HDF5 can hold as one compound attribute " + f"(its type must fit 64 KiB). {advice}" + ) + + +def clustering_columns(root: Any, clusterings: Optional[str]) -> List[str]: + """The obs columns to write as clusterings: those a SCope import + recorded, then those `--clusterings` names (checked by the plan).""" + obs = root["obs"] + wanted = [c for c in _recorded(root)[2] if c in obs] + for column in (clusterings or "").split(","): + column = column.strip() + if column and column not in wanted: + wanted.append(column) + return wanted + + +def plan_export( + root: Any, + *, + obs_columns: Sequence[str], + varm_key: str, + clusterings: Optional[str], + default_embedding: Optional[str], + warnings: List[str], + owners: Optional[Dict[str, Dict[str, str]]] = None, +) -> ScopeExport: + """Everything a SCope loom needs, worked out and checked before writing. + + `owners` holds the attribute names already given away (`claim_attr`): a + table asking for one is left out before anything is derived from it, so + the metadata only ever describes tables that are written. + """ + from adata.elements.read import categories_dataset, dataframe_columns, read_categories, resolve_index + + parsed, recorded_embeddings, recorded_clusterings = _recorded(root) + owners = owners if owners is not None else {"col_attrs": {}, "row_attrs": {}} + + # Compound tables: every dataframe in obsm and varm, numeric columns only. + frames = [] + for axis, where in (("col", "obsm"), ("row", varm_key)): + group = _resolve(root, where) + for key in (group if group is not None else []): + obj = group[key] + if spec.encoding_type(obj) != spec.DATAFRAME: + continue + numeric, other = [], [] + for column in dataframe_columns(obj, resolve_index(obj)[1]): + member = obj[column] + ok = (spec.encoding_type(member) in (None, spec.ARRAY) and getattr(member, "ndim", 0) == 1 + and member.dtype.kind in "biuf") + (numeric if ok else other).append(column) + if other: + warnings.append(f"Not written in {where}/{key}: {', '.join(other)} -- a SCope " + "table holds numbers only.") + # Marker tables claim their names once renumbered, below. + if numeric and (axis == "row" and _MARKERS.match(key) or claim_attr( + owners, "col_attrs" if axis == "col" else "row_attrs", key, f"{where}/{key}", warnings + )): + frames.append((axis, where, key, obj, numeric)) + + # One field name per regulon across AUC and Regulons, so that when two + # names meet once spaces go (`A B`, `A_B`) a regulon's scores and its + # targets still share a field. + is_regulons = lambda key: key.endswith("Regulons") or key.endswith("RegulonsAUC") # noqa: E731 + regulon_fields: Dict[str, str] = {} + taken_regulons: Set[str] = set() + for _, _, key, _, numeric in frames: + if is_regulons(key): + for c in numeric: + if c not in regulon_fields: + regulon_fields[c] = claim(c.replace(" ", "_"), taken_regulons) + + tables: List[ExportTable] = [] + for axis, where, key, obj, numeric in frames: + regulons = is_regulons(key) + if regulons: + fields = [regulon_fields[c] for c in numeric] + else: + taken: Set[str] = set() + fields = [claim(c, taken) for c in numeric] + renamed = [f"{c!r} as {f!r}" for c, f in zip(numeric, fields) if c != f] + if renamed: + warnings.append(f"Renamed in {key}: {', '.join(renamed)} -- SCope's regulon names hold no spaces.") + dtype = np.dtype([(f, np.dtype(obj[c].dtype.str)) for f, c in zip(fields, numeric)]) + if not compound_fits(dtype): + raise ValueError( + f"{where}/{key} has {len(fields):,} columns, more than HDF5 can hold " + "as one compound attribute (its type must fit 64 KiB); pySCENIC's own " + "export stops at the same limit. Keep fewer columns first." + ) + tables.append(ExportTable(axis, key, f"{where}/{key}", numeric, fields, dtype)) + + # Embeddings: obsm arrays of two numeric columns. + obsm = root["obsm"] if "obsm" in root else {} + candidates = [ + key for key in obsm + if spec.encoding_type(obsm[key]) in (None, spec.ARRAY) and getattr(obsm[key], "ndim", 0) == 2 + and obsm[key].shape[1] == 2 and obsm[key].dtype.kind in "biuf" + ] + if default_embedding is not None and default_embedding not in candidates: + raise ValueError( + f"--default-embedding {default_embedding}: not an obsm array of two columns; " + f"there are {', '.join(candidates) or 'none'}." + ) + default = default_embedding or next( + (k for k in candidates if recorded_embeddings.get(k, (None,))[0] == -1), None + ) or next((k for k in ("X_umap", "X_UMAP", "X_tsne", "X_tSNE") if k in candidates), None) or ( + candidates[0] if candidates else None + ) + used = {-1} + ids: Dict[str, int] = {} + for key in candidates: + recorded = recorded_embeddings.get(key, (None,))[0] + if key != default and recorded is not None and recorded not in used: + ids[key] = recorded + used.add(recorded) + next_id = 0 + embeddings: List[ExportEmbedding] = [] + for key in candidates: + if key == default: + id = -1 + elif key in ids: + id = ids[key] + else: + while next_id in used: + next_id += 1 + id = next_id + used.add(id) + name = recorded_embeddings.get(key, (None, ""))[1] or (key[2:] if key.startswith("X_") else key) + embeddings.append(ExportEmbedding(id, name, key)) + embeddings.sort(key=lambda e: e.id) + _check_fits([(str(e.id), np.dtype(obsm[e.key].dtype.str)) for e in embeddings], + f"{len(embeddings):,} embeddings (two-column obsm arrays)", "Keep fewer of them first.") + + # Clusterings: the recorded ones, and those --clusterings names. + obs = root["obs"] + wanted = clustering_columns(root, clusterings) + entries = {str(e.get("id")): e for e in parsed.get("clusterings", []) if isinstance(e, dict)} + found: List[ExportClustering] = [] + renumbered: Dict[int, int] = {} + #: Per clustering written, its recorded cluster ids (as marker fields + #: name them) to the new ones, its categories' positions. + cluster_fields: Dict[int, Dict[str, str]] = {} + #: Per clustering written, the marker metrics its import recorded. + recorded_metrics: Dict[int, List[Dict[str, Any]]] = {} + for position, column in enumerate(wanted): + if column not in obs: + raise ValueError(f"--clusterings {column}: no such obs column.") + if not is_categorical(obs[column]): + raise ValueError( + f"--clusterings {column}: not a categorical obs column, which SCope's " + "clusterings need (cluster ids with descriptions)." + ) + old, name = recorded_clusterings.get(column, (None, column)) + recorded = entries.get(str(old), {}) if old is not None else {} + if old is not None: + renumbered[old] = position + categories = [str(c) for c in read_categories(obs[column], parent_group=obs)] + if old is not None: + # The import labelled the categories by description, in id order. + clusters = _clusters(recorded.get("clusters")) + labels = make_index_unique([label for _, label in clusters]) + position_of = {c: i for i, c in enumerate(categories)} + cluster_fields[position] = { + _label(number): str(position_of[label]) + for (number, _), label in zip(clusters, labels) if label in position_of + } if clusters else { + # Undescribed clusters: the import labelled them by their ids. + label: str(i) for label, i in position_of.items() + } + entry = { + "id": position, + "group": recorded.get("group", ""), + "name": name or column, + "clusters": [{"id": i, "description": c} for i, c in enumerate(categories)], + } + recorded_metrics[position] = [m for m in recorded.get("clusterMarkerMetrics") or [] + if isinstance(m, dict) and m.get("accessor")] + found.append(ExportClustering(position, column, entry)) + _check_fits([(str(c.id), np.int64) for c in found], + f"{len(found):,} clusterings", "Name fewer with --clusterings.") + # A clustering no import recorded owns the marker tables of its new id. + for clustering in found: + if clustering.id not in cluster_fields: + renumbered.setdefault(clustering.id, clustering.id) + + # Cluster markers: only for a clustering written, and with its metrics. + markers = [t for t in tables if t.axis == "row" and _MARKERS.match(t.key)] + keep: List[ExportTable] = [t for t in tables if t not in markers] + by_old: Dict[int, List[ExportTable]] = {} + for table in markers: + by_old.setdefault(int(_MARKERS.match(table.key).group(1)), []).append(table) + for old, group_tables in sorted(by_old.items()): + new = renumbered.get(old) + if new is None: + warnings.append(f"Not written: ClusterMarkers_{old} -- no clustering written has id {old}.") + continue + entry = found[new].entry + for table in group_tables: + accessor = _MARKERS.match(table.key).group(2) + table.key = f"ClusterMarkers_{new}" + (f"_{accessor}" if accessor else "") + named = [t for t in group_tables if claim_attr(owners, "row_attrs", t.key, t.source, warnings)] + base = f"ClusterMarkers_{new}" + if len(named) < len(group_tables) and any(t.key == base for t in group_tables) \ + and not any(t.key == base for t in named): + warnings.append(f"Not written: the rest of {base} -- its table of markers is left out.") + continue + # Fields are cluster ids, renumbered with the clusters; a table left + # with none is not written, so the metrics are worked out after. + fields_map = cluster_fields.get(new) + group_tables = [] + for table in named: + if fields_map is not None: + gone = [f for f in table.fields if f not in fields_map] + if gone: + warnings.append(f"Not written in {table.source}: {', '.join(gone)} -- no " + f"cluster of {found[new].column!r} has that id.") + kept = [(fields_map[f], c, table.dtype[f]) + for f, c in zip(table.fields, table.columns) if f in fields_map] + table.fields = [f for f, _, _ in kept] + table.columns = [c for _, c, _ in kept] + table.dtype = np.dtype([(f, t) for f, _, t in kept]) + if table.fields: + group_tables.append(table) + accessors = [m.group(2) for t in group_tables if (m := _MARKERS.match(t.key)).group(2)] + # The metrics of the tables written: those recorded, named as they + # were, then any table added since, named after its accessor. + metrics = [m for m in recorded_metrics.get(new, []) if m["accessor"] in accessors] + listed = {m["accessor"] for m in metrics} + metrics += [{"accessor": a, "name": a, "description": ""} for a in accessors if a not in listed] + if not metrics or not any(t.key == base for t in group_tables): + warnings.append( + f"Not written: ClusterMarkers_{old} -- without metrics (_avg_logFC, " + "_pval, ...) or markers, SCope cannot show them." + ) + continue + entry["clusterMarkerMetrics"] = metrics + keep.extend(group_tables) + tables = keep + + # Metrics and annotations: the obs columns written; clusterings are SCope's own. + clustered = {c.column for c in found} + metrics, annotations = [], [] + for column in obs_columns: + if found and column == "ClusterID" or column in clustered: + continue # written afresh from the first clustering; a clustering + member = obs[column] + enc = spec.encoding_type(member) + if is_categorical(member): + if len(categories_dataset(member, obs)) > MAX_ANNOTATION_VALUES: + continue # too many to list, and not read to find that out + values = [str(c) for c in read_categories(member, parent_group=obs)] + elif enc in (None, spec.ARRAY, spec.STRING_ARRAY) and getattr(member, "ndim", 0) == 1 \ + and member.dtype.kind in "biuf": + metrics.append({"name": column}) + continue + elif enc in (spec.STRING_ARRAY, spec.NULLABLE_STRING_ARRAY) or ( + getattr(member, "dtype", None) is not None and member.dtype.kind in "OSUT" + ): + values = _distinct_text(member, obs, MAX_ANNOTATION_VALUES) + if values is None: + continue + elif enc in (spec.NULLABLE_INTEGER, spec.NULLABLE_BOOLEAN): + metrics.append({"name": column}) + continue + else: + continue + if len(values) <= MAX_ANNOTATION_VALUES: + annotations.append({"name": column, "values": values}) + + # Regulon thresholds, for the regulons written. + auc_fields = {f for t in tables if t.axis == "col" and t.key.endswith("RegulonsAUC") for f in t.fields} + field_of = lambda regulon: regulon_fields.get(regulon, regulon.replace(" ", "_")) # noqa: E731 + recorded_thresholds: Dict[str, Dict[str, Any]] = {} + for entry in parsed.get("regulonThresholds", []) or []: + if isinstance(entry, dict) and field_of(str(entry.get("regulon", ""))) in auc_fields: + recorded_thresholds[field_of(str(entry["regulon"]))] = entry + table = _resolve(root, "uns/scope/regulon_thresholds") + if table is not None: + # The table is what may have been edited since the import: its rows + # and values win, the metadata adding what the table does not hold + # (allThresholds, and fields a regulon never had). + thresholds = [] + for regulon, row in _frame_rows(table).items(): + name = field_of(regulon) + if name not in auc_fields: + continue + entry = dict(recorded_thresholds.get(name, {}), regulon=name) + for key, value in row.items(): + if key == "defaultThresholdValue": + number = _float(value) + if np.isfinite(number): + entry[key] = number + else: + entry.pop(key, None) # missing in the table + elif value != "" or key in entry: + entry[key] = value + thresholds.append(entry) + else: + thresholds = [dict(e, regulon=name) for name, e in recorded_thresholds.items()] + + metadata = dict(parsed) + metadata.update( + embeddings=[{"id": e.id, "name": e.name} for e in embeddings], + clusterings=[c.entry for c in found], + annotations=annotations, + metrics=metrics, + regulonThresholds=thresholds, + ) + regulons = sum(len(t.fields) for t in tables if t.axis == "col" and t.key.endswith("RegulonsAUC")) or sum( + len(t.fields) for t in tables if t.axis == "row" and t.key.endswith("Regulons") + ) + return ScopeExport( + tables, embeddings, found, metadata, + cluster_id=bool(found), regulons=regulons, + ) + + +def claim_attr(owners: Dict[str, Dict[str, str]], group: str, name: str, what: str, + warnings: List[str]) -> bool: + """Give `name` in a loom's flat `group` (col_attrs, row_attrs) to `what`, + unless something already has it: then `what` is left out, with a warning.""" + owner = owners[group].get(name) + if owner is None: + owners[group][name] = what + return True + warnings.append(f"Not written: {what} -- {group}/{name} is {owner}.") + return False + + +def json_safe(value: Any) -> Any: + """`value` with NaN and infinities as null: SCope's client parses MetaData + as strict JSON, which has no token for them.""" + if isinstance(value, float): + return value if np.isfinite(value) else None + if isinstance(value, dict): + return {k: json_safe(v) for k, v in value.items()} + if isinstance(value, (list, tuple)): + return [json_safe(v) for v in value] + return value + + +def is_categorical(member: Any) -> bool: + """A categorical column, in the current layout or anndata 0.7's: codes + whose `categories` attribute points into `__categories`.""" + if is_group(member): + return spec.encoding_type(member) == spec.CATEGORICAL or "codes" in member + return is_dataset(member) and "categories" in member.attrs + + +def categorical_codes(member: Any) -> Any: + return member["codes"] if is_group(member) else member + + +def _distinct_text(member: Any, obs: Any, limit: int) -> Optional[List[str]]: + """A text column's distinct values, read a block at a time; None as soon + as there are more than `limit`, without reading the rest.""" + from adata.conversion.features import STRING_STEP + from adata.elements.read import element_len, read_str_chunk + from adata.util.names import natural_sorted + + seen: Set[str] = set() + cache: Dict[str, np.ndarray] = {} + n = element_len(member) + for lo in range(0, n, STRING_STEP): + seen.update(read_str_chunk(member, lo, min(lo + STRING_STEP, n), cache, parent_group=obs)) + if len(seen) > limit: + return None + return natural_sorted(seen) + + +def _resolve(root: Any, path: str) -> Any: + node = root + for part in path.split("/"): + if part not in node: + return None + node = node[part] + return node diff --git a/src/adata/elements/read.py b/src/adata/elements/read.py index f2c20b5..49f5d62 100644 --- a/src/adata/elements/read.py +++ b/src/adata/elements/read.py @@ -160,31 +160,28 @@ def _categorical_cache_key(col: Any, parent_group: Any | None = None) -> str: return repr(col) -def read_categories(col: Any, parent_group: Any | None = None) -> np.ndarray: - """Read the category labels of a categorical column, modern or legacy.""" +def categories_dataset(col: Any, parent_group: Any | None = None) -> Any: + """The dataset holding a categorical column's labels, modern or legacy, + unread: its length says how many there are without reading them.""" if is_group(col): - return np.asarray(decode_str_array(col["categories"][...]), dtype=str) + return col["categories"] cats_ref = col.attrs.get("categories", None) if cats_ref is not None and is_hdf5_dataset(col): # Legacy 0.7.x HDF5: an object reference into __categories. - return np.asarray(decode_str_array(col.file[cats_ref][...]), dtype=str) + return col.file[cats_ref] if parent_group is not None and "__categories" in parent_group: col_name = str(getattr(col, "name", "")).split("/")[-1] cats_grp = parent_group["__categories"] if col_name in cats_grp: - return np.asarray(decode_str_array(cats_grp[col_name][...]), dtype=str) + return cats_grp[col_name] if cats_ref is not None and isinstance(cats_ref, (str, bytes)): # Legacy Zarr: `categories` is an absolute path within the store. path = spec.decode_attr(cats_ref).lstrip("/") - root = getattr(col, "store_path", None) if parent_group is not None and path.split("/")[-1] in parent_group: - return np.asarray( - decode_str_array(parent_group[path.split("/")[-1]][...]), dtype=str - ) - del root + return parent_group[path.split("/")[-1]] raise KeyError( f"Cannot find categories for categorical column " @@ -192,6 +189,11 @@ def read_categories(col: Any, parent_group: Any | None = None) -> np.ndarray: ) +def read_categories(col: Any, parent_group: Any | None = None) -> np.ndarray: + """Read the category labels of a categorical column, modern or legacy.""" + return np.asarray(decode_str_array(categories_dataset(col, parent_group)[...]), dtype=str) + + def is_ordered(col: Any) -> bool: """Whether a categorical column declares an order over its categories.""" return bool(spec.decode_attr(col.attrs.get("ordered", False))) diff --git a/tests/test_external_readers.py b/tests/test_external_readers.py index 9328e5e..3a0ba12 100644 --- a/tests/test_external_readers.py +++ b/tests/test_external_readers.py @@ -627,6 +627,19 @@ def test_a_scope_loom_imports_as_the_pyscenic_recipe_reads_it(tmp_path, default_ assert problems == [], problems +def test_a_scope_export_reads_with_the_pyscenic_recipe(tmp_path): + """Our --scope loom, read the pySCENIC tutorial's way with loompy, holds + what the AnnData holds: a SCope loom (fixture) to h5ad, to a SCope loom + of ours, and back, compared through the recipe.""" + truth = write_scope_loom(tmp_path / "scope.loom", compressed=False) + first, back, again = tmp_path / "a.h5ad", tmp_path / "back.loom", tmp_path / "b.h5ad" + _convert(str(truth.path), "-o", str(first)) + _convert(str(first), "-o", str(back), "--scope") + _convert(str(back), "-o", str(again)) + problems = _run(["anndata", "loompy", "pandas"], _COMPARE_SCOPE, str(back), str(again)) + assert problems == [], problems + + _COMPARE_SCOPE = r""" import base64, json, sys, warnings, zlib warnings.filterwarnings("ignore") @@ -635,7 +648,10 @@ def test_a_scope_loom_imports_as_the_pyscenic_recipe_reads_it(tmp_path, default_ got = ad.read_h5ad(sys.argv[2]) problems = [] with lp.connect(sys.argv[1], mode="r", validate=False) as lf: - meta = json.loads(zlib.decompress(base64.b64decode(lf.attrs.MetaData))) + try: + meta = json.loads(lf.attrs.MetaData) + except ValueError: # the tutorial's form: base64 of zlib + meta = json.loads(zlib.decompress(base64.b64decode(lf.attrs.MetaData))) auc = pd.DataFrame(lf.ca.RegulonsAUC, index=lf.ca.CellID) regulons = pd.DataFrame(lf.ra.Regulons, index=lf.ra.Gene) drx = pd.DataFrame(lf.ca.Embeddings_X, index=lf.ca.CellID) diff --git a/tests/test_performance.py b/tests/test_performance.py index ca20ea5..e6e76b7 100644 --- a/tests/test_performance.py +++ b/tests/test_performance.py @@ -676,6 +676,25 @@ def measure(n: int) -> int: assert_grows_linearly(measure, what="convert SCope loom to h5ad", axis="n_cells") +def test_scope_loom_export_work_grows_linearly(tmp_path, monkeypatch): + """Compound attributes are written a block of rows at a time, each field + read once per block.""" + from adata.commands.convert import convert_store + from adata.conversion import Options, loom + from tests.conftest import write_scope_loom + + monkeypatch.setattr(loom, "STRING_STEP", 16) + + def measure(n: int) -> int: + h5ad = tmp_path / f"se{n}.h5ad" + convert_store(write_scope_loom(tmp_path / f"se{n}.loom", n_cells=n).path, h5ad, QUIET) + with count_io() as io: + convert_store(h5ad, tmp_path / f"ose{n}.loom", QUIET, options=Options(scope=True)) + return io.work + + assert_grows_linearly(measure, what="convert h5ad to a SCope loom", axis="n_cells") + + def test_scope_loom_import_writes_tables_in_blocks(tmp_path): """A block of rows is written per field, not a row: the write calls do not depend on the cell count.""" diff --git a/tests/test_scope.py b/tests/test_scope.py index 7ac0df9..f89fb03 100644 --- a/tests/test_scope.py +++ b/tests/test_scope.py @@ -37,6 +37,10 @@ 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 = runner.invoke(app, ["convert", *map(str, args)]) assert result.exit_code == 0, _out(result) @@ -311,3 +315,635 @@ def test_a_two_dimensional_compound_attribute_is_left_out_with_a_warning(tmp_pat result = _ok(truth.path, "-o", out) assert "Not read: col_attrs Pairs" in _out(result) assert "Pairs" not in _read(out).obsm + + +# --------------------------------------------------------------------------- +# export: --scope + + +def _scope_loom(path: Path) -> dict: + """A SCope loom's parts, read as SCope's server reads them.""" + with h5py.File(path, "r") as f: + cols, rows = f["col_attrs"], f["row_attrs"] + meta = json.loads(f["attrs/MetaData"][()]) + return { + "meta": meta, + "cols": {k: cols[k][...] for k in cols}, + "rows": {k: rows[k][...] for k in rows}, + } + + +def test_a_scope_loom_round_trips_through_h5ad(tmp_path): + truth = write_scope_loom(tmp_path / "scope.loom") + first, back, again = tmp_path / "a.h5ad", tmp_path / "back.loom", tmp_path / "b.h5ad" + _ok(truth.path, "-o", first) + result = _ok(first, "-o", back, "--scope") + # write_loom's 2-column `embedding` attribute is an embedding too. + assert "loom 3.0.0 (SCope)" in _out(result) and "3 regulons, 3 embeddings, 2 clusterings" in _out(result) + _ok(back, "-o", again) + a, b = _read(first), _read(again) + # SCope's regulon names hold no spaces: "SOX2 (12g)" comes back as "SOX2_(12g)". + unspaced = lambda frame: frame.rename(columns=lambda c: c.replace(" ", "_")) # noqa: E731 + pd.testing.assert_frame_equal(unspaced(a.obsm["RegulonsAUC"]), b.obsm["RegulonsAUC"]) + pd.testing.assert_frame_equal(unspaced(a.varm["Regulons"]), b.varm["Regulons"]) + for key in ("X_tSNE_default", "X_UMAP"): + assert np.array_equal(a.obsm[key], b.obsm[key]) + assert a.obs["Cell Type"].tolist() == b.obs["Cell Type"].tolist() + assert a.obs["Leiden resolution 0.4"].tolist() == b.obs["Leiden resolution 0.4"].tolist() + assert list(b.uns["scope"]["embeddings"]["name"]) == ["tSNE (default)", "UMAP", "embedding"] + thresholds = a.uns["scope"]["regulon_thresholds"].rename(index=lambda r: r.replace(" ", "_")) + pd.testing.assert_frame_equal(thresholds, b.uns["scope"]["regulon_thresholds"]) + + loom = _scope_loom(back) + assert loom["meta"]["embeddings"] == [ + {"id": -1, "name": "tSNE (default)"}, {"id": 0, "name": "UMAP"}, {"id": 1, "name": "embedding"}, + ], "recorded ids and names kept, a new one numbered after them" + assert [c["id"] for c in loom["meta"]["clusterings"]] == [0, 1], "ids are list positions" + assert loom["cols"]["Embeddings_X"].dtype.names == ("-1", "0", "1") + assert set(loom["cols"]["Embedding"].dtype.names) == {"_X", "_Y"} + assert "X_UMAP" not in loom["cols"], "an embedding is not written twice" + + +def test_a_plain_anndata_becomes_a_scope_loom(tmp_path): + rng = np.random.default_rng(0) + n = 30 + obj = ad.AnnData( + X=rng.poisson(1, (n, 6)).astype(np.float32), + obs=pd.DataFrame({"leiden": pd.Categorical([str(i % 3) for i in range(n)]), + "n_counts": rng.random(n), "batch": pd.Categorical(["a", "b"] * 15)}, + index=[f"c{i}" for i in range(n)]), + var=pd.DataFrame(index=[f"G{j}" for j in range(6)]), + ) + obj.obs.loc["c0", "leiden"] = np.nan + obj.obsm["X_pca"] = rng.normal(size=(n, 5)) + obj.obsm["X_umap"] = rng.normal(size=(n, 2)) + obj.obsm["X_draw"] = rng.normal(size=(n, 2)) + obj.obsm["RegulonsAUC"] = pd.DataFrame(rng.random((n, 2)), index=obj.obs_names, columns=["TP53 (5g)", "MYC (9g)"]) + obj.varm["Regulons"] = pd.DataFrame({"TP53 (5g)": [1, 0, 1, 0, 0, 1], "MYC (9g)": [0] * 6}, index=obj.var_names) + source = tmp_path / "plain.h5ad" + obj.write_h5ad(source) + out = tmp_path / "out.loom" + result = _ok(source, "-o", out, "--scope", "--clusterings", "leiden") + assert "'TP53 (5g)' as 'TP53_(5g)'" in _out(result), "SCope's names hold no spaces" + loom = _scope_loom(out) + meta = loom["meta"] + assert meta["embeddings"] == [{"id": -1, "name": "umap"}, {"id": 0, "name": "draw"}] + assert np.array_equal(loom["cols"]["Embedding"]["_X"], obj.obsm["X_umap"][:, 0]) + assert np.array_equal(loom["cols"]["Embeddings_Y"]["0"], obj.obsm["X_draw"][:, 1]) + assert loom["cols"]["X_pca"].shape == (n, 5), "other obsm arrays stay plain" + assert meta["clusterings"] == [{"id": 0, "group": "", "name": "leiden", + "clusters": [{"id": 0, "description": "0"}, {"id": 1, "description": "1"}, + {"id": 2, "description": "2"}]}] + codes = loom["cols"]["Clusterings"]["0"] + assert codes[0] == -1 and codes[1:4].tolist() == [1, 2, 0] + assert loom["cols"]["ClusterID"].tolist() == codes.tolist() + assert loom["cols"]["RegulonsAUC"].dtype.names == ("TP53_(5g)", "MYC_(9g)") + assert loom["rows"]["Regulons"].dtype.names == ("TP53_(5g)", "MYC_(9g)") + assert {m["name"] for m in meta["metrics"]} == {"n_counts"} + assert {a["name"]: a["values"] for a in meta["annotations"]}["batch"] == ["a", "b"] + + result = _convert(source, "-o", tmp_path / "x.loom", "--scope", "--default-embedding", "X_draw") + assert result.exit_code == 0 + assert _scope_loom(tmp_path / "x.loom")["meta"]["embeddings"][0] == {"id": -1, "name": "draw"} + + +@pytest.mark.parametrize( + "extra,message", + [ + (["--scope", "--clusterings", "n_counts"], "--clusterings n_counts: not a categorical obs column"), + (["--scope", "--clusterings", "nope"], "--clusterings nope: no such obs column"), + (["--scope", "--default-embedding", "X_pca"], "--default-embedding X_pca: not an obsm array of two columns"), + (["--clusterings", "leiden"], "--clusterings shapes a SCope loom: pass --scope too"), + ], +) +def test_scope_export_refusals_leave_nothing(tmp_path, extra, message): + obj = ad.AnnData( + X=np.ones((4, 2), dtype=np.float32), + obs=pd.DataFrame({"leiden": pd.Categorical(["a", "b", "a", "b"]), "n_counts": np.arange(4.0)}, + index=list("wxyz")), + var=pd.DataFrame(index=["g1", "g2"]), + ) + obj.obsm["X_pca"] = np.zeros((4, 3)) + source = tmp_path / "s.h5ad" + obj.write_h5ad(source) + out = tmp_path / "out.loom" + result = _convert(source, "-o", out, *extra) + assert result.exit_code == 1 and message in _out(result) + assert not out.exists() + + +def test_too_many_regulons_for_hdf5_are_refused_before_writing(tmp_path): + n = 4 + names = [f"REGULON_NUMBER_{i:05d}_(+)" for i in range(2500)] + obj = ad.AnnData(X=np.ones((n, 2), dtype=np.float32), obs=pd.DataFrame(index=list("wxyz")), + var=pd.DataFrame(index=["g1", "g2"])) + obj.obsm["RegulonsAUC"] = pd.DataFrame(np.zeros((n, len(names))), index=obj.obs_names, columns=names) + source = tmp_path / "s.h5ad" + obj.write_h5ad(source) + out = tmp_path / "out.loom" + result = _convert(source, "-o", out, "--scope") + assert result.exit_code == 1 + assert "obsm/RegulonsAUC has 2,500 columns, more than HDF5 can hold" in _out(result) + assert not out.exists() + + +def test_cluster_markers_need_their_metrics(tmp_path): + """SCope crashes on ClusterMarkers_ when its clustering lists no + clusterMarkerMetrics: such markers are left out, those with their + _ tables get metrics named after them.""" + truth = write_scope_loom(tmp_path / "scope.loom") + with h5py.File(truth.path, "a") as f: + rows = f["row_attrs"] + marks = np.zeros(10, dtype=[("0", "i8"), ("1", "i8"), ("2", "i8")]) + marks["0"][:3] = 1 + rows.create_dataset("ClusterMarkers_0", data=marks) + rows.create_dataset("ClusterMarkers_0_avg_logFC", data=np.zeros(10, dtype=[("0", "f4"), ("1", "f4"), ("2", "f4")])) + rows.create_dataset("ClusterMarkers_1", data=np.zeros(10, dtype=[("0", "i8"), ("1", "i8")])) + first, back = tmp_path / "a.h5ad", tmp_path / "back.loom" + _ok(truth.path, "-o", first) + result = _ok(first, "-o", back, "--scope") + assert "Not written: ClusterMarkers_1 -- without metrics" in _out(result) + loom = _scope_loom(back) + assert "ClusterMarkers_0" in loom["rows"] and "ClusterMarkers_0_avg_logFC" in loom["rows"] + assert "ClusterMarkers_1" not in loom["rows"] + assert loom["meta"]["clusterings"][0]["clusterMarkerMetrics"] == [ + {"accessor": "avg_logFC", "name": "avg_logFC", "description": ""} + ] + + +def test_a_plain_export_of_a_scope_anndata_points_to_scope(tmp_path): + truth = write_scope_loom(tmp_path / "scope.loom") + first = tmp_path / "a.h5ad" + _ok(truth.path, "-o", first) + result = _ok(first, "-o", tmp_path / "plain.loom") + assert "This came from a SCope loom: pass --scope" in _out(result) + with h5py.File(tmp_path / "plain.loom", "r") as f: + assert "RegulonsAUC" not in f["col_attrs"] + + +def test_scope_export_writes_in_small_blocks_the_same(tmp_path, monkeypatch): + from adata.conversion import loom + + monkeypatch.setattr(loom, "STRING_STEP", 7) + truth = write_scope_loom(tmp_path / "scope.loom", n_cells=50) + first, back = tmp_path / "a.h5ad", tmp_path / "back.loom" + _ok(truth.path, "-o", first) + _ok(first, "-o", back, "--scope") + loom_parts = _scope_loom(back) + got = np.column_stack([loom_parts["cols"]["RegulonsAUC"][r.replace(" ", "_")] for r in truth.regulons]) + assert np.array_equal(got, truth.auc) + assert np.array_equal(loom_parts["cols"]["Embeddings_X"]["0"], truth.embeddings["0"][:, 0]) + assert loom_parts["cols"]["Clusterings"]["1"].tolist() == truth.clusterings["1"].tolist() + + +def test_metadata_from_columns_and_a_thresholds_table(tmp_path): + """Text and nullable obs columns are annotations and metrics too; with no + recorded MetaData, regulon thresholds come from uns["scope"]'s table.""" + n = 6 + obj = ad.AnnData( + X=np.ones((n, 2), dtype=np.float32), + obs=pd.DataFrame({"donor": ["d1", "d2", "d3", "d4", "d5", "d6"], + "n_genes": pd.array([1, None, 3, 4, 5, 6], dtype="Int64")}, + index=[f"c{i}" for i in range(n)]), + var=pd.DataFrame(index=["g1", "g2"]), + ) + obj.obsm["RegulonsAUC"] = pd.DataFrame(np.zeros((n, 2)), index=obj.obs_names, columns=["A_(+)", "B (2g)"]) + obj.uns["scope"] = {"regulon_thresholds": pd.DataFrame( + {"defaultThresholdValue": [0.5, 0.25], "motifData": ["a.png", "b.png"]}, index=["A_(+)", "B (2g)"])} + source = tmp_path / "s.h5ad" + obj.write_h5ad(source) + out = tmp_path / "out.loom" + _ok(source, "-o", out, "--scope") + meta = _scope_loom(out)["meta"] + assert {"name": "donor", "values": ["d1", "d2", "d3", "d4", "d5", "d6"]} in meta["annotations"] + assert {"name": "n_genes"} in meta["metrics"] + assert meta["regulonThresholds"] == [ + {"regulon": "A_(+)", "defaultThresholdValue": 0.5, "motifData": "a.png"}, + {"regulon": "B_(2g)", "defaultThresholdValue": 0.25, "motifData": "b.png"}, + ] + + +def _small_anndata(tmp_path: Path, n: int = 12, **obs): + rng = np.random.default_rng(1) + frame = {"leiden": pd.Categorical([str(i % 3) for i in range(n)]), **obs} + obj = ad.AnnData( + X=rng.poisson(1, (n, 4)).astype(np.float32), + obs=pd.DataFrame(frame, index=[f"c{i}" for i in range(n)]), + var=pd.DataFrame(index=[f"G{j}" for j in range(4)]), + ) + obj.obsm["X_umap"] = rng.normal(size=(n, 2)) + return obj + + +def test_regulons_that_meet_once_spaces_go_keep_one_field_in_both_tables(tmp_path): + """`A B` and `A_B` were renamed table by table: with the columns in + another order in Regulons, one field held one regulon's AUC and the + other's targets.""" + obj = _small_anndata(tmp_path) + obj.obsm["RegulonsAUC"] = pd.DataFrame({"A B": np.full(12, 0.25), "A_B": np.full(12, 0.75)}, + index=obj.obs_names) + obj.varm["Regulons"] = pd.DataFrame({"A_B": [0, 0, 0, 1], "A B": [1, 1, 1, 0]}, index=obj.var_names) + source, out = tmp_path / "in.h5ad", tmp_path / "out.loom" + obj.write_h5ad(source) + _ok(source, "-o", out, "--scope") + loom = _scope_loom(out) + auc, targets = loom["cols"]["RegulonsAUC"], loom["rows"]["Regulons"] + assert sorted(auc.dtype.names) == sorted(targets.dtype.names) + for field in auc.dtype.names: + spaced = auc[field][0] == 0.25 + assert list(targets[field]) == ([1, 1, 1, 0] if spaced else [0, 0, 0, 1]), field + + +def test_cluster_id_is_written_from_the_first_clustering_not_a_stale_column(tmp_path): + """An obs column named ClusterID was written as it was, and SCope's own + ClusterID left out: after the clustering was edited the two disagreed.""" + obj = _small_anndata(tmp_path, ClusterID=np.full(12, 7)) + source, out = tmp_path / "in.h5ad", tmp_path / "out.loom" + obj.write_h5ad(source) + result = _ok(source, "-o", out, "--scope", "--clusterings", "leiden") + assert "'ClusterID' is written afresh from the first clustering, 'leiden'" in _out(result) + loom = _scope_loom(out) + assert list(loom["cols"]["ClusterID"]) == [i % 3 for i in range(12)] + assert "ClusterID" not in [m["name"] for m in loom["meta"]["metrics"]] + + +def test_a_text_column_with_too_many_values_is_no_annotation_and_is_not_read_whole(tmp_path, monkeypatch): + """Annotation values were found by reading the whole column, however many + values it had, to throw them away past 245.""" + import adata.elements.read as read + from adata.conversion import features, scope + + n = 4000 + obj = _small_anndata(tmp_path, n=n, barcode=[f"b{i}" for i in range(n)], + sample=["s1", "s2"] * (n // 2)) + source, out = tmp_path / "in.h5ad", tmp_path / "out.loom" + obj.write_h5ad(source) + _ok(source, "-o", out, "--scope") + names = [a["name"] for a in _scope_loom(out)["meta"]["annotations"]] + assert "barcode" not in names and "sample" in names + + monkeypatch.setattr(features, "STRING_STEP", 10) + real, rows = read.read_str_chunk, [] + + def counted(obj, lo, hi, *args, **kwargs): + rows.append(hi - lo) + return real(obj, lo, hi, *args, **kwargs) + + monkeypatch.setattr(read, "read_str_chunk", counted) + with h5py.File(source) as f: + assert scope._distinct_text(f["obs/barcode"], f["obs"], scope.MAX_ANNOTATION_VALUES) is None + assert sum(rows) == 250, "stops at the first block past 245 values" + assert scope._distinct_text(f["obs/sample"], f["obs"], scope.MAX_ANNOTATION_VALUES) == ["s1", "s2"] + + +def test_marker_fields_are_renumbered_with_their_clusters(tmp_path): + """Clusters numbered 2, 5 and 9 were renumbered 0, 1 and 2 in MetaData, + but the marker tables kept fields 2, 5 and 9: SCope matched no marker + to any cluster.""" + truth = write_scope_loom(tmp_path / "scope.loom", compressed=False) + meta = json.loads(json.dumps(truth.metadata)) + ids = {0: 2, 1: 5, 2: 9} + for cluster in meta["clusterings"][0]["clusters"]: + cluster["id"] = ids[cluster["id"]] + meta["clusterings"][0]["clusterMarkerMetrics"] = [{"accessor": "avg_logFC", "name": "logFC", "description": ""}] + _set_metadata(truth.path, json.dumps(meta)) + with h5py.File(truth.path, "a") as f: + cols, rows = f["col_attrs"], f["row_attrs"] + clusterings = cols["Clusterings"][...] + clusterings["0"] = np.vectorize(ids.get)(clusterings["0"]) + del cols["Clusterings"] + cols.create_dataset("Clusterings", data=clusterings) + marks = np.zeros(10, dtype=[("2", "i8"), ("5", "i8"), ("9", "i8")]) + logfc = np.zeros(10, dtype=[("2", "f4"), ("5", "f4"), ("9", "f4")]) + for k, field in enumerate(("2", "5", "9")): + marks[field][k] = 1 + logfc[field][k] = k + 0.5 + rows.create_dataset("ClusterMarkers_0", data=marks) + rows.create_dataset("ClusterMarkers_0_avg_logFC", data=logfc) + first, back = tmp_path / "a.h5ad", tmp_path / "back.loom" + _ok(truth.path, "-o", first) + _ok(first, "-o", back, "--scope") + loom = _scope_loom(back) + clusters = loom["meta"]["clusterings"][0]["clusters"] + assert [c["description"] for c in clusters] == ["B cell", "T cell", "NK"] + marks, logfc = loom["rows"]["ClusterMarkers_0"], loom["rows"]["ClusterMarkers_0_avg_logFC"] + assert marks.dtype.names == logfc.dtype.names == ("0", "1", "2") + for k in range(3): + assert list(np.flatnonzero(marks[str(k)])) == [k] + assert logfc[str(k)][k] == k + 0.5 + + +def test_a_clustering_named_with_clusterings_keeps_its_marker_tables(tmp_path): + """Marker tables were kept only for a clustering an import recorded, so + a plain AnnData's lost theirs.""" + obj = _small_anndata(tmp_path) + marks = pd.DataFrame({str(k): [int(g == k) for g in range(4)] for k in range(3)}, index=obj.var_names) + obj.varm["ClusterMarkers_0"] = marks + obj.varm["ClusterMarkers_0_avg_logFC"] = marks.astype(np.float32) * 2 + source, out = tmp_path / "in.h5ad", tmp_path / "out.loom" + obj.write_h5ad(source) + _ok(source, "-o", out, "--scope", "--clusterings", "leiden") + loom = _scope_loom(out) + assert loom["rows"]["ClusterMarkers_0"].dtype.names == ("0", "1", "2") + assert "ClusterMarkers_0_avg_logFC" in loom["rows"] + assert loom["meta"]["clusterings"][0]["clusterMarkerMetrics"][0]["accessor"] == "avg_logFC" + + +def test_a_categorical_with_too_many_categories_is_no_annotation_and_its_categories_are_not_read( + tmp_path, monkeypatch +): + import adata.elements.read as read + + n = 300 + obj = _small_anndata(tmp_path, n=n, barcode=pd.Categorical([f"b{i}" for i in range(n)])) + source, out = tmp_path / "in.h5ad", tmp_path / "out.loom" + obj.write_h5ad(source) + real, seen = read.read_categories, [] + + def counted(column, *args, **kwargs): + seen.append(column.name) + return real(column, *args, **kwargs) + + monkeypatch.setattr(read, "read_categories", counted) + from adata.conversion import scope + + with h5py.File(source) as f: + found = scope.plan_export(f, obs_columns=["barcode", "leiden"], varm_key="varm", clusterings=None, + default_embedding=None, warnings=[]) + assert [a["name"] for a in found.metadata["annotations"]] == ["leiden"] + assert "/obs/barcode" not in seen + + +def test_markers_of_undescribed_clusters_are_kept_on_a_round_trip(tmp_path): + """With no clusters in MetaData the import labels them by id, and the + export, finding no descriptions to match, dropped every marker field.""" + truth = write_scope_loom(tmp_path / "scope.loom", compressed=False) + meta = json.loads(json.dumps(truth.metadata)) + meta["clusterings"][0]["clusters"] = [] + meta["clusterings"][0]["clusterMarkerMetrics"] = [{"accessor": "avg_logFC", "name": "logFC", "description": ""}] + _set_metadata(truth.path, json.dumps(meta)) + with h5py.File(truth.path, "a") as f: + logfc = np.zeros(10, dtype=[("0", "f4"), ("1", "f4"), ("2", "f4")]) + marks = np.zeros(10, dtype=[("0", "i8"), ("1", "i8"), ("2", "i8")]) + logfc["1"][4], marks["1"][4] = 2.5, 1 + f["row_attrs"].create_dataset("ClusterMarkers_0", data=marks) + f["row_attrs"].create_dataset("ClusterMarkers_0_avg_logFC", data=logfc) + first, back = tmp_path / "a.h5ad", tmp_path / "back.loom" + _ok(truth.path, "-o", first) + _ok(first, "-o", back, "--scope") + loom = _scope_loom(back) + assert loom["rows"]["ClusterMarkers_0_avg_logFC"].dtype.names == ("0", "1", "2") + assert loom["rows"]["ClusterMarkers_0_avg_logFC"]["1"][4] == 2.5 + + +@pytest.mark.parametrize("scope_flag", [True, False]) +def test_two_elements_asking_for_one_attribute_name_keep_the_first_and_warn(tmp_path, scope_flag): + """An obs column and an obsm entry of one name were both written to + col_attrs: the second failed with the output half written.""" + obj = _small_anndata(tmp_path, scores=np.arange(12.0), Clusterings=pd.Categorical(["x"] * 12)) + obj.obsm["scores"] = pd.DataFrame({"A": np.ones(12)}, index=obj.obs_names) if scope_flag else np.ones((12, 3)) + obj.varm["Gene"] = np.ones((4, 3)) + source, out = tmp_path / "in.h5ad", tmp_path / "out.loom" + obj.write_h5ad(source) + args = ["--scope", "--clusterings", "leiden"] if scope_flag else [] + result = _ok(source, "-o", out, *args) + text = _out(result) + assert "Not written: obsm/scores -- col_attrs/scores is obs column 'scores'." in text + assert "Not written: varm/Gene -- row_attrs/Gene is the gene names." in text + with h5py.File(out) as f: + assert list(f["col_attrs/scores"][...]) == list(np.arange(12.0)) + if scope_flag: + assert "Not written: obs column 'Clusterings' -- col_attrs/Clusterings is SCope's clusterings." in text + assert f["col_attrs/Clusterings"].dtype.names == ("0",) + + +def test_an_edited_thresholds_table_wins_over_the_recorded_metadata(tmp_path): + """The thresholds table was read only when the recorded MetaData had + none, so editing it changed nothing in the export.""" + truth = write_scope_loom(tmp_path / "scope.loom") + first, back = tmp_path / "a.h5ad", tmp_path / "back.loom" + _ok(truth.path, "-o", first) + with h5py.File(first, "a") as f: + f["uns/scope/regulon_thresholds/defaultThresholdValue"][0] = 0.9 + _ok(first, "-o", back, "--scope") + thresholds = _scope_loom(back)["meta"]["regulonThresholds"] + assert thresholds[0]["defaultThresholdValue"] == 0.9 + assert [t["defaultThresholdValue"] for t in thresholds[1:]] == pytest.approx([0.2, 0.3]) + assert thresholds[0]["allThresholds"] == truth.metadata["regulonThresholds"][0]["allThresholds"] + + +def test_marker_metrics_list_the_tables_written(tmp_path): + """clusterMarkerMetrics was copied as recorded, listing a metric whose + table had gone and missing one added since.""" + truth = write_scope_loom(tmp_path / "scope.loom", compressed=False) + meta = json.loads(json.dumps(truth.metadata)) + meta["clusterings"][0]["clusterMarkerMetrics"] = [ + {"accessor": "avg_logFC", "name": "logFC", "description": "log fold change"}, + {"accessor": "pval", "name": "p", "description": ""}, + ] + _set_metadata(truth.path, json.dumps(meta)) + fields = [("0", "f4"), ("1", "f4"), ("2", "f4")] + with h5py.File(truth.path, "a") as f: + rows = f["row_attrs"] + rows.create_dataset("ClusterMarkers_0", data=np.zeros(10, dtype=[(n, "i8") for n, _ in fields])) + rows.create_dataset("ClusterMarkers_0_avg_logFC", data=np.zeros(10, dtype=fields)) + rows.create_dataset("ClusterMarkers_0_score", data=np.zeros(10, dtype=fields)) + first, back = tmp_path / "a.h5ad", tmp_path / "back.loom" + _ok(truth.path, "-o", first) + _ok(first, "-o", back, "--scope") + assert _scope_loom(back)["meta"]["clusterings"][0]["clusterMarkerMetrics"] == [ + {"accessor": "avg_logFC", "name": "logFC", "description": "log fold change"}, + {"accessor": "score", "name": "score", "description": ""}, + ] + + +def test_a_clustering_is_written_once_and_comes_back_under_its_own_name(tmp_path): + """A clustering went out as Clusterings and as a plain attribute too; the + loom read back gave the plain one the name, the clustering `leiden-1`, + and each round trip added another.""" + obj = _small_anndata(tmp_path) + source, out, back = tmp_path / "in.h5ad", tmp_path / "out.loom", tmp_path / "back.h5ad" + obj.write_h5ad(source) + _ok(source, "-o", out, "--scope", "--clusterings", "leiden") + loom = _scope_loom(out) + assert "leiden" not in loom["cols"] + assert "leiden" not in [a["name"] for a in loom["meta"]["annotations"]] + _ok(out, "-o", back) + obs = _read(back).obs + assert [c for c in obs.columns if c.startswith("leiden")] == ["leiden"] + assert list(obs["leiden"].astype(str)) == list(obj.obs["leiden"].astype(str)) + + +def test_markers_whose_table_lost_its_name_are_left_out_with_their_metrics(tmp_path): + """A marker table that lost its name to a var column was left out, but + clusterMarkerMetrics still listed it.""" + obj = _small_anndata(tmp_path) + marks = pd.DataFrame({str(k): [int(g == k) for g in range(4)] for k in range(3)}, index=obj.var_names) + obj.varm["ClusterMarkers_0"] = marks + obj.varm["ClusterMarkers_0_avg_logFC"] = marks.astype(np.float32) + obj.var["ClusterMarkers_0_avg_logFC"] = np.arange(4.0) + source, out = tmp_path / "in.h5ad", tmp_path / "out.loom" + obj.write_h5ad(source) + result = _ok(source, "-o", out, "--scope", "--clusterings", "leiden") + text = _out(result) + assert ("Not written: varm/ClusterMarkers_0_avg_logFC -- row_attrs/ClusterMarkers_0_avg_logFC " + "is var column 'ClusterMarkers_0_avg_logFC'.") in text + assert "Not written: ClusterMarkers_0 -- without metrics" in text + loom = _scope_loom(out) + assert "ClusterMarkers_0" not in loom["rows"] + assert loom["rows"]["ClusterMarkers_0_avg_logFC"].dtype.names is None, "the var column" + assert "clusterMarkerMetrics" not in loom["meta"]["clusterings"][0] + + +def test_an_anndata_07_categorical_can_be_a_clustering(tmp_path, sample_legacy_v010_h5ad): + """--clusterings took only the current categorical layout, and refused + anndata 0.7's codes with a `categories` reference.""" + out = tmp_path / "out.loom" + _ok(sample_legacy_v010_h5ad, "-o", out, "--scope", "--clusterings", "cell_type") + loom = _scope_loom(out) + assert list(loom["cols"]["Clusterings"]["0"]) == [0, 1, 0, 2] + assert [c["description"] for c in loom["meta"]["clusterings"][0]["clusters"]] == ["TypeA", "TypeB", "TypeC"] + + +def test_an_anndata_07_categorical_with_too_many_categories_is_not_read(tmp_path, sample_legacy_v010_h5ad): + """The count of a legacy categorical's categories was checked only for + the current layout: anndata 0.7's were read whole to be thrown away.""" + import adata.elements.read as read + from adata.conversion import scope + + with h5py.File(sample_legacy_v010_h5ad, "a") as f: + obs = f["obs"] + cats = obs["__categories"].create_dataset("barcode", data=np.array([f"b{i}" for i in range(300)], dtype="S")) + codes = obs.create_dataset("barcode", data=np.array([0, 1, 2, 3], dtype=np.int16)) + codes.attrs["categories"] = cats.ref + real, seen = read.read_categories, [] + + def counted(column, *args, **kwargs): + seen.append(column.name) + return real(column, *args, **kwargs) + + with pytest.MonkeyPatch.context() as patch: + patch.setattr(read, "read_categories", counted) + with h5py.File(sample_legacy_v010_h5ad) as f: + found = scope.plan_export(f, obs_columns=["cell_type", "barcode"], varm_key="varm", + clusterings=None, default_embedding=None, warnings=[]) + assert [a["name"] for a in found.metadata["annotations"]] == ["cell_type"] + assert "/obs/barcode" not in seen + + +def test_a_missing_threshold_is_left_out_of_strict_json(tmp_path): + """A missing defaultThresholdValue went into MetaData as NaN, which no + strict JSON reader (SCope's client) parses.""" + truth = write_scope_loom(tmp_path / "scope.loom") + first, back = tmp_path / "a.h5ad", tmp_path / "back.loom" + _ok(truth.path, "-o", first) + with h5py.File(first, "a") as f: + f["uns/scope/regulon_thresholds/defaultThresholdValue"][0] = np.nan + _ok(first, "-o", back, "--scope") + + def refuse(token): + raise ValueError(f"not JSON: {token}") + + with h5py.File(back) as f: + meta = json.loads(f["attrs/MetaData"][()], parse_constant=refuse) + assert "defaultThresholdValue" not in meta["regulonThresholds"][0] + assert meta["regulonThresholds"][1]["defaultThresholdValue"] == pytest.approx(0.2) + + +def test_an_auc_table_that_loses_its_name_takes_its_regulon_metadata_with_it(tmp_path): + """An obs column named RegulonsAUC kept the name, and the AUC table was + left out, but its regulon thresholds stayed in MetaData.""" + truth = write_scope_loom(tmp_path / "scope.loom") + first, edited, back = tmp_path / "a.h5ad", tmp_path / "b.h5ad", tmp_path / "back.loom" + _ok(truth.path, "-o", first) + obj = _read(first) + obj.obs["RegulonsAUC"] = np.arange(obj.n_obs, dtype=float) + obj.write_h5ad(edited) + result = _ok(edited, "-o", back, "--scope") + assert "Not written: obsm/RegulonsAUC -- col_attrs/RegulonsAUC is obs column 'RegulonsAUC'." in _out(result) + loom = _scope_loom(back) + assert loom["cols"]["RegulonsAUC"].dtype.names is None + assert loom["meta"]["regulonThresholds"] == [] + + +def test_a_marker_table_left_with_no_cluster_is_not_listed_as_a_metric(tmp_path): + """A metric table whose fields named no cluster any more was dropped + after clusterMarkerMetrics had listed it.""" + truth = write_scope_loom(tmp_path / "scope.loom", compressed=False) + meta = json.loads(json.dumps(truth.metadata)) + meta["clusterings"][0]["clusterMarkerMetrics"] = [ + {"accessor": "avg_logFC", "name": "logFC", "description": ""}, + {"accessor": "pval", "name": "p", "description": ""}, + ] + _set_metadata(truth.path, json.dumps(meta)) + fields = [("0", "f4"), ("1", "f4"), ("2", "f4")] + with h5py.File(truth.path, "a") as f: + rows = f["row_attrs"] + rows.create_dataset("ClusterMarkers_0", data=np.zeros(10, dtype=[(n, "i8") for n, _ in fields])) + rows.create_dataset("ClusterMarkers_0_avg_logFC", data=np.zeros(10, dtype=fields)) + rows.create_dataset("ClusterMarkers_0_pval", data=np.zeros(10, dtype=[("7", "f4")])) + first, back = tmp_path / "a.h5ad", tmp_path / "back.loom" + _ok(truth.path, "-o", first) + result = _ok(first, "-o", back, "--scope") + assert "Not written in varm/ClusterMarkers_0_pval: 7" in _out(result) + loom = _scope_loom(back) + assert "ClusterMarkers_0_pval" not in loom["rows"] + assert [m["accessor"] for m in loom["meta"]["clusterings"][0]["clusterMarkerMetrics"]] == ["avg_logFC"] + + +def test_a_clustering_leaves_its_name_to_an_obsm_entry(tmp_path): + """A clustering column claimed its attribute name and then went out only + inside Clusterings, so an obsm entry of that name was dropped for + nothing.""" + obj = _small_anndata(tmp_path) + obj.obsm["leiden"] = pd.DataFrame({"A": np.full(12, 0.5)}, index=obj.obs_names) + source, out = tmp_path / "in.h5ad", tmp_path / "out.loom" + obj.write_h5ad(source) + result = _ok(source, "-o", out, "--scope", "--clusterings", "leiden") + assert "Not written: obsm/leiden" not in _out(result) + loom = _scope_loom(out) + assert loom["cols"]["leiden"].dtype.names == ("A",) + assert list(loom["cols"]["Clusterings"]["0"]) == [i % 3 for i in range(12)] + + +def test_metric_tables_without_their_markers_table_are_left_out(tmp_path): + """ClusterMarkers_0_avg_logFC with no ClusterMarkers_0 was written, and + listed as a metric, though SCope finds markers through the latter.""" + obj = _small_anndata(tmp_path) + obj.varm["ClusterMarkers_0_avg_logFC"] = pd.DataFrame( + {str(k): np.ones(4, dtype=np.float32) for k in range(3)}, index=obj.var_names) + source, out = tmp_path / "in.h5ad", tmp_path / "out.loom" + obj.write_h5ad(source) + result = _ok(source, "-o", out, "--scope", "--clusterings", "leiden") + assert "Not written: ClusterMarkers_0" in _out(result) + loom = _scope_loom(out) + assert "ClusterMarkers_0_avg_logFC" not in loom["rows"] + assert "clusterMarkerMetrics" not in loom["meta"]["clusterings"][0] + + +def test_an_obsm_table_named_like_markers_claims_its_name(tmp_path): + """An obsm dataframe named ClusterMarkers_0 skipped the name check meant + for varm's marker tables, and was written over an obs column of that + name with the output half written.""" + obj = _small_anndata(tmp_path, ClusterMarkers_0=np.arange(12.0)) + obj.obsm["ClusterMarkers_0"] = pd.DataFrame({"A": np.ones(12)}, index=obj.obs_names) + source, out = tmp_path / "in.h5ad", tmp_path / "out.loom" + obj.write_h5ad(source) + result = _ok(source, "-o", out, "--scope") + assert "Not written: obsm/ClusterMarkers_0 -- col_attrs/ClusterMarkers_0 is obs column" in _out(result) + assert list(_scope_loom(out)["cols"]["ClusterMarkers_0"]) == list(np.arange(12.0)) + + +def test_too_many_embeddings_for_hdf5_are_refused_before_writing(tmp_path, monkeypatch): + from adata.conversion import scope + + obj = _small_anndata(tmp_path) + obj.obsm["X_tsne"] = np.zeros((12, 2)) + source, out = tmp_path / "in.h5ad", tmp_path / "out.loom" + obj.write_h5ad(source) + # Small enough to test, the same check as at 64 KiB. + monkeypatch.setattr(scope, "compound_fits", lambda dtype: len(dtype.names) < 2) + result = _convert(source, "-o", out, "--scope") + assert result.exit_code == 1 + assert "2 embeddings (two-column obsm arrays) are more than HDF5 can hold" in _out(result) + assert not out.exists()