Align several studies' cubes by sample, gene and cell type

README.md

@@ -16,6 +16,7 @@ what is planned.
 | Provenance and fingerprints | A record of each cube's source (the caller's record, kept unchanged), of the settings that built it and of what each operation kept, merged or pooled, with the fingerprint of the cube it was applied to; `Cube.history()` reads it back. A content fingerprint (SHA-256, documented encoding) of the tables, cell numbers, normalization and layers | a reimplementation of the documented encoding, zarr-python |
 | Cube files | `write_cube`, `read_cube` (a subset of Zarr v3, with orthogonal selections; a read opens only the chunks that hold the selected features; gzip, and zstd on Python 3.14 or newer; stored digests, checked when a file is opened and when a read covers a whole group) | zarr-python |
 | AnnData and MuData | `to_anndata`, `from_anndata`, `to_mudata`, `from_mudata` (optional extras `anndata` and `mudata`): cubes to and from the long AnnData of `scanpy.get.aggregate` and `decoupler.pp.pseudobulk` and the MuData of LIANA's `adata_to_views`, whose layouts are recognised; an import takes sums and counts of cells, refuses means, and gives a cube exported and read back with the same fingerprint | scanpy, decoupler, LIANA, anndata, mudata |
+| Several studies | `CubeCollection`: the cubes of several studies aligned, each kept as it is: samples named `study::sample_id`, with donor ids alike; genes aligned by id, with which study measured which, so that a gene a study did not measure is missing there, not zero; cell types aligned by Cell Ontology id, group id or a mapping for each study, merging those that a study splits; each study's cube on the collection's axes, for analyses study by study; and one cube of the genes and cell types every study has, for a joint model. `align_study_results` lines up the studies' results gene by gene, with the reason a study has none for a gene | anndata, pandas, limma |
 | Normalization and filtering | TMM, upper-quartile and RLE factors, DESeq2 size factors, CPM, log-CPM, `filter_by_expression` | edgeR, DESeq2 |
 | Designs | Model matrices from R formulas, named contrasts, design checks | R, limma, lme4 |
 | Linear models | Ordinary, weighted and generalised least squares for every gene (in parallel when each gene needs its own fit), contrasts (named contrasts are matched to the coefficients by name) | limma |
@@ -55,6 +56,8 @@ These are not implemented yet:
 - Dose-response fits, benchmark doses and synergy
 - Meta-regression (meta-analysis with study-level covariates) and multivariate
   meta-analysis of correlated effects
+- Combined p-values and replication measures across studies, and cell types aligned
+  up the Cell Ontology
 
 ## Known limits
 

cytepool/_columns.py

@@ -11,8 +11,9 @@
 import numpy as np
 
 
-def text_column(values: Sequence[str]) -> np.ndarray:
-    """A text column for a ``Table``, also when it is empty.
+def text_column(values: Sequence[str | None]) -> np.ndarray:
+    """A text column for a ``Table``, also when it is empty; None marks a missing
+    value.
 
     numpy reads an empty list as float64, which would make an empty column of ids or
     reasons a column of numbers; an object array stays a text column.

cytepool/cube/__init__.py

@@ -1,6 +1,7 @@
 """Pseudobulk cubes: additive statistics per sample, group and feature."""
 
+from .collection import CubeCollection
 from .cube import Cube
 from .pool import from_bulk, from_cells
 
-__all__ = ["Cube", "from_bulk", "from_cells"]
+__all__ = ["Cube", "CubeCollection", "from_bulk", "from_cells"]

cytepool/cube/_axes.py

@@ -8,6 +8,8 @@
 - ``checked_axis_table``: an axis table checked, with the position of each id.
 - ``checked_n_cells``: the number of cells of every sample in every group, from an
   integer array or a DataFrame labelled by the ids, as a read-only int64 array.
+- ``selected_positions``: the positions that ids, a boolean mask or positions select
+  along an axis, as ``Cube.select`` and ``CubeCollection.select`` read them.
 - ``axis_codes``: the position on an axis of every cell's label, with the axis table.
 - ``count_dtype``: the dtype of a counts layer that holds pooled counts exactly.
 - ``first_index``: the index of the first entry, in C order, where a mask is true.
@@ -190,6 +192,86 @@ def first_index(mask: np.ndarray) -> tuple[int, ...]:
     return tuple(int(position) for position in np.argwhere(mask)[0])
 
 
+def selected_positions(
+    selector: Any,
+    ids: np.ndarray,
+    index: dict[str, int],
+    axis: str,
+    single: str | None = None,
+    named: bool = False,
+) -> np.ndarray | None:
+    """Convert ids, a boolean mask or positions along one axis to positions.
+
+    A single id or position selects one entry; ``None`` selects every entry and
+    gives None. ``axis`` names the axis in messages, such as ``"samples"``, and
+    ``single`` one of its entries (by default ``axis`` without its last letter);
+    ``named`` says that its ids are names, as those of studies are.
+
+    Raises:
+        KeyError: If an id is unknown.
+        IndexError: If a position is out of range.
+        TypeError: If the selection is not ids, a boolean mask or integers, or is a
+            numpy masked array.
+        ValueError: If a mask has the wrong length or the selection repeats an entry.
+    """
+    if selector is None:
+        return None
+    n = len(index)
+    single = axis[:-1] if single is None else single
+    if isinstance(selector, str) or (
+        isinstance(selector, int | np.integer)
+        and not isinstance(selector, bool | np.bool_)
+    ):
+        selector = [selector]
+    values = unmasked_array(selector, axis)
+    if values.ndim != 1:
+        raise ValueError(f"{axis} must be a 1-D selection, got shape {values.shape}")
+    kind = values.dtype.kind
+    if kind == "b":
+        if values.size != n:
+            raise ValueError(
+                f"a boolean mask of {axis} needs {n} values, got {values.size}"
+            )
+        return np.flatnonzero(values)
+    if values.size == 0:
+        return np.empty(0, dtype=np.intp)
+    if kind in "iu":
+        outside = (values < -n) | (values >= n)
+        if outside.any():
+            raise IndexError(
+                f"{single} position {values[outside][0]} is out of range for {n} {axis}"
+            )
+        positions = values.astype(np.intp) % n
+    elif kind in "UTO":
+        positions = np.empty(values.size, dtype=np.intp)
+        for i, value in enumerate(values.tolist()):
+            if not isinstance(value, str):
+                raise TypeError(
+                    f"{axis} must be ids, a boolean mask or integer positions; got "
+                    f"{type(value).__name__} {value!r} among ids"
+                )
+            position = index.get(value)
+            if position is None:
+                raise KeyError(
+                    f"no {single} {'named' if named else 'with id'} {value!r}"
+                )
+            positions[i] = position
+    else:
+        raise TypeError(
+            f"{axis} must be ids, a boolean mask or integer positions, got dtype "
+            f"{values.dtype}"
+        )
+    seen: set[int] = set()
+    for position in positions.tolist():
+        if position in seen:
+            raise ValueError(
+                f"the selection of {axis} repeats {ids[position]!r}; each {single} can "
+                "be selected once"
+            )
+        seen.add(position)
+    return positions
+
+
 def axis_codes(
     labels: Any,
     name: str,