Round-trip analysis results through AnnData

cytearc/datastore/namespaces.py

@@ -180,6 +180,9 @@ class ReductionAccessor(_Accessor):
     lsi = _bind(_GraphOperationsMixin._reduction_lsi)
     custom = _bind(_GraphOperationsMixin._reduction_custom)
     harmony = _bind(_GraphOperationsMixin._reduction_harmony)
+    import_coordinates = _bind(
+        _PresentationOperationsMixin._reduction_import_coordinates
+    )
 
 
 class RegulonsAccessor(_Accessor):

cytearc/embeddings/imported.py

@@ -1,3 +1,4 @@
+import math
 from collections.abc import Callable, Iterable, Iterator, Mapping, Sequence
 from dataclasses import dataclass
 from typing import Any, cast
@@ -121,6 +122,47 @@ def _validate_source_digest(value: bytes) -> bytes:
     return value
 
 
+def _source_value(value: Any, path: str) -> Any:
+    if value is None or isinstance(value, bool | str):
+        return value
+    if isinstance(value, int | float):
+        if isinstance(value, float) and not math.isfinite(value):
+            raise ValueError(f"{path} must be a finite number")
+        return value
+    if isinstance(value, Mapping):
+        if any(not isinstance(key, str) for key in value):
+            raise TypeError(f"{path} must have text keys")
+        return {
+            key: _source_value(item, f"{path}[{key!r}]") for key, item in value.items()
+        }
+    if isinstance(value, list | tuple):
+        return [
+            _source_value(item, f"{path}[{index}]") for index, item in enumerate(value)
+        ]
+    raise TypeError(
+        f"{path} holds a {type(value).__name__}; source records hold text, "
+        "numbers, booleans, None, and lists or mappings of them"
+    )
+
+
+def imported_source_record(source: Mapping[str, Any]) -> dict[str, Any]:
+    """Return a validated record of the tool that computed imported coordinates.
+
+    The record, such as ``{"tool": "scvi-tools", "version": "1.2.0",
+    "parameters": {"n_latent": 30}}``, is stored in the artifact's
+    provenance, so it is part of its identity and of its lineage.
+
+    Raises:
+        TypeError: If ``source`` is not a mapping with text keys, or holds a
+            value that is not text, a number, a boolean, None, or a list or
+            mapping of them.
+        ValueError: If it holds a non-finite number.
+    """
+    if not isinstance(source, Mapping):
+        raise TypeError("source must be a mapping")
+    return cast(dict[str, Any], _source_value(source, "source"))
+
+
 def _string_source_length(values: ImportedStringSource) -> int:
     shape = getattr(values, "shape", None)
     if shape is not None:
@@ -435,13 +477,26 @@ def write_imported_coordinates(
     stdev_dtype: Any | None = None,
     block_rows: int = 100_000,
     invalidate_cache: bool = False,
+    source: Mapping[str, Any] | None = None,
+    require_writable: Callable[[], None] | None = None,
 ) -> ArtifactRef:
+    """Store, or reuse, imported graph coordinates over an exact cell selection.
+
+    The artifact's identity holds the source key, role, dimensions, which
+    optional payloads are stored, the source digest, the fingerprint of each
+    payload, the fingerprint of the ordered source cell IDs, and the cell
+    selection. ``source``, a record of the tool that computed the
+    coordinates (:func:`imported_source_record`), is stored in the
+    parameters when given. ``require_writable`` is called before anything is
+    written, when no stored artifact matches.
+    """
     storage = ImportedArtifactStorage(root)
     source_digest = _validate_source_digest(source_digest)
     if not isinstance(dimreduc_key, str) or not dimreduc_key:
         raise ValueError("dimreduc_key must be a non-empty string")
     if not isinstance(role, str) or not role or role.lower() in {"umap", "tsne"}:
         raise ValueError("Imported graph coordinates require a non-layout role")
+    source_record = None if source is None else imported_source_record(source)
     block_rows = integer_argument(block_rows, "block_rows", minimum=1)
     data_source = _resolve_source(
         coordinates,
@@ -529,17 +584,20 @@ def write_imported_coordinates(
         )
         for name, shape in shapes.items()
     )
