From 4642e2a3b65354976154d45476e5bd94627edd26 Mon Sep 17 00:00:00 2001 From: Aljes Date: Wed, 21 Oct 2026 20:18:00 +0100 Subject: [PATCH 01/10] convert: write SCope looms with --scope (#61) A loom export wrote a plain loom, so results edited as AnnData could not go back to the SCope viewer. With --scope: - obsm and varm dataframes (RegulonsAUC, Regulons, ClusterMarkers_) become compound attributes of their numeric columns, a block of rows at a time; regulon names lose their spaces, alike in AUC and Regulons, as SCope needs; - two-column obsm arrays become Embeddings_X/_Y (and Embedding for the default, id -1), not again plain attributes; - the categorical columns a SCope import recorded, and --clusterings, become Clusterings and ClusterID, numbered 0..n-1 with their clusters described, as SCope's server indexes them; - MetaData (plain JSON): embeddings, clusterings (markers only with their metrics, else left out: SCope crashes on them), metrics, annotations, regulon thresholds; recorded ids and names are kept, so a SCope loom round-trips through h5ad. A non-categorical clustering, an unknown --default-embedding, or a table too wide for one HDF5 compound type (64 KiB) is refused before writing. --clusterings and --default-embedding need --scope; a plain export of a SCope-born AnnData says --scope would keep its regulons. Checked on the Zenodo pySCENIC loom: loom -> h5ad -> --scope loom equals the original through the pySCENIC recipe, and meets SCope's reading rules. Co-Authored-By: Claude Opus 5.5 --- docs/COMMANDS.md | 47 ++++- docs/TESTING.md | 2 +- src/adata/cli.py | 27 ++- src/adata/conversion/base.py | 5 + src/adata/conversion/loom.py | 106 ++++++++++- src/adata/conversion/scope.py | 328 +++++++++++++++++++++++++++++++++ tests/test_external_readers.py | 18 +- tests/test_performance.py | 19 ++ tests/test_scope.py | 210 +++++++++++++++++++++ 9 files changed, 748 insertions(+), 14 deletions(-) diff --git a/docs/COMMANDS.md b/docs/COMMANDS.md index 0b3310e..fa72c90 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,49 @@ 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, 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, 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 with +their metrics (`ClusterMarkers__avg_logFC` and the like, or those a +SCope import recorded). `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. + +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..3fb5950 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,34 @@ 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) + 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, + ) + # 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." + ) 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 +843,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 +894,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(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(obs[col]["codes"][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: @@ -946,13 +1034,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..f240541 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 @@ -395,3 +396,330 @@ 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 plan_export( + root: Any, + *, + obs_columns: Sequence[str], + varm_key: str, + clusterings: Optional[str], + default_embedding: Optional[str], + warnings: List[str], +) -> ScopeExport: + """Everything a SCope loom needs, worked out and checked before writing.""" + from adata.elements.read import dataframe_columns, read_categories, read_str_all, resolve_index + from adata.util.names import natural_sorted + + parsed, recorded_embeddings, recorded_clusterings = _recorded(root) + + # Compound tables: every dataframe in obsm and varm. + tables: List[ExportTable] = [] + 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.") + if not numeric: + continue + regulons = key.endswith("Regulons") or key.endswith("RegulonsAUC") + taken: Set[str] = set() + fields = [claim(c.replace(" ", "_") if regulons else 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) + + # Clusterings: the recorded ones, and those --clusterings names. + obs = root["obs"] + wanted = [c for c in recorded_clusterings if c in obs] + for column in (clusterings or "").split(","): + column = column.strip() + if column and column not in wanted: + wanted.append(column) + entries = {str(e.get("id")): e for e in parsed.get("clusterings", []) if isinstance(e, dict)} + found: List[ExportClustering] = [] + renumbered: Dict[int, int] = {} + for position, column in enumerate(wanted): + if column not in obs: + raise ValueError(f"--clusterings {column}: no such obs column.") + if spec.encoding_type(obs[column]) != spec.CATEGORICAL: + 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)] + entry = { + "id": position, + "group": recorded.get("group", ""), + "name": name or column, + "clusters": [{"id": i, "description": c} for i, c in enumerate(categories)], + } + if recorded.get("clusterMarkerMetrics"): + entry["clusterMarkerMetrics"] = recorded["clusterMarkerMetrics"] + found.append(ExportClustering(position, column, entry)) + + # 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 + accessors = [m.group(2) for t in group_tables if (m := _MARKERS.match(t.key)).group(2)] + if not entry.get("clusterMarkerMetrics"): + if not accessors: + warnings.append( + f"Not written: ClusterMarkers_{old} -- without metrics (_avg_logFC, " + "_pval, ...) SCope cannot show markers." + ) + continue + entry["clusterMarkerMetrics"] = [ + {"accessor": a, "name": a, "description": ""} for a in accessors + ] + for table in group_tables: + accessor = _MARKERS.match(table.key).group(2) + table.key = f"ClusterMarkers_{new}" + (f"_{accessor}" if accessor else "") + keep.append(table) + tables = keep + + # Metrics and annotations: the obs columns written. + metrics, annotations = [], [] + for column in obs_columns: + member = obs[column] + enc = spec.encoding_type(member) + if enc == spec.CATEGORICAL: + 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 = natural_sorted(set(read_str_all(member, parent_group=obs))) + 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} + thresholds = [] + for entry in parsed.get("regulonThresholds", []) or []: + if isinstance(entry, dict) and str(entry.get("regulon", "")).replace(" ", "_") in auc_fields: + thresholds.append(dict(entry, regulon=str(entry["regulon"]).replace(" ", "_"))) + if not thresholds and "uns" in root and "scope" in root["uns"] and "regulon_thresholds" in root["uns"]["scope"]: + for regulon, row in _frame_rows(root["uns"]["scope"]["regulon_thresholds"]).items(): + name = regulon.replace(" ", "_") + if name in auc_fields: + entry = {"regulon": name} + for key, value in row.items(): + entry[key] = _float(value) if key == "defaultThresholdValue" else value + thresholds.append(entry) + + 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) and "ClusterID" not in obs_columns, regulons=regulons, + ) + + +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/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..4d63998 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,209 @@ 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"}, + ] From fe330212d52ad8b7f8649b1d7382b942ee0ba8ff Mon Sep 17 00:00:00 2001 From: Aljes Date: Thu, 22 Oct 2026 15:51:22 +0100 Subject: [PATCH 02/10] convert --scope: one field per regulon, ClusterID from the clustering, bounded annotation scan - Regulon names that meet once spaces go (`A B`, `A_B`) were given suffixes table by table, so a field could pair one regulon's AUC with another's targets; one map now names them in every table and in the thresholds. - An obs column named ClusterID was written as it was in place of SCope's own; ClusterID now always comes from the first clustering. - A text column's distinct values are read a block at a time and the scan stops once there are more than an annotation holds. Co-Authored-By: Claude Opus 5.5 --- docs/COMMANDS.md | 4 +- src/adata/conversion/loom.py | 6 +++ src/adata/conversion/scope.py | 88 +++++++++++++++++++++++++---------- tests/test_scope.py | 73 +++++++++++++++++++++++++++++ 4 files changed, 146 insertions(+), 25 deletions(-) diff --git a/docs/COMMANDS.md b/docs/COMMANDS.md index fa72c90..ac19f18 100644 --- a/docs/COMMANDS.md +++ b/docs/COMMANDS.md @@ -483,7 +483,9 @@ 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 with their metrics (`ClusterMarkers__avg_logFC` and the like, or those a -SCope import recorded). `MetaData` is plain JSON, which SCope reads as it +SCope import recorded). `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 diff --git a/src/adata/conversion/loom.py b/src/adata/conversion/loom.py index 3fb5950..8cdd265 100644 --- a/src/adata/conversion/loom.py +++ b/src/adata/conversion/loom.py @@ -808,6 +808,12 @@ def plan( varm_key=varm_key, clusterings=options.clusterings, default_embedding=options.default_embedding, warnings=warnings, ) + 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] diff --git a/src/adata/conversion/scope.py b/src/adata/conversion/scope.py index f240541..4588ab0 100644 --- a/src/adata/conversion/scope.py +++ b/src/adata/conversion/scope.py @@ -520,13 +520,12 @@ def plan_export( warnings: List[str], ) -> ScopeExport: """Everything a SCope loom needs, worked out and checked before writing.""" - from adata.elements.read import dataframe_columns, read_categories, read_str_all, resolve_index - from adata.util.names import natural_sorted + from adata.elements.read import dataframe_columns, read_categories, resolve_index parsed, recorded_embeddings, recorded_clusterings = _recorded(root) - # Compound tables: every dataframe in obsm and varm. - tables: List[ExportTable] = [] + # 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 []): @@ -542,22 +541,40 @@ def plan_export( if other: warnings.append(f"Not written in {where}/{key}: {', '.join(other)} -- a SCope " "table holds numbers only.") - if not numeric: - continue - regulons = key.endswith("Regulons") or key.endswith("RegulonsAUC") + if numeric: + 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.replace(" ", "_") if regulons else 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)) + 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 {} @@ -664,6 +681,8 @@ def plan_export( # Metrics and annotations: the obs columns written. metrics, annotations = [], [] for column in obs_columns: + if found and column == "ClusterID": + continue # written afresh from the first clustering member = obs[column] enc = spec.encoding_type(member) if enc == spec.CATEGORICAL: @@ -675,7 +694,9 @@ def plan_export( 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 = natural_sorted(set(read_str_all(member, parent_group=obs))) + 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 @@ -688,11 +709,13 @@ def plan_export( auc_fields = {f for t in tables if t.axis == "col" and t.key.endswith("RegulonsAUC") for f in t.fields} thresholds = [] for entry in parsed.get("regulonThresholds", []) or []: - if isinstance(entry, dict) and str(entry.get("regulon", "")).replace(" ", "_") in auc_fields: - thresholds.append(dict(entry, regulon=str(entry["regulon"]).replace(" ", "_"))) + name = str(entry.get("regulon", "")) if isinstance(entry, dict) else "" + field = regulon_fields.get(name, name.replace(" ", "_")) + if field in auc_fields: + thresholds.append(dict(entry, regulon=field)) if not thresholds and "uns" in root and "scope" in root["uns"] and "regulon_thresholds" in root["uns"]["scope"]: for regulon, row in _frame_rows(root["uns"]["scope"]["regulon_thresholds"]).items(): - name = regulon.replace(" ", "_") + name = regulon_fields.get(regulon, regulon.replace(" ", "_")) if name in auc_fields: entry = {"regulon": name} for key, value in row.items(): @@ -712,10 +735,27 @@ def plan_export( ) return ScopeExport( tables, embeddings, found, metadata, - cluster_id=bool(found) and "ClusterID" not in obs_columns, regulons=regulons, + cluster_id=bool(found), regulons=regulons, ) +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("/"): diff --git a/tests/test_scope.py b/tests/test_scope.py index 4d63998..4ffc6c2 100644 --- a/tests/test_scope.py +++ b/tests/test_scope.py @@ -521,3 +521,76 @@ def test_metadata_from_columns_and_a_thresholds_table(tmp_path): {"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"] From dab51569191ea5b5389240ceb16b39c8faa55b83 Mon Sep 17 00:00:00 2001 From: Aljes Date: Thu, 22 Oct 2026 16:03:27 +0100 Subject: [PATCH 03/10] convert --scope: renumber marker fields, keep new clusterings' markers, bound category scan - Marker tables' fields are cluster ids: renumbered with the clusters, so recorded ids 2, 5, 9 become 0, 1, 2 there as in MetaData. - A clustering named with --clusterings owns ClusterMarkers_; only recorded clusterings' tables were kept. - A categorical with more categories than an annotation lists is skipped by its length, without reading the categories. Co-Authored-By: Claude Opus 5.5 --- docs/COMMANDS.md | 3 +- src/adata/conversion/scope.py | 32 ++++++++++++++ tests/test_scope.py | 78 +++++++++++++++++++++++++++++++++++ 3 files changed, 112 insertions(+), 1 deletion(-) diff --git a/docs/COMMANDS.md b/docs/COMMANDS.md index ac19f18..cca305e 100644 --- a/docs/COMMANDS.md +++ b/docs/COMMANDS.md @@ -483,7 +483,8 @@ 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 with their metrics (`ClusterMarkers__avg_logFC` and the like, or those a -SCope import recorded). `ClusterID` is always written afresh from the first +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 diff --git a/src/adata/conversion/scope.py b/src/adata/conversion/scope.py index 4588ab0..b75709c 100644 --- a/src/adata/conversion/scope.py +++ b/src/adata/conversion/scope.py @@ -626,6 +626,9 @@ def plan_export( 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]] = {} for position, column in enumerate(wanted): if column not in obs: raise ValueError(f"--clusterings {column}: no such obs column.") @@ -639,6 +642,15 @@ def plan_export( 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 + } entry = { "id": position, "group": recorded.get("group", ""), @@ -648,6 +660,10 @@ def plan_export( if recorded.get("clusterMarkerMetrics"): entry["clusterMarkerMetrics"] = recorded["clusterMarkerMetrics"] found.append(ExportClustering(position, column, entry)) + # 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)] @@ -672,9 +688,23 @@ def plan_export( entry["clusterMarkerMetrics"] = [ {"accessor": a, "name": a, "description": ""} for a in accessors ] + fields_map = cluster_fields.get(new) for table in group_tables: accessor = _MARKERS.match(table.key).group(2) table.key = f"ClusterMarkers_{new}" + (f"_{accessor}" if accessor else "") + if fields_map is not None: + # Fields are cluster ids, renumbered with the clusters. + 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 not table.fields: + continue keep.append(table) tables = keep @@ -686,6 +716,8 @@ def plan_export( member = obs[column] enc = spec.encoding_type(member) if enc == spec.CATEGORICAL: + if "categories" in member and len(member["categories"]) > 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": diff --git a/tests/test_scope.py b/tests/test_scope.py index 4ffc6c2..076d005 100644 --- a/tests/test_scope.py +++ b/tests/test_scope.py @@ -594,3 +594,81 @@ def counted(obj, lo, hi, *args, **kwargs): 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 From 24a1b362996d510330a98955152e00d00275a4fe Mon Sep 17 00:00:00 2001 From: Aljes Date: Thu, 22 Oct 2026 16:21:41 +0100 Subject: [PATCH 04/10] convert to loom: one owner per attribute name; markers of undescribed clusters - An obs column and an obsm entry (or a SCope part, or Gene) of one name were both written to col_attrs/row_attrs, and the second write failed with the output half written. The plan now gives each name to the first that asks and leaves the rest out with a warning. - A clustering whose clusters MetaData did not describe is labelled by id on import; on export its marker fields now map through those labels instead of being dropped. Co-Authored-By: Claude Opus 5.5 --- docs/COMMANDS.md | 6 +++- src/adata/conversion/loom.py | 53 +++++++++++++++++++++++++++++++++++ src/adata/conversion/scope.py | 3 ++ tests/test_scope.py | 43 ++++++++++++++++++++++++++++ 4 files changed, 104 insertions(+), 1 deletion(-) diff --git a/docs/COMMANDS.md b/docs/COMMANDS.md index cca305e..a15cd48 100644 --- a/docs/COMMANDS.md +++ b/docs/COMMANDS.md @@ -490,7 +490,11 @@ 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. +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 diff --git a/src/adata/conversion/loom.py b/src/adata/conversion/loom.py index 8cdd265..5e8c6c3 100644 --- a/src/adata/conversion/loom.py +++ b/src/adata/conversion/loom.py @@ -824,6 +824,13 @@ def plan( "leaves them out." ) + obs_columns, obsm, var_columns, varm = _one_owner_per_name( + obs_columns, obsm, var_columns, varm, found, gene, accession, varm_key, warnings + ) + if found is not None: + for part in ("metrics", "annotations"): + found.metadata[part] = [m for m in found.metadata[part] if m["name"] in obs_columns] + dense = n_obs * n_var * max(1, main.dtype.itemsize) * (1 + len(layers)) summary = ( f"{src_format} → loom 3.0.0{' (SCope)' if found is not None else ''}, " @@ -1022,6 +1029,52 @@ def _columns(group: Any, index_name: str, warnings: List[str], axis: str) -> Lis return kept +def _one_owner_per_name( + obs_columns: List[str], obsm: List[str], var_columns: List[str], varm: List[str], + found: Optional[scope.ScopeExport], gene: str, accession: Optional[str], varm_key: str, + warnings: List[str], +) -> Tuple[List[str], List[str], List[str], List[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 the same name, and the second + write failed with the output half written. The first to ask keeps it: + loom's and SCope's own attributes, then columns, SCope's tables, arrays. + 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 found is not None: + own = owners["col_attrs"] + if any(e.id == -1 for e in found.embeddings): + own["Embedding"] = "SCope's default embedding" + if found.embeddings: + own.update(Embeddings_X="SCope's embeddings", Embeddings_Y="SCope's embeddings") + if found.clusterings: + own["Clusterings"] = "SCope's clusterings" + if found.cluster_id: + own["ClusterID"] = "SCope's first clustering" + + def claim_name(group: str, name: str, what: str) -> bool: + 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 + + skip = {gene, accession, "Gene", "Accession"} + obs_columns = [c for c in obs_columns if c == "CellID" or claim_name("col_attrs", c, f"obs column {c!r}")] + var_columns = [c for c in var_columns if c in skip or claim_name("row_attrs", c, f"var column {c!r}")] + if found is not None: + found.tables = [ + t for t in found.tables + if claim_name("col_attrs" if t.axis == "col" else "row_attrs", t.key, t.source) + ] + obsm = [k for k in obsm if claim_name("col_attrs", k, f"obsm/{k}")] + varm = [k for k in varm if claim_name("row_attrs", k, f"{varm_key}/{k}")] + return obs_columns, obsm, var_columns, varm + + def _gene_attrs(columns: List[str]) -> Tuple[str, Optional[str]]: """Where loom's Gene (symbol) and Accession (id) come from. diff --git a/src/adata/conversion/scope.py b/src/adata/conversion/scope.py index b75709c..3ada942 100644 --- a/src/adata/conversion/scope.py +++ b/src/adata/conversion/scope.py @@ -650,6 +650,9 @@ def plan_export( 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, diff --git a/tests/test_scope.py b/tests/test_scope.py index 076d005..7df1cff 100644 --- a/tests/test_scope.py +++ b/tests/test_scope.py @@ -672,3 +672,46 @@ def counted(column, *args, **kwargs): 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",) From f29dc92c14098cf9224a061362ccf1eb81f9c2db Mon Sep 17 00:00:00 2001 From: Aljes Date: Thu, 22 Oct 2026 17:12:05 +0100 Subject: [PATCH 05/10] convert --scope: thresholds from the table as edited, marker metrics from the tables written Co-Authored-By: Claude Opus 5.5 --- docs/COMMANDS.md | 6 ++-- src/adata/conversion/scope.py | 64 +++++++++++++++++++++-------------- tests/test_scope.py | 40 ++++++++++++++++++++++ 3 files changed, 82 insertions(+), 28 deletions(-) diff --git a/docs/COMMANDS.md b/docs/COMMANDS.md index a15cd48..5c7f805 100644 --- a/docs/COMMANDS.md +++ b/docs/COMMANDS.md @@ -468,7 +468,7 @@ adata convert scenic.h5ad -o scenic.loom --scope --default-embedding X_umap | `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, 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, and the recorded ids and names | +| `uns["scope"]` | `MetaData`'s regulon thresholds (from the `regulon_thresholds` table, as edited), and the recorded ids and names | | Flag | Meaning | |---|---| @@ -482,8 +482,8 @@ 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 with -their metrics (`ClusterMarkers__avg_logFC` and the like, or those a -SCope import recorded), their fields renumbered with the clusters; a +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 diff --git a/src/adata/conversion/scope.py b/src/adata/conversion/scope.py index 3ada942..7353dee 100644 --- a/src/adata/conversion/scope.py +++ b/src/adata/conversion/scope.py @@ -629,6 +629,8 @@ def plan_export( #: 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.") @@ -660,8 +662,8 @@ def plan_export( "name": name or column, "clusters": [{"id": i, "description": c} for i, c in enumerate(categories)], } - if recorded.get("clusterMarkerMetrics"): - entry["clusterMarkerMetrics"] = recorded["clusterMarkerMetrics"] + 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)) # A clustering no import recorded owns the marker tables of its new id. for clustering in found: @@ -681,16 +683,18 @@ def plan_export( continue entry = found[new].entry accessors = [m.group(2) for t in group_tables if (m := _MARKERS.match(t.key)).group(2)] - if not entry.get("clusterMarkerMetrics"): - if not accessors: - warnings.append( - f"Not written: ClusterMarkers_{old} -- without metrics (_avg_logFC, " - "_pval, ...) SCope cannot show markers." - ) - continue - entry["clusterMarkerMetrics"] = [ - {"accessor": a, "name": a, "description": ""} for a in accessors - ] + # 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: + warnings.append( + f"Not written: ClusterMarkers_{old} -- without metrics (_avg_logFC, " + "_pval, ...) SCope cannot show markers." + ) + continue + entry["clusterMarkerMetrics"] = metrics fields_map = cluster_fields.get(new) for table in group_tables: accessor = _MARKERS.match(table.key).group(2) @@ -742,20 +746,30 @@ def plan_export( # 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} - thresholds = [] + 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 []: - name = str(entry.get("regulon", "")) if isinstance(entry, dict) else "" - field = regulon_fields.get(name, name.replace(" ", "_")) - if field in auc_fields: - thresholds.append(dict(entry, regulon=field)) - if not thresholds and "uns" in root and "scope" in root["uns"] and "regulon_thresholds" in root["uns"]["scope"]: - for regulon, row in _frame_rows(root["uns"]["scope"]["regulon_thresholds"]).items(): - name = regulon_fields.get(regulon, regulon.replace(" ", "_")) - if name in auc_fields: - entry = {"regulon": name} - for key, value in row.items(): - entry[key] = _float(value) if key == "defaultThresholdValue" else value - thresholds.append(entry) + 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": + entry[key] = _float(value) + 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( diff --git a/tests/test_scope.py b/tests/test_scope.py index 7df1cff..7b6df01 100644 --- a/tests/test_scope.py +++ b/tests/test_scope.py @@ -715,3 +715,43 @@ def test_two_elements_asking_for_one_attribute_name_keep_the_first_and_warn(tmp_ 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": ""}, + ] From 3f4b477fc3b14ffbe619817951c328d6301c5ac5 Mon Sep 17 00:00:00 2001 From: Aljes Date: Thu, 22 Oct 2026 17:19:52 +0100 Subject: [PATCH 06/10] convert --scope: clusterings once, markers pruned after name clashes, anndata 0.7 categoricals - A clustering went out as Clusterings and as a plain attribute; read back, the copy took its name and the clustering became leiden-1. - Marker metrics are pruned again once tables lose a name to another attribute, and markers without their base table or any metric are left out whole. - --clusterings takes anndata 0.7's categoricals (codes with a categories reference), and they are annotations, not metrics. Co-Authored-By: Claude Opus 5.5 --- docs/COMMANDS.md | 2 +- src/adata/conversion/loom.py | 7 +++++- src/adata/conversion/scope.py | 47 +++++++++++++++++++++++++++++------ tests/test_scope.py | 45 +++++++++++++++++++++++++++++++++ 4 files changed, 92 insertions(+), 9 deletions(-) diff --git a/docs/COMMANDS.md b/docs/COMMANDS.md index 5c7f805..a965536 100644 --- a/docs/COMMANDS.md +++ b/docs/COMMANDS.md @@ -466,7 +466,7 @@ adata convert scenic.h5ad -o scenic.loom --scope --default-embedding X_umap |---|---| | 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, with `ClusterID` for the first; their clusters described in `MetaData` | +| 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 | diff --git a/src/adata/conversion/loom.py b/src/adata/conversion/loom.py index 5e8c6c3..1492171 100644 --- a/src/adata/conversion/loom.py +++ b/src/adata/conversion/loom.py @@ -808,6 +808,10 @@ def plan( varm_key=varm_key, clusterings=options.clusterings, default_embedding=options.default_embedding, warnings=warnings, ) + # A clustering is written once, as SCope's: a plain copy too + # would come back from the loom as a second column. + clustered = {c.column for c in found.clusterings} + obs_columns = [c for c in obs_columns if c not in clustered] if found.cluster_id and "ClusterID" in obs_columns: # A stale copy would disagree with Clusterings once the # first clustering is edited. @@ -828,6 +832,7 @@ def plan( obs_columns, obsm, var_columns, varm, found, gene, accession, varm_key, warnings ) if found is not None: + scope.prune_markers(found, warnings) for part in ("metrics", "annotations"): found.metadata[part] = [m for m in found.metadata[part] if m["name"] in obs_columns] @@ -951,7 +956,7 @@ def _export_scope( 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(obs[col]["codes"][lo:hi], dtype=np.int64))) + (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) diff --git a/src/adata/conversion/scope.py b/src/adata/conversion/scope.py index 7353dee..1e4e1b9 100644 --- a/src/adata/conversion/scope.py +++ b/src/adata/conversion/scope.py @@ -43,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" @@ -634,7 +634,7 @@ def plan_export( for position, column in enumerate(wanted): if column not in obs: raise ValueError(f"--clusterings {column}: no such obs column.") - if spec.encoding_type(obs[column]) != spec.CATEGORICAL: + 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)." @@ -715,15 +715,16 @@ def plan_export( keep.append(table) tables = keep - # Metrics and annotations: the obs columns written. + # Metrics and annotations: the obs columns written; clusterings are SCope's own. + clustering_columns = {c.column for c in found} metrics, annotations = [], [] for column in obs_columns: - if found and column == "ClusterID": - continue # written afresh from the first clustering + if found and column == "ClusterID" or column in clustering_columns: + continue # written afresh from the first clustering; a clustering member = obs[column] enc = spec.encoding_type(member) - if enc == spec.CATEGORICAL: - if "categories" in member and len(member["categories"]) > MAX_ANNOTATION_VALUES: + if is_categorical(member): + if is_group(member) and len(member["categories"]) > 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 \ @@ -788,6 +789,38 @@ def plan_export( ) +def prune_markers(found: "ScopeExport", warnings: List[str]) -> None: + """Keep the cluster markers whole once tables have been left out after + planning (a name another attribute took): a clustering's metrics list + only the accessor tables still written, and without its + `ClusterMarkers_` table, or any metric, none of them is written.""" + keys = {t.key for t in found.tables if t.axis == "row"} + for clustering in found.clusterings: + base = f"ClusterMarkers_{clustering.id}" + metrics = [m for m in clustering.entry.get("clusterMarkerMetrics", []) + if f"{base}_{m['accessor']}" in keys] + if base in keys and metrics: + clustering.entry["clusterMarkerMetrics"] = metrics + continue + if base in keys or any(k.startswith(base + "_") for k in keys): + warnings.append(f"Not written: {base} and its metrics -- part of them is left out.") + clustering.entry.pop("clusterMarkerMetrics", None) + found.tables = [t for t in found.tables + if not (t.axis == "row" and (t.key == base or t.key.startswith(base + "_")))] + + +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.""" diff --git a/tests/test_scope.py b/tests/test_scope.py index 7b6df01..25e4f99 100644 --- a/tests/test_scope.py +++ b/tests/test_scope.py @@ -755,3 +755,48 @@ def test_marker_metrics_list_the_tables_written(tmp_path): {"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") + assert "Not written: ClusterMarkers_0 and its metrics" in _out(result) + 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"] From b68045af323ccedeb8c8dd2ffefd8303dfc4bc82 Mon Sep 17 00:00:00 2001 From: Aljes Date: Thu, 22 Oct 2026 17:29:20 +0100 Subject: [PATCH 07/10] convert --scope: settle attribute names before planning; strict JSON; bounded legacy categories - Names are given out before the plan derives anything from the tables: a table that loses its name (an obs column RegulonsAUC) is gone before thresholds, the regulon count or marker metrics are worked out from it, so MetaData only describes tables written. Marker tables claim their names once renumbered. - NaN and infinities in MetaData become null, and a missing defaultThresholdValue is left out: SCope parses strict JSON. - A legacy (anndata 0.7) categorical's categories are counted by the length of the dataset they point to, before any is read (elements.read.categories_dataset). Co-Authored-By: Claude Opus 5.5 --- src/adata/conversion/loom.py | 77 ++++++++++++++--------------------- src/adata/conversion/scope.py | 76 ++++++++++++++++++++++------------ src/adata/elements/read.py | 22 +++++----- tests/test_scope.py | 66 +++++++++++++++++++++++++++++- 4 files changed, 159 insertions(+), 82 deletions(-) diff --git a/src/adata/conversion/loom.py b/src/adata/conversion/loom.py index 1492171..98c9226 100644 --- a/src/adata/conversion/loom.py +++ b/src/adata/conversion/loom.py @@ -801,12 +801,16 @@ def plan( obsp = _graphs(root, "obsp", n_obs, warnings) varp = [] if raw else _graphs(root, "varp", n_var, warnings) uns = _uns_scalars(root, warnings) + obs_columns, var_columns, owners = _column_owners( + obs_columns, 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, ) # A clustering is written once, as SCope's: a plain copy too # would come back from the loom as a second column. @@ -828,13 +832,10 @@ def plan( "leaves them out." ) - obs_columns, obsm, var_columns, varm = _one_owner_per_name( - obs_columns, obsm, var_columns, varm, found, gene, accession, varm_key, warnings - ) - if found is not None: - scope.prune_markers(found, warnings) - for part in ("metrics", "annotations"): - found.metadata[part] = [m for m in found.metadata[part] if m["name"] in obs_columns] + # 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 = ( @@ -916,7 +917,7 @@ def write( _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(detail.scope.metadata)) + _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"): @@ -1034,50 +1035,34 @@ def _columns(group: Any, index_name: str, warnings: List[str], axis: str) -> Lis return kept -def _one_owner_per_name( - obs_columns: List[str], obsm: List[str], var_columns: List[str], varm: List[str], - found: Optional[scope.ScopeExport], gene: str, accession: Optional[str], varm_key: str, - warnings: List[str], -) -> Tuple[List[str], List[str], List[str], List[str]]: +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 the same name, and the second - write failed with the output half written. The first to ask keeps it: - loom's and SCope's own attributes, then columns, SCope's tables, arrays. - The others are left out, each with a warning.""" + 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 found is not None: - own = owners["col_attrs"] - if any(e.id == -1 for e in found.embeddings): - own["Embedding"] = "SCope's default embedding" - if found.embeddings: - own.update(Embeddings_X="SCope's embeddings", Embeddings_Y="SCope's embeddings") - if found.clusterings: - own["Clusterings"] = "SCope's clusterings" - if found.cluster_id: - own["ClusterID"] = "SCope's first clustering" - - def claim_name(group: str, name: str, what: str) -> bool: - 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 - + 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"} - obs_columns = [c for c in obs_columns if c == "CellID" or claim_name("col_attrs", c, f"obs column {c!r}")] - var_columns = [c for c in var_columns if c in skip or claim_name("row_attrs", c, f"var column {c!r}")] - if found is not None: - found.tables = [ - t for t in found.tables - if claim_name("col_attrs" if t.axis == "col" else "row_attrs", t.key, t.source) - ] - obsm = [k for k in obsm if claim_name("col_attrs", k, f"obsm/{k}")] - varm = [k for k in varm if claim_name("row_attrs", k, f"{varm_key}/{k}")] - return obs_columns, obsm, var_columns, varm + 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]]: diff --git a/src/adata/conversion/scope.py b/src/adata/conversion/scope.py index 1e4e1b9..3d6aa36 100644 --- a/src/adata/conversion/scope.py +++ b/src/adata/conversion/scope.py @@ -518,11 +518,18 @@ def plan_export( 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.""" - from adata.elements.read import dataframe_columns, read_categories, resolve_index + """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 = [] @@ -541,7 +548,10 @@ def plan_export( if other: warnings.append(f"Not written in {where}/{key}: {', '.join(other)} -- a SCope " "table holds numbers only.") - if numeric: + # Marker tables claim their names once renumbered, below. + if numeric 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 @@ -682,6 +692,16 @@ def plan_export( 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 + group_tables = named 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. @@ -697,8 +717,6 @@ def plan_export( entry["clusterMarkerMetrics"] = metrics fields_map = cluster_fields.get(new) for table in group_tables: - accessor = _MARKERS.match(table.key).group(2) - table.key = f"ClusterMarkers_{new}" + (f"_{accessor}" if accessor else "") if fields_map is not None: # Fields are cluster ids, renumbered with the clusters. gone = [f for f in table.fields if f not in fields_map] @@ -724,7 +742,7 @@ def plan_export( member = obs[column] enc = spec.encoding_type(member) if is_categorical(member): - if is_group(member) and len(member["categories"]) > MAX_ANNOTATION_VALUES: + 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 \ @@ -765,7 +783,11 @@ def plan_export( entry = dict(recorded_thresholds.get(name, {}), regulon=name) for key, value in row.items(): if key == "defaultThresholdValue": - entry[key] = _float(value) + 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) @@ -789,24 +811,28 @@ def plan_export( ) -def prune_markers(found: "ScopeExport", warnings: List[str]) -> None: - """Keep the cluster markers whole once tables have been left out after - planning (a name another attribute took): a clustering's metrics list - only the accessor tables still written, and without its - `ClusterMarkers_` table, or any metric, none of them is written.""" - keys = {t.key for t in found.tables if t.axis == "row"} - for clustering in found.clusterings: - base = f"ClusterMarkers_{clustering.id}" - metrics = [m for m in clustering.entry.get("clusterMarkerMetrics", []) - if f"{base}_{m['accessor']}" in keys] - if base in keys and metrics: - clustering.entry["clusterMarkerMetrics"] = metrics - continue - if base in keys or any(k.startswith(base + "_") for k in keys): - warnings.append(f"Not written: {base} and its metrics -- part of them is left out.") - clustering.entry.pop("clusterMarkerMetrics", None) - found.tables = [t for t in found.tables - if not (t.axis == "row" and (t.key == base or t.key.startswith(base + "_")))] +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: 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_scope.py b/tests/test_scope.py index 25e4f99..9e48650 100644 --- a/tests/test_scope.py +++ b/tests/test_scope.py @@ -785,7 +785,10 @@ def test_markers_whose_table_lost_its_name_are_left_out_with_their_metrics(tmp_p 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 and its metrics" in _out(result) + 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" @@ -800,3 +803,64 @@ def test_an_anndata_07_categorical_can_be_a_clustering(tmp_path, sample_legacy_v 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"] == [] From 78287f05842ef86292375692291274d14cad3cd0 Mon Sep 17 00:00:00 2001 From: Aljes Date: Thu, 22 Oct 2026 17:38:28 +0100 Subject: [PATCH 08/10] convert --scope: marker metrics from the tables left after renumbering Co-Authored-By: Claude Opus 5.5 --- src/adata/conversion/scope.py | 39 ++++++++++++++++++----------------- tests/test_scope.py | 25 ++++++++++++++++++++++ 2 files changed, 45 insertions(+), 19 deletions(-) diff --git a/src/adata/conversion/scope.py b/src/adata/conversion/scope.py index 3d6aa36..9b13a3d 100644 --- a/src/adata/conversion/scope.py +++ b/src/adata/conversion/scope.py @@ -701,24 +701,12 @@ def plan_export( 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 - group_tables = named - 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: - warnings.append( - f"Not written: ClusterMarkers_{old} -- without metrics (_avg_logFC, " - "_pval, ...) SCope cannot show markers." - ) - continue - entry["clusterMarkerMetrics"] = metrics + # 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) - for table in group_tables: + group_tables = [] + for table in named: if fields_map is not None: - # Fields are cluster ids, renumbered with the clusters. 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 " @@ -728,9 +716,22 @@ def plan_export( 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 not table.fields: - continue - keep.append(table) + 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) and any(t.key == base for t in named): + 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. diff --git a/tests/test_scope.py b/tests/test_scope.py index 9e48650..5dbd9c0 100644 --- a/tests/test_scope.py +++ b/tests/test_scope.py @@ -864,3 +864,28 @@ def test_an_auc_table_that_loses_its_name_takes_its_regulon_metadata_with_it(tmp 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"] From 7ef1d08f8272ce903db18e71ab88f091cc9d5baa Mon Sep 17 00:00:00 2001 From: Aljes Date: Thu, 22 Oct 2026 17:43:50 +0100 Subject: [PATCH 09/10] convert --scope: clusterings claim no attribute name; markers need their base table Co-Authored-By: Claude Opus 5.5 --- docs/COMMANDS.md | 4 ++-- src/adata/conversion/loom.py | 11 ++++++----- src/adata/conversion/scope.py | 24 ++++++++++++++++-------- tests/test_scope.py | 30 ++++++++++++++++++++++++++++++ 4 files changed, 54 insertions(+), 15 deletions(-) diff --git a/docs/COMMANDS.md b/docs/COMMANDS.md index a965536..95137b1 100644 --- a/docs/COMMANDS.md +++ b/docs/COMMANDS.md @@ -481,8 +481,8 @@ 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 with -their metrics (`ClusterMarkers__avg_logFC` and the like, which +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 diff --git a/src/adata/conversion/loom.py b/src/adata/conversion/loom.py index 98c9226..90b724c 100644 --- a/src/adata/conversion/loom.py +++ b/src/adata/conversion/loom.py @@ -801,8 +801,13 @@ def plan( 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( - obs_columns, var_columns, gene, accession, options.scope, warnings + [c for c in obs_columns if c not in clustered], var_columns, gene, accession, + options.scope, warnings, ) found = None if options.scope: @@ -812,10 +817,6 @@ def plan( default_embedding=options.default_embedding, warnings=warnings, owners=owners, ) - # A clustering is written once, as SCope's: a plain copy too - # would come back from the loom as a second column. - clustered = {c.column for c in found.clusterings} - obs_columns = [c for c in obs_columns if c not in clustered] if found.cluster_id and "ClusterID" in obs_columns: # A stale copy would disagree with Clusterings once the # first clustering is edited. diff --git a/src/adata/conversion/scope.py b/src/adata/conversion/scope.py index 9b13a3d..33e9bc8 100644 --- a/src/adata/conversion/scope.py +++ b/src/adata/conversion/scope.py @@ -510,6 +510,18 @@ def _scalar_text(ds: Any) -> str: return value.decode("utf-8") if isinstance(value, bytes) else str(value) +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, *, @@ -628,11 +640,7 @@ def plan_export( # Clusterings: the recorded ones, and those --clusterings names. obs = root["obs"] - wanted = [c for c in recorded_clusterings if c in obs] - for column in (clusterings or "").split(","): - column = column.strip() - if column and column not in wanted: - wanted.append(column) + 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] = {} @@ -724,7 +732,7 @@ def plan_export( 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) and any(t.key == base for t in named): + 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." @@ -735,10 +743,10 @@ def plan_export( tables = keep # Metrics and annotations: the obs columns written; clusterings are SCope's own. - clustering_columns = {c.column for c in found} + clustered = {c.column for c in found} metrics, annotations = [], [] for column in obs_columns: - if found and column == "ClusterID" or column in clustering_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) diff --git a/tests/test_scope.py b/tests/test_scope.py index 5dbd9c0..99dd6fe 100644 --- a/tests/test_scope.py +++ b/tests/test_scope.py @@ -889,3 +889,33 @@ def test_a_marker_table_left_with_no_cluster_is_not_listed_as_a_metric(tmp_path) 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] From a733beff9342d5586f7867cac2402a5bdbd3652c Mon Sep 17 00:00:00 2001 From: Aljes Date: Thu, 22 Oct 2026 18:04:46 +0100 Subject: [PATCH 10/10] convert --scope: obsm tables claim their names; embeddings and clusterings checked against HDF5's limit in the plan Co-Authored-By: Claude Opus 5.5 --- src/adata/conversion/scope.py | 15 ++++++++++++++- tests/test_scope.py | 28 ++++++++++++++++++++++++++++ 2 files changed, 42 insertions(+), 1 deletion(-) diff --git a/src/adata/conversion/scope.py b/src/adata/conversion/scope.py index 33e9bc8..596bff2 100644 --- a/src/adata/conversion/scope.py +++ b/src/adata/conversion/scope.py @@ -510,6 +510,15 @@ def _scalar_text(ds: Any) -> str: 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).""" @@ -561,7 +570,7 @@ def plan_export( 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 (_MARKERS.match(key) or claim_attr( + 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)) @@ -637,6 +646,8 @@ def plan_export( 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"] @@ -683,6 +694,8 @@ def plan_export( 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: diff --git a/tests/test_scope.py b/tests/test_scope.py index 99dd6fe..f89fb03 100644 --- a/tests/test_scope.py +++ b/tests/test_scope.py @@ -919,3 +919,31 @@ def test_metric_tables_without_their_markers_table_are_left_out(tmp_path): 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()