From 65aa56c2fdb2da0f0b7a3fc3a48b9939c33fb1eb Mon Sep 17 00:00:00 2001 From: parashardhapola Date: Mon, 26 Oct 2026 18:26:51 +0100 Subject: [PATCH] Align several studies' cubes by sample, gene and cell type CubeCollection aligns the cubes of several named studies, each kept as it was given: samples keyed study::sample_id (and donor columns the caller names), genes by feature id with a per-study measured mask so that a gene a study did not measure is missing rather than zero, and cell types by Cell Ontology id, group id or a mapping per study. It gives each study's cube on the collection's axes, selections, and one stacked cube for a joint model, which refuses anything a study lacks rather than filling it with zeros. align_study_results lines up the studies' differential expression gene by gene, with the reason a study has no result for a gene. Checked against anndata 0.13.4, pandas 3.1.0 and limma 3.62.1 (validation/cube_collection.py: 15 checks, none failed). --- README.md | 3 + cytepool/_columns.py | 5 +- cytepool/cube/__init__.py | 3 +- cytepool/cube/_axes.py | 82 ++ cytepool/cube/collection.py | 1485 ++++++++++++++++++++ cytepool/cube/cube.py | 74 +- cytepool/cube/provenance.py | 68 +- cytepool/meta_analysis.py | 212 ++- docs/coverage_exemptions.toml | 1 + docs/source/index.md | 4 +- docs/source/reference/api/cubes.md | 5 +- docs/source/reference/api/meta_analysis.md | 5 +- docs/source/reference/bibliography.md | 3 +- docs/source/reference/capabilities.md | 3 + docs/source/reference/cube_files.md | 2 +- docs/source/references.bib | 11 + docs/source/statistics/cubes.md | 33 +- docs/source/statistics/meta_analysis.md | 8 + docs/source/validation.md | 47 +- tests/data/cube_collection/studies.json | 621 ++++++++ tests/data/cube_collection/studies.npz | Bin 0 -> 58047 bytes tests/test_cube_collection.py | 883 ++++++++++++ tests/test_masked_inputs.py | 54 +- validation/cube_collection.py | 686 +++++++++ 24 files changed, 4217 insertions(+), 81 deletions(-) create mode 100644 cytepool/cube/collection.py create mode 100644 tests/data/cube_collection/studies.json create mode 100644 tests/data/cube_collection/studies.npz create mode 100644 tests/test_cube_collection.py create mode 100644 validation/cube_collection.py diff --git a/README.md b/README.md index 2e493dd..1439db5 100644 --- a/README.md +++ b/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 diff --git a/cytepool/_columns.py b/cytepool/_columns.py index 47cd825..dd26540 100644 --- a/cytepool/_columns.py +++ b/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. diff --git a/cytepool/cube/__init__.py b/cytepool/cube/__init__.py index b02af22..d0210aa 100644 --- a/cytepool/cube/__init__.py +++ b/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"] diff --git a/cytepool/cube/_axes.py b/cytepool/cube/_axes.py index 9b81aae..c9945eb 100644 --- a/cytepool/cube/_axes.py +++ b/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, diff --git a/cytepool/cube/collection.py b/cytepool/cube/collection.py new file mode 100644 index 0000000..53d87da --- /dev/null +++ b/cytepool/cube/collection.py @@ -0,0 +1,1485 @@ +"""Several studies' cubes, aligned by sample, gene and group: ``CubeCollection``. + +Studies of one question, such as cohorts of a disease, are compared or combined only +once their samples, genes and cell types are aligned, and that is where errors enter: +donor ids that repeat across studies, genes that a study did not measure read as +zeros. ``CubeCollection`` aligns them and says what each study measured and has. + +References: + Nakatsuka N, Adler D, Jiang L, Hartman A, Cheng E, Klann E, Satija R (2025). + Improving reproducibility of differentially expressed genes in single-cell + transcriptomic studies of neurodegenerative diseases through meta-analysis. + Nature Communications 16, 7436. doi:10.1038/s41467-025-62579-z + + Diehl AD, Meehan TF, Bradford YM, Brush MH, Dahdul WM, Dougall DS, He Y, + Osumi-Sutherland D, Ruttenberg A, Sarntivijai S, Van Slyke CE, Vasilevsky NA, + Haendel MA, Blake JA, Mungall CJ (2016). The Cell Ontology 2016: enhanced content, + modularization, and ontology interoperability. Journal of Biomedical Semantics 7, + 44. doi:10.1186/s13326-016-0088-7 +""" + +import json +import re +from collections.abc import Callable, Mapping, Sequence +from types import MappingProxyType +from typing import Any + +import numpy as np + +from .. import __version__ +from .._columns import text_column, values_equal +from .._masked_arrays import reject_masked +from .._messages import close_matches, listing, quoted +from ..table import Table +from ._axes import INT64_MAX, selected_positions +from .cube import LAYER_NAMES, Cube + +__all__ = ["CubeCollection"] + +_SEPARATOR = "::" # between a study's name and a sample id: study::sample_id +_ONTOLOGY_COLUMN = "cell_ontology_id" +_NORMALIZED_LAYERS = ("log_sum", "log_sum_sq", "linear_sum") +_FLOAT_LIMIT = 2**53 # float64 holds every whole number below it +# An Ensembl version suffix, as in ENSG00000141510.17. +_VERSION = re.compile(r"\.\d+$") +_KINDS = { + "O": "text", + "b": "booleans", + "i": "integers", + "u": "integers", + "f": "floats", + "c": "complex numbers", + "M": "datetimes", + "m": "timedeltas", +} +_NUMBERS = frozenset({"integers", "floats", "complex numbers"}) + + +class CubeCollection: + """Several studies' cubes, aligned by sample, gene and group. + + Most genes found in one study are not found in the next (Nakatsuka et al. + 2025), and the studies are compared or combined only once their samples, genes + and groups (cell types) are aligned. A collection aligns the cubes of several + named studies and keeps each cube as it was given (``cubes``): it reads no layer, + and each analysis of a study reads that study's cube. It is immutable: + ``select`` returns a new collection. + + Samples: + Sample ids mean something only inside their study, so the collection names + every sample ``study::sample_id``: sample ``S1`` of one study is not sample + ``S1`` of another, and study names cannot hold ``::``. The same holds for + donor ids and other columns of ids, which ``within_study`` names: their + values are written ``study::value`` too, so that a design does not take + donor ``D1`` of two studies for one person. ``samples`` has ``sample_id``, + ``study`` and every other column of any study's samples, in the order in + which the studies, taken in order, first have them; a study without a column + has missing values in it. A column keeps its type when every study gives it + the same; integers of different types become int64, integers and floats + float64, datetimes take the finer unit, and integer and boolean columns that + need a missing value become float64 with NaN, as in ``Table``. A column that + holds different kinds of values in different studies, such as text in one + and numbers in another, is refused. + + Genes: + Features are aligned by ``feature_id``, which must hold ids of one kind in + every study, such as Ensembl gene ids. A study measured a feature when its + cube has it and does not mark it False in the feature column ``measured``: + tools that give every study the same genes fill those a study did not + measure with zeros, and ``measured`` tells them apart. The collection's + features are those that at least one study measured, in the order in which + the studies first list them, and ``measured`` says which study measured + which: a feature that a study did not measure is missing there, not zero. + ``features`` has ``feature_id`` and, when a study's features have one, + ``symbol``, the first symbol a study gives the feature, taking the studies + in order. + + Groups: + Groups are aligned by the values of a text column of the group tables, + ``align_groups_by``: by default ``cell_ontology_id``, the Cell Ontology id + (Diehl et al. 2016), when every study's group table has that column, and + ``group_id`` otherwise. ``group_mapping`` gives the aligned id of some of a + study's groups instead, or None to leave a group out. Each aligned group is + named by its aligned id, and groups of one study with the same aligned id are + merged: their layers and cell numbers are summed, exactly, as + ``Cube.regroup`` sums them. ``group_members`` lists the groups of each study + under their aligned ids. ``groups`` has ``group_id``, the aligned id, and + every other column of any study's group table: a group keeps a column's + value where the studies that have the column agree on it for every member, + and has a missing value otherwise. + + Analyses of each study: + ``aligned_cube(study)`` is the study's cube on the collection's axes: its + groups under their aligned ids and the features it measured, in the + collection's order, and its own samples with their own ids, so that a design + built from its sample table applies. Any analysis of one cube, such as + ``limma_voom`` with each study's own design, then gives results keyed by + the aligned groups, which ``align_study_results`` lines up gene by gene and + ``differential_expression_meta_analysis`` combines + (``cytepool.meta_analysis``). + + Stacking: + ``stack()`` puts the studies into one cube, for a joint model with the study + as a covariate: its samples are those of every study, as ``samples`` lists + them, and its features and groups those that every study measured and has. + A cube holds no missing values, so ``stack`` refuses a feature or a group + that a study lacks rather than fill it with zeros. Each study keeps its own + library sizes, and normalized layers are stacked only from studies + normalized alike. + + Provenance: + ``provenance`` records each study's cube by its fingerprint, with its own + provenance, the settings of the alignment with each study's groups under + their aligned ids, and each selection since. A stacked cube records the + studies and the alignment in its provenance (``cytepool.cube.provenance``). + + Not done: + Groups are aligned by equal values only: rolling cell types up to a common + ancestor in the Cell Ontology needs the ontology as a resource, so give such + merges in ``group_mapping``. Feature ids are not converted, so Ensembl ids + with versions (``ENSG00000141510.17``) must have them removed first. + ``stack`` gives a cube in memory. + + Args: + cubes: Study name to cube, at least one, in the order the collection keeps. + Study names are non-empty and hold no ``::``. + align_groups_by: The text column of the group tables whose values align the + groups, such as ``"cell_ontology_id"`` or ``"group_id"``. By default + ``"cell_ontology_id"`` when every study's group table has it, and + ``"group_id"`` otherwise. + group_mapping: Study name to a mapping of its group ids to aligned ids, which + take the place of those groups' values of ``align_groups_by``; None + leaves a group out. Groups of one study given the same aligned id are + merged. + within_study: Sample columns whose values mean something only inside their + study, such as donor ids, written ``study::value`` as sample ids are. + + Raises: + TypeError: If an argument has the wrong type, a column of the sample or + group tables holds different kinds of values in different studies, a + column of ``within_study`` holds neither text nor integers, or + ``group_mapping`` or ``within_study`` is or holds a numpy masked array. + ValueError: If there is no study, a study name is empty or holds ``::``, a + group has no value of ``align_groups_by`` and no aligned id, a study has + no group or measured no feature, or a study shares no feature, or no + group, with any other study. + KeyError: If ``group_mapping`` names a study or a group that is not there, + or ``within_study`` or ``align_groups_by`` a column that no study has. + + Example: + >>> import numpy as np + >>> from cytepool.cube import Cube + >>> def study(genes, donors): + ... return Cube( + ... {"counts": np.ones((1, len(donors), len(genes)), dtype=np.int64)}, + ... samples={"sample_id": donors}, + ... groups={"group_id": ["B"]}, + ... features={"feature_id": genes}, + ... ) + >>> collection = CubeCollection( + ... {"A": study(["g1", "g2"], ["D1", "D2"]), "B": study(["g2", "g3"], ["D1"])} + ... ) + >>> collection.samples["sample_id"].tolist() + ['A::D1', 'A::D2', 'B::D1'] + >>> collection.measured + array([[ True, True, False], + [False, True, True]]) + >>> collection.stack().features["feature_id"].tolist() + ['g2'] + """ + + __slots__ = ( + "_studies", + "_cubes", + "_members", + "_measured", + "_samples", + "_groups", + "_features", + "_settings", + "_operations", + "_study_index", + "_group_index", + "_feature_index", + ) + + _studies: tuple[str, ...] + _cubes: dict[str, Cube] + _members: dict[str, dict[str, tuple[str, ...]]] # study: aligned id: its groups + _measured: np.ndarray # bool, (n_studies, n_features), read-only + _samples: Table + _groups: Table + _features: Table + _settings: str # JSON: the alignment's settings and the groups it aligned + _operations: tuple[str, ...] # JSON texts of the selections, oldest first + _study_index: dict[str, int] + _group_index: dict[str, int] + _feature_index: dict[str, int] + + def __init__( + self, + cubes: Mapping[str, Cube], + *, + align_groups_by: str | None = None, + group_mapping: Mapping[str, Mapping[str, str | None]] | None = None, + within_study: str | Sequence[str] = (), + ) -> None: + studies = _checked_cubes(cubes) + given = {study: cubes[study] for study in studies} + within = _checked_within_study(within_study, given) + mapping = _checked_group_mapping(group_mapping, given) + by = _alignment_column(align_groups_by, given) + members = { + study: _aligned_groups(study, cube, by, mapping.get(study, {})) + for study, cube in given.items() + } + group_ids = list( + dict.fromkeys(group for study in studies for group in members[study]) + ) + feature_ids, measured = _measured_features(given) + for position, study in enumerate(studies): + if not members[study]: + raise ValueError( + f"study {study!r} has no group to align: group_mapping leaves out " + "every group of its cube" + ) + if not measured[position].any(): + raise ValueError( + f"study {study!r} measured no feature: its cube has no features, " + "or marks every one as not measured in its feature column " + "'measured'" + ) + if len(studies) > 1: + _check_shared_features(studies, feature_ids, measured) + _check_shared_groups(studies, members, by) + settings = { + "cytepool_version": __version__, + "align_groups_by": by, + "group_mapping": {study: dict(groups) for study, groups in mapping.items()} + or None, + "within_study": within, + "groups": { + study: {aligned: list(ids) for aligned, ids in members[study].items()} + for study in studies + }, + } + _assemble( + self, + studies, + given, + { + study: {aligned: tuple(ids) for aligned, ids in members[study].items()} + for study in studies + }, + group_ids, + feature_ids, + measured, + json.dumps(settings), + (), + ) + + # Immutability, copying and pickling. + + def __setattr__(self, name: str, value: object) -> None: + raise AttributeError("CubeCollection is immutable; select() returns a new one") + + def __delattr__(self, name: str) -> None: + raise AttributeError("CubeCollection is immutable") + + def __copy__(self) -> "CubeCollection": + return self + + def __deepcopy__(self, memo: dict[int, object]) -> "CubeCollection": + return self + + def __reduce__(self) -> tuple[Callable[..., "CubeCollection"], tuple[Any, ...]]: + return ( + _restore_collection, + ( + self._studies, + self._cubes, + self._members, + self._measured, + self._samples, + self._groups, + self._features, + self._settings, + self._operations, + ), + ) + + # Sizes and tables. + + @property + def studies(self) -> tuple[str, ...]: + """The study names, in order.""" + return self._studies + + @property + def n_studies(self) -> int: + """The number of studies.""" + return len(self._studies) + + @property + def n_samples(self) -> int: + """The number of samples of every study together.""" + return len(self._samples) + + @property + def n_groups(self) -> int: + """The number of aligned groups.""" + return len(self._groups) + + @property + def n_features(self) -> int: + """The number of features that at least one study measured.""" + return len(self._features) + + @property + def cubes(self) -> Mapping[str, Cube]: + """The cube of every study as it was given, by study name (read-only).""" + return MappingProxyType(self._cubes) + + @property + def samples(self) -> Table: + """The samples of every study, with ``sample_id`` (``study::sample_id``), + ``study`` and every other sample column.""" + return self._samples + + @property + def groups(self) -> Table: + """The aligned groups, with ``group_id``, the aligned id, and the group + columns on whose values the studies agree.""" + return self._groups + + @property + def features(self) -> Table: + """The features that at least one study measured, with ``feature_id`` and + ``symbol`` when a study's features have it.""" + return self._features + + @property + def measured(self) -> np.ndarray: + """Whether each study measured each feature, as read-only booleans of shape + ``(n_studies, n_features)``.""" + return self._measured.view() + + @property + def group_members(self) -> Table: + """The groups of each study under their aligned ids, one row per group of a + study that the collection holds, with ``study``, ``group_id`` (the aligned + id) and ``study_group_id`` (the group's id in the study's cube), studies and + aligned groups in order. An aligned group with several rows for one study + merges them; a study without a row for an aligned group does not have it.""" + rows = [ + (study, aligned, member) + for study in self._studies + for aligned, ids in self._members[study].items() + for member in ids + ] + return Table( + { + name: text_column([row[position] for row in rows]) + for position, name in enumerate(("study", "group_id", "study_group_id")) + } + ) + + @property + def provenance(self) -> dict[str, Any]: + """A JSON-serialisable record of how the collection was made, as a new dict. + + It holds ``studies``: each study's cube, by study name in order, with its + ``fingerprint`` and its own ``provenance``; ``alignment``: the version of + CytePool that aligned them, ``align_groups_by``, ``group_mapping``, + ``within_study`` and ``groups``, each study's groups under their aligned + ids when the collection was made; and ``operations``: each selection since, + with the studies, groups and features it kept (null for an axis given as + None). The fingerprints cost what ``Cube.fingerprint`` costs: nothing for a + cube file, which stores them, and a pass over the values of a cube in memory + the first time, at about 1 s per GB. + """ + alignment = json.loads(self._settings) + return { + "studies": { + study: { + "fingerprint": self._cubes[study].fingerprint(), + "provenance": self._cubes[study].provenance, + } + for study in self._studies + }, + "alignment": alignment, + "operations": [json.loads(text) for text in self._operations], + } + + def __repr__(self) -> str: + settings = json.loads(self._settings) + lines = [ + f"CubeCollection of {self.n_studies} stud" + f"{'y' if self.n_studies == 1 else 'ies'}: {self.n_samples:,} samples, " + f"{self.n_groups:,} groups, {self.n_features:,} features", + f" groups aligned by {settings['align_groups_by']!r}", + ] + for position, study in enumerate(self._studies): + members = self._members[study] + merged = sum(len(ids) > 1 for ids in members.values()) + n_samples = self._cubes[study].n_samples + text = ( + f" {study}: {n_samples:,} samples, {len(members):,} groups" + + (f" ({merged} merged)" if merged else "") + + f", {int(self._measured[position].sum()):,} features measured" + ) + lines.append(text) + if self.n_studies > 1: + shared_groups = sum( + all(group in self._members[study] for study in self._studies) + for group in self._groups["group_id"].tolist() + ) + shared_features = int(self._measured.all(axis=0).sum()) + lines.append( + f" in every study: {shared_groups:,} groups, {shared_features:,} " + "features" + ) + return "\n".join(lines) + + # Each study on the collection's axes. + + def aligned_cube(self, study: str) -> Cube: + """Return a study's cube on the collection's axes. + + The cube has the study's groups under their aligned ids, merged where several + share one, and the features it measured, both in the collection's order, and + the study's own samples with their own ids. It is derived from the cube as + given with ``Cube.select`` and ``Cube.regroup``, whose records its provenance + holds; a cube that needs neither is returned as it is. A lazy cube, such as + one read with ``read_cube``, stays lazy; a cube in memory is copied. + + Args: + study: A study name. + + Raises: + KeyError: If there is no such study. + TypeError: If ``study`` is not a string. + OverflowError: If a merged value of a cube in memory does not fit (see + "Sums" in ``Cube``). + """ + position = self._study_position(study) + members = self._members[study] + groups = [group for ids in members.values() for group in ids] + features = [ + feature + for feature, measured in zip( + self._features["feature_id"].tolist(), + self._measured[position].tolist(), + strict=True, + ) + if measured + ] + return _on_axes(self._cubes[study], members, groups, features) + + # Deriving collections and cubes. + + def select( + self, studies: Any = None, groups: Any = None, features: Any = None + ) -> "CubeCollection": + """Return a collection of some studies, groups and features. + + Names, ids and positions also set the order of the result; a single one + selects one entry. Without ``features``, the result keeps the features that + one of its studies measured, and without ``groups`` the groups that one of + them has, in the order they had. The tables of samples, groups and features + are those of the studies kept, as a collection of them would have them. The + provenance records the selection. + + Args: + studies: Study names, a boolean mask or positions; ``None`` keeps every + study. + groups: Aligned group ids, a boolean mask or positions; ``None`` keeps + the groups of the studies kept. + features: Feature ids, a boolean mask or positions; ``None`` keeps the + features the studies kept measured. + + Raises: + KeyError: If a name or id is unknown. + IndexError: If a position is out of range. + TypeError: If a selection is not names or ids, a boolean mask or + integers, or is a numpy masked array. + ValueError: If a mask has the wrong length, a selection repeats an entry + or selects no study, or a group or feature selected is in none of the + studies kept. + """ + study_positions = selected_positions( + studies, + text_column(self._studies), + self._study_index, + "studies", + "study", + named=True, + ) + group_positions = selected_positions( + groups, self._groups["group_id"], self._group_index, "groups" + ) + feature_positions = selected_positions( + features, self._features["feature_id"], self._feature_index, "features" + ) + selections = (study_positions, group_positions, feature_positions) + if all(positions is None for positions in selections): + return self + if study_positions is None: + study_positions = np.arange(self.n_studies) + if study_positions.size == 0: + raise ValueError("the selection of studies is empty; select at least one") + kept = [self._studies[position] for position in study_positions.tolist()] + measured = self._measured[study_positions] + group_ids = self._groups["group_id"].tolist() + has = self._has_groups(kept) + group_positions = _kept_positions(group_positions, has, group_ids, "group") + feature_ids = self._features["feature_id"].tolist() + feature_positions = _kept_positions( + feature_positions, measured, feature_ids, "feature" + ) + kept_groups = [group_ids[position] for position in group_positions.tolist()] + kept_features = [ + feature_ids[position] for position in feature_positions.tolist() + ] + record = { + "operation": "select", + "studies": None if studies is None else kept, + "groups": None if groups is None else kept_groups, + "features": None if features is None else kept_features, + } + collection = object.__new__(CubeCollection) + _assemble( + collection, + tuple(kept), + {study: self._cubes[study] for study in kept}, + { + study: { + group: self._members[study][group] + for group in kept_groups + if group in self._members[study] + } + for study in kept + }, + kept_groups, + kept_features, + measured[:, feature_positions], + self._settings, + (*self._operations, json.dumps(record)), + ) + return collection + + def stack( + self, *, groups: Any = None, features: Any = None, layers: Any = None + ) -> Cube: + """Return the studies stacked into one cube, for a joint model. + + The cube's samples are those of every study, named ``study::sample_id``, + with ``samples`` as its sample table, so that a design can take the study + as a covariate. Its groups are aligned groups that every study has, and its + features are features that every study measured, in the collection's order: + a cube holds no missing values, so a group or feature that a study lacks is + refused rather than filled with zeros. Each study's part is its aligned + cube's (``aligned_cube``): its own counts and library sizes, with merged + groups summed. The cube is in memory; of a lazy study cube, such as one read + with ``read_cube``, only the groups and features stacked are read. + + Layers keep their type when every study gives them the same; integer layers + of different types become int64, and integer counts become float64 when + another study's counts are float64. Cell numbers are stacked when every + study records them. The normalization is that of the studies, which must be + the same for a normalized layer to be stacked, and None for counts and + detected cells alone. + + The provenance records, for every study, its cube's fingerprint and + provenance and the groups stacked from it under their aligned ids, with the + settings of the alignment and the layers (``cytepool.cube.provenance``). + + Args: + groups: Aligned group ids, a boolean mask or positions; by default the + groups every study has. + features: Feature ids, a boolean mask or positions; by default the + features every study measured. + layers: Layer names, ``counts`` among them; by default the layers every + study has. + + Returns: + A cube with every study's samples, in the order of the studies and of + each study's samples. + + Raises: + KeyError: If an id is unknown. + IndexError: If a position is out of range. + TypeError: If a selection is not ids, a boolean mask or integers, or is + or holds a numpy masked array. + ValueError: If a group or feature is not in every study, no group or + feature is in every study, a layer is unknown, not in every study or + not given with ``counts``, some studies record cell numbers and + others do not, the studies' normalizations differ for a normalized + layer, or integer counts of 2**53 or more would become float64. + OverflowError: If a merged value does not fit (see "Sums" in ``Cube``). + """ + cubes = [self._cubes[study] for study in self._studies] + group_ids = self._groups["group_id"].tolist() + has = self._has_groups(self._studies) + chosen_groups = _shared_positions( + groups, self._groups["group_id"], self._group_index, has, self, "group" + ) + chosen_features = _shared_positions( + features, + self._features["feature_id"], + self._feature_index, + self._measured, + self, + "feature", + ) + names = _stacked_layers(layers, self._studies, cubes) + with_cells = [cube.n_cells is not None for cube in cubes] + if any(with_cells) and not all(with_cells): + have = [ + s for s, cells in zip(self._studies, with_cells, strict=True) if cells + ] + lack = [s for s in self._studies if s not in have] + raise ValueError( + f"the cubes of {quoted(have)} record the cells of every sample and " + f"group and those of {quoted(lack)} do not: a cube records cell " + "numbers for all its samples or for none, so these studies are not " + "stacked together; select studies alike (select(studies=...))" + ) + normalization = _shared_normalization(names, self._studies, cubes) + stacked_groups = [group_ids[position] for position in chosen_groups.tolist()] + feature_ids = self._features["feature_id"] + stacked_features = [feature_ids[p] for p in chosen_features.tolist()] + parts = [ + _on_axes( + self._cubes[study], + {group: self._members[study][group] for group in stacked_groups}, + [ + member + for group in stacked_groups + for member in self._members[study][group] + ], + stacked_features, + ) + for study in self._studies + ] + stacked = {name: _stacked_layer(name, parts, self._studies) for name in names} + n_cells = None + if all(with_cells): + n_cells = np.concatenate( + [part.n_cells for part in parts if part.n_cells is not None] + ) + settings = json.loads(self._settings) + construction = { + "function": "CubeCollection.stack", + "cytepool_version": __version__, + "studies": { + study: { + "fingerprint": self._cubes[study].fingerprint(), + "provenance": self._cubes[study].provenance, + "groups": { + group: list(self._members[study][group]) + for group in stacked_groups + }, + } + for study in self._studies + }, + "align_groups_by": settings["align_groups_by"], + "within_study": settings["within_study"], + "layers": names, + } + return Cube( + stacked, + samples=self._samples, + groups=self._groups.take(chosen_groups), + features=self._features.take(chosen_features), + n_cells=n_cells, + normalization=normalization, + provenance={"construction": construction}, + copy=False, + ) + + # Internals. + + def _has_groups(self, studies: Sequence[str]) -> np.ndarray: + """Whether each of these studies has each aligned group, as booleans of + shape ``(len(studies), n_groups)``.""" + group_ids = self._groups["group_id"].tolist() + has = [ + [group in self._members[study] for group in group_ids] for study in studies + ] + return np.array(has, dtype=bool).reshape(len(studies), len(group_ids)) + + def _study_position(self, study: object) -> int: + reject_masked(study, "study") + if not isinstance(study, str): + raise TypeError( + f"study must be a study name (str), got {type(study).__name__}" + ) + position = self._study_index.get(study) + if position is None: + raise KeyError( + f"no study named {study!r}; studies are: {quoted(list(self._studies))}" + + close_matches(study, self._studies) + ) + return position + + +def _assemble( + collection: CubeCollection, + studies: tuple[str, ...], + cubes: dict[str, Cube], + members: dict[str, dict[str, tuple[str, ...]]], + group_ids: list[str], + feature_ids: list[str], + measured: np.ndarray, + settings: str, + operations: tuple[str, ...], +) -> None: + """Build a collection's tables from its aligned parts and set its state.""" + within = json.loads(settings)["within_study"] + measured = np.array(measured, dtype=bool) # a copy the collection owns + measured.flags.writeable = False + # Each study's aligned groups in the collection's order. + members = { + study: { + group: members[study][group] + for group in group_ids + if group in members[study] + } + for study in studies + } + _set_state( + collection, + studies, + cubes, + members, + measured, + _sample_table(studies, cubes, within), + _group_table(studies, cubes, members, group_ids), + _feature_table(cubes, feature_ids), + settings, + operations, + ) + + +def _set_state( + collection: CubeCollection, + studies: tuple[str, ...], + cubes: dict[str, Cube], + members: dict[str, dict[str, tuple[str, ...]]], + measured: np.ndarray, + samples: Table, + groups: Table, + features: Table, + settings: str, + operations: tuple[str, ...], +) -> None: + state = { + "_studies": studies, + "_cubes": cubes, + "_members": members, + "_measured": measured, + "_samples": samples, + "_groups": groups, + "_features": features, + "_settings": settings, + "_operations": operations, + } + for axis, ids in ( + ("study", list(studies)), + ("group", groups["group_id"].tolist()), + ("feature", features["feature_id"].tolist()), + ): + state[f"_{axis}_index"] = {key: position for position, key in enumerate(ids)} + for name, value in state.items(): + object.__setattr__(collection, name, value) + + +def _restore_collection( + studies: tuple[str, ...], + cubes: dict[str, Cube], + members: dict[str, dict[str, tuple[str, ...]]], + measured: np.ndarray, + samples: Table, + groups: Table, + features: Table, + settings: str, + operations: tuple[str, ...], +) -> CubeCollection: + collection = object.__new__(CubeCollection) + measured.flags.writeable = False # unpickled arrays are new and writable + _set_state( + collection, + studies, + cubes, + members, + measured, + samples, + groups, + features, + settings, + operations, + ) + return collection + + +# Checking the arguments. + + +def _checked_cubes(cubes: Any) -> tuple[str, ...]: + if not isinstance(cubes, Mapping): + raise TypeError( + "cubes must be a mapping of study names to cubes, such as " + f"{{'cohort_a': cube_a, 'cohort_b': cube_b}}, got {type(cubes).__name__}" + ) + if not cubes: + raise ValueError("cubes is empty; a collection needs at least one study") + for study, cube in cubes.items(): + if not isinstance(study, str): + raise TypeError( + f"study names must be strings, got {type(study).__name__} {study!r}" + ) + if not study: + raise ValueError("study names must be non-empty") + if _SEPARATOR in study: + raise ValueError( + f"the study name {study!r} holds {_SEPARATOR!r}, which separates " + f"a study from a sample id in 'study{_SEPARATOR}sample_id'; " + "rename the study" + ) + if not isinstance(cube, Cube): + raise TypeError( + f"cubes[{study!r}] must be a Cube, got {type(cube).__name__}" + ) + return tuple(cubes) + + +def _checked_within_study(value: Any, cubes: dict[str, Cube]) -> list[str]: + if isinstance(value, str): + names: list[Any] = [value] + else: + reject_masked(value, "within_study") + try: + names = list(value) + except TypeError: + raise TypeError( + "within_study must be a sample column name or a list of them, got " + f"{type(value).__name__}" + ) from None + columns = list( + dict.fromkeys(name for cube in cubes.values() for name in cube.samples.columns) + ) + for name in names: + if not isinstance(name, str): + raise TypeError( + f"within_study must name sample columns (str), got " + f"{type(name).__name__} {name!r}" + ) + if name in ("sample_id", "study"): + raise ValueError( + f"within_study names {name!r}, which the collection writes itself; " + "name other sample columns, such as donor ids" + ) + if names.count(name) > 1: + raise ValueError(f"within_study names {name!r} more than once") + if name not in columns: + raise KeyError( + f"within_study names {name!r}, which is not a sample column of any " + f"study; their columns are: {listing(columns)}" + + close_matches(name, columns) + ) + for study, cube in cubes.items(): + if name in cube.samples: + column = cube.samples[name] + if column.size and column.dtype.kind not in "Oiu": + raise TypeError( + f"within_study names {name!r}, whose values in study " + f"{study!r} are {column.dtype}; such columns hold ids, as text " + "or integers" + ) + return names + + +def _checked_group_mapping( + value: Any, cubes: dict[str, Cube] +) -> dict[str, dict[str, str | None]]: + if value is None: + return {} + if not isinstance(value, Mapping): + raise TypeError( + "group_mapping must map study names to mappings of group ids to aligned " + f"ids, got {type(value).__name__}" + ) + reject_masked(value, "group_mapping") + checked: dict[str, dict[str, str | None]] = {} + for study, groups in value.items(): + if not isinstance(study, str): + raise TypeError( + f"group_mapping must be keyed by study names (str), got " + f"{type(study).__name__} {study!r}" + ) + if study not in cubes: + raise KeyError( + f"group_mapping names the study {study!r}, which is not in cubes; " + f"studies are: {quoted(list(cubes))}" + + close_matches(study, list(cubes)) + ) + if not isinstance(groups, Mapping): + raise TypeError( + f"group_mapping[{study!r}] must map the study's group ids to aligned " + f"ids, got {type(groups).__name__}" + ) + ids = cubes[study].groups["group_id"].tolist() + known = set(ids) + checked[study] = {} + for group, aligned in groups.items(): + if not isinstance(group, str): + raise TypeError( + f"group_mapping[{study!r}] must be keyed by group ids (str), got " + f"{type(group).__name__} {group!r}" + ) + if group not in known: + raise KeyError( + f"group_mapping[{study!r}] names the group {group!r}, which the " + f"cube of study {study!r} does not have; its groups are: " + f"{quoted(ids)}" + close_matches(group, ids) + ) + if aligned is not None and not isinstance(aligned, str): + raise TypeError( + f"group_mapping[{study!r}][{group!r}] must be an aligned id (str), " + f"or None to leave the group out; got {type(aligned).__name__}" + ) + if aligned == "": + raise ValueError( + f"group_mapping[{study!r}][{group!r}] is empty; aligned ids are " + "non-empty (None leaves the group out)" + ) + checked[study][group] = aligned + return checked + + +def _alignment_column(value: Any, cubes: dict[str, Cube]) -> str: + if value is None: + if all(_ONTOLOGY_COLUMN in cube.groups for cube in cubes.values()): + return _ONTOLOGY_COLUMN + return "group_id" + if not isinstance(value, str): + raise TypeError( + f"align_groups_by must be a group column name, got {type(value).__name__}" + ) + columns = list( + dict.fromkeys(name for cube in cubes.values() for name in cube.groups.columns) + ) + if value not in columns: + raise KeyError( + f"align_groups_by names {value!r}, which is not a group column of any " + f"study; their columns are: {listing(columns)}" + + close_matches(value, columns) + ) + for study, cube in cubes.items(): + if value in cube.groups: + column = cube.groups[value] + if column.size and column.dtype.kind != "O": + raise TypeError( + f"align_groups_by names {value!r}, which holds {column.dtype} in " + f"study {study!r}; groups are aligned by the values of a text " + "column" + ) + return value + + +# Aligning. + + +def _aligned_groups( + study: str, cube: Cube, by: str, mapping: dict[str, str | None] +) -> dict[str, list[str]]: + """The study's groups by aligned id, in the order of their first member.""" + ids = cube.groups["group_id"].tolist() + values = cube.groups[by].tolist() if by in cube.groups else None + aligned: dict[str, list[str]] = {} + for position, group in enumerate(ids): + if group in mapping: + key = mapping[group] + if key is None: + continue + elif values is None: + raise ValueError( + f"the groups of study {study!r} are aligned by {by!r}, which its group " + f"table does not have; give the aligned ids of its groups in " + f"group_mapping[{study!r}], or align the groups by a column every " + "study has (align_groups_by)" + ) + else: + key = values[position] + if key is None or key == "": + raise ValueError( + f"group {group!r} of study {study!r} has no {by}, by which the " + f"groups are aligned; give its aligned id in " + f"group_mapping[{study!r}], or None to leave it out" + ) + aligned.setdefault(key, []).append(group) + return aligned + + +def _measured_features(cubes: dict[str, Cube]) -> tuple[list[str], np.ndarray]: + """The features that a study measured, in order of first appearance, and + whether each study measured each.""" + index: dict[str, int] = {} + measured_ids = [] + for cube in cubes.values(): + ids = cube.features["feature_id"].tolist() + if "measured" in cube.features: + flags = cube.features["measured"].tolist() + ids = [feature for feature, flag in zip(ids, flags, strict=True) if flag] + measured_ids.append(ids) + for feature in ids: + index.setdefault(feature, len(index)) + measured = np.zeros((len(cubes), len(index)), dtype=bool) + for row, ids in enumerate(measured_ids): + measured[row, [index[feature] for feature in ids]] = True + return list(index), measured + + +def _check_shared_features( + studies: tuple[str, ...], feature_ids: list[str], measured: np.ndarray +) -> None: + counts = measured.sum(axis=0) + for row, study in enumerate(studies): + if (measured[row] & (counts > 1)).any(): + continue + own = [f for f, flag in zip(feature_ids, measured[row], strict=True) if flag] + elsewhere = np.delete(measured, row, axis=0).any(axis=0) + others = [f for f, flag in zip(feature_ids, elsewhere, strict=True) if flag] + hint = "" + if {_VERSION.sub("", f) for f in own} & {_VERSION.sub("", f) for f in others}: + hint = ( + " Some ids differ only in a version, such as '.17' in " + "'ENSG00000141510.17': remove the versions from the feature ids first." + ) + raise ValueError( + f"study {study!r} shares no feature with the other studies: none of the " + f"{len(own):,} features it measured, such as {quoted(own[:3])}, was " + f"measured by another study, which measured features such as " + f"{quoted(others[:3])}. Features are aligned by feature_id, which must " + "hold ids of one kind in every study, such as Ensembl gene ids." + hint + ) + + +def _check_shared_groups( + studies: tuple[str, ...], members: dict[str, dict[str, list[str]]], by: str +) -> None: + for study in studies: + others = { + group for other in studies if other != study for group in members[other] + } + if set(members[study]) & others: + continue + raise ValueError( + f"study {study!r} shares no group with the other studies: aligned by " + f"{by!r}, its groups are {quoted(list(members[study]))}, and theirs " + f"{quoted(sorted(others))}. Align the groups by another column " + "(align_groups_by), or give their aligned ids in group_mapping" + ) + + +def _on_axes( + cube: Cube, + members: Mapping[str, Sequence[str]], + groups: list[str], + features: list[str], +) -> Cube: + """A study's cube with the groups ``groups``, merged and renamed by ``members`` + (aligned id to the groups it sums, in the order of ``groups``), and the features + ``features``; selecting an axis only where it changes.""" + selected = cube.select( + groups=None if groups == cube.groups["group_id"].tolist() else groups, + features=None if features == cube.features["feature_id"].tolist() else features, + ) + renamed = { + aligned: list(ids) for aligned, ids in members.items() if list(ids) != [aligned] + } + return selected.regroup(renamed) if renamed else selected + + +# The tables of a collection. + + +def _sample_table( + studies: tuple[str, ...], cubes: dict[str, Cube], within: list[str] +) -> Table: + ids: list[str] = [] + labels: list[str] = [] + for study in studies: + sample_ids = cubes[study].samples["sample_id"].tolist() + ids += [f"{study}{_SEPARATOR}{sample}" for sample in sample_ids] + labels += [study] * len(sample_ids) + if "study" in cubes[study].samples: + found = cubes[study].samples["study"].tolist() + if any(value != study for value in found): + other = next(value for value in found if value != study) + raise ValueError( + f"the cube of study {study!r} has a sample column 'study' with the " + f"value {other!r}, where the collection writes the study's name; " + "rename that column" + ) + names = list( + dict.fromkeys( + name + for study in studies + for name in cubes[study].samples.columns + if name not in ("sample_id", "study") + ) + ) + columns: dict[str, np.ndarray] = { + "sample_id": text_column(ids), + "study": text_column(labels), + } + for name in names: + parts = [] + for study in studies: + table = cubes[study].samples + column = table[name] if name in table else None + if column is not None and name in within: + column = text_column( + [ + None if value is None else f"{study}{_SEPARATOR}{value}" + for value in column.tolist() + ] + ) + parts.append((study, column, len(table))) + columns[name] = _stacked_column(name, "sample", parts) + return Table(columns) + + +def _group_table( + studies: tuple[str, ...], + cubes: dict[str, Cube], + members: dict[str, dict[str, tuple[str, ...]]], + group_ids: list[str], +) -> Table: + names = list( + dict.fromkeys( + name + for study in studies + for name in cubes[study].groups.columns + if name != "group_id" + ) + ) + columns: dict[str, np.ndarray] = {"group_id": text_column(group_ids)} + # The rows of each study's group table that each aligned group sums. + rows = { + study: { + group: _positions_of(cubes[study].groups["group_id"], ids) + for group, ids in members[study].items() + } + for study in studies + } + for name in names: + present = [ + (study, cubes[study].groups[name]) + for study in studies + if name in cubes[study].groups + ] + dtype = _common_dtype(name, "group", present) + values: list[Any] = [] + for group in group_ids: + found = [ + column[rows[study][group]].astype(dtype) + for study, column in present + if group in rows[study] + ] + if not found: + values.append(_MISSING) + continue + agreed = np.concatenate(found) + same = values_equal(agreed, np.repeat(agreed[:1], agreed.size)).all() + values.append(agreed[0] if same else _MISSING) + columns[name] = _column_from_values(values, dtype) + return Table(columns) + + +def _feature_table(cubes: dict[str, Cube], feature_ids: list[str]) -> Table: + columns = {"feature_id": text_column(feature_ids)} + if any("symbol" in cube.features for cube in cubes.values()): + index = {feature: position for position, feature in enumerate(feature_ids)} + symbols: list[str | None] = [None] * len(feature_ids) + for cube in cubes.values(): + if "symbol" not in cube.features: + continue + for feature, symbol in zip( + cube.features["feature_id"].tolist(), + cube.features["symbol"].tolist(), + strict=True, + ): + position = index.get(feature) + if position is not None and symbols[position] is None: + symbols[position] = symbol + columns["symbol"] = text_column(symbols) + return Table(columns) + + +def _positions_of(ids: np.ndarray, chosen: Sequence[str]) -> np.ndarray: + where = {identifier: position for position, identifier in enumerate(ids.tolist())} + return np.array([where[identifier] for identifier in chosen], dtype=np.intp) + + +class _Missing: + """A missing value of a column that is built value by value.""" + + +_MISSING = _Missing() + + +def _common_dtype( + name: str, axis: str, present: list[tuple[str, np.ndarray]] +) -> np.dtype: + """The dtype that holds a table column of every study that has it. + + Columns without rows do not count. Raises TypeError for columns whose kinds of + values differ, such as text in one study and numbers in another. + """ + filled = [(study, column) for study, column in present if column.size] + if not filled: + return present[0][1].dtype if present else np.dtype(object) + dtypes = [column.dtype for _, column in filled] + if all(dtype == dtypes[0] for dtype in dtypes): + return dtypes[0] + kinds = {_KINDS[dtype.kind] for dtype in dtypes} + if kinds == {"integers"}: + for study, column in filled: + if column.dtype == np.uint64 and column.max() > INT64_MAX: + raise TypeError( + f"the {axis} column {name!r} of study {study!r} holds integers " + "above 2**63 - 1, which int64 does not hold, and other studies " + "hold signed integers in it" + ) + return np.dtype(np.int64) + if kinds <= {"integers", "floats"}: + return np.dtype(np.float64) + if kinds <= _NUMBERS: + return np.dtype(np.complex128) + if kinds in ({"datetimes"}, {"timedeltas"}): + return np.result_type(*dtypes) + first_study, first = filled[0] + other_study, other = next( + (study, column) + for study, column in filled + if _KINDS[column.dtype.kind] != _KINDS[first.dtype.kind] + ) + raise TypeError( + f"the {axis} column {name!r} holds {_KINDS[first.dtype.kind]} " + f"({first.dtype}) in study {first_study!r} and {_KINDS[other.dtype.kind]} " + f"({other.dtype}) in study {other_study!r}; a column must hold one kind of " + "values in every study: convert it in one of the cubes' tables, or rename it " + "in one of them" + ) + + +def _with_missing_dtype(dtype: np.dtype) -> np.dtype: + """The dtype of a column that needs a missing value: integers and booleans + become float64 with NaN.""" + return np.dtype(np.float64) if dtype.kind in "biu" else dtype + + +def _missing_value(dtype: np.dtype) -> Any: + if dtype.kind == "O": + return None + if dtype.kind in "mM": + return np.array("NaT", dtype=dtype) + return np.nan + + +def _stacked_column( + name: str, axis: str, parts: list[tuple[str, np.ndarray | None, int]] +) -> np.ndarray: + """A table column of every study, one after another; missing values where a + study does not have it.""" + present = [(study, column) for study, column, _ in parts if column is not None] + dtype = _common_dtype(name, axis, present) + if any(column is None and rows for _, column, rows in parts): + dtype = _with_missing_dtype(dtype) + result = np.empty(sum(rows for _, _, rows in parts), dtype=dtype) + start = 0 + for _, column, rows in parts: + if column is None: + result[start : start + rows] = _missing_value(dtype) + elif rows: + result[start : start + rows] = column.astype(dtype) + start += rows + return result + + +def _column_from_values(values: list[Any], dtype: np.dtype) -> np.ndarray: + """A column of the values given, missing where a value is ``_MISSING``.""" + missing = [value is _MISSING for value in values] + if any(missing): + dtype = _with_missing_dtype(dtype) + result = np.empty(len(values), dtype=dtype) + for position, value in enumerate(values): + result[position] = _missing_value(dtype) if missing[position] else value + return result + + +# Stacking. + + +def _kept_positions( + positions: np.ndarray | None, present: np.ndarray, ids: list[str], single: str +) -> np.ndarray: + """The groups or features a selection keeps: those given, each of which a kept + study must have, or by default those a kept study has, in order.""" + if positions is None: + return np.flatnonzero(present.any(axis=0)) + lacking = [ids[p] for p in positions.tolist() if not present[:, p].any()] + if lacking: + word = "measured" if single == "feature" else "has" + raise ValueError( + f"no study selected {word} the {single}" + f"{'' if len(lacking) == 1 else 's'} {quoted(lacking)}; select studies " + f"that do, or other {single}s" + ) + return positions + + +def _shared_positions( + selector: Any, + ids: np.ndarray, + index: dict[str, int], + present: np.ndarray, + collection: CubeCollection, + single: str, +) -> np.ndarray: + """The groups or features to stack: those given, each of which every study must + have, or by default those every study has, in the collection's order.""" + word = "measured" if single == "feature" else "has" + positions = selected_positions(selector, ids, index, f"{single}s") + if positions is None: + shared = np.flatnonzero(present.all(axis=0)) + if shared.size == 0: + raise ValueError( + f"no {single} is in every study, so there is nothing to stack: a cube " + f"holds no missing values, so stack() takes only the {single}s that " + f"every study {word}; select studies that share {single}s " + "(select(studies=...))" + ) + return shared + for position in positions.tolist(): + lacking = [ + study + for row, study in enumerate(collection.studies) + if not present[row, position] + ] + if lacking: + raise ValueError( + f"the {single} {ids[position]!r} is not in " + f"{'study' if len(lacking) == 1 else 'studies'} {quoted(lacking)}: a " + f"cube holds no missing values, so stack() takes only {single}s that " + f"every study {word}, rather than fill the others with zeros; leave " + f"out the {single} or those studies (select(studies=...))" + ) + return positions + + +def _stacked_layers( + value: Any, studies: tuple[str, ...], cubes: list[Cube] +) -> list[str]: + if value is None: + return [ + name + for name in LAYER_NAMES + if all(name in cube.layer_names for cube in cubes) + ] + reject_masked(value, "layers") + names: list[Any] = [value] if isinstance(value, str) else list(value) + for name in names: + if not isinstance(name, str): + raise TypeError( + f"layers must be layer names (str), got {type(name).__name__} {name!r}" + ) + if name not in LAYER_NAMES: + raise ValueError( + f"unknown layer {name!r}; layers are named: {', '.join(LAYER_NAMES)}" + ) + if names.count(name) > 1: + raise ValueError(f"layers names {name!r} more than once") + lacking = [ + study + for study, cube in zip(studies, cubes, strict=True) + if name not in cube.layer_names + ] + if lacking: + shared = [ + layer + for layer in LAYER_NAMES + if all(layer in cube.layer_names for cube in cubes) + ] + raise ValueError( + f"the cube of {'study' if len(lacking) == 1 else 'studies'} " + f"{quoted(lacking)} has no {name!r} layer; every study has " + f"{quoted(shared)}" + ) + if "counts" not in names: + raise ValueError("layers must include 'counts', which every cube holds") + return [name for name in LAYER_NAMES if name in names] + + +def _shared_normalization( + names: list[str], studies: tuple[str, ...], cubes: list[Cube] +) -> dict[str, Any] | None: + normalized = [name for name in names if name in _NORMALIZED_LAYERS] + if not normalized: + return None + first = cubes[0].normalization + for study, cube in zip(studies[1:], cubes[1:], strict=True): + if cube.normalization != first: + raise ValueError( + f"the normalized layers {quoted(normalized)} of study {study!r} were " + f"computed with the normalization {cube.normalization!r}, and those of " + f"study {studies[0]!r} with {first!r}: stacked, they would mix two " + "normalizations. Stack counts and detected cells alone " + "(layers=['counts', 'detected']), or studies normalized alike" + ) + return first + + +def _stacked_layer( + name: str, parts: list[Cube], studies: tuple[str, ...] +) -> np.ndarray: + """A layer of every study's part, with the samples of each after another.""" + n_groups, n_features = parts[0].n_groups, parts[0].n_features + n_samples = sum(part.n_samples for part in parts) + dtype: np.dtype + if name in _NORMALIZED_LAYERS: + dtype = np.dtype(np.float64) + else: + dtypes = ( + [part.layer(name, 0).dtype for part in parts] + if n_groups + else [np.dtype(np.int64)] + ) + if any(dtype.kind == "f" for dtype in dtypes): + dtype = np.dtype(np.float64) + elif all(dtype == dtypes[0] for dtype in dtypes): + dtype = dtypes[0] + else: + dtype = np.dtype(np.int64) + result = np.empty((n_groups, n_samples, n_features), dtype=dtype) + start = 0 + for study, part in zip(studies, parts, strict=True): + block = part.layer(name) + if ( + dtype.kind == "f" + and block.dtype.kind in "iu" + and block.size + and block.max() >= _FLOAT_LIMIT + ): + raise ValueError( + f"the integer {name} of study {study!r} reach {block.max()}, and float64, " + "which another study's counts need, holds whole numbers exactly only " + "below 2**53" + ) + result[:, start : start + part.n_samples] = block + start += part.n_samples + del block + return result diff --git a/cytepool/cube/cube.py b/cytepool/cube/cube.py index 73c1ad5..24742c9 100644 --- a/cytepool/cube/cube.py +++ b/cytepool/cube/cube.py @@ -10,7 +10,13 @@ from .._columns import values_equal, with_missing from .._masked_arrays import reject_masked, reject_masked_table, unmasked_array from ..table import Table -from ._axes import INT64_MAX, checked_axis_table, checked_n_cells, first_index +from ._axes import ( + INT64_MAX, + checked_axis_table, + checked_n_cells, + first_index, + selected_positions, +) from ._provenance_records import ( add_record, plain_json, @@ -720,13 +726,13 @@ def select( ValueError: If a mask has the wrong length or a selection repeats an entry. """ - sample_positions = _positions( + sample_positions = selected_positions( samples, self._samples["sample_id"], self._sample_index, "samples" ) - group_positions = _positions( + group_positions = selected_positions( groups, self._groups["group_id"], self._group_index, "groups" ) - feature_positions = _positions( + feature_positions = selected_positions( features, self._features["feature_id"], self._feature_index, "features" ) selections = (sample_positions, group_positions, feature_positions) @@ -2266,66 +2272,6 @@ def _check_moment_layers( _check_moments(numerator, squares[group], group, ids) -def _positions( - selector: Any, ids: np.ndarray, index: dict[str, int], axis: str -) -> np.ndarray | None: - """Convert ids, a boolean mask or positions along one axis to positions.""" - if selector is None: - return None - n = len(index) - single = axis[:-1] - 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} 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 _take_rows(table: Table, positions: np.ndarray | None) -> Table: return table if positions is None else table.take(positions) diff --git a/cytepool/cube/provenance.py b/cytepool/cube/provenance.py index 28d378d..da81dd1 100644 --- a/cytepool/cube/provenance.py +++ b/cytepool/cube/provenance.py @@ -18,8 +18,8 @@ ``construction``: How the cube was built: ``function`` (``"from_cells"``, ``"from_bulk"``, - ``"from_anndata"`` or ``"from_mudata"``), ``cytepool_version``, and the - function's settings. ``from_cells`` records + ``"from_anndata"``, ``"from_mudata"`` or ``"CubeCollection.stack"``), + ``cytepool_version``, and the function's settings. ``from_cells`` records - ``n_cells``: the number of cells pooled (the cells of each sample and group are the cube's ``n_cells``); @@ -61,6 +61,18 @@ ``sample_column`` and ``group_column``; ``from_mudata`` lists the columns of the MuData's ``obs`` placed in the sample table. + ``CubeCollection.stack`` (``cytepool.cube.collection``), which stacks the cubes + of several studies, records + + - ``studies``: for each study, in the order stacked, the cube that the + collection was given, by its ``fingerprint``, with that cube's own + ``provenance``, and ``groups``: each group stacked, by its aligned id, with + the ids of the study's groups that it sums; + - ``align_groups_by``: the group column whose values aligned the groups; + - ``within_study``: the sample columns whose values were written + ``study::value``, as the sample ids are; + - ``layers``: the layers stacked. + ``operations``: The operations that derived the cube, oldest first. Each is an object with ``operation``, the method, and ``input``, the fingerprint of the cube it was @@ -174,6 +186,7 @@ def describe_provenance(provenance: Mapping[str, Any]) -> str: references = isinstance(record["source"], dict) if "construction" in record: lines += _construction_lines(record["construction"]) + references |= _stacked_studies(record["construction"]) is not None for number, operation in enumerate(record.get("operations", []), start=1): lines += _operation_lines(number, operation) references |= isinstance(operation, dict) and operation.get("input") is not None @@ -220,7 +233,24 @@ def _construction_lines(construction: Any) -> list[str]: if layout == "given" else f", in the layout of {layout}" ) + studies = _stacked_studies(construction) + if studies is not None: + header += f" from {len(studies):,} stud{'y' if len(studies) == 1 else 'ies'}" lines = [header + ":"] + if studies is not None: + del settings["studies"] + lines += _study_lines(studies) + by = settings.get("align_groups_by") + if isinstance(by, str): + del settings["align_groups_by"] + lines.append(f" groups aligned by the values of {by!r}") + within = settings.get("within_study") + if _is_text_list(within): + del settings["within_study"] + if within: + lines.append( + f" sample columns written study::value: {', '.join(within)}" + ) for axis in ("samples", "groups"): value = settings.get(axis) if ( @@ -264,6 +294,40 @@ def _construction_lines(construction: Any) -> list[str]: return lines +def _stacked_studies(construction: Any) -> dict[str, Any] | None: + """The studies of a cube that ``CubeCollection.stack`` built, else None.""" + if not isinstance(construction, dict): + return None + studies = construction.get("studies") + stacked = construction.get("function") == "CubeCollection.stack" + return studies if stacked and isinstance(studies, dict) else None + + +def _study_lines(studies: dict[str, Any]) -> list[str]: + """The studies of a stacked cube: each cube by its fingerprint, and the groups + that were merged or renamed.""" + lines = [] + for name, study in studies.items(): + if not isinstance(study, dict): + lines.append(f" study {name!r}: {_json(study)}") + continue + found = study.get("fingerprint") + named = found if isinstance(found, str) else _json(found) + lines.append(f" study {name!r}: the cube {named}, with its own record") + groups = study.get("groups") + if isinstance(groups, dict) and _is_membership(groups): + changed = { + aligned: members + for aligned, members in groups.items() + if members != [aligned] + } + for aligned, members in list(changed.items())[:_SHOWN_MEMBERS]: + lines.append(f" {aligned}: {_listing(members, ' + ')}") + if len(changed) > _SHOWN_MEMBERS: + lines.append(f" ... ({len(changed):,} in all)") + return lines + + def _operation_lines(number: int, operation: Any) -> list[str]: if not isinstance(operation, dict): return [f"{number}. {_json(operation)}"] diff --git a/cytepool/meta_analysis.py b/cytepool/meta_analysis.py index 64f2adb..4309ef5 100644 --- a/cytepool/meta_analysis.py +++ b/cytepool/meta_analysis.py @@ -6,7 +6,11 @@ ``rma.uni`` does for one feature without moderators (Viechtbauer 2010). :func:`differential_expression_meta_analysis` does it for the differential-expression results of several cubes, such as the same cell types in several cohorts, tissues or -platforms, matched by group, contrast and feature id. +platforms, matched by group, contrast and feature id. :func:`align_study_results` +lines up the results of the studies of a :class:`~cytepool.cube.CubeCollection` for +one group and contrast, feature by feature, with the reason a study has no result for +a feature: it did not measure it, did not test it, or did not analyse the group or +contrast. The model: Study i estimates the effect of a feature by y_i, with standard error s_i, and its @@ -163,6 +167,7 @@ from ._masked_arrays import read_only, reject_masked, unmasked_array from ._messages import PACKAGE_DIRECTORY, listing, quoted from ._options import checked_choice, index_integer, real_number +from .cube import CubeCollection from .differential_expression import DifferentialExpressionResult from .stats.distributions import ( chi2_sf_scalar, @@ -177,6 +182,7 @@ __all__ = [ "MetaAnalysisFit", "MetaAnalysisResult", + "align_study_results", "differential_expression_meta_analysis", "meta_analysis", ] @@ -211,6 +217,11 @@ # of tau^2, to bracket every maximum of the restricted likelihood. _GRID_PER_DECADE = 8 _ANALYSED = "analysed" # the status of a group and contrast that was meta-analysed +# Why a study has no result for a feature, as align_study_results gives it. +_NOT_TESTED = ( + "the study measured the feature and tested the group and contrast, but not the " + "feature, as when its expression filter left the feature out" +) @dataclass(frozen=True, eq=False, repr=False) @@ -1195,6 +1206,205 @@ def differential_expression_meta_analysis( ) +def align_study_results( + collection: CubeCollection, + results: Mapping[str, DifferentialExpressionResult], + group: str, + contrast: str, +) -> Table: + """Line up the differential expression of every study, feature by feature. + + Each study's result is the analysis of its aligned cube + (:meth:`~cytepool.cube.CubeCollection.aligned_cube`), such as + :func:`~cytepool.differential_expression.limma_voom` with the study's own design, + so that its groups are the collection's aligned groups and its features those + the study measured. For one group and contrast, the table has a row for every + study and every feature of the collection: the studies in order, and the + features in the collection's order within each study, so that + ``table[column].reshape(collection.n_studies, collection.n_features)`` gives a + column with a row per study. ``status`` says whether the study tested the + feature, and why not: + + - ``"tested"``: the study's table of the group and contrast has the feature; + - ``"not_tested"``: the study measured the feature and has a table of the group + and contrast, without the feature, as when its expression filter left it out; + - ``"not_measured"``: the study did not measure the feature + (:attr:`~cytepool.cube.CubeCollection.measured`), so no analysis of it could + test it; + - ``"not_analysed"``: the study measured the feature but has no table of the + group and contrast: it does not have the group, its result leaves the group + out, the group was not fitted, or the contrast cannot be estimated from its + samples; ``reason`` says which, in the words of the study's result. + + The columns are ``study``, ``feature_id``, ``symbol`` (when the collection's + features have it), ``status``, ``reason`` (empty where tested), and then the + numeric columns of the studies' tables, such as ``log2_fold_change``, + ``standard_error``, ``t``, ``p_value`` and ``adjusted_p_value``, in the order in + which the studies first have them, NaN where a study did not test the feature. + A missing feature is thus never a zero, and a feature that a study did not + measure is told apart from one that it measured and left out. + + Args: + collection: The collection whose aligned cubes the studies' results + analysed. + results: The differential-expression result of every study of the + collection, by study name, from + :func:`~cytepool.differential_expression.limma_voom`, + :func:`~cytepool.differential_expression.limma_trend` or + :func:`~cytepool.differential_expression.dream_voom`. + group: An aligned group id of the collection. + contrast: A contrast name, or the name of an F-test of + :func:`~cytepool.differential_expression.dream_voom`, that a study + tested. + + Returns: + A table of ``collection.n_studies * collection.n_features`` rows. + + Raises: + TypeError: If ``collection`` is not a ``CubeCollection``, ``results`` not a + mapping of study names to ``DifferentialExpressionResult``, ``group`` or + ``contrast`` not a string, or ``group`` or ``contrast`` is a numpy masked + array. + ValueError: If ``results`` does not give the result of every study of the + collection and no other, or a study's result analyses a group that the + study does not have in the collection or tests a feature that it did not + measure there: it did not analyse the study's aligned cube. + KeyError: If ``group`` is not a group of the collection, or no study tested + ``contrast``. + """ + if not isinstance(collection, CubeCollection): + raise TypeError( + f"collection must be a CubeCollection, got {type(collection).__name__}" + ) + studies = collection.studies + if not isinstance(results, Mapping): + raise TypeError( + "results must be a mapping of study names to DifferentialExpressionResult, " + f"one for every study of the collection, got {type(results).__name__}" + ) + missing = [study for study in studies if study not in results] + extra = [name for name in results if name not in studies] + if missing or extra: + raise ValueError( + "results must give the result of every study of the collection, and no " + "other: " + + "; ".join( + part + for part in ( + f"missing {quoted(missing)}" if missing else "", + f"not in the collection: {quoted(extra)}" if extra else "", + ) + if part + ) + ) + for study in studies: + if not isinstance(results[study], DifferentialExpressionResult): + raise TypeError( + f"results[{study!r}] must be a DifferentialExpressionResult, from " + f"limma_voom, limma_trend or dream_voom, got " + f"{type(results[study]).__name__}" + ) + for name, value in (("group", group), ("contrast", contrast)): + reject_masked(value, name) + if not isinstance(value, str): + raise TypeError(f"{name} must be a string, got {type(value).__name__}") + group_ids = collection.groups["group_id"].tolist() + if group not in group_ids: + raise KeyError( + f"no group {group!r} in the collection; its groups are: " + f"{listing(group_ids) or '(none)'}" + ) + tested_names = list( + dict.fromkeys( + name + for study in studies + for name in (*results[study].contrasts.names, *results[study].f_tests) + ) + ) + if contrast not in tested_names: + raise KeyError( + f"no study tested the contrast {contrast!r}; the studies tested: " + f"{listing(tested_names) or '(none)'}" + ) + members = collection.group_members + held: dict[str, set[str]] = {study: set() for study in studies} + for study, aligned in zip( + members["study"].tolist(), members["group_id"].tolist(), strict=True + ): + held[study].add(aligned) + feature_ids = collection.features["feature_id"].tolist() + index = {feature: position for position, feature in enumerate(feature_ids)} + n_features = len(feature_ids) + measured = collection.measured + statuses: list[np.ndarray] = [] + reasons: list[list[str]] = [] + values: dict[str, list[tuple[int, np.ndarray, np.ndarray]]] = {} + for row, study in enumerate(studies): + result = results[study] + foreign = [name for name in result.groups if name not in held[study]] + if foreign: + raise ValueError( + f"the result of study {study!r} analyses the group" + f"{'' if len(foreign) == 1 else 's'} {quoted(foreign)}, which the " + "study does not have in this collection; analyse the study's aligned " + f"cube, collection.aligned_cube({study!r})" + ) + table = None + if group not in held[study]: + absence = f"study {study!r} has no group {group!r}" + else: + table = result.tables.get((group, contrast)) + absence = "" if table is not None else _absence(result, group, contrast) + status = np.where(measured[row], "not_analysed", "not_measured").astype(object) + reason = [ + absence if flag else f"study {study!r} did not measure the feature" + for flag in measured[row].tolist() + ] + if table is not None: + positions = np.array( + [index.get(feature, -1) for feature in table["feature_id"].tolist()], + dtype=np.intp, + ) + outside = (positions < 0) | ~measured[row, np.maximum(positions, 0)] + if outside.any(): + feature = table["feature_id"][int(np.flatnonzero(outside)[0])] + raise ValueError( + f"the result of study {study!r} tests the feature {feature!r}, " + "which the study did not measure in this collection; analyse the " + f"study's aligned cube, collection.aligned_cube({study!r})" + ) + status[measured[row]] = "not_tested" + status[positions] = "tested" + for position in np.flatnonzero(measured[row]).tolist(): + reason[position] = _NOT_TESTED + for position in positions.tolist(): + reason[position] = "" + for name in table.columns: + column = table[name] + if name not in ("feature_id", "symbol") and column.dtype.kind in "biuf": + values.setdefault(name, []).append((row, positions, column)) + statuses.append(status) + reasons.append(reason) + columns: dict[str, Any] = { + "study": text_column([study for study in studies for _ in range(n_features)]), + "feature_id": text_column(feature_ids * len(studies)), + } + if "symbol" in collection.features: + columns["symbol"] = text_column( + collection.features["symbol"].tolist() * len(studies) + ) + columns["status"] = text_column( + [status for array in statuses for status in array.tolist()] + ) + columns["reason"] = text_column([text for part in reasons for text in part]) + for name, parts in values.items(): + column = np.full((len(studies), n_features), np.nan) + for row, positions, found in parts: + column[row, positions] = found + columns[name] = column.reshape(-1) + return Table(columns) + + def _meta_analysed_table( tested: Mapping[str, Table], method: str, test: str, level: float, fewest: int ) -> tuple[Table, int, list[str]]: diff --git a/docs/coverage_exemptions.toml b/docs/coverage_exemptions.toml index 297ff55..b6cd21f 100644 --- a/docs/coverage_exemptions.toml +++ b/docs/coverage_exemptions.toml @@ -69,6 +69,7 @@ "cytepool.plotting.plot_composition" = "guides/composition: cell-type proportions per sample" "cytepool.plotting.plot_dispersions" = "guides/bulk_rna_seq: DESeq2's dispersion estimates" "cytepool.plotting.plot_forest" = "guides/combine_studies: one gene across studies" +"cytepool.cube.CubeCollection.select" = "guides/combine_studies: the studies, cell types and genes of an analysis" "cytepool.plotting.plot_gene_sets" = "guides/gene_sets: the results of camera, fry and mroast" "cytepool.plotting.plot_md" = "guides/genes_per_cell_type: MD plot of a contrast" "cytepool.plotting.plot_mds" = "quickstart and guides/genes_per_cell_type: samples on an MDS plot" diff --git a/docs/source/index.md b/docs/source/index.md index 7dbedb7..590fcf6 100644 --- a/docs/source/index.md +++ b/docs/source/index.md @@ -36,12 +36,12 @@ what CytePool offers for each today. | Which cell types expand or contract? | propeller. Not yet: precision-weighted mixed models or count models of composition | [Composition tests](statistics/composition.md) | | Which pathways and regulators respond? | camera, fry, roast, mroast and romer on any contrast; gene set enrichment analysis (GSEA) of ranked genes; over-representation analysis of the genes that change; gene sets read from GMT files. Not yet: activity scores of pathways and transcription factors | [Gene-set tests](statistics/gene_set_tests.md) | | What drives the variation between donors, and what must be adjusted for? | Design checks for confounded factors, sample quality weights, mixed models that estimate a variance per random effect for every gene, and variance partitioning: the share of each gene's variance that each variable explains, in every cell type | [Designs, linear models and contrasts](statistics/linear_models.md), [Mixed models, REML and dream](statistics/mixed_models.md) | -| Does it replicate, and how should studies be combined? | Meta-analysis of every gene across studies, with the between-study variance, Q and I². Not yet: replication measures and shrinkage across studies | [Meta-analysis](statistics/meta_analysis.md) | +| Does it replicate, and how should studies be combined? | Studies aligned by sample, gene and cell type, each analysed on its own and lined up gene by gene; meta-analysis of every gene across studies, with the between-study variance, Q and I². Not yet: replication measures and shrinkage across studies | [Meta-analysis](statistics/meta_analysis.md), [The cube and its additive layers](statistics/cubes.md) | | Can patients be stratified from their profiles? | Not yet; cubes convert to AnnData and MuData for other tools | [Capability matrix](reference/capabilities.md) | | Which cell-cell interactions change? | Not yet | [Related software](reference/related_software.md) | | Do effects depend on genotype, time, sex or context? | Interactions and F-tests in designs, likelihood ratio tests, random intercepts for repeated measures, with F-tests and small-sample (Kenward-Roger) tests in mixed models. Not yet: random slopes. QTL mapping is left to dedicated tools | [Designs, linear models and contrasts](statistics/linear_models.md), [Count models](statistics/count_models.md), [Mixed models, REML and dream](statistics/mixed_models.md) | | How do cells respond to drugs and doses? | A dose as a factor or a covariate of a design. Not yet: dose-response fits, benchmark doses, synergy | [Capability matrix](reference/capabilities.md) | -| Can public studies be pooled into one cube? | Not yet; cubes import the pseudobulk layouts of scanpy, decoupler and LIANA, and meta-analysis combines results across cubes | [Export a cube from another tool](guides/export_cubes.ipynb), [Meta-analysis](statistics/meta_analysis.md) | +| Can public studies be pooled into one cube? | The cubes of several studies stacked into one, of the genes and cell types that every study has, with sample and donor ids kept apart by study; cubes import the pseudobulk layouts of scanpy, decoupler and LIANA | [The cube and its additive layers](statistics/cubes.md), [Export a cube from another tool](guides/export_cubes.ipynb) | The [Literature](reference/literature.md) page collects studies behind these questions, and the [Capability matrix](reference/capabilities.md) lists every method with the diff --git a/docs/source/reference/api/cubes.md b/docs/source/reference/api/cubes.md index 29eb1f9..0ef3606 100644 --- a/docs/source/reference/api/cubes.md +++ b/docs/source/reference/api/cubes.md @@ -1,7 +1,7 @@ --- myst: html_meta: - description: Build pseudobulk cubes from cell-level or bulk counts with CytePool, keep sample, group and feature metadata in tables, and record what each cube pooled with a provenance record and a content fingerprint. + description: Build pseudobulk cubes from cell-level or bulk counts with CytePool, keep sample, group and feature metadata in tables, record what each cube pooled with a provenance record and a content fingerprint, and align the cubes of several studies by sample, gene and cell type. --- # Pseudobulk cubes and tables @@ -13,6 +13,9 @@ samples pooled without going back to the cells. {func}`~cytepool.cube.from_cells pools a cell-level count matrix in one pass, and {func}`~cytepool.cube.from_bulk` wraps a samples x features count matrix as a cube with one group. Sample, group and feature metadata are {class}`~cytepool.table.Table` objects. +{class}`~cytepool.cube.CubeCollection` aligns the cubes of several studies by sample, +gene and cell type, with what each study measured, for analyses of each study and for +one cube of them all. ```{eval-rst} .. automodule:: cytepool.cube diff --git a/docs/source/reference/api/meta_analysis.md b/docs/source/reference/api/meta_analysis.md index 28bda48..07a9785 100644 --- a/docs/source/reference/api/meta_analysis.md +++ b/docs/source/reference/api/meta_analysis.md @@ -1,7 +1,7 @@ --- myst: html_meta: - description: Meta-analysis of the effects of many genes across studies in CytePool, by fixed and random effects with the between-study variance by DerSimonian-Laird, REML or Paule-Mandel, Cochran's Q, I-squared, H-squared and the Knapp-Hartung adjustment, as metafor computes them, and of the differential expression of several cubes. + description: Meta-analysis of the effects of many genes across studies in CytePool, by fixed and random effects with the between-study variance by DerSimonian-Laird, REML or Paule-Mandel, Cochran's Q, I-squared, H-squared and the Knapp-Hartung adjustment, as metafor computes them, and of the differential expression of several cubes, with each study's results lined up gene by gene. --- # Meta-analysis @@ -20,6 +20,9 @@ same cell types in several cohorts: it matches groups, contrasts and genes by th reads each table's log2 fold changes and moderated standard errors, and returns a table per cell type and contrast with p-values adjusted over its genes, the studies used and left out with the reason, and the provenance of every study. +{func}`~cytepool.meta_analysis.align_study_results` lines up the results of the studies +of a {class}`~cytepool.cube.CubeCollection` gene by gene, with the reason a study has no +result for a gene. ```{eval-rst} .. automodule:: cytepool.meta_analysis diff --git a/docs/source/reference/bibliography.md b/docs/source/reference/bibliography.md index 93d26af..a2816c9 100644 --- a/docs/source/reference/bibliography.md +++ b/docs/source/reference/bibliography.md @@ -79,4 +79,5 @@ densities, and the check of the colour palette for colour-vision deficiency [@go [@virshup2024; @bredikhin2022; @wolf2018; @heumos2026; @miles2025; @zarr2023; @maitinshepard2025; @deutsch1996a; @deutsch1996b; @collet2021; @nist2015; @yergeau2003; -@bray2017; @ieee2019; @iso2019] +@bray2017; @ieee2019; @iso2019]; the Cell Ontology, whose ids align the cell types of +several studies [@diehl2016] diff --git a/docs/source/reference/capabilities.md b/docs/source/reference/capabilities.md index 0657884..076f56e 100644 --- a/docs/source/reference/capabilities.md +++ b/docs/source/reference/capabilities.md @@ -23,6 +23,7 @@ size; its sections give every check. | Reads of selected features from cube files | {api}`read_cube`, {api}`Cube.select` | zarr-python 3.4.1 | exact | [`cube_reads.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/cube_reads.py) | | Fingerprints of cubes | {api}`Cube.fingerprint`, {api}`Cube.digests` | the documented encoding, with zarr-python 3.4.1 | exact | [`cube_fingerprints.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/cube_fingerprints.py) | | AnnData and MuData layouts of scanpy, decoupler and LIANA | {api}`to_anndata`, {api}`from_anndata`, {api}`to_mudata`, {api}`from_mudata` | scanpy 1.12.4, decoupler 2.2.0, LIANA 1.10.0, anndata 0.13.4, mudata 0.4.1 | exact | [`anndata_conversion.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/anndata_conversion.py) | +| Several studies' cubes aligned by sample, gene and cell type, with what each study measured; stacked for a joint model; each study's results lined up gene by gene | {api}`CubeCollection`, {api}`align_study_results` | anndata 0.13.4, pandas 3.1.0, limma 3.62.1, edgeR 4.4.0 | exact | [`cube_collection.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/cube_collection.py) | | TMM, upper-quartile and RLE factors; CPM and log-CPM; `filterByExpr` | {api}`tmm_factors`, {api}`upper_quartile_factors`, {api}`rle_factors`, {api}`cpm`, {api}`log_cpm`, {api}`filter_by_expression` | edgeR 4.4.0 | exact | [`normalization.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/normalization.py) | | DESeq2's median-of-ratios size factors | {api}`median_ratio_size_factors` | DESeq2 1.46.0 | exact | [`normalization.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/normalization.py) | | Model matrices from R formulas, named contrasts, design checks | {api}`Design`, {api}`model_matrix`, {api}`make_contrasts`, {api}`Design.check` | R 4.4.3, limma 3.62.1, lme4 2.0.6 | exact | [`design.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/design.py) | @@ -65,5 +66,7 @@ have no CytePool function yet: - dose-response fits, benchmark doses and drug synergy; - factorization of samples across cell types; - shrinkage of effects across conditions (mashr) and meta-regression; +- combined p-values and replication measures across studies, and cell types aligned + up the Cell Ontology; - edgeR's quasi-likelihood F-tests and DESeq2's tests on every cell type of a cube, and DESeq2's shrunken fold changes. diff --git a/docs/source/reference/cube_files.md b/docs/source/reference/cube_files.md index 4dc9f2e..fedad8b 100644 --- a/docs/source/reference/cube_files.md +++ b/docs/source/reference/cube_files.md @@ -99,7 +99,7 @@ independent implementation ([Validation](../validation.md#cube-fingerprints)). | Key | What it records | | --- | --- | | `source` | Where the data came from, as the caller gives it to {api}`from_cells`, {api}`from_bulk`, {api}`from_anndata`, {api}`from_mudata` or `Cube(provenance={"source": ...})`. Stored as given and never interpreted; its references are reported as unresolved, since following them needs the tool that wrote them. | -| `construction` | How the cube was built: the function, the CytePool version and its settings, such as the layers computed, the target sum of the normalization and where the sample, group and feature ids came from. | +| `construction` | How the cube was built: the function, the CytePool version and its settings, such as the layers computed, the target sum of the normalization and where the sample, group and feature ids came from. A cube that {api}`CubeCollection.stack` built records each study's cube by its fingerprint, with that cube's provenance and its groups under their aligned ids. | | `operations` | Every operation that derived the cube, oldest first: `select`, `regroup` or `collapse_samples`, with what it kept, merged or pooled, and `input`, the fingerprint of the cube it was applied to. | Because each operation names its input by fingerprint, a chain can be checked link by diff --git a/docs/source/references.bib b/docs/source/references.bib index c70473b..63f498d 100644 --- a/docs/source/references.bib +++ b/docs/source/references.bib @@ -393,6 +393,17 @@ @techreport{deutsch1996b doi = {10.17487/RFC1952}, } +@article{diehl2016, + author = {Diehl, Alexander D. and Meehan, Terrence F. and Bradford, Yvonne M. and Brush, Matthew H. and Dahdul, Wasila M. and Dougall, David S. and He, Yongqun and Osumi-Sutherland, David and Ruttenberg, Alan and Sarntivijai, Sirarat and Van Slyke, Ceri E. and Vasilevsky, Nicole A. and Haendel, Melissa A. and Blake, Judith A. and Mungall, Christopher J.}, + title = {The Cell Ontology 2016: enhanced content, modularization, and ontology interoperability}, + journal = {Journal of Biomedical Semantics}, + year = {2016}, + volume = {7}, + number = {1}, + pages = {44}, + doi = {10.1186/s13326-016-0088-7}, +} + @article{dimitrov2022, author = {Dimitrov, Daniel and Türei, Dénes and Garrido-Rodriguez, Martin and Burmedi, Paul L. and Nagai, James S. and Boys, Charlotte and Ramirez Flores, Ricardo O. and Kim, Hyojin and Szalai, Bence and Costa, Ivan G. and Valdeolivas, Alberto and Dugourd, Aurélien and Saez-Rodriguez, Julio}, title = {Comparison of methods and resources for cell-cell communication inference from single-cell RNA-Seq data}, diff --git a/docs/source/statistics/cubes.md b/docs/source/statistics/cubes.md index 6c65b4f..d943190 100644 --- a/docs/source/statistics/cubes.md +++ b/docs/source/statistics/cubes.md @@ -1,5 +1,5 @@ --- -description: What a CytePool cube holds, why every layer is a sum over cells, and how merging cell types and pooling samples stay exact. +description: What a CytePool cube holds, why every layer is a sum over cells, how merging cell types and pooling samples stay exact, and how the cubes of several studies are aligned. --- # The cube and its additive layers @@ -84,3 +84,34 @@ lazy layer at once. {api}`to_anndata` and {api}`to_mudata` convert a cube to the pseudobulk layouts of scverse, and {api}`from_anndata` and {api}`from_mudata` read them back. + +## Several studies + +Results of one study often do not replicate in the next [@nakatsuka2025], and studies +are compared or combined only once their samples, genes and cell types are aligned. +Aligned by hand, they go wrong in two ways: donor ids repeat across studies, and genes +that a study did not measure are read as zeros. {api}`CubeCollection` aligns the cubes +of several studies and keeps each as it was given: + +- **Samples** are named `study::sample_id`, and so are the values of the columns of ids + that `within_study` names, such as donors: donor `D1` of one study is not donor `D1` + of another. +- **Genes** are aligned by `feature_id`, such as Ensembl ids, and + {api}`CubeCollection.measured` says which study measured which: a gene that a study's + cube lacks, or marks as not measured in its feature column `measured`, is missing in + that study, not zero. +- **Cell types** are aligned by their Cell Ontology ids [@diehl2016] when every cube + has them, by group id otherwise, or by a mapping for each study. Groups of one study + that share an id are merged, exactly, since every layer is a sum. + +Each study is then analysed on its own: {api}`CubeCollection.aligned_cube` gives its +cube with the aligned cell types and the genes it measured, and with its own samples and +library sizes, and an analysis such as {api}`limma_voom` with the study's own design +gives results keyed by the aligned cell types. {api}`align_study_results` lines them up +gene by gene, and {api}`differential_expression_meta_analysis` combines them +([Meta-analysis](meta_analysis.md)). For a joint model with the study as a covariate, +{api}`CubeCollection.stack` puts the studies into one cube with the genes that every +study measured and the cell types that every study has: a cube holds no missing values, +so it refuses the others rather than fill them with zeros. The collection's provenance +names each study's cube by its fingerprint, with the alignment, and so does the +provenance of a stacked cube. diff --git a/docs/source/statistics/meta_analysis.md b/docs/source/statistics/meta_analysis.md index f0445fd..56c1721 100644 --- a/docs/source/statistics/meta_analysis.md +++ b/docs/source/statistics/meta_analysis.md @@ -64,6 +64,14 @@ tested by at least `min_studies` studies (two by default, since one study says n about $\tau^2$). It takes limma's moderated standard errors as the known sampling errors. +A gene missing from a study's results was not measured there, or was measured and left +out, as by the expression filter, which are different findings. For the studies of a +{api}`CubeCollection`, each analysed on its aligned cube, +{api}`align_study_results` lines up the results of one cell type and contrast with a +row for every study and gene, which says whether the study tested the gene and, where +it did not, why: it did not measure the gene, measured it and left it out, or did not +analyse the cell type or the contrast. + ## Limits The studies are taken as independent, with known sampling variances. Meta-regression on diff --git a/docs/source/validation.md b/docs/source/validation.md index cfa1633..0e4ae10 100644 --- a/docs/source/validation.md +++ b/docs/source/validation.md @@ -16,8 +16,8 @@ exits with an error when any check falls outside its tolerance. This page summarises the tables the scripts printed when they were run against the package in this revision of the repository: on **2026-10-26** for `plotting.py`, -`variance_partitioning.py`, `gene_set_files_and_ora.py`, `kenward_roger_and_f_tests.py` -and `gsea.py`, and for `dream.py`, which printed the same table again after changes to +`variance_partitioning.py`, `gene_set_files_and_ora.py`, `kenward_roger_and_f_tests.py`, +`gsea.py` and `cube_collection.py`, and for `dream.py`, which printed the same table again after changes to the code it checks; on **2026-10-25** for the six scripts on cubes and their files (`cube.py`, `pool_cells.py`, `cube_files.py`, `cube_reads.py`, @@ -43,6 +43,7 @@ covers and the largest difference among them. Run a script to see every check. | [`cube_reads.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/cube_reads.py) | Reads of selected features | zarr-python 3.4.1 | 5 | 5 pass | 1 min 46 s | 500 MiB | | [`cube_fingerprints.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/cube_fingerprints.py) | Cube fingerprints | the documented encoding, with zarr-python 3.4.1 | 7 | 7 pass | 4.4 s | 122 MiB | | [`anndata_conversion.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/anndata_conversion.py) | AnnData and MuData | scanpy 1.12.4, decoupler 2.2.0, LIANA 1.10.0, anndata 0.13.4, mudata 0.4.1 | 18 | 18 pass | 1 min 5 s | 1,483 MiB | +| [`cube_collection.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/cube_collection.py) | Several studies | anndata 0.13.4, pandas 3.1.0, limma 3.62.1, edgeR 4.4.0 | 15 | 14 pass, 1 control | 26.7 s | 495 MiB | | [`normalization.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/normalization.py) | Normalization and filtering | edgeR 4.4.0, DESeq2 1.46.0 | 50 | 46 pass, 3 deliberate differences, 1 control | 72.9 s | 745 MiB | | [`design.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/design.py) | Designs and contrasts | R 4.4.3, limma 3.62.1, lme4 2.0.6 | 27 | 26 pass, 1 information | 4.3 s | 221 MiB | | [`linear_models.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/linear_models.py) | Linear models | limma 3.62.1 | 218 | 200 pass, 18 information | 4.5 s | 177 MiB | @@ -69,7 +70,7 @@ covers and the largest difference among them. Run a script to see every check. | [`plotting.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/plotting.py) | Figures | limma 3.62.1, edgeR 4.4.0, DESeq2 1.46.0, metafor 5.0.1, speckle 1.6.0 | 26 | 23 pass, 2 information, 1 deliberate difference | 25.4 s | 810 MiB | | [`distributions.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/distributions.py) | Distributions and special functions | mpmath 1.3.0 | 46 | 46 pass | 6 min 39 s | 221 MiB | | [`smoothing.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/smoothing.py) | LOWESS and spline bases | R 4.4.3 | 90 | 87 pass, 1 information, 2 deliberate differences | 1.5 s | 175 MiB | -| **All** | | | **2,174** | **1,874 pass, none failed; 172 information, 127 deliberate differences, 1 control** | | | +| **All** | | | **2,189** | **1,888 pass, none failed; 172 information, 127 deliberate differences, 2 controls** | | | Run time is the wall time of the whole script and peak memory the largest resident set of any of its processes, reference programs included, both measured with GNU @@ -280,6 +281,45 @@ a missing profile: the import refuses them until cell numbers are given. decoupl gives its row without cells the metadata of another row (found by probing); CytePool reads a column as sample or cell-type metadata from the rows with cells. +## Several studies + +[`validation/cube_collection.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/cube_collection.py) +simulates the cells of three studies of one stimulation (seed 20261026) and pools each +with `from_cells`: `lupus`, 8 samples of 4 donors, each control and stimulated, 4 cell +types with Cell Ontology ids and 700 genes; `covid`, 10 samples of 5 donors whose +sample and donor ids repeat those of `lupus`, 5 cell types of which two share an +ontology id and merge, and 750 genes, 40 of them filled with zeros and marked as not +measured; and `atlas`, 6 unpaired samples with integer donor ids, 5 cell types without +ontology ids, aligned by a mapping that leaves one out, and 950 genes. The collection +aligns the cell types by ontology id and writes donor ids `study::donor`. anndata and +pandas run as a separate program: each study's cube is exported with `to_anndata`, the +genes marked as not measured are dropped, the rows of each sample and aligned cell type +are summed, and the studies are joined with `anndata.concat` (outer and inner joins, +labelled by study); pandas joins the sample tables. R runs limma-voom on each study's +own counts, summing the cell types that share an aligned id, and on the stacked counts +of `lupus` and `covid` with `~ donor + condition`, by the steps of +`differential_expression.py`. Reference: anndata 0.13.4, pandas 3.1.0, numpy 2.5.4, +Python 3.13.12; R 4.4.3, limma 3.62.1, edgeR 4.4.0. CytePool ran with numpy 2.5.4 and +numba 0.68.0. + +| Check family | Reference | Largest difference | Tolerance | Outcome | +| --- | --- | --- | --- | --- | +| The sample table: ids `study::sample_id`, donors `study::donor`, the union of the columns, ages that are integers in two studies and floats in one as float64, and missing values where a study lacks a column | `pandas.concat` | 0 mismatches | 0 | 1 pass | +| The genes, as a set and in the order in which the studies first list them, and which study measured which | `anndata.concat`, outer join (NaN where a study lacks a gene), `pandas.unique` | identical; 0 mismatches | 0 | 2 pass | +| Every layer of each study's aligned cube: merged cell types, measured genes | the outer join's values | 0 | 1e-10 | 1 pass | +| `stack()`: genes in order, cell types and samples; every layer; the cells of each entry, and an entry without cells | `anndata.concat`, inner join, kept to the cell types of every study | 0 | 0 and 1e-10 | 3 pass | +| `limma_voom` on each study's aligned cube, lined up by `align_study_results`: logFC, AveExpr, t, P.Value, adj.P.Val, B and the standard error, in 12 cell types of 3 studies | limma, study by study | 3.8e-13 | 1e-10 | 3 pass | +| The status of every gene of every study: tested, measured and left out, or not measured | limma's tables and each study's measured genes | 0 mismatches | 0 | 3 pass | +| `limma_voom` on the stack of `lupus` and `covid`, `~ donor + condition` with donors `study::donor` | limma on the stacked counts | 1.4e-13 | 1e-10 | 1 pass | +| The same with the donor ids as the studies give them | limma on the stacked counts, donors `study::donor` | 1.9 (log2 fold changes) | must exceed 1e-10 | 1 control | + +limma ran on each study's own counts, with the genes filled with zeros, and its tables +equal those of the aligned cubes, which leave those genes out: genes without counts +change no library size, expression filter or normalization factor, so a gene that a +study did not measure can be marked missing without changing any other gene's result. +The control shows why donor ids are written with their study: taken for one person, +donor D1 of two studies changes the fit of every gene. + ## Normalization and filtering [`validation/normalization.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/normalization.py) @@ -2168,6 +2208,7 @@ python validation/cube_fingerprints.py SCANPY_PYTHON=$REF/python-scanpy/bin/python DECOUPLER_PYTHON=$REF/python/bin/python \ LIANA_PYTHON=$REF/python-liana/bin/python \ $REF/python-anndata/bin/python validation/anndata_conversion.py +REFERENCE_PYTHON=$REF/python-scanpy/bin/python python validation/cube_collection.py ``` Here `python` is `$REF/python/bin/python`, the environment below. With diff --git a/tests/data/cube_collection/studies.json b/tests/data/cube_collection/studies.json new file mode 100644 index 0000000..70f7c71 --- /dev/null +++ b/tests/data/cube_collection/studies.json @@ -0,0 +1,621 @@ +{ + "studies": { + "lupus": { + "samples": { + "sample_id": [ + "S1", + "S2", + "S3", + "S4", + "S5", + "S6", + "S7", + "S8" + ], + "donor": [ + "D1", + "D1", + "D2", + "D2", + "D3", + "D3", + "D4", + "D4" + ], + "condition": [ + "ctrl", + "stim", + "ctrl", + "stim", + "ctrl", + "stim", + "ctrl", + "stim" + ], + "age": [ + 30, + 33, + 36, + 39, + 42, + 45, + 48, + 51 + ], + "sex": [ + "female", + "female", + "male", + "male", + "female", + "female", + "male", + "male" + ] + }, + "groups": { + "group_id": [ + "B cell", + "CD4 T cell", + "CD8 T cell", + "monocyte" + ], + "cell_ontology_id": [ + "CL:0000236", + "CL:0000624", + "CL:0000625", + "CL:0000860" + ] + }, + "features": { + "feature_id": [ + "ENSG00000000001", + "ENSG00000000002", + "ENSG00000000003", + "ENSG00000000004", + "ENSG00000000005", + "ENSG00000000006", + "ENSG00000000007", + "ENSG00000000008", + "ENSG00000000009", + "ENSG00000000010", + "ENSG00000000011", + "ENSG00000000012", + "ENSG00000000013", + "ENSG00000000014", + "ENSG00000000015", + "ENSG00000000016", + "ENSG00000000017", + "ENSG00000000018", + "ENSG00000000019", + "ENSG00000000020", + "ENSG00000000021", + "ENSG00000000022", + "ENSG00000000023", + "ENSG00000000024", + "ENSG00000000025", + "ENSG00000000026", + "ENSG00000000027", + "ENSG00000000028" + ], + "measured": [ + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true + ] + } + }, + "covid": { + "samples": { + "sample_id": [ + "S1", + "S2", + "S3", + "S4", + "S5", + "S6", + "S7", + "S8", + "S9", + "S10" + ], + "donor": [ + "D1", + "D1", + "D2", + "D2", + "D3", + "D3", + "D4", + "D4", + "D5", + "D5" + ], + "condition": [ + "ctrl", + "stim", + "ctrl", + "stim", + "ctrl", + "stim", + "ctrl", + "stim", + "ctrl", + "stim" + ], + "age": [ + 30.0, + 33.0, + 36.0, + NaN, + 42.0, + 45.0, + 48.0, + 51.0, + 54.0, + 57.0 + ], + "batch": [ + "b1", + "b1", + "b1", + "b1", + "b2", + "b2", + "b2", + "b2", + "b2", + "b2" + ] + }, + "groups": { + "group_id": [ + "B", + "T4 naive", + "T4 memory", + "T8", + "NK" + ], + "cell_ontology_id": [ + "CL:0000236", + "CL:0000624", + "CL:0000624", + "CL:0000625", + "CL:0000623" + ] + }, + "features": { + "feature_id": [ + "ENSG00000000005", + "ENSG00000000006", + "ENSG00000000007", + "ENSG00000000008", + "ENSG00000000009", + "ENSG00000000010", + "ENSG00000000011", + "ENSG00000000012", + "ENSG00000000013", + "ENSG00000000014", + "ENSG00000000015", + "ENSG00000000016", + "ENSG00000000017", + "ENSG00000000018", + "ENSG00000000019", + "ENSG00000000020", + "ENSG00000000021", + "ENSG00000000022", + "ENSG00000000023", + "ENSG00000000024", + "ENSG00000000025", + "ENSG00000000026", + "ENSG00000000027", + "ENSG00000000028", + "ENSG00000000029", + "ENSG00000000030", + "ENSG00000000031", + "ENSG00000000032", + "ENSG00000000033", + "ENSG00000000034" + ], + "measured": [ + true, + false, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + false, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true + ] + } + }, + "atlas": { + "samples": { + "sample_id": [ + "S1", + "S2", + "S3", + "S4", + "S5", + "S6" + ], + "donor": [ + 101, + 102, + 103, + 104, + 105, + 106 + ], + "condition": [ + "ctrl", + "stim", + "ctrl", + "stim", + "ctrl", + "stim" + ], + "age": [ + 40, + 41, + 42, + 43, + 44, + 45 + ] + }, + "groups": { + "group_id": [ + "Bcells", + "CD4", + "CD8", + "Mono", + "Doublets" + ] + }, + "features": { + "feature_id": [ + "ENSG00000000003", + "ENSG00000000004", + "ENSG00000000005", + "ENSG00000000006", + "ENSG00000000007", + "ENSG00000000008", + "ENSG00000000009", + "ENSG00000000010", + "ENSG00000000011", + "ENSG00000000012", + "ENSG00000000013", + "ENSG00000000014", + "ENSG00000000015", + "ENSG00000000016", + "ENSG00000000017", + "ENSG00000000018", + "ENSG00000000019", + "ENSG00000000020", + "ENSG00000000021", + "ENSG00000000022", + "ENSG00000000023", + "ENSG00000000024", + "ENSG00000000025", + "ENSG00000000026", + "ENSG00000000027", + "ENSG00000000028", + "ENSG00000000029", + "ENSG00000000030", + "ENSG00000000031", + "ENSG00000000032", + "ENSG00000000033", + "ENSG00000000034", + "ENSG00000000035", + "ENSG00000000036", + "ENSG00000000037", + "ENSG00000000038", + "ENSG00000000039", + "ENSG00000000040" + ], + "measured": [ + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true, + true + ] + } + } + }, + "samples": { + "sample_id": { + "dtype": "str", + "values": [ + "lupus::S1", + "lupus::S2", + "lupus::S3", + "lupus::S4", + "lupus::S5", + "lupus::S6", + "lupus::S7", + "lupus::S8", + "covid::S1", + "covid::S2", + "covid::S3", + "covid::S4", + "covid::S5", + "covid::S6", + "covid::S7", + "covid::S8", + "covid::S9", + "covid::S10", + "atlas::S1", + "atlas::S2", + "atlas::S3", + "atlas::S4", + "atlas::S5", + "atlas::S6" + ] + }, + "donor": { + "dtype": "str", + "values": [ + "lupus::D1", + "lupus::D1", + "lupus::D2", + "lupus::D2", + "lupus::D3", + "lupus::D3", + "lupus::D4", + "lupus::D4", + "covid::D1", + "covid::D1", + "covid::D2", + "covid::D2", + "covid::D3", + "covid::D3", + "covid::D4", + "covid::D4", + "covid::D5", + "covid::D5", + "atlas::101", + "atlas::102", + "atlas::103", + "atlas::104", + "atlas::105", + "atlas::106" + ] + }, + "condition": { + "dtype": "str", + "values": [ + "ctrl", + "stim", + "ctrl", + "stim", + "ctrl", + "stim", + "ctrl", + "stim", + "ctrl", + "stim", + "ctrl", + "stim", + "ctrl", + "stim", + "ctrl", + "stim", + "ctrl", + "stim", + "ctrl", + "stim", + "ctrl", + "stim", + "ctrl", + "stim" + ] + }, + "age": { + "dtype": "float64", + "values": [ + 30.0, + 33.0, + 36.0, + 39.0, + 42.0, + 45.0, + 48.0, + 51.0, + 30.0, + 33.0, + 36.0, + null, + 42.0, + 45.0, + 48.0, + 51.0, + 54.0, + 57.0, + 40.0, + 41.0, + 42.0, + 43.0, + 44.0, + 45.0 + ] + }, + "sex": { + "dtype": "str", + "values": [ + "female", + "female", + "male", + "male", + "female", + "female", + "male", + "male", + null, + null, + null, + null, + null, + null, + null, + null, + null, + null, + null, + null, + null, + null, + null, + null + ] + }, + "study": { + "dtype": "str", + "values": [ + "lupus", + "lupus", + "lupus", + "lupus", + "lupus", + "lupus", + "lupus", + "lupus", + "covid", + "covid", + "covid", + "covid", + "covid", + "covid", + "covid", + "covid", + "covid", + "covid", + "atlas", + "atlas", + "atlas", + "atlas", + "atlas", + "atlas" + ] + }, + "batch": { + "dtype": "str", + "values": [ + null, + null, + null, + null, + null, + null, + null, + null, + "b1", + "b1", + "b1", + "b1", + "b2", + "b2", + "b2", + "b2", + "b2", + "b2", + null, + null, + null, + null, + null, + null + ] + } + }, + "mapping": { + "atlas": { + "Bcells": "CL:0000236", + "CD4": "CL:0000624", + "CD8": "CL:0000625", + "Mono": "CL:0000860", + "Doublets": null + } + } +} diff --git a/tests/data/cube_collection/studies.npz b/tests/data/cube_collection/studies.npz new file mode 100644 index 0000000000000000000000000000000000000000..c51470f82be65c2aad26bc043f575755ce2ac42b GIT binary patch literal 58047 zcmbrkcTiJX7d}i8P(eVYh%_5bq$?ewfFL5Eaut*k0Vztap$7y-q>1zz5D=6SkuH#c zKtNjPog@%?3lKs`|9pSPpK;?qpZpp5nV~^L0etO9 zT-$aysOG!x?^4U@cM?kW3BF$|B~zPYoMx?a`kjUG*XnvDy0MS>igf&Bhoc>`1_+5vnxPcV)VwW4Qb4usWDDC~2m z)T@0Ec@ch=qJwEFyr;AL{^dz;DeqJ=8_R7D*M?8{+p-u7jLLCi6Im{*`gY|mn=XqP zGl<0*(*j3sw{-GE3kDm=NAc*DT@4iv)6na>Lp-TOcx!lflGpP3Rv>SGq~`i63nwNx zh2K#{s7+wx%FJHDiOJNa&4Ax$Le8WhJF0IT3m?9_4BS;}yPG3?+*XJ_kDqGy3VwaE zCX=A)pnyk@)&=~c;#q`2P-r$Q{ zQ!{HaG{dBcevLgYVwn9)o+ys?ErT1WH*xLQy@SZdY3DjdwlmK}FwN}986djf-<-rD zB5D^CnQ@VxJt6Kf6Z#-`sTWgX)Q=R}&7QUPINqKpfnxxtS#)D+K>x{#iKf?I8kuHR zldj|oKlH1FbeL-dy&T>Q;`HuyQ>87`7L4YpKLtYc;4ES%baL9<9IeHhY#+FA9dJDJ za)N)w#O;0Hb;dv?4>@~divbyfyV2G0RJ*;$dwj4CX~@K^!4&!ggkYFQ9;QVRMdr3{ zK1^zr{ca4G^Jo5%w|d0`6ELC~1>&1VlF!!s6k#j)6`{XAQZ^FLVR;9CYFQ``#k$KL zKHx@3AG*jtb7q_Ky@5}IVP0FdRTa)$(iu(w8666<+&Q8>rr6-C2wdTjg)XDuyXxPD zlY7t{XEUUB2;J=wVX8(%`ygm^BNoSU`~$i9g*y<|vhLzg{#!G_RQ_Q!q4srWNWqcd ze#r#y@0W~M$Mfqo5z|Zn;D972uy)b>DP%8Xv#}i_NEXC6CtAZz1UCg*$ipu4LQht$ zlI!o_-1dyN$EttcgGh%&J@3B@5NP8*9%5#^Z{=MG!l?fU+zWB7?(J4y?`Mu*3?2&p z#8ms8)8V`UTwa~$f`V#vgg|cYU~dOx+Kk=YNUR8NNzt#2t!q&0L~}okP(D@uZd7B3 z=Ve{f4FhPUTyt#onYO_in90M?8ljMZm8nYz%CJ_y*xYoqW1sQ5?sM~9V*gwCxl3WB z-}!Add5d=;gF8=pXB7uR^}wYnoK5WGhD4_-R*We|e*H9c*{>j4H@9e(M>H%DByzT* zChTzZrnv)XCiWuFVHlLul1<3bau*nWI)zhi^BM<<7G@$!{2qt$bDZC6PA0TC5Iz%> zM5t9+P=`I6dukm*=Q*8@BN=I14_qL&)yHz>E) zhMl~b8zvSmG$>n1vb)2&v{b>h>LBx)wzZHQnlcB zUiFGWy*l14!ks3DwJq*AdY2i_XkG5CAr7v(<_gdDTh9139w|48&{z^XXb}=(6~A5I zM=kpZZvAGipgj}Y`W@e!WHcLI5Jw9Z5JydZKD_*s)>>NXKpYw=e{_q z8K6byjZF6%#o8hqR9+TNIp5)HRY^JW9I=x-8#+hmK%z19-V3?PIaJtSw_uV|Na|$$ z(7WePvEbKoK^_@pqp9_d+Y6<)N-YbICPdaNlG;pB=R(0^0-?me$nKLX%_?Own*AF6 zLcX`M!sA7Ai+~6#Pl!cM7R<yxk(Iu%#ZaMZdsD zi}?`xkorXB=+xmA1b*=BL^}{+&^YJt$bXpFH0}468Z6FavKSs)2Z-$dAzz z4KslHi6_sh0GUs0pCVh1iWgdSc4s~O_JCe?MAJ1i+#QuWqzjx}`6jPmDHuSq5F_*> z?tc~G<@~9jbZPVC+Y4Hg-EiKMq7GqM%wNdC(ukxLsz%M4sMJLY{4&$5()8A6x5F9n zg;@Jv$HO6dJpnDdU>5c9`!~Aoy?kJ1P|}kyqK>}T>vA+hv@`nESf-mJ<$}`EWA;EC zf8o}ktcsn~UQ2kJr4scZxh$dKMcCVqUe`aaE4kFwL@P)f87xNVj$Re)$rlysb=qV- z%PeOf#_CX(khm$qbWsg7eBSN#lTDtp!>+h>Jn2x7t(?zZ5jFz@f) zhlqz&1)rYh*_3#e7J8<>Dyl9nEXV-<^*Aj{5CTm5bov~R%ym;Wn=1k-_tN#Tp~? z%qkqt@;n!Tyoo$`L5`e=e;L4UWg`8L0so=Qj$VPUUjBBfb^6L@zNh6*k7Dk% zl#glX$&a|!S9}a`tL0u1nf(q;Cb~>AC|YZz@B$n$SQoh#5e(~BXetTm-IE^R^1I(N zVJ2>EIa^$ESh@8x<7rx266CpBvFAYAh{YCSVwUnvDYb0sAcEYbH1;~!6BKgSxGafZ z!d5fmN!oxc2^^7V#zPj3*9#KO>Q}en{UsR3`ysp0`KAlCE5ttHq$tgHpb6pQ<;kZf zi}Bgu)nw(6)Ia9<9^>P6J;u7d4T296^*sHiN~G~){eHV-!o@xKDn?}W8CS+HSFGIC z7}fT$9+|+V^tPFHH^-`UCby`;3l>F|?aVk$Xq$-3`%b_$2|uY-TmXjaSUV{=)KCa} z|0BEFRI4oO*-rhWmZUb1sJn?0xEIH;gxUw0TF>U%sq%Z)dsOY{mEl1jsb5o<%SLE1 zyeQ-MXq5}^gW@)7V?m>}M1030K~%E@UPt0D@L z4+}q|om%tC>LjQZ{yvxk8||a!BL16m@g+Fc)&3opTmRT%K;eb+4R0|P8pm^d7LOCE zfOP5<$t?DsZ?gA1WYML#A5Ie4?$GU}_5p6ZWtq%$Hh~PWtC& zdUi5CnI*uwOMk!KpF=O}q41_XH4cY8mx95#rk##vNUef1>;k31t)%vPw13ta)N{UD zeYC!1T49vuyLJ9+-D|>oJ{&H#z0@O-ZiKp_3BLGV4eZ+;OKnhc7dAT2Lck5y8Q^iK zE6246aRwgc^YJ+z&AnavmIy@QoJ1U`!kSva#LK-pk8C_2x1uicsf=BuerU^1kR1RP z&SjG72Y_6oYDN6dedvDenQww{>rFbGFQ-s@&)HC{DFxjl8xnKrlo@AqzfdsgIa+u**FzL!^jw)wZ%r5T2ytRLoGXaSsmfeB zLQWi`2}g9wCe39NcLwAKfVaJxYYtnbd)JT-I~wb-!LWNJ)N!;A&#H{`fS0)v&{>-a z7akS*V)qdjrFo5Fexah%iyjuEB-6`*z7plBM2|j*B(+UIoiLM33Y^~(?B0fx2P$GYe_(t3K!Z%};VuKI<=$TLQ(U30Z5T=23GUY|vQL zn5|Jlg9@z}>6wR~4C{73ofmeonI#zXg&#hV=O^zC$L0W5?tb%d9GPnwpeO@{)6Q2# z1iabpY?qpc@h53r)(@uHuiJ4~316Ma7LUrGsmmc+k9}&C{M4E&N(9~s?xZ>hg%q5V z@vNGB7V#rz-><>9wY#+2K)gqE{GNZ=GcINjtJ)n|7I}rq)lsYc<7)g|VMTPcUqv#~ zZ>TxtfWdaBzQ;H0+}-dMnEgQkwh8UjhtyRo5R!`0641mc0fgb-Jg?=e`OpUyU_SF< z+>C2axwm5qG9%{{)@PO*T5nan4cmv_AvMQR5CaY~AOKP}O1~;o4{Cl9?w;2hY@|}+ z(qC5ARJ9G#__=dpp=cnPR{`RKC@hTV+Z^#UoO}8Wh&(|$eo)p@W7QhIP5zPN;)8ct zr4)i9wQs6hgs6&IbKleS6rkIs{*~PHw(H059WP%i_vgFx%fN#)gQYg`jrr!_ZxE24 zm#(+=n)c{DG78bM*gLK(hx*xsEXm)!TInRokAxHEN76K&<$)CxgHB zq^$VjJI`u2zL~a_Z`la`I3}8>5DbY+BSP?qpuHeh?+M-rH zORZ@k@uKZ&DcD!Ixhgef#_Wp!+MDxb<3E_Rq6f!-eBx1A%u^_vY z!||lGt;qB0Ychrsl#E~g3Mb~lt-J0|734B4wQMW&#Y1!eee~GWPh&3?t)CanT9$~^ zP2s|8ZBYl3Ki3wPZ&jb9cGm1(vn>R(nyFnkB7Xi!p`oUnv`;=&|62n|8d>t=LeSVYSKVn z;k-Ry3+vq7moHm5cjKviU8_6g3k%2Uhbg<7!}PU-g8Z*R{RBFNcS~ivnuW@_Co>tk zEQ=US!bkhY(khV}Oelj|@pm5EbQg`^g%hm^Q-31WvClc=+2?RH*4~U4J1y1duD8{vP54199&Sv2-KeTQ$2>_lFMPAZnqa zqDT!V3{r_w!Z*eN$NFK;W;HxslfAVs_aew?RI7E@)w6 zfaKcV>ZYGm_ht%x)!%&JW4Jg^AE$NU2tPIKsI7DM`CfzCQw7)k$)eSKwNL_fOue{%LF3JVQ)gqLQEl3JyB}f#SuW@lheg{LbmP4n0}tuh{D=8v zmwULoxNB!rrjDhG_#^|WNoCzNE58@^@}h7_H?27ctw$BeC)x^qACsX>Ak)Wach47K z%hs+{_Q?s`K&CC;%izfiYpxn3PEwN6(|u9bI2+BYKanl)aRK?9k#UnzfI zOSb|Q$h5CWAEQdB_}#*J*wlu66cp*Wf-bt&l33$VmATb@;KZ%j{D#PRb!AF6Z`X_g zdj!bt+1$1ohtAir zD^FTdk;WULe_0;23i&tI0)sOS?uF~=Jh`=ZZQ0uE=$38!cn4SaTLLuDG8nxPxZ28E zvug`KOwl!Yv()K76}u*=zIAm-L-!M7FyQ9|owgZWH&n0XYP7mFpu6`x7dmT??w>r? zsU6ygNdM>weF+<;$i_datsxC4BxL06iv)pr`eIr0c57>v)^(u8fr83_A9 zh_8NYwYwc&Ayz0n11lAYef<#9h{QcAr&@;7bRUO@e#FZs@|9sP9&QtwhPR&O#)S2U zRDrBQg*d9sKQu1g3xB{zlZ29L3GT!+W1gpg?OZ4c9iQpPQo#COX<}vS zb^Z=lKrU@)ceU6Q+00fV@*$FQY2ZxPTUM#RMlNx`|IF8PDS{1r;pucAcA~&E?r>Qs z83NM5X78^~Rj4Td;W{a}14qBMmT)*)+uExD9QtFXT20*Z2uQSqpcl5GplEvcMA_Ve zIBln=LZJmf<1laBGRmAIVEbd|2uK2hn)94ae|`7V&V36127BIXbd$z(q1=&5`gxD| z`#jNuam9k4Qe${JWxPYuRX8j}z3xo8!iQmza(&spBfyHEHglhUJSnK(5z^SW9xm@I z`$_Ee_YTkz3UqsI;I2?tztWZ0>FH#?^n6%yy(}R!YH0PU)NWcytl`SukjaWXzi3_P zcOzRUK4JbGB9r-fx&Q~l+QOYx@}_uVr(8bi!%}@P-#-4^xbE7Aox20F{n?UxZ#6Q1 z+lv#pw`r$Gd45_f1zyvchTX^uMzns{a9iHZd&aetaJQYh9MndX6QzwWCKbJN@opt{ zgk!t&9eTX)cz1Z=G~}8xV#tR0>_I zU64Z%+&`<(ro&cgZlrk@wR*1^9^W<<#B31z;xoj&(lRD_O9X0WJfV0c!W3I{Vs2c= zH5})O4ZR#7e!(}Zv)Olz`JHxFd0zC5_j3|q7v2SVRP`KrrIkHuQ+g+`eU}>~_xaH? zO-Bu%t#*!INBhHe(x9R&{W+A|Zo*pEPP?QYPvBQVz1GXkgnhDV6wirjp()Fgb+;Du zoXqA`_oA{=aME^KW6Gn!`$xq0Q>)h3U!+CZpfw6AIp7q~cO}%~O9=|8guXC1>sWFg8S4?G`dafpwdn z_-C^f8B+uKY_K4PYM+P5a#Hok~Op z*)#upXYGELCWqTn`Ja`RYl=LKKoxEe$JpJYW6i;5V%A5V4;tBO(4o4!Av-Iuov0^}&>(Y9d z>P@`*pBI86{>AQ_4Z+(M$TOvEBg=mlw7S+^$cgaBHxF+cS@tA;;g*z7kmSMh474eN z;y(Z8{JhWnzs3=o?;8GN+<)pS_Wthn|87M-JomRzuBZN<2eaJLGM@1FAn~oZ0WFYo z6PJWE>q%DS6DKc7Uu%45oN`dq>y{_?d z0$++EqvYNyUKx0^G>OA8Dkq4#gz6bLNG(!>GUI|qR*I<~zPqC)wYX((Pl)-cfpxJ2 zK6|p@XqUN}^KCOu#Tl7zs};VTFmry>00-0xFcUr0NOg0WRPq|D|FZcZ=iAW2?Bri^ z-^h@NyKDL~qv!c1iSwG#3k-H5-#CDjpPcfOq>t8s%%NZur6y;~z0d>|ZX9J(q z?xHJ1!a1r=wn`|(WMCjl-~?KgkOD^|D>yG}ioe}<;xb<(M?8l&kDEBNnj+o)`;D?} zI1@bg8tZZ3J{tepBoK=EBujOOdF{ToixyHjuuK01{ z=D>_ryI(f?96b1uct84lIg-rL9$a`H58Nzhcbw%xdnv;-g@4(3ztc(8MS?CsSx_gA zq%wSNe9&N}pT_UzRXK#6cwX~Tp3`0`$W`}8VXaJpbxoO0lW<*2M`z>@Scb{Lu!=tjPj*?{ciB7NzR9xHth(!#N~Wf+yKx;0w`21>PQcCx9dqXTp|lcz zvrqOpsSMTgo3937=d&CNb+~y(|6QoXTdAK>~*f%KhMwG^ddxaNy0034= zK5D z?ay|{t~s`4L_*^mwyfQs85xED2G>7r7iEuxuuPWFk7$(=Kcti!nFMRsU>Bs^qQkw< z;LjeYz0e0hEat;c`@aAOFCEXGC1{E~6MLb3{y1h|@r4LlI9gQr?(+znOa<+fHC$zM zM5K=j0&W_{Fkauf*>rAE{-ubTg@gMzE1b>l)C0+ff_;egbKP#QpZ#Qt&*ou4pQ7he zpKH3%F0W>?HVSCCB$(|)AGC<~mt}$F4^q)vL&=gu=dXTTd9^nhuh(=v;KpxVY&-ZQ zBeBO?_QBOUz4J0;#$nmKx7YVw@e_tDRUbhr$5I_imjOHM^Zntmw?nVj=sey-cU`Ub z6g%0V@`a0VjuP&)S*yP{sJvG?{8VDdgH5U6q0{cDIpPUvhhNN3ypwOVjD{H1L{+GeRH8Hts@|S}(=Crh-)J&Ysj&M(-bmZEDSD z&l3Ho?p*u&C@;J77>(hGf~$zJiEw2+t7&5kVf%g_CZ~=*Z37)vno99Ws)tA>G1!FG zMhOY%-f8$dqh?Qgx>8Vw#IbBxrb>i@Pqh#;x~IeKW%RCCJo3~^^$Vcq>3ROG39BRG zgQ{pv2XFaf>rs~?w?HABO9!gvpF>2YPEp!CKZfPee1AmLEt;{`jvyZT$ z7U8arR6mqA7L|o6>UN<>LTfnB_I4ILFehLD*MVkqz>*Rdz&9Mdo*^`C5&57fOFadH z$f}M2&r{SgW5%#`l!9rS%k|^CVZjX!kr9{xZ_XfFniJhjTm zR)N~5(mIN`9KZXV!PF;!;$*_8CSK-))*g9Wr^m!; zS9DSjV7;y-Zr1u`TFx1*aE@Zg@xrDc+)~wW#r}H8+?8Oz>LT(Hp$3qX<70CjiF>{K zTng59u3t7kfNN+g`Hh-)wbwRocEtu)+T(WG)}>UOAZ9TkUYhiXl(C6X5u&6Y->a^w zM6*@{k;EXiX~AV+zQq(^Bm|qLW(a62AOv9E#gR6J>y?tM)0&6=y3=Ck!(8yz#(YaY zj?$a4FYm+;BWf){GzIS|Y%-HofL+HyaTLDp-VHqp@!sykPgy-a2+KR6!ddC{o=@O8 zOrVnLM+}x(K6SO;hz|?sNRx{#Qu*h$%BeQN;ou^-`TG6J5J`gJCMZyHC_WZ%3wWY$ z=2T_90(D%H3wH@Kt3fC2=f>etCf&xwCH(Z)&vQojq zY(tnjxgLQl8FGTF& zl)em+?Ly6r067~rMs-)1xGT{pbZhrR@|W}s3)Mc8g^HW;v_sFm{=K%d&D@4S>&Y%R zsMCUMg{9UB#^>qc))c=Rn`cI>j+W=~9l=mgibu%x0N8qA-Y1-PptXbe4llME^p?Wk zrjwRBaaBOF+Gu4Pkn$V{Dnt2w8>byWi%dj5QjPHy;(l4Jkd?JZ2l{F=fzumlEe1^gpP54S54_fu;_PqR~tcgxPfoR%g)iR z2ewiPk^Np^$qG>{NSlRU?4Bk$c%yX-N9%p<#}i-AQ@VFAkC(-xi$}_!@*w$qZg>3% zgOg}us+eX@0Jt_;XAMt#r=EIj>gzWdUXS{aY2!w#7^1V4o8%+UUQHPK3Hf+Nbpm9% zo=TRxWp{Xy_{ode+*#zRfFv`iJ2!fU7O%h{H4Gx!E1G#i%Q5drr4;nATjGnw_K*@KY{qJzaAWlPSxtImc~9 z>Qj{9x2T;2%H^@FN_hC?<3`mZ{>?{6s)W$nx1Kbc(AtxF@e(uAP2&A-u0XoW9>wSY zW{J62rh6Bcc)13w`D9xJ%8;lutx}*V+@VcHS@{u{w~PDU&RNVo3a{QSZ`{B{96cgE zqIBV037QAc#F|~&*$ba)y$X^;&caF4%Ma^n-4%ogNqL#^BwW~WCeotrnThd z-fzcGpPzLrBTtx(CxlwQmMdOe+7QTNhfbq;a|f-&k`=z=-97(&SCFTVF zK80`19<>>FRIM^*4TX!nddmp;sq5rLS0locON#PQRjTfD zy!xdNfii+!&~9B(a2L`T==qbIgGD6#OoCVX;5rc$&ur zgnMTS2y0tl%zA;tU<>H&JbD;~_C0i59OePu6@#E#(T8QZ^qtsLQ${0_oNXwx@{l#N z9m%}2@0Y9_Y|uqHD+*AZk{xRPu{~78_E&=MlV3lJ2;X(9A|)M&;Uj1! zL7`@01#z{q#C+Lx9YG2_d7-L^ZgQs*TKcC<5jsHS|6u%8N3G&vE0fg~ij~yn!{i3D z_tu67mDw1rU`TU| z4>;BnI}$reQ$Pzs9a7{elN6lY>dIBSb0mJpGIGXpE{mF8rL%jJmsAVPH(ZfJaT;=d%3_kY{T z;r~A+k@G)_=ubP@ll>%bwV78Ut>wE|LP_Wgm+i`v;Y679NV|Yu=UP#OxR8%@K@I<4ND$TEqkk7#wzHyP`K$X7!MUHP?0CXHw-<7 zUgw~`1!KQ#Sd1|s;L;Gp5?Ct|Sk|>zRT9dd*a^45dh_gg_8B8yTi}969bMN>Vdi8; zsE@`&(+=qlBN0aF37t6OG7IQPwolW~$0891^s{7evVZqEpPeuShw z&3&`83^;@0VJ@BGEC=yw4@kg+<6HsV7m(XSK23aj9`}M}5umgXFuw^my`YEz*&u3@ z+v2IoS-@~aV}ifiCaXhF7URZf^w&n_u>AYa$rHV-|B!=Yv21G~bWXcxsd3Y(IbwHjX4N(4Ak6y?r*wYlT{1Uun z@Zw!e+1+CNy5R829pc$xJkGQaClM~RVi}$JZJe z4ERS+GluvbX?zMp^38WqWnH?RGi+CFW!0(O7E^O-hr0xYEm3;|hZr-M54FD9U{d!0FiM%BZRmpd)N8Pb9X@qPAego^lY8RVA&%;5B7TDVrj(0^i0|L=^rg$FHIga^|je5x(R)v=jE#E}kAC3qZ zbxbt?PbujVN}wVN8;I>xbk@|(Q9mK^mqI3%u*+24=)jaIU(I)A=pdRQ4XxX#~qk3(yyu~uiW)iSmE*RgMVF?k#(qKdX+AoW8>6e{Bc3@n&m_FFQe6DBUZ8{L+sNvTBF;I z51DThEEqka364S(dkuoXqiV*^ulC(B3q}bm+ z=&O$=Qe`>WIJ-+}1M;EVpnHoox( zmpFZ5a4?M4PJj41LrU4sNk<1lnHJyBlt_&RkkSKDJz24VUr5pyWHD+Lzb$cN5l=+5 z6{F*!a=-T4P1ftsg*zml6{G9gcXep+(0f~6 z8uwu#=;_<$U(L0GhbJx{bzeegWz)x_(WalMSI|~37?pUb1maz;Zi&Lqt@1UEB6^KU zOG{V^4j5TIZesZ8Py1ZH7`1a~of1atAduv(b$DVL*udFoIIPyxiovJDM)M753EgkK z=-m{y2EJIM^-~c&y`lA3sPgfMDhA!oA8%1J5{9xKC%Tpv%ybm z%=GVxvw_c)H|FDhswQ@NwI|aRu=h9jX}|gZxCvV>&v6~SiU%PGyzhU&@V*RPL=*aw zN9)_bFj|PLaRrel+B4}Vz;UGRaLR}xyVY#fa&70jZRAAkHb~$Q3jcBO*yGOg zP=@E|#%=5s&o96o$4%4GmJI_e)Nr^ua9{TK6gX`#ECD%O4b8XHxeLZCt(%ezz(2uD z=m>&~3^ZBc9L-Ggso9jLjzHKG2d2K3Y70ss`p@?+!@F`f()?ZLR|2W$u+O{*Pf2+Z zd;LPk-RZ@oSyr87?bsdS+1UVqr3?5zSgZtig@G2A}t*G+jaJ zTVKPwNsSt}xEt9cIBAzZT5mJqr%x@Gq@u5tm%$w1vt^}QJS?;({;k=pBg#j3A{G3J zVZ_TQv{}(LzDch>GyCUCu=WF3#C!m4$y*64p?wozHr15Xjad}rvLHlhA2b=iC7ZzJg}<`%zkrV-}}j^s=)*+%2!I0U-AZ~^`derwqTqQEP-7xAJoRad)A zcAw^s9h_vI88jmrwPPVZ2O1l{31ISr4p-F*6Au>vmwP z;6p2_m*cSKBgC49WWg?$=~({bTvxJhgHM|W3snX#ZRF}w2PE7i5=U=GWSfz6WJYk1 zXwQkQ$>Asp@iI;MZ47+A&~{~MFsdy=A*Shvm45-*Xhi{G@9q%o#yAr3g#bkl=Wdydc>H|h?G}eX&>Sgdt7i^88Zu2_ZkX&wWc-stn39g5pTd^EPsuKWvM3p!1kJ>pOd5f{{ zOWa92ECWo2>ck;{i7&)7cpnQRptk-OIcKw|nZMEe@ky2%SI%qoU89w>`)^ipk!jqq zYnS@Q2pkq(T*q6i^vNc(llKyR#!HoWu?=?2N3z$bG=YI)(FD*%F2~a=gpQs&NK1qt zUpX9qdc9PPBbJ!;^dnC9srh&Gl5_}}f^!nTur?b_vyF>>!G|qRXg}toshCBVQ{jb=Z=rLH_uoT?@raA3`30d z4iE-?a!nK3$xM%kQMxF8sCax`FI2b?uD|{k0dM1`1>8iv#~YBxfR-}r3Xj07ZvtyK z;5zsPg!)B?_ZIL8SQ&jOm@2bk{^O78M<0vt$alVpFscP7f|Si>Lr<>-;kZIl3&f*5 z$1w3N?JC=n=%laV$|J|2R|?^B*JsJu_r5qDPcX8L=nfdy7hVLx^`-l(B>soZEq`3+ zwSGPRAhsiVehMV z-Zx@fh?xMV;J;F+DM-qrY~gEz;-;8mE?#}of#y01SE1I*aJo^y zzYo-t#H$S_dy%qH$H&(R@b z{yu5;QKHJVwjA8tuSc6xcrRsP*JURvu}xtm7(v%z%H#RNvq^Th zpI%;6P;i|JJ=C`-J{D&Cr7Q%GI7^G`3z+p&?N~36U<9SCgVZ>v(6C6jQb+h(o_eD< zWu+)BXY}=+$`@Bl@c9Rp_W9x-(ShlWeZk2e>wK(pvojA&#EwT%_BN9Sdvd_-Um(Qq zL(X+`Koi}8=C_1bfr&GX15r~=>PJ1XDEzV(>FZX8D>ad8c0!DNek)Yy6`GM#w(yDC z!pg`vRE<210$U;{abhqAM&}#}eg5_eTCl@JmaN7dE{m=zgqu#EJp@Vf$RPim$jY1G zL#I6QsoD(?|H*>|Ny1fJZ7-vB>d(FQ2?yA$y$f?^Q?5}8Ff9~#YH8S9HKLc2yG-2E z)_<6-5AVaHu5w}r&Jl=ik#Ntw)U9ZWGJzLU@x^crM=IAYr#G-_-ox75LV?c{GT87l zr_hge^Eg?w%QA<(>@ZWKNblXZKILdzYP09epA2Oc3FqEKSJobtuqpEHjEq7KzT#fJ zQpoDlRw&KCn%1w77QfTTNzj>cNmY#)4Lu7D5s#x*7%Mn*l7~{)vI7Bx-9L#I%LImX z9^Jw|2Ja+ZS)uW+Tqiss#**iahxadG`w3ScnFDfZo~M?00+qRuwGwHR>Lg#M>l-fW z(jxGKHagqNdjqUO9V_TFmVG%C*$7~-TZ$67L34_PsAR8kuUi3rgwznV^n)E(6U-^U zY+!RnKhJJ>)K%nvmLXYVkQ1+Z@Pi<6u|KJW|I{8)ih4Dg%7@kXaiQ!$7duj21p5l; zBladFzrL#lTRT6iqLS5kKrowJWF1{FhOmcMT|pC_??)f&E9$O9EgiaPn}6OX-id?@ z3``MiGv*YyuM!^P6LZLY^_(*A>Q2@s#R1+oaG!09GWY z&=&waEOjk9Wbp<0Sr2L7&0C>|m|S;N@*>zT8jrn%ea0ikFooNGH`+P1l=~xcSBf?} zn=AGbMvh8f8}Hv%yZW7nELLzNj%E-!4*`FE9^c9g=k{%v zD^DyUiN* zKg2HZ$`0-IhMz)&74}rR;AS7AJwA|~c6+l2cD@{PE&kAT$YLW@F3haOPy-uB8X~xT zMW-Da^rdUz2@V&Ex9q8=Z&{Xf>f$e2mM_@q^n|aputr~i<01=s=zHO|w0KLO!l2}b z^lJJFWo#)x`ov1O+0J*5H+WiPp!H~leJiRHFP26ALQW*cr9+=vnA5jtdp~$fj?bet z*f7YZEu6*Y9k(96#DV0NEWYD_1bw`%8J1Vz4#y~*GJ14qGz-T~zk-a}In7rbKz&{# zECCP!4-gj{tn9OhVm>`#>aguGHoZQS%EbDZhgZA^I0Zl!TSg)HL&WYU!|>PIgnFvl zu-E!la+Ccqzt!qO-;3k+Tex98<&5soKArM1utI+*XUW$0f&&H8F^!oXrUq>O3vRFa zim+(t9z+Nxjk~q;f>B!gI<7sJ4|N6|D}iTlwx59S>JD-#br#;}zMy9t1`)d;SdYDk zzSuC4sVh1}3SM3s?L;Jk1VHWfy$~W>t{nl5|eI5jX%$-QA?CRgSy)hMfWwvJp{N#fHrR@Z4 z5$AYpmyTmd6EAQ`z$Xj9lE--JV>^X=5kBxi&RIOFf2u^QW@r|lAbb2F&cj7K7`@6o z+GeaU0XBTDQ}}_J5V+j3Il&pNo~=+P^~v)jBe)=h_AD~$P=R3y3AFXRwgUs@hh7P~ z09mJmg9Rc1w5PMt1Rj$7`z_3$n)?#s5sBN@cE<|w%$L5p--a_0X6)c3gX~c43KCIr)@P(|S+2YRSc`v(lt?@(d+)ru)P z_amgmd_k&!fwItuyaE+n@S6c$)uXO3cF#VfHj9Tb#zE)mq2Cl1n6Ey@O~{7`omVZe z4{R4t03B*ShCMXXo!I|k*dt;48n>VCioLIV`2n#Q$GBi53Al*@f|Gn-vnS{^n`nIT zWj=J7UwgVKqoDS}hC>F}t_gJlHt~qtUzvU8dCQKI*$4LjkoT5Rab;WEXmEE3?jb;M zceex$p5X58P)Klx1b3I<8l1u%q|MuV}BaQYj;!qG0}BK91X zAaZ&n=m#I+;=BAe9jMtiKnq&XS1R3P$Jz5`@AIqVuliT% zGX;;EWHuC$ZOn0&70Vh4#zPVSpb$-iSYBA|E*t zMrlH{z1p8Wn8vgCSSF9Nc)vW_`oscKB?rty@FQbn9@+`? zJ2`p?1Qm{2=Npm8ogmq7?PHs_naMjTJAS$^^Bz|P82d?Y0UR4$*(VJ5V8jLdnQYzS z+CzaoF_Dr-+AYkzH|sEcC9UuJ)RBs)D&)%Xw;}>E+MvW(A@%9-d5G(x2(R`tP0Oq$c zFhrN{d?wykQ<0|PVW+j?9>=??=P`ldL8uth@9{DBJk~`y*^yu_LCKsk$7$ZVo`Dox zIVph&Jj4Aqx4=>7|48;77ClPvJi0FFDgBA^DiH&pb(CQ(9xYb{WX^GHDZGJo61^s{ zg73vgvk-Hl?g_e1hW5>t1pq*z%x~3gl@6H90swktxqW!D#F)UyTdLF@cPF3O zi82O5?fRbj7%bkAHfsTD9k(BWCEj*49_OU^TML$OhoCun(6dV&u2FH|>6eFi&-;~z zlaFBEk-gOl+6_GW#}&PXBbY?|oAJtXjmOr&1Z#xz_O=JmdnhGc|JGcQqiK5;d!Ykn zKd*ErKN(9h5V?7_VsX}uj^E2Czu%&a!NL=>{~l#DvvhKH(YLX5{pz2} zg_8ov-J%3iYoJ;60^(rY?KK?jFDkG#)p0pkC~R?I(M`xahdP7Mo5P)^uS2Q?KqJ2G zLN14<_U_E-sz+v@(of=ZBdOg>_RV#dV%2?^uc{gC6!xBY=GK+(PxOStz>b&Ur+du9 zHXDr>+IMJ=I_omFLdq1vr~-1t|uth{P<&Q zdPjws%qn`09;W-_7L&f;o`2!gX{ld=f~cRellOeZebGf zYc_=6|Gsx0sju;PoXpZi|iITzY8Dd@INE}BkkM#hw$VS4bVm~ zoHnq%{h~JTz4`mEHO=AGt=UUss9cNk_XkzXq13HeOJf*a1KYg`E%)JDNQl&B7SFBJ zNL{B1V@ru~_u&Da{93{oqYnZLiUm5PzfEo^lN{yO;`}i=OFO$CDG>iqyw7ld&Z_+9 z1XoxVJL?FN zFk|O;I1#Vc5Z(O7&GNAHJlM5AtM&5r)zNez_qU!R;P&9noeGFL8-cPZm%|>$55k3& z*al-iL=0&T2wh0+9>sHm(s}~t&obD&ov-Ou#C>inaj0fEQQ4cP8!pv~9Q0$0Z!bZz z5b~}bt;T}bXo9;p&o(TP8~D%MsS`7uzw$l)p6owV2o?U%S-Sq0mHe+Nv>o!lsE~JW zoIfvl*w6=FoG?}CKIy*C6iJA^G^igr5sdI-Ba|88e_JA75c>CDOZ12S{JBK`p+D0< zYl{Ax^8Z);`8Ub`i~hudyyZAsQKgRx2laHe-{vLb@;v)Zc-Jg@p(nKddv3b3=)6#Y|Bjd$%s{d-kF_x_FO9O<7uQRf%qi5mE=kiRYE~f7b`(@ zayXO1xeC2W`!r(NdwZtZ#JGAXm4%sEWa7jGHRB59dAE{Jwl{aaI~Sf=A~!c9a7Orhkttd;^O&J~Li?rdeY0?liZo8SP1?sX{XSv(t|#B@M5MKKOxPST-36Y& zSOw%Y{?+XGyi2xA!1D)(@0g?h!rK9uFD>l7T@kaZ$hu+1l{iiiR%c!X{pbRA-sy77 zuHb+@k9~`Ns(U{N(?L8>;z3QO0;&FR6FRuXK zjV^&Q&D9Dv@QfZ1#$LV$%RqHj6?G-hS0Il59CUOL#7cVavbp#y zD-Z~#1}Wbl?DF}9dM1AWNr8n1m0^4}pJ~G?3vUqmpL!e5H5$nEkZJAOZH z*KF<15^TI<^DM4bgy?W>sgSKE73_NZA=nu_M0Y9-viU#uOo@A-9r1Vo37+5_L@ry; zrlLRWd_r(WE~56(JXNeGH(OCX-9sYf?gDiBY)QPsYS;S;#GUl>JHX`r;YW}9;Sh24 z$9p#hg&P<80l=l>Th-(DvgxK*Wz~06F%Xn_sf_|8^;GiDhMP6XVkZ^9m6>%$x@5eK zh1a<9Rzq8(UgBU8o6EK}S>6cPy%ih0gM}2o;(@eGRS};TjXz=-$Ss&JVRqvUJHu_n zLS5~YwN){ruQ84uS@f$H>GK~Dpg+=_E10ZEcTz`LOLQ==vQ11jrTG z{N0#oeQfaK3B~|PjO3q+sY`^$2o8a%m+Wk>ir2z1#%*ib}9ny8}=LfdHH&@~ul7Z8hy>_Alw z@K*NZ^m>V7X5Mf@8DSl;``xV293hXg0pZW|hTe(v2p`$O-*ftwhYK}A->f|13Svou z0ift-JS&9I+3WK{N9_l~sPIWS<|k+y?tYawBkwYVg!vXk40Vv~AmUHl=~ti4E590G z+<-j}t6#c6__y&`fu@wFFI+HPF%chuwtAKkeFFBv3m zH5YV?&tUgaBUW26;Ao5vYP0X^ZnD5vJlw+%O-LOD;bO*kFdkWyaGg`DeNr4 zQ1(D^+ch%fDZ6(H%@`78eEtEV>_L4x#b`m`A$$+P^1R4R?|A6-1tM0i=+ittVR~W~ zt6}sz&Jc>=2fAc0V_cO9>WjYqFegTSJdixBDeP9c0+*EEE;a01=8mxjf6~PcVBq!}`Bz8~He~nNidZVXOX-bl^aTx2RO18m!pti^6d1byL za(_oweLZ8Ks(f!@XNq4eH0sf1qbbI^GFB)er7M4>Q#3d2@CL|6=(SAe$AH6yz`3EXWV>fTX=tw+D*dn=Tl>|zQ!4i<)f+<4~> z*c1xRP9JVs=#Wg@?irVZK~HnNS0Dd^veDwCGe66eT-`An z>t)C?^EI%a^~qoXWunS8U*R-WPSon*@KFgW?~&p=gsQ`nE^O3ZcMa=AY;*8~Fpkf@ zJw|#gRr0ku1|Lx_NPdN(n77d5S)OWuFCo5=iFfhTP(r@A%!}e#-zuC|6s* z&TawTN8~G%WA-d12ml{;35MmIUZLcH}uyUZ%-_V$_3M%ZPT zD0Xc|9aUB*@pj*xjW0JxH_$I-9lP0`^Z42})spns^+7N_@=Us4JN8DeM)3<2n8VI4 zrVXW>0mi(pGG>*S(0+NGx18;=BBSy7B`t(`qz;09QR}L4-+FQb6&=qVXixibmh8dT z2>O)MDPyC})BZuUQU0*&ga^+;M67+#qyL0{!=L2B-~}WSk-(`F1W#1_MY?h<{2)5t zT0*I&GLWR-aaKxD6XOHz)EL->-mNy!vA^#=QtwA^3m_hnvHC&sWlwbSivAt!q;oEK z`77A&g_EGYFB^wNXd1D|_^cLVjF4fMzaQFMMBQ96!X4xAxK zhP`a#1Igc$=je*q>6UXmL>EsQ-F&!u@;>&L$r+328_6fQA@r(7KfaL}JMro_HIC)_ zc;mYeeesfGwC$9B^Qtx7Gr(2&#B?K#jBe6n3QP;-TjPioH4MJW?r$Eb&?jmP#Hb;? zrP(M~@{8)yFviTYV9N|Wnk=IK;6 z9+kei14w_fM);_OVWYh6F4KVL;$1PTy~g7hAuOn0`8 zcF^?BGk!r-J8z@RzT2Vun(J0{i_BB5*{L*7h~j~lAkUle#*}P^Mv!h7(qh2gi{aVG z;j$|XA>*qm_tYorQ#iJc))yMd^8vteV-LNAp&Rz9$uaLdoZG6dkrQ-fR@Sjcp=x1o zVKK6!e(=o6&59psr}KF0juG{v8h9FXKwfT-aVfxZNjy&u+961eGRfx=pJUl=6jq3+ zG-i#MA}<0Db8=`p04apNdh@6W1%QRua=a~hnD-1a zG1|zyIF51eeM?|{rp+?kd?}k;U*5E1E5xBaYKPy*`6LYyt(sbq3f#RWm zFSO*9ezO0ZAQ?u+QTWn%mR<;90qnKP4lz1x+f%waxUUn4>iF?i9UiL)>^neR%mMpO zyMO4@oX1jMUbITTx(25DuC(6ckL90y9qxIs)6I578v&nT1V`B; zEw{v-kdM?I%aB`lfg<}#AjaqDmzmg2j#)Rm2Af}>?Tgg%YW9EG4u=~AGd6i~u5e+z zj``kxx5IpUpIHw38KQ9kc9+)!CEwan5PG&F>=?f7N2{m^W&PxMIh>N28Vc~OXe89n zV!=t2Xx9!sBTc*uyy;cm1T66*yR88$L;OLl+bJh&Ipq2J_PfBU3&;WfX(s9=Se50n z?CE+r@odTN0eo}H_UOIDPjb;d+I%V?L%WXNVLxZS!Vg3VTA_AhZb;Qn0rl_0#U@98 zJ2N=UOgzE$Azy)NIk_HgoV6OLx0mYZX9I4qfq3i$GjK_ZK_##v2O2*#Ru_SgAwtSu3LpicWWyW-AK^%O37-bABs8zl&pu zvV{7BRGTkm$DuD|jpa?t^NaI0&%0C_h*)U?!P3;H#b zvoL2dAAnppt1@0Z`3*vyout7S`1O#557dIPl(XShK3(psh16*rNj)but`%6?~Gs13# z4XM_n0E#}^@|pJy`(@vmvEgd6LGs z4+6tD%oK%cF&()_3#mrm@m{T*WsFe!V8v!;r^N|rQmwb9fzT0EvCK30t_IshNmpob zyfm0EeRG!(0Fo(TyvA$OfUoE4{$U(Rias_gvGP(GH7@(@KLmU62)$nZgkli8%f4%kpe35 zCBzn4k-2L81{#mOLmeWo@ng?5=8hX-V2nA^Le6BPZK0f3&Fk?h@g~hCpM$*WFoUqd zYFD7l1Nt5$!pd%XveJ(_qB_m)EDrfmcB*uC73Ofx#$DIx*L>CouJ*AeC>@1j=2W+$ zF_6e%`zzYVr0xjk`#BuRH=E<#AXB346W27Oog(gf#e2Zd3!xlZA_kscsv^sg$gw)E zI3Bar?qKo+Xnp&e;|GY10A)^iXJcX6XEpJ&cDEvux@g?j#dk`b8kaE`+l}{0;-MQb zG1Mz({zovQ;h>cZ6 z3I0NuUew!fY@TF~%suC(f>lp+cgnmc&$Gra;HE>9Fofpbwp3o{JAypvF8h@Y1QVL8 z(~ryOb*CV@*SuDB4N3#>V@-M!dSY`4lgvpK3SLdKm_kPicPTZH|PePSJSj6sn$2d#pM_^Tet zIJ$?t)7N~I4J=s*?mjUaG{-V^_iG+4WoIY@s>5}T8H-#91~=iO_mtOwj}2iL{`aE= z^b#G(mRkC5%uTjND?_7L#vR(0X|yG6p<7}rSK8VcAEG^xkIieH7LDMzXY-#C4|$I9 zE#9NkUN5;UUcjwj>=<+FKm67#cnOia;MaTMkIKG3-wVH0_BH)1F75h1I`00bMzOzz z8-DriepL4TUpI>FihMa5@Nz6rq8*cnya;_;n4XRaq0lskbJ!30troN_!yg5InyAqs zcPM5ojNS%oo#7IsJZby{Y7Zw0Ubzc8Tu zfq{*st*s%;zYb~8*O!pjXZ`VE=ip}fhqHY8r~kP1CnklHHGeqEgVzrBRsV376IZA^ zs?$bNkY*Ij|KTj}lk!H^pC>tm+Fa_+eX+LqsPx{}Jv2*j%QSsFGf$VGwY+wT-+0W{ zyZ7GB#u(?cy0S|5_04i2Fw}+Ywi_+ellcCiO0qzHvDkJ}+K*YzxI*(^4sveI`I7GE z$O6>1*WkJ&TbLMIm><`F9~x8r8TER%m2IcGo@=cE#;pOS|2CvFe`|0BnN<<~g$C9C z0Sz}p8`r;sb{0S3xcV7G$c(k*od)0zSjHF1bJEX6RixF7aK0{a@hw^U7FUV*o;f6? zmXms!0Kf2~>8nJlBD?yplhX8%r&M`V_U6~KSkX3IiIgTf})>*iCDM!Ll z=6}~)AhscnG7JRTMHoIpyOqckOsJK7Hhtf@NXI9!DH<1h#mO!fEqUXjCy4hc*?#4T zZh=5y!ni86Y3d^g@eL#^oDI2k{@4T5{fv)cXJRE%;p4o0NN5o?-r5(N5w8F&qt8yU zD_mS+o@s%3sT_jzTccHRqaHGt4QQ8Gx;(0#|n8zxH|7&oP08GAeACGe`W{60yq?PJKh z>!rE9TEI9U1NzrM180gGZ} z@)r-7(X}Hyb&VFOtVk<@?b>Fs@lc|1k8i}n7 ztyxPzjMuJzb0{tLsbcsGLfBnIZBZ-7r#gkQ)YmU!We1qXl<$g~Oh?a;0M5eQ(>m)G zniIGWEtA2Y&A8z61nSi6*riw-Ps-|TM=+aSa%u^V?IJ`v?g)ifVk0^;jn(T3!s5vt z44YW87^cH{=4t9e$EDE3x4!FH9d-$yA0vB)57YhCZ36S1O}5ln;z{s0KQ9yNPd)A)fKZ#7!@9j^pEe!)RC@J#?OE zSn=X?6a?%8IhT`pQs>cz60-YY*_FgcNF?GNU4;(~ujQMp8mkCnSG;=8orPF0uxC)s zwFkZ;~<}KTce<9^X`hFcT$qnglEC zQCl^=4>(S<^}9PNs3ZrJ8sG7#XbS6`d@95h}#yoEJ+_*uepcp=<_P)gxeq4Cn2iW zJ8-$wyby7tr*c!$x#pGWelLMwM}|bZ7A1q8%wG1Ui7@gCS-;9Z9-ivu#^meUB$o@- znPy3Ehn_nhN}5Oq+D4?&FWxWD3FS;)+0S;vjJaTy%;B$ne|LQBDm(Qw-mcOJQ9)Sw ziejEZir1x)HBX|4Xy_^@UZlaBuG#;l4}#ZLy!jPHGg< zTj$c|zUt=zT{v0MS5f>PnfsKccIEIBw%-vHgE4Jr%q$XA;0oXCdo?eQtrkwr!x!0; zGI}I@oZ-@di^#2G>Q)Uh7^xQ)L6`{)!r3I=XLxpp@@i--YP$O)XQYFc5^SV?xvQOR z8$tvqpvt$JCDh_u&n54c*0{&Vj0R0bMB-&CjHg*GxYP^MSK;fjx4)%&WPi1$`$Z&w z)9Syvsr`3aZTUs33IC#1sq8ZSPfW<9hml%%ny$5olU;8jqIts-)pUBRPD2$sF+Kud zl25TmT`*2mH5Qk1X_Y*`_q!Jx;P*3^c^h&3;8l;HAG|8DF*VQvbQ1=;C*L1$2#!+` zn{C^!PU(lgCajD7-W&Be`~uz-B_I>)ki_Hy6kW^L^{R_MfE0&uU)rb`HqtBc-82K< z5*p3DdY=BjUBFT|AvjCXy?KBhq~k=KrLta{s^22gp{_HuzrgsLRsR!=|Cv>Rf3d13 z@3be$fNPjq7!&)3MK^ChVY82K(x7_2#P>t0{!akug*Y@Cq6BP}vUjx>lGqGeVCQIS z66#IFLPwN7+0bI!r*ZPD0pT&<$2-B3xt$ANkud?w6G}ko)GOdyAZsl`^qL|q@wZRW zLRVjL`D)E6go2KruVTNx$$<Sdxu7+-JgPx4`5&hw04YLQPQ+Vp6)i66p5d0Xf=F3uWYI)WiP zn2+@KFplIKI;a$Wqc72xlbKTB+FJLSa7|pqnr=#zeVoPmfSdl)^h@;g+osH~tTfwgNWn4{} z=Af9_wpB`)<*W{b4bMih^0U#B98~^-;3bOnM8Dl|5U-q9E>HJQhJ7~fVt)o{QM zkNdF9)QG7#sJ`Pu;t-NAAaS4RY$+P;KcHtY7#pp7Wph(z{LOPMzxJVYJ{BEsM&cr3 z9tnyza^IRLCAb7$TmBs$gWm@2yR!Tc%c9TDSkXa_Ats@WoNCo|l}296Bz)ynWqjs~ zQ}PLb*^wS?ujfZPF_Ao0+R?DbeTFZ&_FzQWs-!73{$(2Z4GHqEf^@vV;K^PGv-o6Xg}^^eq4e5=TcJNJ3NK6D$^YXUwMRp>nU!V)5!xn-T0vdWlNbUB~2@y?j=+AV9z}h z>Wfm7XcKoB=UP_-PHu8t5k6f9_f!mL)YxcOy^H{`oAvXN(Com840VdNW%DKhypqpM zLMWI1PFrd_%uTPd-m4_Xqqe?FZ=GT=&k^tgj(ja8mgRm+c=^coQG2M3)ZG&nUEqiS zjq;1Aj$mqjEbmK&mxYfLQl!<75W>69BUNHStHlFAA*8F4Yq_;d9oA^0OVCFm1YmD5 zd^>{lhn>6CrTBBhd0jF*97;WeO8J9e>_)O%PJY3iT@iRO%1^j8)55!GnqQ+YEkpHK z)G3{>*a(z{G)$M$q1M)4rf$xPe0WSh9uAHo^H|DuQZ)=h5E%wmhf{9(99^~eg>LN^ zAo21X&*)iAoRKyHR7ILz-d)YLcgOp&r3>4l){Rgd^?1&w zurfy$o#p}BKn{-c^l%br=BEW3DSl((${sl-pS`k3BtgCc`Og%qNkyY8&_MRv{BKtF zC3TVyvw$A4h#=fYXZ<;gEg7mT)k*v2R#hUNYIVW>Hx#Sf!To83=qk6A1Nj=>d9dCB z9d^)`T10-#*mMr8DPgpkYAtr#uY#?m7~*0ole)&t@6?~qz21=FmwS-r6=?Qh0 zJuC}u*jq>>f#N3jE2{Hs8?XBAiSzsqKa{{0TPwTKKu)a_V>65uyqM@aMR`3Qy=L zF}aQa9Vh~dI)x(J3IbD2bq8Bx(y)<$lzo^?0tMQGzK}DrdA1tH*Sxdn40(RW_j)H? zn`n!0>j&Nxc>E6s@eq{PJiG~x?DgagC3PG3Vd^eeb}H?_?4|DPp@v6ts{JX)b`FEWlD*I%}6UACykuYS)|Ib07u=E5+1$V_yBUGYu7)+$!$T#wbZ*~ zqZbdu82T#g_8C9Hu9_fWn4YmHE;lji!rvjmYxk8Vx$cs|G5x$vn8u|cOt-movP57v5Zdmb@iw#iO zCuHDHPJk42ZoqS?Ia45kps(AO79T8>svutc$YaNiP7RaFC`>1YO0*IQ(EW8ve4Nwl z-)Fw>3c}3z5#(i6$HwScFT^D^i;|PN@a~D)aZ2VVz3Z8NoY@U92m}Z#MNO!reotEf zSWl21>)P8KD^Ii+>B)e}2=n{>6KKQ3Lo6CX9?tSag*e+jLX@w&p9v}Q?@s87OJ!n! zJj*(I5N9qTXuh%$@n^0vHWA56WxE%UIANvbj;Jy(3r@Upq2q_tF9k+eh1WO$uDEqi7pGnB8|*jLdiy|D;NHlIt3w#*%h55r0oa^h(ELu$ z;ua%*KQ)!MuLwESMU9)AE(lA?U|c-D^&v)WbCa{|mI8lRCH|i-VTdK zusX%0ZX=e8J*IOZ`C4tkw{;vM@$n~v29-z58)hbif*F{C#Ba$%zpc^OH(mK)&KvBj zJc=P+z0$^?nZ6W@vQ-9`eloPku)~s6M_HramRY{2M4@tDMz#|?C2Fe0DXjCM^0<5K zx<{$!Z>nbFeq}2Ot()IwCa#QjivApUL}wl9En;RcTS(>sZ0@EPqsK5=NDUR?#pp^r zmoR2O$DB{V(Gf&V%dr@DW^*(Zx(&Bd3XYV*#LJ|_3zP1;N0^#jQWW0xGV86SEnszg z;hUSWh>?Hx^{kJLpcxm@eCNBu|5iZ`o>r=T_%@l83ivM@n|IhG!}>IJ8nm)x|yq*K{|& zka&RSejAf9t`5w#pbZ78-wfLxc9#0MpHK9Ll+@j0kKabe%trDoFNO(PAZYU2@Q~eP z6&Uo*(2Fn2qJ{f<1pLHsfOZ;> z?kPOY8WL0U*1cd3)u{9b59!t#xl!M)y(!v%n+DJP_BA8Za=th69m->$IR}*!Z!;#CuBZrVfcRz}!;s(^ZVz_x1wv^RS(L;U9i~(T6CvUxLFb( zmI)C8K7LTT3F=RNhJbu*g#jMqC9B(~< z+f3W}3N8cW*6-pROeQJ5D^7A~=-tPYr_pWLGeUXxMN{eZc37@UZ3DA`rvS;+>!vnP zB}~KIDA5R=tAv;|ihZU#!b8SRMOv>-V8?H1xU|olLVxl7|7IO=vvd8^C;zKq`SLI8 zX!3XKNb7&Ij#jZ2P8T0RwF~5m)Sh3p8!PSG^yu3wCv|9L_4%JJuOLfzvulRhATRvK zo4Nb7?wj7{UpbY@rLxbYCBIc0&okRx@7NGn*CD3M%zFmuO3qVtGrvhk=dqqbw5DQQ z#QTr&hgn4Jttow5xjJ8I@cs8ueyZD5tNVVfwYJ+GZs2$0DiQ7GxA-oJ=c(Z@y#6=) z=nuU9X&;UMw2z{H+DGjltvRQ`L<6C+bcG$Z)~c;oX2YoQlX}dSS$;gq%!Sy?1da6w zJ@YgnF0Yc8uYxk8iMD2S`LV1nz|>BuX|FAZqlQxA9KSL}P@9iNI_n*}@08b#UTn7l zL0__V!0(kx7s=*GSF8$*UkGziuuf*ne4$k(4X`PBZeoR*e1j58y4|w{2m8q3ZW1a9 zUwS@h;O3F@Eq(Fb#|m4^=bOCeegX^q-D|B*1Y-KrZmdFrWd1Q03i%*O^f+#7MQaXg zEFM<_iN^$zGyHv)cmXb|l@xb!l(Eo|<=hgqq$*6~T%X77tevpLk7_5_i_tnpksg7y z+#VWn(gOAx(AtB@vYrnB7f3yiN4X`8y>HOe!lXuZOLNxc)hVloJ<-x?4pn4tX zb~T`=X3RuDDbs)glq+JEv_semYa@;_w+q zD*o6B_wTaBC2XbWQDWu z4ekI8lB>M;MDF|I04p59yOsU2`gm0DWpC+`DQR|JH%nzt?Uy*@%Z{vEwjyZRpFiGB z>T9iGpu|p$;ihenL25K{ja_sWH~G1^W}$C}9*+82J5tumdt36`)0xiB7_987Xw=$H z6C3@Dnb0~u;3cw~Yi?-k6BS!!`)i6(Lg>N)rh1U65Rxm%Nei@@O4bqk>(G-RZl}rhOadrZM1;Of#JPK21>!gB1W> z+vr)ckwkpJ%+d8ZP13k-wh|g;!>-E5Rw`&J_u}~w3p}wBskc{L>G+?; zqbC!T@32y%Fk}OR)eS88RVf!ip%M1ID3!w}QZ~tM4dATzuTS1_7t(m35}AD3#oHm3 zZ4w3B##EDQp=WD}E7RoyH^M@acVJ5|-1!!L8VL8oT+TVizIM6mUlx|KY`yJ}bMCD? z=c5Wzp8-^1P?ItW_tbB|2y(hTR}+2{`j#9jf5qV0u9OXvEr`-WqsUTc`4(zaEp(fY zXf9BdRF2!K@R?almW)&{MblkzM_8SP=$ftZIqiE42L@-QjA)CHX1KEb{-kiRG(0NP zCVE50H{UaKHDEj+F>~!3?zQ(cUA6-5FB)Oe^SgJ^fg-l=BKUMibirt^Qcl;_4I+Pd8r&-s90Th*cbkYR!?aR#! zB-GL+i{ED_2qr??vN)TKdyS1_`*>muC39@dM`m#dl-)%l`gV+ZITkjfG#I5#mNYs0 z$$G;O4ka6S38;h2Oq8w(qrzhJ${QVtr3*aNM_tIe!d9lR7C9tXjX)I(JJTgJz0eqj zlmhV7W(14Wc7#T~3zk=&p1uT*=d`5zv9Ds>>Snx>m5h4h$E%}0utxnji<*4Dc)|Y7 zP!b^wm{vZjq?sZ=$v)XS;)+da@Qr(Dqvin-*cvS=41-MqY(iZ9V&8p?=rIOG% zRG8Gq1n12BwtZp6>uGhKpI!;Cv^SZS0=HK`#jKb_t?DQldp!kU0cf1Pd%W55w5LYD z>dtfweCen2kcCI_o#zHl7UQHX$n#Y(53L252Z$!Wfbx-&C%ehTNzC3{KmiL`UGs>6Gz(BK+ru6w}A=MYTU=rw>Lr&^sc-;JQ=0P7N71~@M<<`@tQTDDVm#I)be)% zZFW*%o|rsBthcJ(O~Iivy9TT|@DUskd(FR>TOf$wakF#& z(;)i)PDS7P?ne24+Bg3>{eRjwZ-2&3KtD@PGh+ap)1CzVGI`@Ewoa*vjh1}`(h4Om zFHUZ3NQpxfN=(Z?MWt%ZS-2?E&sw&;RhFw0=hjLyWw2%IKhchtqkl0JN0`RV5%q2^ zOxJqs`$p{Y(>3Je$ML|PE40gi#yNC9Ar%t(tS*dbv~WY-DMID|1=Q!6lwlEkxHu$) zxOhDnkVSVC^jNgaGvgZsnb1)SZQAEP%lt1ReFGl-#+rD%LL>0?(Rhnz82SpK7U+a{ zaA;ZUi(Xy{^kY>Q*u})THzB9-aHqfh*85SKpF=@o2d79y2QXMf-l=9dcm|LVJYY4! z;7!It@x3U%G}#Jh3@%=hoL0o`m?cSQr5DYl@$DyMgi#*rYok_GH^^oW8(1;JvUa;K zH`f+INHz-8*K)G=ub7nAdOxLf%iorn(#V45t63w#g0M{US$G+KZ3O!itjOj-wFITG z3tD03FwZByM82dhdQX?5?4;;7^~$B7k7DAQj~L0N4Cbq=e8Q#Z8mORN${n)q%Xny> zce^C4r3eGVH@k=Ck+4WztM|y|&W~bSRr{rs%rLLcsEls``uVFCrmOejte9VoUpK$_ z7;JY|6&E2pms6o{(cGiUIPH1M>NEEdM%H@sx#MukU|(YK42Ch~MqCYTStRLOL%WSa zQ6<2oct%3E96{(@w{dy_JMulycAr|E6W*g@tQ3ID_2NWxs16Mk>&pp=t(pcDWpap& z@b_`~TaE0>A~$By+jOBX8rS4uw$qVH`YpR*_8qTg;U&4cl5jBmj&y|a$?ZAhT&)N~ z>P6Etf_3@FnE-^WHS*+yZ(2P9544eV;b~VlGwFGGRc`N{J`&y5hVR~+U@NA^%y-gS z<6q&iWjho|<+_>#3fk>|>So|~?UdG~^;9n0P%DfaO z&021)eZ&3Y%`R9;H~=u}n(2PS$g*6MAd*yiYhGoQwD^sndtRr{=9u16CubMH1s|=a#ugXO89A8_e(rKADzI2+#tnvILwbT=AHqX%d1z#p@qd_2hTllUA@j_tE zebv=w=xQa=bH&pA zRry+M$Jr#&iJ8+V%omqwUpUNgKNH0hd6cM8Kwty-Q@ zfjR=%1js~}cDB<$UDRp7ww)@nbAL_@{9$7@iaX%XhRBRrspmvn^uukj5Dl!yaF zfwUsskNafL&X0TfUo{*SPO%hhk`j>=8knZGdkdz<{9nEaDQa(*VaJ^Tv8 zt6HiJ0|@LsJ;!<}t2*_0V%S2fKH~B*o~u-u%36yF$u@5$PL2Hq4PnV=g7w{# znN|@R6Njqe>q0XvwFoD)mg8m_@{ElkpANwt&rz`O3E101HFG&LMcc#_eMnNp2O{nr zP$D;|rv22X2A{g`$I8n&fAC?a%Rs-w?@C<{Eb-vGO9^O0zk{rDTEl<0VapVYdHDoV z`6S>~fgtDi?4X&WX}_G^pwe38a47|Yrylk$1XTBdlQy`mxCXY8s6ACZgBg}7aJd1l zhU-*ItH{u$@|hp_lU_}?IGHKC?t2wG;^UVN%D$?ff(BB*WAl5Hui+QT>`|XHWr(&` zHHuW@P6v?U#qMZytuSfu*R|E2B~976L{@9fuheh6bs{HE$+Qf~i;?x7DhZQRa%Cs^X}i%paJ6cdXxV4;h|?&~;P_IDXO;GX5vOHr zDQH()Nt5gXx0|)k3oX_t46oW}x=5J4f1cUX4u_Z^)p7m})j_5!i()Hx$pa zY~IsJ($90sY`SL*oSKiK&d0TxD=KD6dw@$R_RHVeuR^K$PIQelRoFj&Zn1T_qAoq( z$TBNc+U_FHCc_1^2#a_l7x5LP1(kK33VF>ep3SVK?#9cd*ZWX@fj9R!6{`r~t+jo{ z|EeLfMC0N|m4MG)?q}?ir_}WJ`vq^)Q=gC&>A9QG{R3ItM*{A{Z3%B@;Y6vjXGq8 z1WAOVjXaN?#?-b%_oMNoi6MUi^kT9b2nac&x%G8sySRxmFgd=_Xe#ICE*V+SpEMile z5J+Kr>I?H#3v~f0UHA#AHIv>H+z4{H`~0RQV{+ovX*&9W|c#lrbo_?tt5>d6X%K9(M?b8k1?Pr+C_=gqg(ajf};B9TJDnubr)S`u#Tg-k9$Wvyyab? zMT-X8xRp2iu_!Y5_@?j7sKM~kl76f?z;*z)2W*jTOu>bFm%2KWe*|U0jI^m~vuxeJ zC!xAhAjZ4ftzzV>tnH@Awdp~1fivv($FF0E_!?mc&QT);dVs2-jJ=o^lANU%|B+!> z;8ectxnv4%?^dxl;0FJAA{V1Vtmk?ki;cxn7p* z`TQiX}vUeeFx9J!p8D**YHlF%GyU5NV=^JlH)WL zWNew(4W2&Gt|LIy4eBN7hV^t1X!E1xGkPDwRJXd}oki96J2DWIq2rpZo+H0o*VlpczdNCShK?iu-O4_Pj?2UAA^P&9;>wt(k(OquenZB6CKVOy zekjJp?EsXC@dU}_M-Tur7!F)!?oAh-Mqr?QhbJs8$c;$Wcg3R&>6;WKpb{tyJj2cm zg$OZrlXm8Im~u6LRkJmJB{63=JTV}oAdqSqbIE7U4-=**mQ_+SNlzHv$bmx4mb^NAR49x!M!i z$47`0*9K}C*Y6FzrHVXLgmGbcSDlbT|Hr^w3iR3x`%jnhRlAAKt5d*83H_MHq*22g;YI=A^r z)@~7fF}HRruR4k+f4^TJQC@zfYksirsk_=sZJu7$PgHP4DNkSybAeqr((toe&)H4d9TBKQI|^hDjYha)xGEg0hK8ZVwUV z!I#k@`CN7b8*e!?=qA%bo^?pl+#SRd2O|xs3t_ZKLI+;M#TwiS;m>mhH~J&*J4%^T1%XR1iygFMcTDI+mIhI!?N{K;kk2MR)sw216|-JZ3Eg zD@nWY*~TLtt@Iw?x2`tplvhI_6m@FUpj@^AUTja;K+Zx?}5@>LIrt&u=P9S-9=7NB70q&79ktG zF{v3(eOEM7I~tDJe7Dq5gD77vYKkiFy19Owy(u(25e1qm+@VsRF%mAfTuhYaN*7v$ zwm>aFsKx*IobFFKjx*c%jTfU%qMbW2HO!-l;utddnkglUVElYZ%p1~DrYdP39RTK(*?fq0dXT3)zDvQ~T zrWQY@DY+egt+WtEq#X^+YxY6aO3Vn`hh8^I&vH&Zen|Tx5 zSKDb8gSUFuH;x4xmJfo`j$F&gn|;dzy_(Esw?t!$UbjS%W_3^&9Lz%JRBzHxxI5(h zzQ0pz;!KBR$tq+c(balVYQRJPP*G#QDG@pGnFl@~Rw30vBkG+TcwMkqc1&A#`*X%XwG#J?wt z<*T1%&XYoAjZ#(i-gAh~&ZG0G{%GT6<(Nqr$Q%qrdmDpI-TirqRtuM+$=NDF?=3AG z>7;sz5s9NQ`{sg^G}^;%Dc7=i4lAn5u^LE76Y@ksX@wPox_2mY_X#JQ#JT?VSfq)e zSoo3S=N!rGx3X(P7Cl@!0lYfj7!fbX`RNj_D967~@vHF41nCVI%R)H`PqD;MyR0dPP zgB!k!{G?{Mpp%DYNX*zn<5KPGA!t@tL5U(QH1@+UtKeesa}0m<+NBK;Yn~?ST)-9j zd6H35I_oE1Ul{|(Q&NK{92c$h1}lhnb6P=RYU*v;oPjv{nXpJ~pfZCs`jPy6NF5@4 zc*;aI){8C;fbY_KGn*z%TS)^IW357M)=c$sv+P5Gu~@S*TaVDa8}%8c;<_1`ABcw( zglbXL28W)RP%hir=%kGH;DM(ZxPDxc7eTR7^%E<+L->qA^<+Brk>*ht+PO(>15FWX z@IDYwx{#o|dvzOQG%j`DVPC?Il_w(_Y5vcW3x+%A5u+H!93k zwE6ejK`9>kev_&85^UE>n{U$$keUtj2Sa>_G@4;^+Wk>X%~G{J4p(XQpoj1BW@c2J zGooo|^i=lMqz^QD78}w5V=T0$6?-|3SciGOP;+6k=kG%cKb2+>@5=zuBxobPfl+1nfftUl5%M~=@4$7&~csAY2qdxJZg zs)+_VuDFTSq#8K}^c{qv#y&{aCjElBJ^(|3L(}BPJzw&hl^rUOux2xT``-3-&TcKv zkfL(eYpm*37wzBrc@+=*8}}RKs-e7 z8EB)|Qi2rB;F!e$KbluIUzxCf4FOJ)4e7_BwqM+0hF!(RR^&LXDWF0!FZos+0ulr< zzG?sGSsW{H8f!2ST~kGZ)xCl?a9s$0Ha3e!1*&NAA9Xj3&I~?~-80h@J+l6Oj2$J3 zk_2~5r;Ev@nGwo6hi_Pt(sbK1_QiC24_HE}~lp9kh( zy^oQ_5rUQ-V0R%|-_?R8v>|P8Q?>S#gm1>|F?JjUhrGjo>o~8%P8<=OP6>*WKY1Nf zpZX1=T&Ci&PKN!RSrghlLoU$b8Q4|EN07LP{`8QFt7V;(v*z(CXJVDMK1jzd5Sy$B zM%zdL_MM@G*ys8m*c{Sdb{|B|;})?9(O23?6C{0l(y~&`itLL*ZJ_-ZvfOJ$*oP6V zz75nkAE&1r4M-3+=d;IXz$(1g@XMr}=0E1SR8cb;K%>MbGHsGPjB1IMdut&mh+KGk z-Pvef1Z)td<$LfJ2N~gt$X-yXIpr5u@xK%Y{imz=_37Gg8PNYdhl@51Z6Uwg&EIeT zwO9D(;o|Ji!^IFc_Oa*WtcfgDR{Rh>WFu`ftpLpEMXZPEC|u|&>2xsYtnPDdq?MJz z+DMp*kf=eC9X~<(yz>XtfO}Q==ufdEpW&h9@HY-0m&59w4j+cjp0|{rj;6@kysAv< zmdFb-RB04ja7LU!6!S@^}F_9&HW4)>^w>pTDs892Y!u2{;AiD6E6VbPC zzpS)k@gW`V)I-d42^|sD2&c-cPos!MconWle~W72|MIqW0i?EaJGH^nD!7;5?p-TG zb_pUeO-&Z{z$jqPty@`MX%$Zd@y_IeoHhQtQWy$6x?z&!k3n8z;kNakm$K;5tWOX2 zjT%uZyQ?WuPtp-vl^bscIWfv_ zEhGObg5Kra;J1!87~+xY!&pRQk9;OM(AReFr1@n@@lZGoPB376a4lX?zfCftnDY#e z5nnPl(fo&RV=&z;)lOJ5Np}uMWHkU2HoWV1t4YoBBPzs(3hB?M1V~a1x1(Vb=uA(f z@KM@!Zf0MmT+uuzBgqO=)Kl8Z?&H)K4r{MpJ6#>(7*u@Lg?H6;)EUFWSqHGToBwj zA`;QUC=&pNQ(hk%su50wLzY*yG>E?~Inl5U48KQOt6LhDKsYAk*&BLyw0}D z-brR&;;u>8+gGGv==QY4ybw#kn(&a)Xt30^&)!s4fsgJcYun&j04BM~ z=zMHJv|8g~Vo!m5nnCM@xYMKbtj~C z39q9vbSE87W!X!u?_GpxP~1*-!hG{W*OW7C0PLt2ftRY+9%9OBiU_HV3yL3(9+OzxiQjt>Ir4h02v!tlTu*Q6qzzr#v-13NB$+KGUl-EHU*zk)h2D=k}R|?2Ie#K4D~eBygF0(zr(ttXp-91deVQ z*ElF3IFu;16F3D4zvSQz$x4-F8G}p?bCgwCfP#4xd|b2Hz-_I!bntoVGF~<_vBvsJ z<~QFRKZqe<$CGcG!w^S;zwd~a{Ps%aX>6G4fDUpwL;NNrbm>DpeSZSw&qOKJq*Vml)d;V1Q$0 zWGqVOm=_qJ&XzRLi&l0bNR-nRb#PG&HQk9bCK=KS);rL3)XQ4waF^=U+92by?6C-O zBn?Q%ZJKM%Dv;^0en>~O>W3SFc8goK43@S_%=1Pu3M^gykW|#j#EaO4C3b2=F%o4h z$)7WXz8}rsBrry@`M)N7c?i9=TbTx*>O8O`?Wli#E-3& zs3e*?@KEO>dO3GtJm|=bfJ~6JDd( zTI__o#rz;hia}Fjq zPsfq_B8oLOgxC=rOkYIE^^AY>LrvM?; zeBHPL`xriA#-0GW_qF5jZBY}+dl?n%>@XU=vohAYDrNc9IHmYv+n7(E#^rxxRJhu0 z)g?cb2Nu8*IZSFZB_uqURSG#~?1xJ3Y2DN6L-FyaYc4VONg_U%vh?W*1)JTTFi}B1;*=d`xMc5y)6;qK52)-gmppwv(J0=6YL_w5-{fGTfTw^+TW>xG*A3{!lJ|q!VW!~?NovdlPve1&@g1bSMs4YYNy#t@iQ7QFj*8mPw zg;k4fR%YT+yDu%kOvLhidnM+OQb>!u&iaHyf8-fzo)_KU!jH2aTdwIReftgwdlp%` zY6M7sl?jQnrdsVRby4#u28cfwSFOJ90%`CT)`Zbbi7yR%BCsRRdos{cFvjVF@NL#X zmBj9HBX5RObpy<9Itrd8N`+>c`Moc5EcCC5!v9Gs%y0Jge{m`dt{vw;XcYcwj$fw= zzx3B@bNn+EX8vcYu;g>8jXTUquq6CWD%n|`coeWAu_mxdDd_czPO zTs-9C##I7~K7EJw(khfJIcCZjo*;Q7Le4?f=;Dx~c7AxwnsWs4xgFR`1JodT*krN6 zQ+M%-f3N=V<>Kk8#{?|T?tRKOF0xmerAeSqdy(t{U(2{8vC0^`YV=o7pM^sF-327+ z6_XDWl@;$pUAvdWzzW|+DC=0LmA(CO_tu$lrqR9VTi&I*p)NeKq*28jxhF|A`-uJ& z*p*NJS&^puJ0t759m|F{#&2g&rmu>4@Qbu|B~d-zl*XhvfLxf9pD6I}wK3vy72#bU z1YzRtnzx&&Vt({wi4cunukgy7ioXuam}t*Zc`(sir`MX-<26p2;TP9O6DU&DCk)Om z5A=K*;bC9YQktxXqSl$S6`Ui3f3XGx6iAQ#}vX8_hy+ z0njfEI=Oq&mS*wK`N9$*d}33-wVQzKP0w_(qkdOCW@4eSm%E8agryhXUbEeGW;` z7c`MByI@F{VdR|EKEW_7v<3tBEo#25QZxP4WJf;+RuD(~>5t=QQVO;9COjcq$%N|8 zN1jHmum`u*jb&E|UnoA9A>0ZIOUa!mS$Fm3BD&n~_RC5+chHZ-x-d3stwLEee76;Z zmvHy8l_tjYcHiSKw;x$-YfZUUTKc#TJXT0xvh_ImWC5V&pUY>A1`t(Bb%uCuFJ{;lCPB_B&~J6J(||BQDPgcFe@xYLF6%vy?Krqc z^8HkV-~rX}#dzL6Lb%T^22&}S04Tz$Ab9I7nBw_EuJ|Z2`+UNZt)m*~%1OR6m)D)Y za!mK9o47GP;~qYY`c+t7yJ>uNHl4&_M&Ya)Rw$nLLN`(B{4^GYH_M7cmto#Ta)H(O zGU5o{M_fk~9N-p=tYlOL=!D&wo*QP{3fPnOT_poYZ+QqO_Ds0Zwgr-K$?D;o!58AG zbOmXckS%8fi|aS?$rzw8gfHmd6ze$0n{ZI*Xq{Yr4{0{ zC?rn20jKQI`e2G=f0b}E&vf9^&7Trr9*X=UZF@grMbG{Urzo218emq?^UM%e1^djn z#Z}P}mmR~dCT&Q_9LM=nzKAdKWT5;(w%00l?jcGBf6<%w2F$FYXo|d5136l<#w6`T zwRCXe-vAS?bCRO9;B|swmu$aX@CG-xx<%f>ASvzmlRW>!t%r@SioJ2V5HWhRq8N7XdR0U$yV%|Ht#j&X`yG0`*rBa#|AOC(Fq;&ap4MfvT$9%XB@2QtVee5aeVAz=}7n# zj-w;dmksjrX#@0p>kG)a%M%q(4@VCr*wauVFae856Y!Igtu8leq_zVLHL#@o*$5L3GE1#-qEZRTSc6F4N3FUIm1H%_PwPf zy?r`H<=frd>u>0zLFg)aDA!+Pig}o?dnET=%{w?qF_z|CoS_5q1@+z4svY%O6A^C8 z*R2cY29wYc9*Apdk^Z(PVA+8I7ih-8b9eFHmsUG5jPQ~q?YW)~!yt(0wm+B~{u!k2 zx|SLi)`Bn2w-_=50Tkbll>M1h#$CfidvQ{}Of==Xv3(t%?y1KUbWrFs{nFIJF|+xZ zeJ4>(KglvTVL6{>1k2K#`_0RiUt#*5u3se3LTeLu$lYl5?%ZI2L@oOxLLLNe&S8>v zNMD(KCN$fnX)k}W6srK@Q;p0o;;0hoW+K~QC4%FtN$}m9KswzfRa+L+RP28=2xzWUEYIBRMb=+UD+3Qbc}BsD(5vRLb|=8%*WuAY ziDf4*%Y01q+Mn^Ctgn8{todKOzUm86tNh*F@p^sr(qI36eKqrQeI;%5yuL!z+hSy> zS*PXuAw#20iAz$IzZ9N1rcNyL%|t3az>+L56DH6&LR};cSGSiy_V$Sm&L;|H^V|H} z8F|axg+u3IhvdPsro@X4mqj}cqlP2ZTr_niW82tRh-fz6hYj|;sm{R zG!Q0&IC6mngJnfBZ*r1a&uf|D8e9s^wof$Fd*-$=#05}5eQ`_&2!ZYbjBL_*37}|3 zBr3fX3d^U#(pjUpZJ1^X7t9R!Ahmv9T9~&u5GJl^mL^1L47@cr3Sq$P%dC=UNz8Bv zd^GQFJdpWiz&q|Y4vIJBY|P+y2yAFr>5PwIz8JDQ>X zun>A%48|z}`6P2!bSGRJVum$FFXbXC2Mp@)n3|lzLO&XCMTwLK2M^d)l{_Cuy|M6) z469BByYK$0MP>>G-MlBc)vKZCtZ=G_JB*|089`uRTWxGA!nNrpCe{$}6PGiI2G{-E zsyj7}8>0l<@(>x8A3Bi3tr$&uh>Bqz?4GyCMy#``Jr&ivB0_Wuw&sPS+r&p|`L#-5 z5G=irV_>d9GMz~I-t_ItMnRoi0D=4>TF} zqLKs;@jG9qbH3+#8$*+_7eTPS!a8@$QCa#z!LzqEnDpgFVzW3Yc)rhuzLsQG{4MIZ z`NN)o3aB3EP-^5BW1vG6UFq6ZA2x}X(!zKo5XWVsOQHuPFlk;tDJlBPj2irEcP+?KPz^wsD{*yEo7=W4zaab?&b z6C%8YKvv>%Z_lJl7ZZ;Z%rwD75P~lFl%0oRHs=BoZ934RGcw@(e#47!^xyZ zbG`vAj^xI)A;XEri3K1dgCSZ>f@ph&mJ1V`K!ews-xVe!$-eKdR6qtupvlKJEy% zNnt8~b23`ck0PP|HO{vBG}MwKe~T4XVh&fB#l7d?qB)`drMkutqv>!5uCd#IWysDU zKDn3K`_~%-f~_#Kj%u>~NRN_>?J9)v zFe6h=b(9A53a#4COgD>>las5(7Jjs4D+~=~YJRjcBbcy4vSx32WG21{m>$k`e9*v$ z4-(7{5dr1$%3y61Fr~1oHvCu_=%iXR!kg&KuTIokr@nrEq*^*T*p49iNG60}4dS%S z6G%|RUil_tUDSiowZr0j)MBBB-4JX@0GhBpcWAh@3HOKCALkbI6+P7&+XfrQ{1W#A zCYI+KgeX|h3Gx(`_$$+wF?7Dr!qqD{MYKz%81pPvdafgaw3s##GVNGfX zDe{%OcV_aS_KUY>xII_rQ49>EWRv+C7CeyBQRgEBL?+CGg2KprDSRV^J#y4H{I_Mi z5xi(pK!p?Z8QtIC`?IwnIPD{~;Hj&Dt}gw6(32@C-o4O`0e>gacV_YpRPl;RjY#grSPXcUO7@&j5`Ap|%sM zW-!YHwdX*Xjd`X71oaZa(zSNSiSyUja>fRMZ}Di9>oOF#$ZW1uGcjMvI;=Zt zvtu>rdGwl`47YM9OlrKK?V5t0+ry3{gJf+dAoNwB!N3~1^c zL)?kYKgxTZtum`_YH%4NLf@xEY#y=c1PWK(DrPP!`=7>^bf7}V5$CV2CfCpC2@k*& zX=kCN;pB`ne&8ZBG_nIjnKlVV1}`@&n(w$nP2}4g&RZ}y9%raz zltR)g9DMBWTW=SO_cQwFN-bPadDh8w?mjs`5g!|N#Y>}T!#;M#M!Yq#XflL~bI-_9 zyWLmoPuDwXtPtl8fr8URqiwb_+FnTmRsjl}teR|820QcC^%e+T zs)!+VZQ_8Ep-uNcg1>qS3YcVnGX}VHHe(K1SOfwM{xXxpR^dEZk%)9kng%^KNzSYY zWk%jxtnS9v-M-?tR~^PpcGl@nk%RH0GxEdkL_sr4w3d+yg)b550~NZNH4o`Ul6H(b zw766Q>bT*)ni)INf;g+4ob~?54@n==!_?V; z3Zr5Nx{Q~YygZHM_al3>Q(SmkzMTd!&KrXdVpd_O)Hh5SYw*1 zl2Yb676n*2gjG^!rr10{h1A(fTQu3)a0A7AZp>cR_@@U6q@MtXo&brD^v#j=k9G}) z3$fm_e`ns(UwLN`%c&Tu2I5PRrD+xJfg7FgCxJJQ>84)f7uL)6As&)68Exq}|iY6Z|YRjl%SS7*pe;BGkUP{bPodbrpZ2)S(UKwt{jp0a~g0N)psZpTdl z!5k-M7xRKNhv?fZvJwMmG_HZ(rsqhK}vfgM}>KC*69|Y|EmUs5w6R@`m#p~&%{&98$GUM`uHK4h>5wTMzXR z%8khyPBu1@+xXI6o*8dy#wB!qi21`uOJ@7>TWVJ)YJ=j{@U%SL&M^yRcyZl~d`?_z zNtw2NleE7roDm7%^K=Cn=(?2$Wf~-t)#Ev1_nAVI%9!KfA{zHx{)&vJx*U7}cS`y& zslf;Q(&szX1U+A=L0`w9e-{kBF1;6$RBLNexOc9v%=vKlol(!d>bbVf;iS$f#JgVd z)}%#RZ#~(?_DAE70jD!Y6X!`6(k9~!tm$(=m=N!GQS~xXN%M`}DteI~%u|sSk9EJ& zNXz6BrSNpsCO`Cc_F#*l2TD?z_^7K|c!->1j8C-aYB*n^;%B&gwQe(DXbr5cT$&ks ze==LhMNC^RC~BGQU?+wwJm2eEbr&`g&S3|Z8a|C}HobP+(P-=6f@#rkiTL5OvQHF^ z_i`yll9GkD(;Oe^GjnY#CE z;87Adhr`JY$>q?r!@T2J?;z@p!&xP-X6o^5yA^CdL0yrtwiYrIg#%N=koywKGl`Dy zz&rCNrM;~`kx*AlVHY%M-|u8AppB~h&OaC|o1S!(PPbHe4>2#CKd1g49){DP;Y(sh zSk!VTxQYeq^y(g)E6(hB>Zys0s_rs-J1q`(V)l58cUJ}nKM~t7gxpx(J1jqU{%*Pb zhI_=?BR_=|%TDaGU;>e_Xm^bS#Ie#!{bgP5>Cz|&=KkG}pU0Q+DA=zDQ*;RueeO+< z7(fs}P^os0jtww<(A*{?0|h9|54jh{SJp_Wq4Sj^vf03Gey}xe)N3fu9mV!nxq@f8 zOhAWqJOQFi_(C~S-W1O{*Eo&NbBY(1MwR<+{qm#MABkhP4>g$7;TC6)7ak4 zl+kaWmHgn33F-8v&%y|D13#i?YqR=l^}Q=IAqzgnL(}eC zRCrSVfuqG$ez#9-fu{+zvyl+S{OvR56^a1b25`Rr4J`#xivNI}!Q2M==TMfEs`_(m zek?MnQrS$$gwH_v66(H89bfhRMlN0W^W(8mZ0m+08yXQG-Z+GDWAjqn4>7eCq{u3U zS9IJxsH+}Il}a_y)hV?AVl22n{VoTrm1|v1 z;8J)+;}Ca!`%=-rKP2H?S|lx%G*w*8>D@sK^tCHMQ5AtWf}W3(kCIufzAgb-jkKh;LzoTkD)6c4uZ z3ULWb$ruh~;|Zw91*S>#iKgAU9)T)-^jTO%rX%AzGEY`vurwN?dnnnO;kvs&bKpFB zHZCrUDP;+5KwfynsG#)k%WSRWI5$h?0+yfr19LMIkjktq)VlXKvrnK9AXKKCe z#Wx?kYrQU|2D6`lN<)Ka0=C%}5=;WKTUKRfa{z($Vxl@Xo#aDlm!+?;$-6P1B`xe$ zEP9K?#Z>JOKO0Qdab8L5^60)^vB&#Koy&V~os1u=sJE@Gk-p&ABZql8Ym;Ng2n4aR zGOHpLaacG=G4Ck3GkcqeZsDsaj6l&$T2&URnJFR+RT}Uo$bBu$cqM2WHiZ=CW6b5i ztqwAvMpbdfb!5`<-Xm&)zH%mYHE+6z`~;P=U3WS%TaJlAUcbIhbjqkX?3U{asBwBS zYleyqwlfvA%s%k5kY(0%SV>=-#%9ci8_%_;2TP&$=n1*YgD)9*g;1!qH1zVA3T|ZE zbWZZeA98+lkHGnSpAs0CGlAujJ&1zq6Z|3%-SjEFQlJLbgyk)Raa2tp$^?529EXZm zov4J#SE=iCl1DKo1WY*Aop5cx>i2DdKRPbEzf+_y>6EdkauqtBE*DZbo|%-SZheDC z4QYTgi!=;F&8#q?AR?D3`BQ<|X;i-8AitPp7GyY(f*4<6tL^0rvzXz@ zv+t@I{3v7q)mVUjbqVneSiyw?D_;Ye|WlgD3b$e8HZw z33zkzGmVYi@;jIuaGQ9h8J{K#}LtjAIZr&2pN$f`m_xRtJ#C>Mn+1!>B~bc$B4Q6^;; z<^p)o^I;beBX$#QUYov+$z|DrU_lwg`>>2l7BnFW z5R}boE)9EmoAc=qYh^1u8>{<{ zv#z`cS8=UTOX!MDr72w>el!N@-=J2oOpZ1|cAZ+HC72$G#^TO2F;xh3J%C48cmr7{q<~~~@4xYg}LT6hv z;fR;0=zcmJf*hIu)1LmW3RuWH;t$`Y69f_8)3Apne)&qoBiHC6%erH6#lcnBl=u;9 zf{Y3(9X3!S*$F3whm??bBP?)x-LJX9gh7@5R8yp?E~|+p;irlRE8PrbPSX&0}W9; zi0>QyYL=`4LGSA6wYiBRW6#Ia^_j{g%dQ%Hr>sW@8*Qj*oLF#FTgsyytSwIlI<>G4 z2$%N*hIJPZfT-9aO=1aw@>&B?a!6^Mm&etqR|(y3>;;|4^I?OlW80unz<4H%r_U`* z-ibn{Kw@Iu!4@Wv%|Yq&4v(z(=Up=L-p1|h8zRwcXB$hSY(lhP~L;BQ+N+7E{=YB{b^*Sy0rs zYu*skuJGC6;DZHO&36DTM@Juqq%rVT2_=U$ro01(+V|UP8sIoJ7L?H=^?fN~4G2W> zhT=KM{iU)p4!}i1ViWoEukP;w0AT(e`u_6a4*2r(e?r8mw9RZxt*t0c%^iM~p3UIn zZGQgG7iqW`fxlnwfTqtn|0eyf6u!U8Z{AhxyFSajvso&L^RV&Z9GP5wZhZOq_Fv!n zmH+^No`I#6^(SKkD{Ebye>dpJxAnc>t8I*Y|Q^%>n{aX^Ph!ZydA#vNFZ{zT&(L2ehF{O}1kIoFe=$*B z>My-8@V=B(!g>DudAS3oe13kz{ojaue;vLr+_x_bwl5`#aDVE0{&WW{(E78we+Ihz zs_%sY@rB9ocX*fjKlQyn_j#3N`>VPa62TWLz~AA$G5oW-Kf`;Wpn74FdSR7%-FqL* zf9iW3*YlU~UO1Xwn3i5xmtOZCpT(cm{Tbd%VAczj!^>;@H9R_pf3NSA%k$UXd!a~p zAv$=8^nMLb+wsrp{tWM5Ul_09@q7N%_p3|tU+R(k3h$-D@k^cJ7cb>&czj-eR`(b1 zo=;v|o_p_Q%l;akt^ZGbugBi2VxYgmd#RlHQvdR0zxx_qV8EZ%{Tbe$XYOp+Pkpab z=>GEIzLXw&DIoUy)Y$&|XLWyu_i}Fi?|V-^_NTsId+#qE?n`Bsm-;Hd_nv;-pVj>Z zyniuPuODt&@=txQ$KGGg+?NFSmo)a@$6jj6pVj>t-ivSP#YOb{*c;6Jsqgi}eU-rS z>)3lquzX3A{Jr;;PU!pF5hc{UHXLWyu_vfi`_x#NN)_cEBjeiO7_!ZtuG{j3x!|x9l zy7|xQ{sP_~r-p3jPkpba#$UczUXGex4vc=E8ZupfR`+Lkf1VoA{Xg}+o*IAoa9^C3 jFRsVmdoO$7&+7ga9?0{=0sx?Y{wjH11zwN6{PlkTc2$ tuple[dict[str, np.ndarray], dict[str, Any]]: + arrays = dict(np.load(DATA / "studies.npz")) + record = json.loads((DATA / "studies.json").read_text()) + return arrays, record + + +@cache +def fixture_cubes() -> dict[str, Cube]: + arrays, record = fixture() + return { + name: Cube( + {layer: arrays[f"{name}/{layer}"] for layer in ("counts", "detected")}, + samples=tables["samples"], + groups=tables["groups"], + features=tables["features"], + n_cells=arrays[f"{name}/n_cells"], + ) + for name, tables in record["studies"].items() + } + + +@cache +def fixture_collection() -> CubeCollection: + return CubeCollection( + fixture_cubes(), + align_groups_by="cell_ontology_id", + group_mapping=fixture()[1]["mapping"], + within_study=["donor"], + ) + + +@cache +def fixture_results() -> dict[str, DifferentialExpressionResult]: + collection = fixture_collection() + results = {} + for study in collection.studies: + cube = collection.aligned_cube(study) + design = Design(cube.samples, FORMULAS[study], levels=LEVELS) + results[study] = limma_voom(cube, design, contrasts={"stim": "conditionstim"}) + return results + + +def scaled_difference(ours: np.ndarray, reference: np.ndarray, p_value: bool) -> float: + """|ours - reference| over max(|reference|, 1), or for p-values over + max(p, 1e-300) max(1, |log p|), as validation/cube_collection.py scales them.""" + if p_value: + scale = np.maximum(reference, 1e-300) + scale = scale * np.maximum(1.0, np.abs(np.log(scale))) + else: + scale = np.maximum(np.abs(reference), 1.0) + return float(np.max(np.abs(ours - reference) / scale, initial=0.0)) + + +def limma_table(key: str) -> tuple[list[str], np.ndarray]: + arrays, _ = fixture() + return arrays[f"limma/{key}/genes"].tolist(), arrays[f"limma/{key}/values"] + + +# The references: pandas, anndata and limma. + + +def test_samples_are_those_of_pandas_concat() -> None: + samples = fixture_collection().samples + reference = fixture()[1]["samples"] + assert set(samples.columns) == set(reference) + assert samples.columns[:2] == ("sample_id", "study") + for column, expected in reference.items(): + ours = samples[column].tolist() + missing = [v is None or (isinstance(v, float) and np.isnan(v)) for v in ours] + assert missing == [value is None for value in expected["values"]], column + present = [v for v, gap in zip(ours, missing, strict=True) if not gap] + assert present == [v for v in expected["values"] if v is not None], column + # Integer ages in two studies and floats in one become float64, as in pandas. + assert samples["age"].dtype == np.float64 + assert reference["age"]["dtype"] == "float64" + # Sample and donor ids repeat across studies and are kept apart by the study. + assert len(set(samples["sample_id"].tolist())) == len(samples) + assert samples["donor"][0] == "lupus::D1" and "covid::D1" in samples["donor"] + + +def test_features_and_measured_are_anndata_outer_join() -> None: + collection = fixture_collection() + arrays, _ = fixture() + genes = collection.features["feature_id"].tolist() + assert genes == arrays["first_listed"].tolist() + assert sorted(genes) == sorted(arrays["outer_genes"].tolist()) + where = {gene: k for k, gene in enumerate(arrays["outer_genes"].tolist())} + counts = arrays["outer_counts"][:, [where[gene] for gene in genes]] + for row, study in enumerate(collection.studies): + mine = counts[arrays["outer_study"] == study] + # The outer join fills the genes a study lacks with NaN, stored as -1. + assert_array_equal((mine >= 0).all(axis=0), collection.measured[row]) + assert_array_equal((mine < 0).all(axis=0), ~collection.measured[row]) + cube = collection.aligned_cube(study) + groups = cube.groups["group_id"].tolist() + samples = cube.samples["sample_id"].tolist() + positions = ( + [groups.index(g) for g in arrays["outer_group"][arrays["outer_study"] == study]], + [samples.index(s) for s in arrays["outer_sample_id"][arrays["outer_study"] == study]], + ) # fmt: skip + assert_array_equal( + cube.layer("counts")[positions], mine[:, collection.measured[row]] + ) + + +def test_stack_is_anndata_inner_join() -> None: + collection = fixture_collection() + arrays, _ = fixture() + stacked = collection.stack() + assert stacked.features["feature_id"].tolist() == arrays["inner_genes"].tolist() + assert set(stacked.groups["group_id"]) == set(arrays["inner_group"].tolist()) + assert stacked.samples.equals(collection.samples) + groups = [ + stacked.groups["group_id"].tolist().index(g) for g in arrays["inner_group"] + ] + ids = stacked.samples["sample_id"].tolist() + keys = zip(arrays["inner_study"], arrays["inner_sample_id"], strict=True) + samples = [ids.index(f"{study}::{sample}") for study, sample in keys] + for layer in ("counts", "detected"): + assert_array_equal( + stacked.layer(layer)[groups, samples], arrays[f"inner_{layer}"] + ) + cells = stacked.n_cells + assert cells is not None + assert_array_equal(cells[samples, groups], arrays["inner_n_cells"]) + # An entry without a row in the join has no cells, and zeros in every layer. + empty = np.ones(stacked.layer("counts").shape[:2], dtype=bool) + empty[groups, samples] = False + assert empty.any() and not cells.T[empty].any() + assert not stacked.layer("counts")[empty].any() + + +def test_studies_analysed_alone_and_lined_up_give_limma() -> None: + collection = fixture_collection() + results = fixture_results() + arrays, _ = fixture() + for group in collection.groups["group_id"].tolist(): + table = align_study_results(collection, results, group, "stim") + assert len(table) == collection.n_studies * collection.n_features + for row, study in enumerate(collection.studies): + mine = table.filter(table["study"] == study) + key = f"{study}__{group.replace(':', '_')}" + if f"limma/{key}/genes" not in arrays: # the study lacks the group + assert set(mine["status"].tolist()) <= {"not_analysed", "not_measured"} + assert np.isnan(mine["t"]).all() + continue + genes, values = limma_table(key) + where = {gene: k for k, gene in enumerate(mine["feature_id"].tolist())} + positions = [where[gene] for gene in genes] + expected = np.where(collection.measured[row], "not_tested", "not_measured") + expected = expected.astype(object) + expected[positions] = "tested" + assert_array_equal(mine["status"], expected) + for k, column in enumerate(STATISTICS): + difference = scaled_difference( + mine[column][positions], values[k], "p_value" in column + ) + assert difference <= 1e-10, (study, group, column, difference) + + +def test_the_stack_of_two_studies_fits_their_joint_model_as_limma() -> None: + stacked = fixture_collection().select(studies=["lupus", "covid"]).stack() + design = Design(stacked.samples, "~ donor + condition", levels=LEVELS) + result = limma_voom(stacked, design, contrasts={"stim": "conditionstim"}) + for group in stacked.groups["group_id"].tolist(): + genes, values = limma_table(f"joint__{group.replace(':', '_')}") + table = result.table(group, "stim") + assert table["feature_id"].tolist() == genes + for k, column in enumerate(STATISTICS): + difference = scaled_difference( + table[column], values[k], "p_value" in column + ) + assert difference <= 1e-10, (group, column, difference) + + +def test_lined_up_results_feed_the_meta_analysis() -> None: + results = fixture_results() + meta = differential_expression_meta_analysis(results) + collection = fixture_collection() + group = collection.groups["group_id"][0] + table = align_study_results(collection, results, group, "stim") + lined_up = meta.table(group, "stim") + tested = table["status"] == "tested" + counts = { + gene: int(np.count_nonzero(tested & (table["feature_id"] == gene))) + for gene in lined_up["feature_id"].tolist() + } + assert lined_up["n_studies"].tolist() == list(counts.values()) + + +# A small example worked out by hand. + + +def study_a() -> Cube: + return Cube( + {"counts": np.arange(1, 13, dtype=np.int64).reshape(2, 2, 3)}, + samples={ + "sample_id": ["S1", "S2"], + "donor": ["D1", "D1"], + "condition": ["ctrl", "stim"], + "age": [30, 31], + }, + groups={"group_id": ["B cell", "T cell"], "cell_ontology_id": ["CL:1", "CL:2"]}, + features={"feature_id": ["g1", "g2", "g3"], "symbol": ["G1", "G2", None]}, + n_cells=np.array([[5, 6], [7, 8]]), + provenance={"source": {"study": "A"}}, + ) + + +def study_b() -> Cube: + # Features g2, g4, g3 and g5, which is not measured and filled with zeros. + counts = np.array( + [ + [[1, 0, 2, 0], [3, 1, 4, 0]], + [[5, 1, 6, 0], [7, 2, 8, 0]], + [[9, 3, 10, 0], [11, 4, 12, 0]], + [[13, 5, 14, 0], [15, 6, 16, 0]], + ], + dtype=np.int64, + ) + return Cube( + {"counts": counts}, + samples={ + "sample_id": ["S1", "S3"], + "donor": ["D1", "D2"], + "condition": ["ctrl", "stim"], + "age": [40.5, None], + "batch": ["x", "y"], + }, + groups={ + "group_id": ["B", "CD4", "CD8", "NK"], + "cell_ontology_id": ["CL:1", "CL:2", "CL:2", "CL:3"], + }, + features={ + "feature_id": ["g2", "g4", "g3", "g5"], + "symbol": ["G2b", "G4", "G3", "G5"], + "measured": [True, True, True, False], + }, + n_cells=np.array([[3, 4, 5, 6], [7, 8, 9, 10]]), + ) + + +def test_a_small_collection_worked_out_by_hand() -> None: + a, b = study_a(), study_b() + collection = CubeCollection({"A": a, "B": b}, within_study="donor") + assert collection.studies == ("A", "B") + assert collection.samples.to_dict()["sample_id"].tolist() == [ + "A::S1", "A::S2", "B::S1", "B::S3" + ] # fmt: skip + assert collection.samples["donor"].tolist() == ["A::D1", "A::D1", "B::D1", "B::D2"] + assert collection.samples["study"].tolist() == ["A", "A", "B", "B"] + assert_array_equal(collection.samples["age"], [30.0, 31.0, 40.5, np.nan]) + assert collection.samples["batch"].tolist() == [None, None, "x", "y"] + # Genes in the order the studies first list them; g5 is measured by none. + assert collection.features["feature_id"].tolist() == ["g1", "g2", "g3", "g4"] + assert collection.features["symbol"].tolist() == ["G1", "G2", "G3", "G4"] + assert_array_equal( + collection.measured, [[True, True, True, False], [False, True, True, True]] + ) + # Every cube has cell_ontology_id, which aligns the groups; CD4 and CD8 merge. + assert collection.groups["group_id"].tolist() == ["CL:1", "CL:2", "CL:3"] + assert collection.groups["cell_ontology_id"].tolist() == ["CL:1", "CL:2", "CL:3"] + members = collection.group_members + assert list(zip(*members.to_dict().values(), strict=True)) == [ + ("A", "CL:1", "B cell"), + ("A", "CL:2", "T cell"), + ("B", "CL:1", "B"), + ("B", "CL:2", "CD4"), + ("B", "CL:2", "CD8"), + ("B", "CL:3", "NK"), + ] + aligned = collection.aligned_cube("B") + assert aligned.groups["group_id"].tolist() == ["CL:1", "CL:2", "CL:3"] + assert aligned.features["feature_id"].tolist() == ["g2", "g3", "g4"] + assert aligned.samples["sample_id"].tolist() == ["S1", "S3"] + assert_array_equal( + aligned.layer("counts"), + [ + [[1, 2, 0], [3, 4, 1]], + [[5 + 9, 6 + 10, 1 + 3], [7 + 11, 8 + 12, 2 + 4]], + [[13, 14, 5], [15, 16, 6]], + ], + ) + assert_array_equal(aligned.n_cells, [[3, 4 + 5, 6], [7, 8 + 9, 10]]) + stacked = collection.stack() + assert stacked.groups["group_id"].tolist() == ["CL:1", "CL:2"] + assert stacked.features["feature_id"].tolist() == ["g2", "g3"] + assert stacked.samples.equals(collection.samples) + assert_array_equal( + stacked.layer("counts"), + [ + [[2, 3], [5, 6], [1, 2], [3, 4]], + [[8, 9], [11, 12], [14, 16], [18, 20]], + ], + ) + assert_array_equal(stacked.n_cells, [[5, 6], [7, 8], [3, 9], [7, 17]]) + assert stacked.normalization is None + + +def test_groups_by_id_with_a_mapping_that_renames_merges_and_leaves_out() -> None: + mapping = {"B": {"B": "B cell", "CD4": "T cell", "CD8": "T cell", "NK": None}} + collection = CubeCollection( + {"A": study_a(), "B": study_b()}, + align_groups_by="group_id", + group_mapping=mapping, + ) + assert collection.groups["group_id"].tolist() == ["B cell", "T cell"] + # The studies agree on the ontology ids of both groups. + assert collection.groups["cell_ontology_id"].tolist() == ["CL:1", "CL:2"] + aligned = collection.aligned_cube("B") + assert aligned.groups["group_id"].tolist() == ["B cell", "T cell"] + assert_array_equal(aligned.layer("counts", "T cell"), [[14, 16, 4], [18, 20, 6]]) + # The cube of A already has the aligned ids and the collection's features. + assert collection.aligned_cube("A") is collection.cubes["A"] + + +def test_a_group_keeps_a_column_where_the_studies_that_have_it_agree() -> None: + a = study_a() + lineage = Cube( + {"counts": a.layer("counts")}, + samples=a.samples, + groups=a.groups.with_columns(lineage=["lymphoid", "lymphoid"]), + features=a.features, + n_cells=a.n_cells, + ) + # T cell merges CD4, CD8 (CL:2) and NK (CL:3) of study B. + mapping = {"B": {"B": "B cell", "CD4": "T cell", "CD8": "T cell", "NK": "T cell"}} + collection = CubeCollection( + {"A": lineage, "B": study_b()}, + align_groups_by="group_id", + group_mapping=mapping, + ) + # Study B has no lineage column, so A's values stand; the members of T cell + # disagree on their ontology id, which is then missing. + assert collection.groups["lineage"].tolist() == ["lymphoid", "lymphoid"] + assert collection.groups["cell_ontology_id"].tolist() == ["CL:1", None] + assert_array_equal( + collection.aligned_cube("B").n_cells, [[3, 4 + 5 + 6], [7, 8 + 9 + 10]] + ) + + +def test_groups_without_an_ontology_id_take_one_from_the_mapping() -> None: + a = study_a() + c = Cube( + {"counts": np.ones((2, 1, 2), dtype=np.int64)}, + samples={"sample_id": ["S1"]}, + groups={"group_id": ["Bcells", "junk"], "cell_ontology_id": ["CL:1", None]}, + features={"feature_id": ["g1", "g2"]}, + ) + with pytest.raises(ValueError, match=r"group 'junk' of study 'C' has no cell_onto"): + CubeCollection({"A": a, "C": c}) + collection = CubeCollection({"A": a, "C": c}, group_mapping={"C": {"junk": None}}) + assert collection.group_members["study_group_id"].tolist() == [ + "B cell", "T cell", "Bcells" + ] # fmt: skip + # A group column that only some studies have aligns only with a mapping for + # the others. + plain = Cube( + {"counts": np.ones((1, 1, 1), dtype=np.int64)}, + samples={"sample_id": ["S1"]}, + groups={"group_id": ["B"]}, + features={"feature_id": ["g1"]}, + ) + with pytest.raises(ValueError, match=r"aligned by 'cell_ontology_id', which its"): + CubeCollection({"A": a, "P": plain}, align_groups_by="cell_ontology_id") + by_mapping = CubeCollection( + {"A": a, "P": plain}, + align_groups_by="cell_ontology_id", + group_mapping={"P": {"B": "CL:1"}}, + ) + assert by_mapping.groups["group_id"].tolist() == ["CL:1", "CL:2"] + # Without the column in every cube, groups are aligned by group_id. + assert CubeCollection({"A": a, "P": plain}, group_mapping={"P": {"B": "B cell"}}) + + +# Invariants. + + +def test_the_order_of_studies_changes_no_value() -> None: + a, b = study_a(), study_b() + forward = CubeCollection({"A": a, "B": b}) + backward = CubeCollection({"B": b, "A": a}) + assert backward.features["feature_id"].tolist() == ["g2", "g4", "g3", "g1"] + for study in ("A", "B"): + ours = forward.aligned_cube(study) + theirs = backward.aligned_cube(study) + ids = theirs.features["feature_id"].tolist() + assert ours.select(features=ids).fingerprint() == theirs.fingerprint() + row = forward.studies.index(study) + measured = dict( + zip(forward.features["feature_id"].tolist(), forward.measured[row].tolist()) + ) + genes = backward.features["feature_id"].tolist() + other = backward.studies.index(study) + assert [measured[g] for g in genes] == backward.measured[other].tolist() + stacked = backward.stack() + reordered = forward.stack().select(samples=stacked.samples["sample_id"]) + assert_array_equal(reordered.layer("counts"), stacked.layer("counts")) + assert_array_equal(reordered.n_cells, stacked.n_cells) + assert reordered.samples.equals(stacked.samples) + assert reordered.groups.equals(stacked.groups) + # Only the symbol of g2 differs: the first study that gives one names it. + assert reordered.features["symbol"].tolist() == ["G2", "G3"] + assert stacked.features["symbol"].tolist() == ["G2b", "G3"] + + +def test_the_order_of_genes_within_a_cube_changes_no_value() -> None: + a, b = study_a(), study_b() + order = [3, 2, 0, 1] + shuffled = b.select(features=order) + first = CubeCollection({"A": a, "B": b}) + second = CubeCollection({"A": a, "B": shuffled}) + # The shared genes come in A's order, so the stacked cube is the same. + assert first.stack().fingerprint() == second.stack().fingerprint() + assert ( + first.aligned_cube("B").fingerprint() == second.aligned_cube("B").fingerprint() + ) + assert_array_equal(first.measured, second.measured) + + +def test_selecting_then_stacking_is_stacking_then_selecting() -> None: + collection = CubeCollection({"A": study_a(), "B": study_b()}, within_study="donor") + chosen = {"groups": ["CL:2"], "features": ["g3", "g2"]} + selected = collection.select(**chosen).stack() + assert selected.fingerprint() == collection.stack().select(**chosen).fingerprint() + assert collection.stack(**chosen).fingerprint() == selected.fingerprint() + # The order of a selection is that of every study's aligned cube. + reordered = collection.select(groups=["CL:2", "CL:1"]).aligned_cube("B") + assert reordered.groups["group_id"].tolist() == ["CL:2", "CL:1"] + assert_array_equal( + reordered.layer("counts"), + collection.aligned_cube("B").select(groups=["CL:2", "CL:1"]).layer("counts"), + ) + # A study selected alone stacks to its own aligned cube with keyed samples. + alone = collection.select(studies="B") + assert alone.studies == ("B",) + assert alone.features["feature_id"].tolist() == ["g2", "g3", "g4"] + assert alone.groups["group_id"].tolist() == ["CL:1", "CL:2", "CL:3"] + assert "age" in alone.samples and "batch" in alone.samples + stacked = alone.stack() + assert_array_equal( + stacked.layer("counts"), collection.aligned_cube("B").layer("counts") + ) + assert stacked.samples["sample_id"].tolist() == ["B::S1", "B::S3"] + + +def test_a_single_study_is_its_own_cube() -> None: + a = study_a() + collection = CubeCollection({"A": a}, align_groups_by="group_id") + assert collection.aligned_cube("A") is a + assert_array_equal(collection.measured, [[True, True, True]]) + stacked = collection.stack() + assert_array_equal(stacked.layer("counts"), a.layer("counts")) + assert_array_equal(stacked.n_cells, a.n_cells) + assert stacked.groups.equals(a.groups) + assert stacked.features.equals(a.features) + assert stacked.samples["sample_id"].tolist() == ["A::S1", "A::S2"] + + +def test_cube_files_stay_lazy_and_align_as_cubes_in_memory(tmp_path: Path) -> None: + a, b = study_a(), study_b() + write_cube(a, tmp_path / "a") + write_cube(b, tmp_path / "b") + lazy = CubeCollection( + {"A": read_cube(tmp_path / "a"), "B": read_cube(tmp_path / "b")} + ) + eager = CubeCollection({"A": a, "B": b}) + aligned = lazy.aligned_cube("B") + assert "lazy" in repr(aligned) + assert aligned.fingerprint() == eager.aligned_cube("B").fingerprint() + assert lazy.stack().fingerprint() == eager.stack().fingerprint() + # A cube file's fingerprint is the one stored with it: no layer is read. + assert lazy.provenance["studies"]["B"]["fingerprint"] == b.fingerprint() + + +def test_stacked_layers_take_a_type_that_holds_every_study() -> None: + def cube(counts: np.ndarray) -> Cube: + return Cube( + {"counts": counts}, + samples={"sample_id": ["S1"]}, + groups={"group_id": ["B"]}, + features={"feature_id": ["g1", "g2"]}, + ) + + small = cube(np.array([[[1, 2]]], dtype=np.uint32)) + wide = cube(np.array([[[3, 4]]], dtype=np.uint64)) + real = cube(np.array([[[0.5, 2.0]]])) + assert CubeCollection({"A": small, "B": small}).stack().layer("counts").dtype == ( + np.uint32 + ) + assert CubeCollection({"A": small, "B": wide}).stack().layer("counts").dtype == ( + np.int64 + ) + mixed = CubeCollection({"A": small, "C": real}).stack().layer("counts") + assert mixed.dtype == np.float64 + assert_array_equal(mixed, [[[1.0, 2.0], [0.5, 2.0]]]) + with pytest.raises(ValueError, match="holds whole numbers exactly only below"): + CubeCollection( + {"A": cube(np.array([[[2**53, 1]]], np.uint64)), "C": real} + ).stack() + + +def test_normalized_layers_stack_only_from_studies_normalized_alike() -> None: + def cube(target: float) -> Cube: + return Cube( + { + "counts": np.ones((1, 1, 2), dtype=np.int64), + "log_sum": np.full((1, 1, 2), 0.5), + }, + samples={"sample_id": ["S1"]}, + groups={"group_id": ["B"]}, + features={"feature_id": ["g1", "g2"]}, + n_cells=np.array([[2]]), + normalization={"method": "library_size", "target_sum": target}, + ) + + same = CubeCollection({"A": cube(1e4), "B": cube(1e4)}).stack() + assert same.layer_names == ("counts", "log_sum") + assert same.normalization == {"method": "library_size", "target_sum": 1e4} + differing = CubeCollection({"A": cube(1e4), "B": cube(1e6)}) + with pytest.raises(ValueError, match=r"would mix two normalizations"): + differing.stack() + counts = differing.stack(layers="counts") + assert counts.layer_names == ("counts",) and counts.normalization is None + + +# Provenance, pickling and immutability. + + +def test_provenance_records_each_cube_and_the_alignment() -> None: + a, b = study_a(), study_b() + collection = CubeCollection({"A": a, "B": b}, within_study=["donor"]) + record = collection.provenance + assert json.loads(json.dumps(record)) == record + assert record["studies"] == { + "A": {"fingerprint": a.fingerprint(), "provenance": {"source": {"study": "A"}}}, + "B": {"fingerprint": b.fingerprint(), "provenance": {}}, + } + alignment = record["alignment"] + assert alignment["align_groups_by"] == "cell_ontology_id" + assert alignment["group_mapping"] is None + assert alignment["within_study"] == ["donor"] + assert alignment["groups"]["B"] == { + "CL:1": ["B"], "CL:2": ["CD4", "CD8"], "CL:3": ["NK"] + } # fmt: skip + assert record["operations"] == [] + selected = collection.select(studies="B", features=["g4", "g2"]) + assert selected.provenance["operations"] == [ + { + "operation": "select", + "studies": ["B"], + "groups": None, + "features": ["g4", "g2"], + } + ] + assert selected.provenance["alignment"] == alignment + # The aligned cube records how it was derived from the cube as given. + operations = collection.aligned_cube("B").provenance["operations"] + assert [op["operation"] for op in operations] == ["select", "regroup"] + assert operations[0]["input"] == b.fingerprint() + assert operations[0]["features"] == ["g2", "g3", "g4"] + assert operations[1]["membership"] == {"CL:2": ["CD4", "CD8"]} | { + "CL:1": ["B"], "CL:3": ["NK"] + } # fmt: skip + + +def test_a_stacked_cube_records_the_studies_it_stacks() -> None: + a, b = study_a(), study_b() + stacked = CubeCollection({"A": a, "B": b}, within_study="donor").stack() + construction = stacked.provenance["construction"] + assert construction["function"] == "CubeCollection.stack" + assert construction["studies"]["B"] == { + "fingerprint": b.fingerprint(), + "provenance": {}, + "groups": {"CL:1": ["B"], "CL:2": ["CD4", "CD8"]}, + } + assert construction["studies"]["A"]["provenance"] == a.provenance + assert construction["align_groups_by"] == "cell_ontology_id" + assert construction["within_study"] == ["donor"] + assert construction["layers"] == ["counts"] + history = stacked.history() + assert "from 2 studies" in history + assert f"study 'B': the cube {b.fingerprint()}" in history + assert " CL:2: CD4 + CD8" in history + assert "sample columns written study::value: donor" in history + + +def test_pickling_copying_and_immutability() -> None: + collection = CubeCollection({"A": study_a(), "B": study_b()}).select( + features=["g2", "g3"] + ) + restored = pickle.loads(pickle.dumps(collection)) + assert restored.studies == collection.studies + for name in ("samples", "groups", "features", "group_members"): + assert getattr(restored, name).equals(getattr(collection, name)) + assert_array_equal(restored.measured, collection.measured) + assert restored.provenance == collection.provenance + assert repr(restored) == repr(collection) + assert restored.stack().fingerprint() == collection.stack().fingerprint() + for measured in (collection.measured, restored.measured): + with pytest.raises(ValueError, match="read-only"): + measured[0, 0] = False + with pytest.raises(AttributeError, match="immutable"): + collection.studies = () # type: ignore[misc] + with pytest.raises(AttributeError, match="immutable"): + del collection._studies + assert copy.copy(collection) is collection + assert copy.deepcopy(collection) is collection + with pytest.raises(TypeError): + collection.cubes["A"] = study_a() # type: ignore[index] + + +def test_repr_reports_the_studies_and_what_they_share() -> None: + lines = repr(CubeCollection({"A": study_a(), "B": study_b()})).splitlines() + assert lines == [ + "CubeCollection of 2 studies: 4 samples, 3 groups, 4 features", + " groups aligned by 'cell_ontology_id'", + " A: 2 samples, 2 groups, 3 features measured", + " B: 2 samples, 3 groups (1 merged), 3 features measured", + " in every study: 2 groups, 2 features", + ] + + +# What goes wrong in practice. + + +@pytest.mark.parametrize( + ("cubes", "error", "message"), + [ + ([study_a()], TypeError, "cubes must be a mapping of study names"), + ({}, ValueError, "cubes is empty"), + ({"a::b": study_a()}, ValueError, r"'a::b' holds '::', which separates"), + ({"": study_a()}, ValueError, "study names must be non-empty"), + ({1: study_a()}, TypeError, "study names must be strings"), + ({"A": "cube"}, TypeError, r"cubes\['A'\] must be a Cube, got str"), + ], +) +def test_invalid_studies_are_refused(cubes: Any, error: type, message: str) -> None: + with pytest.raises(error, match=message): + CubeCollection(cubes) + + +def test_studies_that_share_no_gene_or_group_are_refused() -> None: + a = study_a() + versioned = Cube( + {"counts": np.ones((2, 1, 2), dtype=np.int64)}, + samples={"sample_id": ["S1"]}, + groups=a.groups, + features={"feature_id": ["g1.7", "g2.7"]}, + ) + with pytest.raises( + ValueError, match=r"study 'A' shares no feature .* differ only in a version" + ): + CubeCollection({"A": a, "V": versioned}) + other = Cube( + {"counts": np.ones((1, 1, 2), dtype=np.int64)}, + samples={"sample_id": ["S1"]}, + groups={"group_id": ["neuron"], "cell_ontology_id": ["CL:0000540"]}, + features={"feature_id": ["g1", "g2"]}, + ) + with pytest.raises( + ValueError, match=r"study 'A' shares no group .* theirs 'CL:0000540'" + ): + CubeCollection({"A": a, "N": other}) + unmeasured = Cube( + {"counts": np.zeros((2, 1, 1), dtype=np.int64)}, + samples={"sample_id": ["S1"]}, + groups=a.groups, + features={"feature_id": ["g1"], "measured": [False]}, + ) + with pytest.raises(ValueError, match="study 'U' measured no feature"): + CubeCollection({"A": a, "U": unmeasured}) + + +@pytest.mark.parametrize( + ("mapping", "error", "message"), + [ + ({"Z": {}}, KeyError, r"names the study 'Z', which is not in cubes"), + ({"B": {"CD5": "T"}}, KeyError, r"group 'CD5', which the cube of study 'B'.*did you mean 'CD8' or 'CD4'"), + ({"B": {"CD4": 4}}, TypeError, r"group_mapping\['B'\]\['CD4'\] must be an aligned id"), + ({"B": {"CD4": ""}}, ValueError, "is empty; aligned ids are non-empty"), + ({"B": ["CD4"]}, TypeError, r"group_mapping\['B'\] must map"), + (["B"], TypeError, "group_mapping must map study names"), + ({"B": {name: None for name in ("B", "CD4", "CD8", "NK")}}, ValueError, "study 'B' has no group"), + ], +) # fmt: skip +def test_invalid_group_mappings_are_refused( + mapping: Any, error: type, message: str +) -> None: + with pytest.raises(error, match=message): + CubeCollection({"A": study_a(), "B": study_b()}, group_mapping=mapping) + + +def test_sample_columns_must_hold_one_kind_of_values() -> None: + a = study_a() + b = study_b() + text_age = Cube( + {"counts": b.layer("counts")}, + samples=b.samples.with_columns(age=["forty", None]), + groups=b.groups, + features=b.features, + n_cells=b.n_cells, + ) + with pytest.raises( + TypeError, + match=r"sample column 'age' holds integers \(int64\) in study 'A' and text", + ): + CubeCollection({"A": a, "B": text_age}) + named = Cube( + {"counts": b.layer("counts")}, + samples=b.samples.with_columns(study=["other", "other"]), + groups=b.groups, + features=b.features, + n_cells=b.n_cells, + ) + with pytest.raises(ValueError, match=r"sample column 'study' with the value 'oth"): + CubeCollection({"A": a, "B": named}) + + +@pytest.mark.parametrize( + ("within", "error", "message"), + [ + (["donr"], KeyError, r"'donr', which is not a sample column.*'donor'"), + (["sample_id"], ValueError, "which the collection writes itself"), + (["age"], TypeError, r"'age', whose values in study 'B' are float64"), + (["donor", "donor"], ValueError, "names 'donor' more than once"), + ([1], TypeError, "within_study must name sample columns"), + (3, TypeError, "within_study must be a sample column name or a list"), + ], +) +def test_invalid_within_study_columns_are_refused( + within: Any, error: type, message: str +) -> None: + with pytest.raises(error, match=message): + CubeCollection({"A": study_a(), "B": study_b()}, within_study=within) + + +def test_invalid_alignment_columns_are_refused() -> None: + cubes = {"A": study_a(), "B": study_b()} + with pytest.raises(KeyError, match=r"'cell_type', which is not a group column"): + CubeCollection(cubes, align_groups_by="cell_type") + with pytest.raises(TypeError, match="align_groups_by must be a group column name"): + CubeCollection(cubes, align_groups_by=1) # type: ignore[arg-type] + + +def test_stacking_refuses_what_a_study_lacks() -> None: + collection = CubeCollection({"A": study_a(), "B": study_b()}) + with pytest.raises(ValueError, match=r"feature 'g4' is not in study 'A'"): + collection.stack(features=["g2", "g4"]) + with pytest.raises(ValueError, match=r"group 'CL:3' is not in study 'A'"): + collection.stack(groups=["CL:3"]) + with pytest.raises(ValueError, match="no group is in every study"): + collection.select(groups=["CL:3"]).stack() + with pytest.raises(ValueError, match="layers must include 'counts'"): + collection.stack(layers=[]) + with pytest.raises(ValueError, match="unknown layer 'means'"): + collection.stack(layers=["counts", "means"]) + with pytest.raises(ValueError, match=r"study 'A' has no 'detected' layer"): + CubeCollection({"A": study_a(), "B": _with_detected(study_b())}).stack( + layers=["counts", "detected"] + ) + bulk = Cube( + {"counts": np.ones((1, 2, 3), dtype=np.int64)}, + samples={"sample_id": ["S1", "S2"]}, + groups={"group_id": ["B cell"], "cell_ontology_id": ["CL:1"]}, + features={"feature_id": ["g1", "g2", "g3"]}, + ) + with pytest.raises(ValueError, match=r"of 'A' record the cells .* those of 'bulk'"): + CubeCollection({"A": study_a(), "bulk": bulk}).stack() + + +def _with_detected(cube: Cube) -> Cube: + counts = cube.layer("counts") + return Cube( + {"counts": counts, "detected": np.minimum(counts, 1)}, + samples=cube.samples, + groups=cube.groups, + features=cube.features, + n_cells=cube.n_cells, + ) + + +def test_selections_are_checked() -> None: + collection = CubeCollection({"A": study_a(), "B": study_b()}) + assert collection.select() is collection + with pytest.raises(KeyError, match="no study named 'C'"): + collection.select(studies=["C"]) + with pytest.raises(ValueError, match="the selection of studies is empty"): + collection.select(studies=[]) + with pytest.raises( + ValueError, match=r"no study selected measured the feature 'g1'" + ): + collection.select(studies="B", features=["g1"]) + with pytest.raises(ValueError, match=r"no study selected has the group 'CL:3'"): + collection.select(studies="A", groups="CL:3") + with pytest.raises(ValueError, match="repeats 'g2'"): + collection.select(features=["g2", "g2"]) + assert collection.select(studies=[False, True]).studies == ("B",) + with pytest.raises(KeyError, match=r"no study named 'b'; studies are: 'A' and 'B'"): + collection.aligned_cube("b") + with pytest.raises(TypeError, match=r"study must be a study name \(str\)"): + collection.aligned_cube(0) # type: ignore[arg-type] + + +def test_lining_up_results_refuses_results_of_other_cubes() -> None: + collection = fixture_collection() + results = fixture_results() + group = collection.groups["group_id"][0] + with pytest.raises(TypeError, match="collection must be a CubeCollection"): + align_study_results(fixture_cubes(), results, group, "stim") # type: ignore[arg-type] + with pytest.raises(ValueError, match=r"missing 'atlas'"): + align_study_results( + collection, {s: results[s] for s in ("lupus", "covid")}, group, "stim" + ) + with pytest.raises(KeyError, match=r"no group 'T cell' in the collection"): + align_study_results(collection, results, "T cell", "stim") + with pytest.raises(KeyError, match="no study tested the contrast 'ctrl'"): + align_study_results(collection, results, group, "ctrl") + # The covid study's own cube has its own group ids, which the collection merges. + own = fixture_cubes()["covid"] + design = Design(own.samples, FORMULAS["covid"], levels=LEVELS) + unaligned = results | { + "covid": limma_voom(own, design, contrasts={"stim": "conditionstim"}) + } + with pytest.raises(ValueError, match=r"analyses the groups 'B', 'T4 naive'"): + align_study_results(collection, unaligned, group, "stim") + # Its unmeasured genes are filled with zeros, which a test without the + # expression filter tests. + alone = CubeCollection({"covid": own}, align_groups_by="group_id") + tested_all = limma_voom( + own, design, contrasts={"stim": "conditionstim"}, expression_filter=False + ) + with pytest.raises(ValueError, match="which the study did not measure"): + align_study_results(alone, {"covid": tested_all}, "B", "stim") diff --git a/tests/test_masked_inputs.py b/tests/test_masked_inputs.py index 8bab538..55a2067 100644 --- a/tests/test_masked_inputs.py +++ b/tests/test_masked_inputs.py @@ -46,7 +46,7 @@ import cytepool from cytepool._masked_arrays import reject_masked from cytepool.composition import propeller -from cytepool.cube import Cube, from_bulk, from_cells +from cytepool.cube import Cube, CubeCollection, from_bulk, from_cells from cytepool.design import Design, make_contrasts from cytepool.differential_expression import limma_voom from cytepool.gene_sets import GeneSets @@ -226,6 +226,30 @@ def _cube_arguments(context: Context) -> dict[str, Any]: } +def _collection_arguments(context: Context) -> dict[str, Any]: + """Two studies' cubes, with a donor column, and a value for every option.""" + cube = Cube( + **_cube_arguments(context) + | {"samples": {"sample_id": ["S1", "S2"], "donor": ["D1", "D2"]}} + ) + return { + "cubes": {"A": cube, "B": cube}, + "align_groups_by": "group_id", + "group_mapping": {"B": {"T": "T"}}, + "within_study": ["donor"], + } + + +@functools.cache +def _collection() -> CubeCollection: + cube = Cube( + {"counts": np.arange(12, dtype=np.uint32).reshape(2, 2, 3)}, + **CUBE_AXES, + n_cells=np.ones((2, 2), dtype=np.int64), + ) + return CubeCollection({"A": cube, "B": cube}) + + def _from_cells_arguments(context: Context) -> dict[str, Any]: return { "counts": CELL_COUNTS, @@ -551,6 +575,25 @@ def _forest_arguments(context: Context) -> dict[str, Any]: instance=lambda c: c.cube, keys={"mapping": "all"}, ), + # The cubes were checked when they were built; the mapping of each study's + # groups and the column names are read value by value. + "cytepool.cube.collection.CubeCollection": Entry( + _collection_arguments, + {"group_mapping": None, "within_study": None}, + options=("align_groups_by",), + objects=("cubes",), + keys={"group_mapping": "B"}, + ), + "cytepool.cube.collection.CubeCollection.select": Entry( + lambda c: {"studies": ["B"], "groups": [True, False], "features": [2, 0]}, + {"studies": None, "groups": None, "features": None}, + instance=lambda c: _collection(), + ), + "cytepool.cube.collection.CubeCollection.stack": Entry( + lambda c: {"groups": ["T"], "features": [1, 0], "layers": ["counts"]}, + {"groups": None, "features": None, "layers": None}, + instance=lambda c: _collection(), + ), "cytepool.cube.pool.from_cells": Entry( _from_cells_arguments, { @@ -2032,6 +2075,10 @@ def _forest_arguments(context: Context) -> dict[str, Any]: "cytepool.table.Table.equals": (("other",), "a table, compared as it is"), "cytepool.table.Table.rename": (("mapping",), "column names, which are str"), "cytepool.cube.cube.Cube.collapse_samples": (("by",), "a column name, a str"), + "cytepool.cube.collection.CubeCollection.aligned_cube": ( + ("study",), + "a study name, looked up", + ), "cytepool.design.contrasts.Contrasts.__getitem__": (("name",), "a contrast name"), "cytepool.design.contrasts.Contrasts.__contains__": (("name",), "a contrast name"), "cytepool.design.contrasts.Contrasts.get": ( @@ -2113,6 +2160,11 @@ def _forest_arguments(context: Context) -> dict[str, Any]: ("group", "contrast"), "a group id and a contrast name, looked up", ), + "cytepool.meta_analysis.align_study_results": ( + ("collection", "results", "group", "contrast"), + "a CubeCollection and results, which were checked when they were built, and a " + "group id and a contrast name, looked up", + ), "cytepool.gene_sets.GeneSets.__getitem__": (("name",), "a set name, looked up"), "cytepool.gene_sets.GeneSets.__contains__": (("name",), "a set name, looked up"), "cytepool.gene_sets.read_gmt": ( diff --git a/validation/cube_collection.py b/validation/cube_collection.py new file mode 100644 index 0000000..e1f2daa --- /dev/null +++ b/validation/cube_collection.py @@ -0,0 +1,686 @@ +"""Check CubeCollection against anndata's concat, pandas and limma run study by study. + +The script simulates the cells of three studies of one stimulation (seed ``SEED``) +and pools each with ``from_cells``: + +- ``lupus``: 8 samples (donors D1 to D4, each control and stimulated), 4 cell types + with Cell Ontology ids, 700 genes; +- ``covid``: 10 samples (donors D1 to D5, paired: the donor and sample ids of + ``lupus`` again), 5 cell types, two of which share a Cell Ontology id and are + merged, 750 genes, of which 40 are filled with zeros and marked not measured; +- ``atlas``: 6 samples of 6 donors with integer ids, unpaired, 5 cell types without + ontology ids, aligned by a mapping that leaves one out, 950 genes. + +The collection aligns the groups by ``cell_ontology_id``, with a group mapping for +``atlas``, and writes donor ids ``study::donor`` (``within_study=["donor"]``). Some +sample and cell type entries have fewer than 10 cells, and one none. + +References, run as separate programs on files in a temporary directory: + +- anndata and pandas: each study's cube is exported with ``to_anndata``; the + reference drops the genes that ``var["measured"]`` marks as not measured, gives + each row its aligned group (the ontology id, or the mapping), sums the rows of each + sample and aligned group, and joins the studies with ``anndata.concat`` + (``label="study"``). The outer join must give the collection's genes (as a set, and + in the order in which the studies first list them, from ``pandas.unique``), + ``measured`` (where its dense layers are not NaN) and the values of every aligned + cube; the inner join, kept to the cell types of every study, must give ``stack()``: + genes, groups, samples, every layer and the cell numbers. ``pandas.concat`` of the + sample tables, with sample and donor ids written ``study::id``, must give + ``samples``. +- limma 3.62.1 and edgeR 4.4.0 in R: limma-voom on each study's own counts, summing + the groups that share an aligned id, and on the stacked counts of ``lupus`` and + ``covid`` with ``~ donor + condition``, the steps of + ``validation/differential_expression.py`` (``min_cells=10``, ``filterByExpr``, + TMM, ``voom`` with ``adaptive.span = TRUE``, ``lmFit``, ``eBayes``, ``topTable``). + ``limma_voom`` on every aligned cube, lined up by ``align_study_results``, must give + limma's tables, ``"not_measured"`` for every gene the study did not measure and + ``"not_tested"`` for the others limma left out; ``limma_voom`` on the stack must + give the joint model. A control fits the joint model with the donor ids as the + studies give them, so that donor D1 of both studies is taken for one person: its + results must differ. + +Tolerances: identical ids, orders, masks and integers; a relative 1e-10 for float +layers and for limma's statistics, which differences are scaled as in +``validation/differential_expression.py`` (by max(|reference|, 1), and p-values by +max(p, 1e-300) max(1, |log p|)). + +Requirements: + This script: numpy, pandas, anndata and h5py, with cytepool from the source tree. + Reference program: anndata and pandas, in REFERENCE_PYTHON (default: this + interpreter). R 4.4 with limma, edgeR and jsonlite: set RSCRIPT to its Rscript + when it is not on PATH. R runs with its address space limited to 8 GB. + +Usage, from the repository root: + PYTHONPATH=. python validation/cube_collection.py + PYTHONPATH=. python validation/cube_collection.py --write-fixtures +""" + +import argparse +import json +import os +import resource +import subprocess +import sys +import tempfile +import textwrap +from pathlib import Path + +import numpy as np + +from cytepool.cube import CubeCollection, from_cells +from cytepool.design import Design +from cytepool.differential_expression import limma_voom +from cytepool.io.anndata import to_anndata +from cytepool.meta_analysis import align_study_results + +SEED = 20261026 +EXACT = 1e-10 +MEMORY = 8_000_000 * 1024 # bytes of address space for R +FIXTURE_DIR = ( + Path(__file__).resolve().parent.parent / "tests" / "data" / "cube_collection" +) +LAYERS = ("counts", "detected", "log_sum", "log_sum_sq", "linear_sum") +TABLE = ("logFC", "AveExpr", "t", "P.Value", "adj.P.Val", "B", "SE") +COLUMNS = ( + "log2_fold_change", + "average_log2_cpm", + "t", + "p_value", + "adjusted_p_value", + "b_statistic", + "standard_error", +) +ONTOLOGY = {"B": "CL:0000236", "CD4": "CL:0000624", "CD8": "CL:0000625"} +MONOCYTE, NK = "CL:0000860", "CL:0000623" +MAPPING = { + "atlas": { + "Bcells": ONTOLOGY["B"], + "CD4": ONTOLOGY["CD4"], + "CD8": ONTOLOGY["CD8"], + "Mono": MONOCYTE, + "Doublets": None, + } +} +FORMULAS = {"lupus": "~ donor + condition", "covid": "~ donor + condition"} +FORMULAS["atlas"] = "~ condition" + +REFERENCE_PROGRAM = textwrap.dedent( + ''' + """anndata's concat of the aligned studies, and pandas' concat of their samples.""" + + import json + import sys + import warnings + from importlib.metadata import version + + import anndata + import numpy as np + import pandas as pd + + warnings.filterwarnings("ignore") + root = sys.argv[1] + spec = json.load(open(f"{root}/spec.json")) + sample_tables = json.load(open(f"{root}/samples.json")) # NaN marks missing values + layers = ("counts", "detected", "log_sum", "log_sum_sq", "linear_sum") + aligned = {} + for name in spec["studies"]: + adata = anndata.read_h5ad(f"{root}/{name}.h5ad") + if "measured" in adata.var: + adata = adata[:, adata.var["measured"].to_numpy(bool)] + mapping = spec["mapping"].get(name, {}) + groups = adata.obs["group_id"].astype(str).tolist() + ontology = ( + adata.obs["cell_ontology_id"].astype(object).tolist() + if "cell_ontology_id" in adata.obs + else [None] * len(groups) + ) + keys = [mapping[g] if g in mapping else o for g, o in zip(groups, ontology)] + rows = np.array([key is not None for key in keys]) + pairs = pd.MultiIndex.from_arrays( + [adata.obs["sample_id"].astype(str)[rows], np.array(keys, object)[rows]] + ) + codes, unique = pd.factorize(pairs) + summed = {} + for layer in layers: + values = np.asarray(adata.layers[layer], dtype=np.float64)[rows] + summed[layer] = np.zeros((len(unique), adata.n_vars)) + np.add.at(summed[layer], codes, values) + cells = np.bincount(codes, weights=adata.obs["n_cells"].to_numpy()[rows]) + obs = pd.DataFrame( + { + "sample_id": unique.get_level_values(0).astype(str), + "group": unique.get_level_values(1).astype(str), + "n_cells": cells, + }, + index=[f"{s}|{g}" for s, g in unique], + ) + aligned[name] = anndata.AnnData( + obs=obs, var=pd.DataFrame(index=adata.var_names), layers=summed + ) + out = {} + for join in ("outer", "inner"): + joined = anndata.concat( + aligned, join=join, label="study", index_unique="::" + ) + if join == "inner": + every = set.intersection(*(set(a.obs["group"]) for a in aligned.values())) + joined = joined[joined.obs["group"].isin(every).to_numpy()] + out[f"{join}_genes"] = joined.var_names.to_numpy().astype(str) + for key in ("study", "sample_id", "group"): + out[f"{join}_{key}"] = np.array(joined.obs[key].astype(str).tolist()) + out[f"{join}_n_cells"] = joined.obs["n_cells"].to_numpy() + for layer in layers: + out[f"{join}_{layer}"] = np.asarray(joined.layers[layer], np.float64) + out["first_listed"] = pd.unique( + np.concatenate([a.var_names.to_numpy() for a in aligned.values()]) + ).astype(str) + np.savez(f"{root}/reference.npz", **out) + frames = [] + for name in spec["studies"]: + columns = sample_tables[name] + frame = pd.DataFrame( + {c: np.array(v["values"], dtype=v["dtype"]) for c, v in columns.items()} + ) + for column in ["sample_id", *spec["within_study"]]: + if column in frame: + frame[column] = name + "::" + frame[column].astype(str) + frames.append(frame.assign(study=name)) + samples = pd.concat(frames, join="outer", ignore_index=True) + table = { + column: { + "dtype": str(samples[column].dtype), + "values": [ + None if pd.isna(value) else value + for value in samples[column].astype(object).tolist() + ], + } + for column in samples.columns + } + versions = {p: version(p) for p in ("anndata", "pandas", "numpy")} + versions["python"] = sys.version.split()[0] + json.dump( + {"samples": table, "versions": versions}, + open(f"{root}/reference.json", "w"), + default=float, + ) + ''' +) + +R_SCRIPT = r""" +suppressPackageStartupMessages({library(limma); library(edgeR); library(jsonlite)}) +root <- commandArgs(trailingOnly = TRUE)[[1]] +spec <- fromJSON(file.path(root, "spec.json"), simplifyVector = FALSE) +rd <- function(path) readBin(path, "double", n = file.size(path) / 8, size = 8) +wr <- function(x, path) writeBin(as.double(x), path, size = 8) +study <- function(name) { + dir <- file.path(root, "r", name) + genes <- readLines(file.path(dir, "genes.txt")) + groups <- readLines(file.path(dir, "groups.txt")) + samples <- read.delim(file.path(dir, "samples.tsv"), colClasses = "character") + # The aligned id of a group: the study's mapping, else its Cell Ontology id. + ontology <- readLines(file.path(dir, "ontology.txt")) + mapping <- spec$mapping[[name]] + aligned <- vapply(seq_along(groups), function(k) { + if (!is.null(mapping) && groups[k] %in% names(mapping)) { + m <- mapping[[groups[k]]]; if (is.null(m)) NA_character_ else m + } else ontology[k] + }, "") + list(name = name, genes = genes, samples = samples, aligned = aligned, + measured = readLines(file.path(dir, "measured.txt")) == "TRUE", + counts = array(rd(file.path(dir, "counts.bin")), + dim = c(length(genes), nrow(samples), length(groups))), + cells = matrix(rd(file.path(dir, "n_cells.bin")), nrow = nrow(samples))) +} +for (a in spec$analyses) { + data <- lapply(unlist(a$studies), study) + # One study: all its genes and aligned groups; several: those every study has. + genes <- if (length(data) == 1) data[[1]]$genes else + Reduce(intersect, lapply(data, function(d) d$genes[d$measured])) + groups <- Reduce(intersect, lapply(data, function(d) unique(na.omit(d$aligned)))) + for (g in groups) { + y <- NULL; cells <- NULL; samples <- NULL + for (d in data) { + members <- which(d$aligned == g) + part <- d$counts[match(genes, d$genes), , members, drop = FALSE] + y <- cbind(y, apply(part, c(1, 2), sum)) + cells <- c(cells, rowSums(d$cells[, members, drop = FALSE])) + s <- d$samples + s$donor <- paste0(d$name, "::", s$donor) + samples <- rbind(samples, s) + } + samples$donor <- factor(samples$donor, levels = sort(unique(samples$donor))) + samples$condition <- factor(samples$condition, levels = c("ctrl", "stim")) + used <- cells >= 10 & colSums(y) > 0 + X <- model.matrix(as.formula(a$formula), samples[used, , drop = FALSE]) + dge <- DGEList(y[, used, drop = FALSE]) + rownames(dge) <- genes + dge <- dge[filterByExpr(dge, design = X), , keep.lib.sizes = FALSE] + dge <- calcNormFactors(dge) + fit <- eBayes(lmFit(voom(dge, X, adaptive.span = TRUE), X)) + tt <- topTable(fit, coef = "conditionstim", number = Inf, sort.by = "none") + se <- sqrt(fit$s2.post) * fit$stdev.unscaled[, "conditionstim"] + out <- file.path(root, "r", paste0(a$name, "__", gsub(":", "_", g))) + dir.create(out) + writeLines(rownames(tt), file.path(out, "genes.txt")) + wr(cbind(as.matrix(tt[, c("logFC", "AveExpr", "t", "P.Value", "adj.P.Val", + "B")]), se), file.path(out, "table.bin")) + } +} +writeLines(c(paste("R", getRversion()), paste("limma", packageVersion("limma")), + paste("edgeR", packageVersion("edgeR"))), file.path(root, "versions.txt")) +""" + + +def simulate(rng, name, genes, groups, samples, cells, unmeasured=()): + """Pool simulated cells of one study: a profile per cell type, a shift per + donor, and a stimulation that changes a tenth of the genes.""" + n = len(genes) + stim = np.where(rng.random(n) < 0.1, rng.normal(0.0, 1.2, n), 0.0) + shift = {d: rng.normal(0.0, 0.3, n) for d in dict.fromkeys(samples["donor"])} + base = rng.normal(0.0, 1.5, n) + blocks, sample_labels, group_labels = [], [], [] + for g, group in enumerate(groups["group_id"]): + profile = base + rng.normal(0.0, 0.7, n) + for s, sample in enumerate(samples["sample_id"]): + log_mean = profile + shift[samples["donor"][s]] + log_mean = log_mean + stim * (samples["condition"][s] == "stim") + share = np.exp(log_mean) / np.exp(log_mean).sum() + depth = rng.uniform(1000.0, 4000.0, (cells[s, g], 1)) + counts = rng.poisson(rng.gamma(4.0, 0.25, (cells[s, g], n)) * depth * share) + counts[:, list(unmeasured)] = 0 + empty = np.flatnonzero(counts.sum(axis=1) == 0) + counts[empty, n - 1] = 1 + blocks.append(counts) + sample_labels += [sample] * cells[s, g] + group_labels += [group] * cells[s, g] + measured = np.ones(n, dtype=bool) + measured[list(unmeasured)] = False + return from_cells( + np.vstack(blocks), + samples=sample_labels, + groups=group_labels, + sample_table=samples, + group_table=groups, + features={"feature_id": genes, "measured": measured}, + source={"study": name, "simulation": "validation/cube_collection.py"}, + ) + + +def build(seed, sizes): + """The three studies' cubes, from ``sizes`` = (genes, cells per entry).""" + rng = np.random.default_rng(seed) + n, per_entry = sizes + universe = [f"ENSG{k:011d}" for k in range(1, n + 1)] + tenth = n // 10 + paired = { + "lupus": {"donors": 4, "genes": (0, 7 * tenth), "age": np.int64}, + "covid": {"donors": 5, "genes": (tenth, 17 * tenth // 2), "age": np.float64}, + } + cubes = {} + for name, spec in paired.items(): + donors = [f"D{k + 1}" for k in range(spec["donors"]) for _ in range(2)] + samples = { + "sample_id": [f"S{k + 1}" for k in range(len(donors))], + "donor": donors, + "condition": ["ctrl", "stim"] * spec["donors"], + "age": np.array([30 + 3 * k for k in range(len(donors))], spec["age"]), + } + if name == "lupus": + samples["sex"] = ["female", "female", "male", "male"] * 2 + groups = { + "group_id": ["B cell", "CD4 T cell", "CD8 T cell", "monocyte"], + "cell_ontology_id": [*ONTOLOGY.values(), MONOCYTE], + } + else: + samples["age"][3] = np.nan + samples["batch"] = ["b1"] * 4 + ["b2"] * 6 + # Two kinds of CD4 T cells, which share an ontology id and merge. + groups = { + "group_id": ["B", "T4 naive", "T4 memory", "T8", "NK"], + "cell_ontology_id": [ + ONTOLOGY["B"], ONTOLOGY["CD4"], ONTOLOGY["CD4"], ONTOLOGY["CD8"], NK + ], + } # fmt: skip + cells = rng.poisson(per_entry, (len(donors), len(groups["group_id"]))) + cells[1, 0] = 4 # too few cells for limma_voom's min_cells + start, stop = spec["genes"] + unmeasured = range(1, 6 * tenth - 1, 15) if name == "covid" else () + cubes[name] = simulate( + rng, name, universe[start:stop], groups, samples, cells, unmeasured + ) + samples = { + "sample_id": [f"S{k + 1}" for k in range(6)], + "donor": np.arange(101, 107), + "condition": ["ctrl", "stim"] * 3, + "age": np.arange(40, 46), + } + groups = {"group_id": list(MAPPING["atlas"])} + cells = rng.poisson(per_entry, (6, 5)) + cells[2, 1] = 0 # an entry without cells, of a cell type every study has + cubes["atlas"] = simulate( + rng, "atlas", universe[tenth // 2 :], groups, samples, cells + ) + collection = CubeCollection( + cubes, + align_groups_by="cell_ontology_id", + group_mapping=MAPPING, + within_study=["donor"], + ) + return cubes, collection + + +def run_references(directory, cubes): + spec = { + "studies": list(cubes), + "mapping": MAPPING, + "within_study": ["donor"], + "analyses": [ + {"name": name, "studies": [name], "formula": FORMULAS[name]} + for name in cubes + ] + + [ + { + "name": "joint", + "studies": ["lupus", "covid"], + "formula": FORMULAS["lupus"], + } + ], + } + samples = {} + for name, cube in cubes.items(): + to_anndata(cube).write_h5ad(directory / f"{name}.h5ad") + samples[name] = { + column: {"dtype": str(values.dtype), "values": values.tolist()} + for column, values in cube.samples.to_dict().items() + } + folder = directory / "r" / name + folder.mkdir(parents=True) + groups = cube.groups + ontology = groups["cell_ontology_id"] if "cell_ontology_id" in groups else None + texts = { + "genes.txt": cube.features["feature_id"].tolist(), + "groups.txt": groups["group_id"].tolist(), + "ontology.txt": ["NA"] * cube.n_groups if ontology is None else ontology, + "measured.txt": [str(bool(m)).upper() for m in cube.features["measured"]], + } + for file, lines in texts.items(): + (folder / file).write_text("\n".join(lines) + "\n") + table = cube.samples + rows = zip(table["sample_id"], table["donor"], table["condition"], strict=True) + lines = ["sample_id\tdonor\tcondition"] + ["\t".join(map(str, r)) for r in rows] + (folder / "samples.tsv").write_text("\n".join(lines) + "\n") + (folder / "counts.bin").write_bytes( + cube.layer("counts").astype("= 0]] = "tested" + wrong += int((mine["status"] != expected).sum()) + positions.count(-1) + for k, column in enumerate(COLUMNS): + ours = mine[column][positions] + largest = max(largest, difference(ours, values[k], "P" in TABLE[k])) + check, *rest = CHECKS["limma"] + rows.append((f"{study}: {check}, {fitted} cell types", *rest, largest, EXACT)) + check, *rest = CHECKS["status"] + rows.append((f"{study}: {check}", *rest, wrong, 0)) + joint = collection.select(studies=["lupus", "covid"]).stack() + own_donors = [d.split("::")[1] for d in joint.samples["donor"]] + for key, samples in ( + ("joint", joint.samples), + ("control", joint.samples.with_columns(donor=own_donors)), + ): + design = Design(samples, FORMULAS["lupus"], levels=levels) + result = limma_voom(joint, design, contrasts=contrast) + # The control compares log2 fold changes, whose difference reads plainly, + # over the genes both tested: its design filters other genes. + columns = COLUMNS[:1] if key == "control" else COLUMNS + largest, other_genes = 0.0, 0 + for group in joint.groups["group_id"].tolist(): + genes, values = reference["limma"][f"joint__{group.replace(':', '_')}"] + table = result.tables[(group, "stim")] + ours = {gene: k for k, gene in enumerate(table["feature_id"].tolist())} + both = [k for k, gene in enumerate(genes) if gene in ours] + other_genes += len(genes) + len(ours) - 2 * len(both) + for k, column in enumerate(columns): + mine = table[column][[ours[genes[j]] for j in both]] + p_value = "P" in TABLE[k] + largest = max(largest, difference(mine, values[k][both], p_value)) + if key == "joint" and other_genes: + largest = np.inf # the genes tested must be the same + rows.append((*CHECKS[key], largest, EXACT if key == "joint" else "control")) + return rows + + +def print_table(rows): + print("| check | reference | metric | value | tolerance | result |") + print("|---|---|---|---|---|---|") + passed = True + for check, reference, metric, value, tolerance in rows: + shown = f"{value:.2g}" if isinstance(value, float) else str(value) + if tolerance == "control": # must exceed the exact tier + ok, limit = bool(value > EXACT), f"> {EXACT:.0e}" + elif isinstance(tolerance, bool): + ok, limit = bool(value == tolerance), str(tolerance) + else: + ok, limit = bool(value <= tolerance), f"{tolerance:.2g}" + result = ("control" if tolerance == "control" else "pass") if ok else "FAIL" + passed &= ok + print(f"| {check} | {reference} | {metric} | {shown} | {limit} | {result} |") + return passed + + +def write_fixtures(directory): + """Small studies, and what anndata, pandas and limma make of them, for tests: + the integer layers and cell numbers, with counts -1 where a study lacks a gene + (NaN in the outer join).""" + cubes, _ = build(SEED + 1, (40, 14)) + reference, _ = run_references(directory, cubes) + FIXTURE_DIR.mkdir(parents=True, exist_ok=True) + arrays = {} + for name, cube in cubes.items(): + for layer in ("counts", "detected"): + arrays[f"{name}/{layer}"] = cube.layer(layer) + arrays[f"{name}/n_cells"] = cube.n_cells + outer = reference["outer_counts"] + arrays["outer_counts"] = np.where(np.isnan(outer), -1, outer).astype(np.int64) + for key in ("first_listed", "outer_genes", "outer_study", "outer_sample_id"): + arrays[key] = reference[key] + arrays["outer_group"] = reference["outer_group"] + for key in ("genes", "study", "sample_id", "group", "counts", "detected"): + arrays[f"inner_{key}"] = reference[f"inner_{key}"] + arrays["inner_n_cells"] = reference["inner_n_cells"].astype(np.int64) + for key in ("counts", "detected"): + arrays[f"inner_{key}"] = arrays[f"inner_{key}"].astype(np.int64) + for key, (genes, values) in reference["limma"].items(): + arrays[f"limma/{key}/genes"] = np.array(genes) + arrays[f"limma/{key}/values"] = values + np.savez_compressed(FIXTURE_DIR / "studies.npz", **arrays) + tables = { + name: { + "samples": {c: v.tolist() for c, v in cube.samples.to_dict().items()}, + "groups": {c: v.tolist() for c, v in cube.groups.to_dict().items()}, + "features": {c: v.tolist() for c, v in cube.features.to_dict().items()}, + } + for name, cube in cubes.items() + } + record = {"studies": tables, "samples": reference["samples"], "mapping": MAPPING} + (FIXTURE_DIR / "studies.json").write_text(json.dumps(record, indent=1) + "\n") + print(f"wrote {FIXTURE_DIR / 'studies.npz'} and {FIXTURE_DIR / 'studies.json'}") + + +def main(): + parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + parser.add_argument( + "--write-fixtures", action="store_true", help="write test fixtures" + ) + arguments = parser.parse_args() + with tempfile.TemporaryDirectory() as temporary: + cubes, collection = build(SEED, (1000, 30)) + reference, r_versions = run_references(Path(temporary), cubes) + if arguments.write_fixtures: + with tempfile.TemporaryDirectory() as small: + write_fixtures(Path(small)) + versions = ", ".join(f"{k} {v}" for k, v in reference["versions"].items()) + print(f"Reference versions: {versions}; {', '.join(r_versions)}") + print(f"CytePool side: numpy {np.__version__}") + print(repr(collection)) + print() + rows = check_alignment(collection, reference) + check_limma(collection, reference) + return 0 if print_table(rows) else 1 + + +if __name__ == "__main__": + sys.exit(main())