+    parameters: dict[str, Any] = {
+        "dimreduc_key": dimreduc_key,
+        "role": role.lower(),
+        "dims": dims,
+        "loadings_stored": loading_source is not None,
+        "feature_ids_stored": feature_source is not None,
+        "stdev_stored": stdev_source is not None,
+    }
+    if source_record is not None:
+        parameters["source"] = source_record
     planned = storage.plan(
         assay=assay,
         kind="imported_coordinates",
-        parameters={
-            "dimreduc_key": dimreduc_key,
-            "role": role.lower(),
-            "dims": dims,
-            "loadings_stored": loading_source is not None,
-            "feature_ids_stored": feature_source is not None,
-            "stdev_stored": stdev_source is not None,
-        },
+        parameters=parameters,
         inputs={
             "source_digest": source_digest,
             "payload_fingerprints": fingerprints,
@@ -557,6 +615,8 @@ def write_imported_coordinates(
         ),
     )
     if not planned.reused:
+        if require_writable is not None:
+            require_writable()
         with storage.transaction(planned) as group:
             _write_numeric_payload(
                 storage,

cytearc/readers/_h5ad_columns.py

@@ -300,6 +300,122 @@ def read_table_column(
     return values, missing
 
 
+@dataclass(frozen=True, slots=True)
+class H5adCategories:
+    """The declared categories of an AnnData categorical column.
+
+    Attributes:
+        values: Every category in its stored order, unused ones included,
+            decoded as the column's values decode: byte strings as text and
+            NumPy scalars as Python scalars.
+        ordered: Whether the categories are ordered.
+    """
+
+    values: tuple[Any, ...]
+    ordered: bool
+
+
+def _category_value(value: Any) -> Any:
+    if isinstance(value, bytes | np.bytes_):
+        return as_text(value)
+    return value.item() if isinstance(value, np.generic) else value
+
+
+def column_categories(
+    h5: h5py.File,
+    table: H5adNode,
+    key: str,
+    category_groups: tuple[str, ...],
+) -> H5adCategories | None:
+    """Return the categories of a categorical column, or None for another column.
+
+    A group-encoded categorical keeps its categories beside its codes and
+    records ``ordered`` as an attribute. AnnData 0.7 and 0.6 kept the
+    categories of integer codes apart from them, as
+    :func:`legacy_categories` finds them, and such categories are unordered
+    unless their dataset records ``ordered``.
+    """
+    if isinstance(table, h5py.Group):
+        node = table[key]
+        if is_categorical(node):
+            return H5adCategories(
+                tuple(_category_value(value) for value in node["categories"][:]),
+                bool(node.attrs.get("ordered", False)),
+            )
+        if not isinstance(node, h5py.Dataset):
+            return None
+        stored = node.dtype
+    else:
+        stored = np.dtype(table.dtype.fields[key][0])
+    if stored.kind not in "iu":
+        return None
+    categories = legacy_categories(h5, table, key, category_groups)
+    if categories is None:
+        return None
+    ordered = False
+    if isinstance(table, h5py.Group):
+        for name in category_groups:
+            group = table.get(name)
+            if isinstance(group, h5py.Group) and key in group:
+                ordered = bool(group[key].attrs.get("ordered", False))
+                break
+    return H5adCategories(
+        tuple(_category_value(value) for value in categories), ordered
+    )
+
+
+_HEX_DIGITS = frozenset("0123456789abcdefABCDEF")
+
+
+def hex_color(value: str) -> str | None:
+    """Return a color as ``#rrggbb``, or None when it is not a color.
+
+    ``#rgb``, ``#rrggbb``, and ``#rrggbbaa`` are read directly, dropping any
+    alpha. Other names, such as ``red`` or ``tab:blue``, are resolved with
+    Matplotlib when it is installed.
+    """
+    text = value.strip()
+    if text.startswith("#") and set(text[1:]) <= _HEX_DIGITS:
+        digits = text[1:]
+        if len(digits) in {3, 4}:
+            digits = "".join(digit * 2 for digit in digits[:3])
+        if len(digits) in {6, 8}:
+            return f"#{digits[:6].lower()}"
+        return None
+    try:
+        from matplotlib.colors import to_hex
+    except ImportError:
+        return None
+    try:
+        return str(to_hex(text, keep_alpha=False))
+    except ValueError:
+        return None
+
+
+def category_colors(h5: h5py.File, key: str, count: int) -> tuple[str, ...] | None:
+    """Return the ``uns["<key>_colors"]`` of a categorical with ``count`` categories.
+
+    AnnData and Scanpy keep one color per category, in category order. Like
+    Scanpy, a list with more colors than categories gives each category the
+    color at its position. A list with fewer colors, or with a value that is
+    not a color, is not used, and None is returned.
+    """
+    node = h5.get(f"uns/{key}_colors")
+    if not isinstance(node, h5py.Dataset) or node.ndim != 1 or count < 1:
+        return None
+    if int(node.shape[0]) < count:
+        return None
+    colors = []
+    for value in node[:count]:
+        if not isinstance(value, str | bytes | np.str_ | np.bytes_):
+            return None
+        color = hex_color(as_text(value))
+        if color is None:
+            return None
+        colors.append(color)
+    return tuple(colors)
+
+
 def table_column_dtype(
     h5: h5py.File,
     table: H5adNode,