Add API conveniences and eight fixes, with every existing call unchanged

cytepool/_analysis_records.py

@@ -6,8 +6,10 @@
 it fitted, with the contrasts it tested; and the versions of the software.
 """
 
+import json
 import platform
 from collections import Counter
+from collections.abc import Mapping
 from typing import Any
 
 import numba
@@ -71,6 +73,73 @@ def mixed_design_record(design: Design, contrasts: Contrasts) -> dict[str, Any]:
     return record
 
 
+def keep_provenance(
+    record: Any, method: str, cube: Cube, settings: Mapping[str, Any]
+) -> None:
+    """Give the record of scores of a cube's group, such as ``PathwayScores``, its
+    provenance: the method, the software, the cube's fingerprint, digests and
+    provenance, the layer read and the settings. ``provenance_of`` reads it back.
+
+    The record keeps it as JSON text in its own attributes, outside its fields, so
+    that a record built from fields alone, or copied with ``dataclasses.replace``,
+    has none.
+    """
+    text = json.dumps(
+        {
+            "method": method,
+            "versions": software_versions(),
+            "cube": cube_record(cube),
+            "layer": "counts",
+            "settings": dict(settings),
+        },
+        allow_nan=False,
+    )
+    object.__setattr__(record, "_provenance", text)
+
+
+def provenance_of(record: Any) -> dict[str, Any] | None:
+    """The provenance that ``keep_provenance`` gave a record, as a new dict, or
+    None."""
+    text = vars(record).get("_provenance")
+    if text is None:
+        return None
+    found: dict[str, Any] = json.loads(text)
+    return found
+
+
+def provenance_text(record: Any) -> str | None:
+    """The JSON text of the provenance that ``keep_provenance`` gave a record, or
+    None: what a pickled record carries, for ``restore_provenance``."""
+    text: str | None = vars(record).get("_provenance")
+    return text
+
+
+def restore_provenance[R](record: R, text: str | None) -> R:
+    """Give a record the provenance it was pickled with, if it had one."""
+    if text is not None:
+        object.__setattr__(record, "_provenance", text)
+    return record
+
+
+def only_key(value: Any, available: list[str], argument: str, what: str) -> Any:
+    """The key of a result's table: ``value`` when given, else the only one.
+
+    Results that keep one table per group, contrast or both let a key be left out
+    when they have one. ``what`` names the keys in the message, such as "groups".
+
+    Raises:
+        ValueError: If ``value`` is None and there is not exactly one key.
+    """
+    if value is not None:
+        return value
+    if len(available) == 1:
+        return available[0]
+    raise ValueError(
+        f"the result has {len(available)} {what}; name one with {argument}"
+        + (f", one of {quoted(available)}" if available else "")
+    )
+
+
 def software_versions() -> dict[str, str]:
     """The versions of CytePool, numpy, numba and Python that run an analysis."""
     return {

cytepool/_expression_values.py

@@ -3,9 +3,11 @@
 
 A cube's group is read as limma_voom reads it: the samples with at least ``min_cells``
 cells and some counts, over the genes that edgeR's ``filterByExpr`` keeps
-(:func:`~cytepool.normalization.filter_by_expression`), as log2 counts per million
-after TMM normalization with a prior count of 2, edgeR's ``cpm(log = TRUE)`` (Robinson
-and Oshlack 2010; Chen, Lun and Smyth 2016), or as counts.
+(:func:`~cytepool.normalization.filter_by_expression`), with the samples as one group
+or grouped by a column of the cube's samples (``{"sample_groups": "condition"}``), as
+log2 counts per million after TMM normalization with a prior count of 2, edgeR's
+``cpm(log = TRUE)`` (Robinson and Oshlack 2010; Chen, Lun and Smyth 2016), or as
+counts.
 
 References:
     Robinson MD, Oshlack A (2010). A scaling normalization method for differential
@@ -17,13 +19,20 @@
     quasi-likelihood pipeline. F1000Research 5, 1438. doi:10.12688/f1000research.8987.2
 """
 
+from collections.abc import Mapping
 from typing import Any, NamedTuple
 
 import numpy as np
 import numpy.typing as npt
 
 from ._columns import text_column
 from ._gene_ids import is_missing
+from ._group_fits import (
+    FILTER_DEFAULTS,
+    SAMPLE_GROUPS,
+    check_sample_groups,
+    filter_setting,
+)
 from ._masked_arrays import Missing, real_array, reject_masked, unmasked_array
 from ._messages import quoted
 from ._options import checked_choice, checked_flag, integral_number
@@ -47,6 +56,7 @@ class ExpressionValues(NamedTuple):
     samples: Table
     gene_ids: list[str | None]  # one per column; None for a gene without an id
     settings: dict[str, Any]  # how a cube's values were taken
+    cube: Cube | None = None  # the cube whose group was read
 
 
 def expression_values(
@@ -73,7 +83,8 @@ def expression_values(
         values: For a cube, ``"log_cpm"`` or ``"counts"``.
         min_cells: For a cube, the fewest cells a sample needs in the group.
         expression_filter: For a cube, whether to keep only the genes that
-            filter_by_expression keeps.
+            filter_by_expression keeps: True or False, or a mapping of its settings,
+            which may name a column of the cube's samples as ``sample_groups``.
         missing: How the values mark missing ones, for messages.
         matched_to: What the gene ids are matched to, for messages, such as "the
             gene sets".
@@ -200,7 +211,13 @@ def _cube_values(
     fewest = integral_number(min_cells, "min_cells")
     if fewest < 0:
         raise ValueError(f"min_cells must be at least 0, got {fewest}")
-    filtered = checked_flag(expression_filter, "expression_filter")
+    if isinstance(expression_filter, Mapping):
+        filter_settings = filter_setting(expression_filter)
+        check_sample_groups({"expression_filter": filter_settings}, cube.samples)
+        filtered: Any = dict(filter_settings or {})
+    else:
+        filtered = checked_flag(expression_filter, "expression_filter")
+        filter_settings = dict(FILTER_DEFAULTS) if filtered else None
     if group is None:
         if cube.n_groups != 1:
             raise ValueError(
@@ -223,8 +240,17 @@ def _cube_values(
         )
     kept = counts[rows]
     genes = np.ones(kept.shape[1], dtype=bool)
-    if filtered:
-        genes = filter_by_expression(kept, group=np.zeros(rows.size))
+    if filter_settings is not None:
+        numbers = {
+            name: value
+            for name, value in filter_settings.items()
+            if name != SAMPLE_GROUPS
+        }
+        if SAMPLE_GROUPS in filter_settings:
+            labels = cube.samples[filter_settings[SAMPLE_GROUPS]][rows]
+        else:
+            labels = np.zeros(rows.size)
+        genes = filter_by_expression(kept, group=labels, **numbers)
         if not genes.any():
             raise ValueError(
                 f"no gene of group {group!r} passes filter_by_expression; pass "
@@ -265,6 +291,7 @@ def _cube_values(
         cube.samples.take(rows),
         ids,
         settings,
+        cube,
     )
 
 

cytepool/_group_fits.py

@@ -14,8 +14,9 @@
    column must be estimable. With random effects, those of the design are checked
    on the samples used.
 3. Features. Features are filtered by expression with the group's design, by edgeR's
-   ``filterByExpr`` (Chen, Lun and Smyth 2016), and the library sizes are taken over
-   the features kept.
+   ``filterByExpr`` (Chen, Lun and Smyth 2016), or with the groups of samples that a
+   column of the design's samples gives (``filterByExpr(y, group = ...)``), and the
+   library sizes are taken over the features kept.
 4. Normalization. TMM normalization factors (Robinson and Oshlack 2010), when asked.
 
 This module also checks the options these workflows share, and builds the table of
@@ -63,6 +64,10 @@
 }
 """edgeR's ``filterByExpr`` defaults."""
 
+SAMPLE_GROUPS = "sample_groups"
+"""The setting of ``expression_filter`` that names a column of the samples whose
+values group them for ``filterByExpr``, in place of the design."""
+
 
 @dataclass
 class GroupRecord:
@@ -211,6 +216,14 @@ def prepare_group(
     filter_settings = settings["expression_filter"]
     if filter_settings is None:
         keep = np.ones(y.shape[1], dtype=bool)
+    elif SAMPLE_GROUPS in filter_settings:
+        numbers = {
+            name: value
+            for name, value in filter_settings.items()
+            if name != SAMPLE_GROUPS
+        }
+        labels = design.samples[filter_settings[SAMPLE_GROUPS]]
+        keep = filter_by_expression(y, group=labels[inputs.model_rows[used]], **numbers)
     else:
         keep = filter_by_expression(y, design=x_kept, **filter_settings)
     fit.tested = np.flatnonzero(keep)
@@ -361,22 +374,35 @@ def threads_setting(value: Any) -> int | None:
     return number
 
 
-def filter_setting(value: Any) -> dict[str, float] | None:
-    """The settings of ``filterByExpr``: its defaults, changed ones, or None."""
+def filter_setting(value: Any) -> dict[str, Any] | None:
+    """The settings of ``filterByExpr``: its defaults, changed ones, or None.
+
+    A mapping may also name, as ``sample_groups``, a column of the samples whose
+    values group them; ``check_sample_groups`` checks it against the samples.
+    """
     reject_masked(value, "expression_filter")
     if isinstance(value, bool | np.bool_):
         return dict(FILTER_DEFAULTS) if value else None
     if not isinstance(value, Mapping):
         raise TypeError(
             "expression_filter must be True, False or a mapping of filterByExpr "
-            f"settings such as {{'min_count': 5}}, got {type(value).__name__}"
+            f"settings such as {{'min_count': 5}} or {{'{SAMPLE_GROUPS}': "
+            f"'condition'}}, got {type(value).__name__}"
         )
-    settings = dict(FILTER_DEFAULTS)
+    settings: dict[str, Any] = dict(FILTER_DEFAULTS)
     for name, number in value.items():
+        if name == SAMPLE_GROUPS:
+            if not isinstance(number, str):
+                raise TypeError(
+                    f"expression_filter[{SAMPLE_GROUPS!r}] must name a column of the "
+                    f"samples (str), got {type(number).__name__}"
+                )
+            settings[name] = number
+            continue
         if name not in FILTER_DEFAULTS:
             raise ValueError(
                 f"expression_filter has the unknown setting {name!r}; the settings are "
-                f"{', '.join(FILTER_DEFAULTS)}"
+                f"{', '.join([*FILTER_DEFAULTS, SAMPLE_GROUPS])}"
             )
         settings[name] = non_negative_number(number, f"expression_filter[{name!r}]")
     if settings["min_prop"] > 1.0:
@@ -386,6 +412,41 @@ def filter_setting(value: Any) -> dict[str, float] | None:
     return settings
 
 
+def check_sample_groups(settings: Mapping[str, Any], samples: Table) -> None:
+    """Refuse a column of ``sample_groups`` that the samples do not have, or in which
+    a sample has no value.
+
+    Raises:
+        KeyError: If ``samples`` has no such column.
+        ValueError: If a value of the column is missing.
+    """
+    filter_settings = settings["expression_filter"]
+    if filter_settings is None or SAMPLE_GROUPS not in filter_settings:
+        return
+    name = filter_settings[SAMPLE_GROUPS]
+    if name not in samples:
+        raise KeyError(
+            f"expression_filter[{SAMPLE_GROUPS!r}] names {name!r}, which is not a "
+            f"column of the samples; their columns are: {listing(list(samples.columns))}"
+        )
+    values = samples[name].tolist()
+    missing = [
+        position
+        for position, value in enumerate(values)
+        if value is None or value != value
+    ]
+    if missing:
+        where = (
+            f"sample {samples['sample_id'][missing[0]]!r}"
+            if "sample_id" in samples
+            else f"the sample in row {missing[0]}"
+        )
+        raise ValueError(
+            f"expression_filter[{SAMPLE_GROUPS!r}] names {name!r}, which is missing for "
+            f"{where}; every sample needs a group for the filter"
+        )
+
+
 def normalization_setting(value: Any) -> str | None:
     """The normalization: "tmm" or None."""
     if value is None or (isinstance(value, str) and value == "tmm"):

cytepool/_messages.py

@@ -1,7 +1,10 @@
-"""Text in messages and reprs: lists of ids and names, hints, ranges of values."""
+"""Text in messages and reprs: lists of ids and names, hints, ranges of values, and
+lines that fit the width of a page."""
 
 import difflib
 import os
+import re
+import textwrap
 from collections.abc import Sequence
 from typing import Any
 
@@ -10,6 +13,17 @@
 LISTED = 8
 """The most items a message names before it says how many more there are."""
 
+WIDTH = 80
+"""The most characters in a line of a repr or a printed record, as pages show them."""
+
+KEEP_TOGETHER = "\ue000"
+"""A space at which :func:`fitted_lines` does not break a line, such as the one in
+``counts (uint32)``; the text it returns has a plain space there."""
+
+_NUMBERED = re.compile(r"\d+\. ")
+# A line that begins with a short key and a colon, such as "  layers: ...".
+_KEYED = re.compile(r" *[^\s:][^:]{0,29}: ")
+
 PACKAGE_DIRECTORY = os.path.join(os.path.dirname(os.path.abspath(__file__)), "")
 """The package's directory, ending in a separator: warnings skip the files under it,
 so that a warning is attributed to the first caller outside the package."""
@@ -39,6 +53,54 @@ def close_matches(name: str, candidates: tuple[str, ...] | list[str]) -> str:
     return " (did you mean " + " or ".join(repr(match) for match in found) + "?)"
 
 
+def fitted_lines(text: str, width: int = WIDTH) -> str:
+    """Text whose lines are at most ``width`` characters, broken between words.
+
+    A broken line continues below the text after its key, as in ``layers: ...``,
+    or after the number of a numbered line (``1. ...``), and otherwise, or where its
+    longest word would not fit there, four spaces deeper than it began. A word
+    longer than a line stays whole, and so do words joined by ``KEEP_TOGETHER``; a
+    fingerprint, 71 characters, fits on a line of its own.
+    """
+    lines = []
+    for line in text.split("\n"):
+        if len(line) <= width:
+            lines.append(line.replace(KEEP_TOGETHER, " "))
+            continue
+        body = line.lstrip(" ")
+        indent = line[: len(line) - len(body)]
+        numbered = _NUMBERED.match(body)
+        keyed = _KEYED.match(line)
+        longest = max(len(word) for word in line.split())
+        if numbered:
+            hanging = indent + " " * numbered.end()
+        elif keyed and keyed.end() <= width // 2:
+            hanging = " " * keyed.end()
+        else:
+            hanging = indent + "    "
+        if len(hanging) + longest > width:
+            hanging = indent + "    "
+        lines.extend(
+            part.replace(KEEP_TOGETHER, " ")
+            for part in textwrap.wrap(
+                line,
+                width=width,
+                subsequent_indent=hanging,
+                break_long_words=False,
+                break_on_hyphens=False,
+            )
+        )
+    return "\n".join(lines)
+
+
+def fitted_message(text: str, width: int = WIDTH) -> str:
+    """A message, such as a warning's, broken into lines of at most ``width``
+    characters between words."""
+    return "\n".join(
+        textwrap.wrap(text, width=width, break_long_words=False, break_on_hyphens=False)
+    )
+
+
 def value_range(values: float | np.ndarray) -> str:
     """A value, or the range of the finite values of an array, as reprs show them."""
     if isinstance(values, float):