diff --git a/README.md b/README.md index 38c7d18..a6b2163 100644 --- a/README.md +++ b/README.md @@ -38,6 +38,7 @@ what is planned. | Replication across studies | `sumrank`: SumRank (Nakatsuka et al. 2025), which ranks every gene within each study by its signed -log10 p-value and sums its best normalized ranks across studies, for genes up and down, with p-values from the Irwin-Hall distribution and from a null pooled over permutations; `differential_expression_sumrank`: the same for one cell type and contrast of the studies of a `CubeCollection`, with permutations that re-run each study's limma analysis from a seeded generator; `held_out_prediction`: each study left out in turn, genes chosen from the others by SumRank and by meta-analysis, its cases scored by UCell's rank score, and the AUC by study and method; with SumRank's Manhattan plot and the held-out AUC of every study | the authors' code, mpmath, scikit-learn | | Composition tests | `propeller`: differences in cell-type proportions between conditions, from the cell numbers of a cube or a table of them, with the logit or arcsine square-root transform, a design matched by sample id, t-tests of its contrasts and the F-test of a factor, robust empirical Bayes moderation, speckle's columns (baseline and mean proportions, proportion ratio, statistic, p-value, FDR) and a provenance record | speckle, limma | | Figures | `cytepool.plotting` (optional extra `plot`, matplotlib): the plots of limma, edgeR, DESeq2, speckle and metafor from CytePool's results, in one call each: MDS plots of pseudobulk samples (`plotMDS`), the cells pooled into each sample and cell type, cell-type composition and `propeller`'s results, voom's mean-variance trend, residual variances and their prior (`plotSA`), sample quality weights, MD and volcano plots, p-value histograms, the genes up and down in each cell type, barcode plots and summaries of gene-set tests, DESeq2's dispersions (`plotDispEsts`), edgeR's quasi-dispersions (`plotQLDisp`), forest plots of meta-analyses, the fractions of variance of every variable over the genes (`plotVarPart`) and of chosen genes (`plotPercentBars`), and heatmaps of pathway scores | limma, edgeR, DESeq2, metafor, speckle, variancePartition, GSVA | +| Samples and cells at a glance | `principal_components` and `plot_pca`: samples on the principal components of their most variable genes, from a matrix or a cell type of a cube, as DESeq2's `plotPCA` computes and draws them; `plot_dotplot`: the fraction of cells that express each gene and their mean expression in each cell type, from a cube's sums and so without the cells, as scanpy's `dotplot` draws them | DESeq2, edgeR, scanpy | | Distributions and smoothers | Tails of the normal, t, F, chi-squared, gamma, beta and hypergeometric distributions, and quantiles of the first four; LOWESS, and limma's rule for its span; natural spline bases | mpmath, R, limma | ## Planned diff --git a/THIRD_PARTY_NOTICES.md b/THIRD_PARTY_NOTICES.md index d2776bb..473d9f9 100644 --- a/THIRD_PARTY_NOTICES.md +++ b/THIRD_PARTY_NOTICES.md @@ -187,6 +187,49 @@ of neurodegenerative diseases through meta-analysis. Nature Communications 16, 7 doi:10.1038/s41467-025-62579-z. The code is distributed under the Apache License, Version 2.0 (Copyright 2024 Nathan Nakatsuka; http://www.apache.org/licenses/LICENSE-2.0). +### scanpy + +`plot_dotplot` in `cytepool/plotting.py` follows three conventions of scanpy 1.12.4's +dot plot (https://github.com/scverse/scanpy, `src/scanpy/plotting/_dotplot.py`) that its +documentation leaves open: the largest fraction rounded up to a tenth, the steps of the +key of the dots' areas (0.2, 0.1 or 0.05), and a scaled row of equal values drawn as 0. +No code is copied. scanpy is described in Wolf FA, Angerer P, Theis FJ (2018). SCANPY: +large-scale single-cell gene expression data analysis. Genome Biology 19, 15. +doi:10.1186/s13059-017-1382-0. It is distributed under the BSD 3-Clause License: + +```text +BSD 3-Clause License + +Copyright (c) 2025 scverse® +Copyright (c) 2017 F. Alexander Wolf, P. Angerer, Theis Lab +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +* Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + +* Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + +* Neither the name of the copyright holder nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +``` + ### numba `cytepool/_jit.py` places the cache of its compiled routines with subclasses of numba's diff --git a/cytepool/plotting.py b/cytepool/plotting.py index 91f8468..3cefa75 100644 --- a/cytepool/plotting.py +++ b/cytepool/plotting.py @@ -17,8 +17,12 @@ The figures: - Samples: :func:`plot_mds`, the multidimensional scaling plot of limma's ``plotMDS`` - (Ritchie et al. 2015), whose distances are leading log2 fold changes. -- Cubes: :func:`plot_cell_numbers`, the cells pooled into each sample and group. + (Ritchie et al. 2015), whose distances are leading log2 fold changes; + :func:`plot_pca`, the principal component plot of DESeq2's ``plotPCA`` (Love et al. + 2014; Pearson 1901). +- Cubes: :func:`plot_cell_numbers`, the cells pooled into each sample and group; + :func:`plot_dotplot`, scanpy's dot plot of the fraction of cells that express each + gene and its mean expression in each group (Wolf et al. 2018), from the cube's sums. - Composition: :func:`plot_composition`, the cell-type proportions of every sample, as speckle's ``plotCellTypeProps`` draws them; :func:`plot_proportions`, the proportions of every cell type by condition, with the adjusted p-values of ``propeller`` @@ -183,6 +187,13 @@ 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 + + Pearson K (1901). On lines and planes of closest fit to systems of points in space. + The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 2, + 559-572. doi:10.1080/14786440109462720 + + Wolf FA, Angerer P, Theis FJ (2018). SCANPY: large-scale single-cell gene expression + data analysis. Genome Biology 19, 15. doi:10.1186/s13059-017-1382-0 """ import importlib @@ -220,6 +231,7 @@ from .meta_analysis import MetaAnalysisResult from .normalization import filter_by_expression, log_cpm, tmm_factors from .pathway_scores import PathwayScores +from .principal_components import PrincipalComponents, principal_components from .replication import ReplicationResult from .stats.deseq2 import DESeq2Dispersions, DESeq2Fit, DESeq2Results from .stats.distributions import f_cdf, f_sf, norm_isf @@ -244,6 +256,7 @@ "plot_cell_numbers", "plot_composition", "plot_dispersions", + "plot_dotplot", "plot_enrichment", "plot_forest", "plot_gene_sets", @@ -256,6 +269,7 @@ "plot_over_representation", "plot_p_values", "plot_pathway_scores", + "plot_pca", "plot_proportions", "plot_quasi_dispersions", "plot_replicated_genes", @@ -4369,3 +4383,594 @@ def plot_held_out_auc( ax.set_xlabel("AUC of the held-out study's cases against its controls") _legend(ax, handles) return ax + + +# Principal components and dot plots. + +# DESeq2's plotPCA takes the 500 genes of largest variance. +_PCA_TOP = 500 +# scanpy's dot plots: the area of the dot of the largest fraction, in points squared, +# and the power of the fraction that sets a dot's area. +_DOT_AREA = 200.0 +_DOT_EXPONENT = 1.5 +_DOT_SCALES = ("gene", "group") + + +def plot_pca( + data: Any, + *, + group: str | None = None, + color_by: str | Sequence[Any] | None = None, + marker_by: str | Sequence[Any] | None = None, + labels: bool | Sequence[Any] = False, + top: int | None = None, + components: Sequence[int] = (1, 2), + min_cells: int | None = None, + ax: "Axes | None" = None, +) -> "Axes": + """Plot samples on two principal components of their most variable genes. + + DESeq2's ``plotPCA`` (Love et al. 2014): the samples' coordinates on two components + of :func:`~cytepool.principal_components.principal_components`, which takes the + ``top`` genes of largest variance over the samples, centres each on its mean, and + decomposes them by singular values (Pearson 1901); the sign of each component makes + its largest coordinate positive. The axes are labelled as ``plotPCA`` labels them, + with each component's share of the variance rounded to a whole percent, as in + "PC1: 45% variance", and drawn to the same scale, as ``plotPCA`` draws them, so that + the distances between samples are drawn true. The coordinates are those of + ``principal_components``; pass its result as ``data`` to draw the components it + holds. + + Draws the samples as points (gid ``"samples"``), coloured by ``color_by`` and + marked by ``marker_by`` as :func:`plot_mds` draws them, with a legend of their + levels, and labels the samples (gid ``"sample labels"``) when asked. + + Args: + data: A :class:`~cytepool.principal_components.PrincipalComponents`, or what + ``principal_components`` takes: a :class:`~cytepool.cube.Cube`, a + :class:`~cytepool.stats.voom.VoomResult`, or ``(n_samples, n_genes)`` + log-expression values, finite. + group: The group (cell type) of a cube to plot; it may be left out when the + cube has one group. Only cubes take it. + color_by: The variable whose levels colour the points: a column of the + samples' table (a cube's, or that of a PrincipalComponents), or one label + per sample (of the cube, for a cube). With more than three levels and no + ``marker_by``, each level also gets a marker. + marker_by: The variable whose levels mark the points, as ``color_by``; at most + 12 levels. + labels: True to write each sample's id next to it (numbers from 1 for an + array), or one label per sample. + top: The number of genes of largest variance used, at least 1; 500 by + default, as DESeq2's ``ntop``. Not for a PrincipalComponents. + components: The two components to plot, counted from 1 in decreasing order + of the variance they explain. + min_cells: The fewest cells a sample needs in the group to be plotted, 10 by + default as for ``limma_voom``; for a cube that records cell numbers. + ax: The Axes to draw into, or None for a new figure. + + Returns: + The Axes. + + Raises: + TypeError: If an argument has the wrong type, or is or holds a numpy masked + array. + ValueError: As ``principal_components``, or if the data have fewer components + than asked for, ``group``, ``top`` or ``min_cells`` is given with a + PrincipalComponents, or a sample variable has the wrong length, a missing + label or too many levels. + KeyError: If the group, or a variable named, is not in the data. + """ + _check_axes(ax) + first, second = _component_pair(components) + _check_variable(color_by, "color_by") + _check_variable(marker_by, "marker_by") + if not isinstance(labels, bool | np.bool_): + reject_masked(labels, "labels") + if isinstance(data, PrincipalComponents): + for name, value in (("group", group), ("top", top), ("min_cells", min_cells)): + if value is not None: + reject_masked(value, name) + raise ValueError( + f"{name} applies to the values that principal_components takes; " + "data holds its components already" + ) + result = data + else: + result = principal_components( + data, + group=group, + top=_PCA_TOP if top is None else top, + min_cells=min_cells, + ) + if result.n_components < max(first, second): + raise ValueError( + f"components {first} and {second} need {max(first, second)} components; " + f"the data have {result.n_components}, one fewer than the samples and at " + "most the genes used" + ) + samples = _component_samples(data, result) + x = np.asarray(result.coordinates[:, first - 1], dtype=np.float64) + y = np.asarray(result.coordinates[:, second - 1], dtype=np.float64) + ax = _axes(ax, (6.5, 4.5)) + _draw_samples(ax, x, y, samples, color_by, marker_by, labels) + share = result.variance_explained + ax.set_xlabel(f"PC{first}: {100 * share[first - 1]:.0f}% variance") + ax.set_ylabel(f"PC{second}: {100 * share[second - 1]:.0f}% variance") + ax.set_aspect("equal", adjustable="datalim") + return ax + + +def _component_pair(value: Any) -> tuple[int, int]: + reject_masked(value, "components") + if isinstance(value, str) or not isinstance(value, Sequence) or len(value) != 2: + raise TypeError( + f"components must be a pair of components, such as (1, 2), got {value!r}" + ) + first = positive_integer(value[0], "components[0]") + second = positive_integer(value[1], "components[1]") + if first == second: + raise ValueError(f"components must be two different components, got {value!r}") + return first, second + + +def _component_samples(data: Any, result: PrincipalComponents) -> _Samples: + """The samples of principal components, with the variables of their table: a + cube's sample table for a cube, so that labels follow its samples.""" + table = result.samples + if isinstance(data, Cube): + position = { + value: i for i, value in enumerate(data.samples["sample_id"].tolist()) + } + rows = [position[value] for value in table["sample_id"].tolist()] + return _cube_samples(data, np.array(rows, dtype=np.intp)) + n = len(table) + ids = [str(value) for value in table["sample_id"].tolist()] + if not isinstance(data, PrincipalComponents): + return _Samples(ids, list(range(n)), n) + return _Samples( + ids=ids, + rows=list(range(n)), + n_rows=n, + variables={name: table[name] for name in table.columns}, + source="the samples' table", + ) + + +def _draw_samples( + ax: "Axes", + x: FloatArray, + y: FloatArray, + samples: _Samples, + color_by: Any, + marker_by: Any, + labels: Any, +) -> None: + """Samples as points, coloured and marked by their variables, as plot_mds draws + them: a hue per level of color_by, a marker per level of marker_by, or of color_by + when it has more levels than scatter plots keep apart; labels when asked; and a + legend of the levels.""" + n_samples = x.size + colors = None if color_by is None else samples.labels(color_by, "color_by") + markers = None if marker_by is None else samples.labels(marker_by, "marker_by") + remedy = "give it as marker_by, or as labels" + hues = {} if colors is None else _hues(colors[1], "color_by", remedy) + shapes_follow_colors = ( + markers is None and colors is not None and len(colors[1]) > _SCATTER_HUES + ) + if markers is not None: + shapes, by_shape = _shapes(markers[1], "marker_by"), markers[0] + elif colors is not None and shapes_follow_colors: + shapes, by_shape = _shapes(colors[1], "color_by"), colors[0] + else: + shapes, by_shape = {}, None + if isinstance(labels, bool | np.bool_): + texts = list(samples.ids) if labels else None + else: + texts = samples.labels(labels, "labels")[0] + color_of = [hues[level] for level in colors[0]] if colors else [_BLUE] * n_samples + marker_of = [shapes[level] for level in by_shape] if by_shape else None + marker_of = marker_of or [_MARKERS[0]] * n_samples + for marker in dict.fromkeys(marker_of): + rows = [i for i in range(n_samples) if marker_of[i] == marker] + ax.scatter( + x[rows], + y[rows], + c=[color_of[i] for i in rows], + marker=marker, + s=_SPARSE, + edgecolors="white", + linewidths=0.6, + zorder=3, + gid="samples", + ) + if texts is not None: + for i, text in enumerate(texts): + ax.annotate( + text, + (x[i], y[i]), + xytext=(4, 4), + textcoords="offset points", + fontsize=8, + gid="sample labels", + ) + handles: list[Any] = [] + if colors is not None: + if isinstance(color_by, str): + handles.append(_key(color_by, marker="")) + for level in colors[1]: + shape = shapes[level] if shapes_follow_colors else _MARKERS[0] + handles.append(_key(level, marker=shape, color=hues[level])) + if markers is not None: + if isinstance(marker_by, str): + handles.append(_key(marker_by, marker="")) + for level in markers[1]: + handles.append(_key(level, marker=shapes[level], color=_GREY)) + if handles: + _legend(ax, handles) + + +def plot_dotplot( + cube: Cube, + genes: Any, + *, + groups: Sequence[str] | None = None, + samples: Any = None, + standard_scale: str | None = None, + mean_only_expressed: bool = False, + ax: "Axes | None" = None, +) -> "Axes": + """Plot how many cells of each group express each gene, and how much, as dots. + + scanpy's dot plot (``scanpy.pl.dotplot``; Wolf et al. 2018), drawn from the sums + that a cube keeps instead of from the cells: for each group (cell type), a row, + and each gene, a column, the fraction of the group's cells with a positive count, + ``detected`` over ``n_cells``, sets the area of the dot, and the mean of the cells' + log-normalized values, ``log_sum`` over ``n_cells``, its colour. Both are pooled + over the samples, every sample's sums added before dividing, so that they are the + fractions and means scanpy computes from the cells the cube was pooled from (with + its default ``expression_cutoff`` of 0); ``samples`` pools some of them only. + A dot's area is 200 (f / f_max)^1.5 points squared, f its fraction and f_max the + largest fraction shown rounded up to a tenth, as scanpy sizes its dots, and the + key of the areas shows fractions in scanpy's steps. The colours run from the + smallest to the largest value shown. + + ``standard_scale`` scales the means of each gene over the groups (``"gene"``, + scanpy's ``"var"``), or of each group over the genes (``"group"``), to run from 0 + to 1, (m - min) / (max - min), 0 where they are all equal, as scanpy does. + ``mean_only_expressed`` takes the mean over the cells that express the gene, + ``log_sum`` over ``detected`` (0 where none does), as scanpy's option of that name; + it needs a cube whose transform maps a zero count to zero, one whose normalization + records a pseudocount of 1, as :func:`~cytepool.cube.from_cells` records it. Genes + given as a mapping are + drawn in groups, under a bracket labelled with the group's name, as scanpy draws + ``var_names`` given as a dict. + + Draws the dots (gid ``"dots"``: their sizes are the areas, and their array the + values that the colours show) and, for genes in groups, the brackets (gid + ``"gene groups"``) and their names (gid ``"gene group names"``). + + Args: + cube: A cube with cell numbers and the layers ``detected`` and ``log_sum``, + as :func:`~cytepool.cube.from_cells` builds it. + genes: The genes, the columns from left to right: feature ids or symbols of + the cube's features; or a mapping of names to such lists, for groups of + genes. + groups: The groups, the rows from the top; by default every group, in the + cube's order. + samples: The samples whose cells are pooled: sample ids, or a boolean mask + over the cube's samples; by default every sample. + standard_scale: ``"gene"``, ``"group"`` or None. + mean_only_expressed: Whether the mean is taken over the cells that express + the gene only. + ax: The Axes to draw into, or None for a new figure. + + Returns: + The Axes. + + Raises: + TypeError: If an argument has the wrong type, or is or holds a numpy masked + array. + ValueError: If the cube lacks cell numbers, ``detected`` or ``log_sum``, a + group has no cells in the samples pooled, ``genes``, ``groups`` or + ``samples`` is empty or names one twice, a symbol names several genes, + ``mean_only_expressed`` is asked of a cube without a pseudocount of 1, or + an option is out of range. + KeyError: If a gene, group or sample is not in the cube. + """ + _check_axes(ax) + if not isinstance(cube, Cube): + raise TypeError(f"cube must be a Cube, got {type(cube).__name__}") + scale = ( + None + if standard_scale is None + else checked_choice(standard_scale, "standard_scale", _DOT_SCALES) + ) + expressed_only = checked_flag(mean_only_expressed, "mean_only_expressed") + lacking = [name for name in ("detected", "log_sum") if name not in cube.layer_names] + if cube.n_cells is None or lacking: + what = ( + "cell numbers" + if cube.n_cells is None + else f"the layer(s) {quoted(lacking)}" + ) + raise ValueError( + f"the cube lacks {what}; a dot plot needs cell numbers and the layers " + "'detected' and 'log_sum', as from_cells computes them" + ) + columns, names, gene_groups = _dot_genes(cube, genes) + rows = _chosen_ids(groups, cube.groups["group_id"].tolist(), "groups", "group") + pooled = _pooled_samples(cube, samples) + if expressed_only: + normalization = cube.normalization or {} + if normalization.get("pseudocount") != 1: + raise ValueError( + "mean_only_expressed needs a cube whose transform maps a zero count to " + "zero, one whose normalization records a pseudocount of 1, as from_cells " + f"records it; this cube's normalization is {normalization or None!r}" + ) + fractions, means = _dot_values(cube, rows, columns, pooled, expressed_only) + shown = ( + means if scale is None else _min_max_scaled(means, 0 if scale == "gene" else 1) + ) + + largest = float(fractions.max()) + f_max = math.ceil(largest * 10) / 10 if largest > 0.0 else 1.0 + areas = _DOT_AREA * (fractions / f_max) ** _DOT_EXPONENT + n_groups, n_genes = fractions.shape + colors = _module("matplotlib.colors") + ramp = colors.LinearSegmentedColormap.from_list("cytepool_blues", _BLUES) + low, high = float(shown.min()), float(shown.max()) + norm = colors.Normalize(vmin=low, vmax=high if high > low else low + 1.0) + width = min(14.0, 3.0 + 0.4 * n_genes) + ax = _axes(ax, (width, min(12.0, max(4.0, 1.6 + 0.4 * n_groups)))) + y, x = np.indices((n_groups, n_genes)) + ax.scatter( + x.ravel(), + y.ravel(), + s=areas.ravel(), + c=shown.ravel(), + cmap=ramp, + norm=norm, + edgecolors=_INK, + linewidths=0.2, + gid="dots", + ) + group_ids = cube.groups["group_id"].tolist() + ax.set_yticks(range(n_groups), [group_ids[g] for g in rows]) + ax.set_xticks(range(n_genes), names, rotation=90) + ax.set_xlim(-0.6, n_genes - 0.4) + ax.set_ylim(n_groups - 0.4, -0.6) + ax.tick_params(length=0) + if gene_groups: + _draw_gene_groups(ax, gene_groups) + if scale is not None: + label = f"Mean expression\nin group, scaled\nper {scale}" + elif expressed_only: + label = "Mean expression in\nexpressing cells" + else: + label = "Mean expression\nin group" + mappable = _module("matplotlib.cm").ScalarMappable(norm=norm, cmap=ramp) + # The bar starts below the top of the plot, to leave room for its label above + # it, where scanpy writes it; a label along the bar ran into the key below it. + bar = ax.figure.colorbar(mappable, ax=ax, shrink=0.4, anchor=(0.0, 0.85)) + handles = [ + _key( + f"{round(100 * share):d}", + color=_GREY, + markeredgecolor=_INK, + markeredgewidth=0.2, + markersize=math.sqrt(_DOT_AREA * (share / f_max) ** _DOT_EXPONENT), + ) + for share in _fraction_keys(f_max) + ] + key = bar.ax.legend( + handles=handles, + title="Fraction of cells\nin group (%)", + loc="upper left", + bbox_to_anchor=(0.0, -0.06), + borderaxespad=0.0, + frameon=False, + alignment="left", + ) + key.get_title().set_fontweight("bold") + bar.ax.set_title( + label, + loc="left", + fontsize=key.get_title().get_fontsize(), + fontweight="bold", + ) + return ax + + +def _dot_genes( + cube: Cube, genes: Any +) -> tuple[list[int], list[str], list[tuple[str, int, int]]]: + """The cube's features of the genes asked for, the names shown, and the groups of + genes as (name, first column, last column).""" + reject_masked(genes, "genes") + if isinstance(genes, str): + named: dict[str, Any] = {"": [genes]} + elif isinstance(genes, Mapping): + named = dict(genes) + else: + named = {"": genes} + features = cube.features + ids = features["feature_id"].tolist() + index = {value: j for j, value in enumerate(ids)} + symbols = features["symbol"].tolist() if "symbol" in features else [] + by_symbol: dict[Any, list[int]] = {} + for j, symbol in enumerate(symbols): + by_symbol.setdefault(symbol, []).append(j) + columns: list[int] = [] + names: list[str] = [] + gene_groups: list[tuple[str, int, int]] = [] + for label, members in named.items(): + if not isinstance(label, str): + raise TypeError( + f"the names of groups of genes must be text, got {type(label).__name__}" + ) + if isinstance(members, str): + members = [members] + if not isinstance(members, Sequence | np.ndarray): + raise TypeError( + "genes must be feature ids or symbols, or a mapping of names to them, " + f"got {type(members).__name__}" + ) + start = len(columns) + for gene in list(members): + if not isinstance(gene, str): + raise TypeError( + f"genes must be feature ids or symbols (str), got " + f"{type(gene).__name__} {gene!r}" + ) + if gene in index: + columns.append(index[gene]) + elif gene in by_symbol: + found = by_symbol[gene] + if len(found) > 1: + raise ValueError( + f"the symbol {gene!r} names {len(found)} genes, " + f"{quoted([ids[j] for j in found])}; give their feature ids" + ) + columns.append(found[0]) + else: + raise KeyError( + f"no gene {gene!r} among the cube's {len(ids):,} features, by " + "feature_id or symbol" + ) + names.append(gene) + if len(columns) == start: + raise ValueError( + "genes is empty; give at least one gene" + if not label + else f"the group of genes {label!r} is empty" + ) + if label: + gene_groups.append((label, start, len(columns) - 1)) + return columns, names, gene_groups + + +def _chosen_ids(value: Any, ids: list[str], name: str, kind: str) -> list[int]: + """The positions of the ids chosen, in the order given; every id for None.""" + if value is None: + return list(range(len(ids))) + reject_masked(value, name) + chosen = [value] if isinstance(value, str) else value + if not isinstance(chosen, Sequence | np.ndarray): + raise TypeError( + f"{name} must be a list of {kind} ids, got {type(value).__name__}" + ) + index = {identifier: i for i, identifier in enumerate(ids)} + positions = [] + for identifier in list(chosen): + if not isinstance(identifier, str): + raise TypeError( + f"{name} must hold {kind} ids (str), got {type(identifier).__name__} " + f"{identifier!r}" + ) + if identifier not in index: + raise KeyError(f"no {kind} {identifier!r}; {kind}s are: {quoted(ids)}") + positions.append(index[identifier]) + if not positions: + raise ValueError(f"{name} is empty; give at least one {kind}") + if len(set(positions)) != len(positions): + raise ValueError(f"{name} names a {kind} more than once") + return positions + + +def _pooled_samples(cube: Cube, samples: Any) -> list[int]: + """The positions of the samples pooled: ids, a boolean mask, or every sample.""" + ids = cube.samples["sample_id"].tolist() + if samples is not None and not isinstance(samples, str): + reject_masked(samples, "samples") + mask = np.asarray(samples) + if mask.dtype == bool: + if mask.shape != (len(ids),): + raise ValueError( + f"a boolean mask of samples needs one entry per sample " + f"({len(ids)}), got shape {mask.shape}" + ) + if not mask.any(): + raise ValueError("samples is empty; the mask selects no sample") + return np.flatnonzero(mask).tolist() + return _chosen_ids(samples, ids, "samples", "sample") + + +def _dot_values( + cube: Cube, + rows: list[int], + columns: list[int], + pooled: list[int], + expressed_only: bool, +) -> tuple[FloatArray, FloatArray]: + """The fraction of cells detected and the mean log value, groups x genes, each + from the sums over the samples pooled.""" + n_cells = np.asarray(cube.n_cells, dtype=np.int64) + fractions = np.empty((len(rows), len(columns))) + means = np.empty((len(rows), len(columns))) + group_ids = cube.groups["group_id"].tolist() + for i, g in enumerate(rows): + cells = int(n_cells[pooled, g].sum()) + if cells == 0: + raise ValueError( + f"group {group_ids[g]!r} has no cells in the samples pooled, so its " + "fractions and means are undefined; leave it out with groups" + ) + detected = np.asarray(cube.layer("detected", g))[pooled][:, columns] + log_sum = np.asarray(cube.layer("log_sum", g), dtype=np.float64)[pooled] + expressing = detected.sum(axis=0, dtype=np.int64).astype(np.float64) + total = log_sum[:, columns].sum(axis=0) + fractions[i] = expressing / cells + if expressed_only: + means[i] = np.divide( + total, expressing, out=np.zeros_like(total), where=expressing > 0 + ) + else: + means[i] = total / cells + return fractions, means + + +def _min_max_scaled(values: FloatArray, axis: int) -> FloatArray: + """(v - min) / (max - min) along an axis, 0 where every value is equal.""" + shifted = values - values.min(axis=axis, keepdims=True) + spread = shifted.max(axis=axis, keepdims=True) + scaled: FloatArray = np.divide( + shifted, spread, out=np.zeros_like(shifted), where=spread > 0.0 + ) + return scaled + + +def _fraction_keys(f_max: float) -> list[float]: + """The fractions of the key of a dot plot's areas, in scanpy's steps: 0.2 apart + from f_max down, 0.1 when f_max is at most 0.6 and 0.05 when at most 0.3.""" + step = 0.2 if f_max > 0.6 else (0.1 if f_max > 0.3 else 0.05) + return [float(value) for value in np.arange(f_max, 0.0, -step)[::-1]] + + +def _draw_gene_groups(ax: "Axes", gene_groups: list[tuple[str, int, int]]) -> None: + """Brackets above the columns of each group of genes, labelled with its name; + the names turned upright when any is longer than four characters.""" + transform = _text_transform(ax) + rotation = 90 if max(len(name) for name, _, _ in gene_groups) > 4 else 0 + for name, start, end in gene_groups: + left, right = start - 0.3, end + 0.3 + ax.plot( + [left, left, right, right], + [1.01, 1.04, 1.04, 1.01], + transform=transform, + color=_INK, + linewidth=0.8, + clip_on=False, + gid="gene groups", + ) + ax.text( + (left + right) / 2.0, + 1.06, + name, + transform=transform, + ha="center", + va="bottom", + rotation=rotation, + fontsize=9, + gid="gene group names", + ) diff --git a/cytepool/principal_components.py b/cytepool/principal_components.py new file mode 100644 index 0000000..0794c0f --- /dev/null +++ b/cytepool/principal_components.py @@ -0,0 +1,293 @@ +"""Principal components of samples, as DESeq2's ``plotPCA`` computes them. + +:func:`principal_components` places the samples of a matrix of log-expression values, or +of one group (cell type) of a cube, on their principal components (Pearson 1901), as +DESeq2's ``plotPCA`` does for its transformed values (Love et al. 2014): from the +``top`` genes of largest variance over the samples, each centred on its mean, by the +singular value decomposition of that samples x genes matrix, X = U S V'. The +coordinates of the samples are U S, the loadings of the genes are V, and the share of +the variance that each component explains is its squared singular value over their +sum, DESeq2's ``percentVar``. :func:`~cytepool.plotting.plot_pca` draws them. + +What DESeq2's help page leaves open was found by running ``plotPCA`` 1.46.0 as a black +box: it ranks the genes by their variance over the samples, keeps genes of equal +variance in their order, takes every gene when ``ntop`` exceeds their number, and +returns the coordinates of R's ``prcomp``, whose signs, its help page says, are +arbitrary. Here the sign of each component makes the sample with its largest absolute +coordinate positive (of near ties, within a relative 1e-9, the first), as +:func:`~cytepool.plotting.plot_mds` orients its dimensions, so that the result does +not depend on the order of the genes or the sign of the values. + +The values: + A matrix holds one row per sample and one column per gene, finite, such as + log-CPM values or the variance-stabilized values that DESeq2's guide takes. A + :class:`~cytepool.stats.voom.VoomResult` gives its log-CPM values. A cube gives + one group's samples prepared as :func:`~cytepool.plotting.plot_mds` prepares them, + as the edgeR User's Guide prepares counts for ``plotMDS``: the samples with fewer + than ``min_cells`` cells or no counts are left out, as + :func:`~cytepool.differential_expression.limma_voom` leaves them out; the genes are + kept by :func:`~cytepool.normalization.filter_by_expression` (Chen, Lun and Smyth + 2016) with the samples as one group; and the values are their log2 counts per + million with TMM normalization factors and an average count of 2 added + (:func:`~cytepool.normalization.log_cpm`). + +References: + Pearson K (1901). On lines and planes of closest fit to systems of points in space. + The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 2, + 559-572. doi:10.1080/14786440109462720 + + Love MI, Huber W, Anders S (2014). Moderated estimation of fold change and + dispersion for RNA-seq data with DESeq2. Genome Biology 15, 550. + doi:10.1186/s13059-014-0550-8 + + Chen Y, Lun ATL, Smyth GK (2016). From reads to genes to pathways: differential + expression analysis of RNA-Seq experiments using Rsubread and the edgeR + quasi-likelihood pipeline. F1000Research 5, 1438. doi:10.12688/f1000research.8987.2 +""" + +from collections.abc import Mapping +from dataclasses import dataclass +from types import MappingProxyType +from typing import Any + +import numpy as np +import numpy.typing as npt + +from ._columns import text_column +from ._masked_arrays import read_only, real_array, reject_masked, unmasked_array +from ._messages import quoted +from ._options import integral_number, positive_integer +from .cube import Cube +from .normalization import filter_by_expression, log_cpm, tmm_factors +from .stats.voom import VoomResult +from .table import Table + +__all__ = ["PrincipalComponents", "principal_components"] + +FloatArray = npt.NDArray[np.float64] +IntArray = npt.NDArray[np.intp] + +# DESeq2's plotPCA takes the 500 genes of largest variance. +_TOP = 500 +# A cube's samples need this many cells, as limma_voom leaves samples out, and its +# log-CPM values add an average count of 2, as edgeR's plotMDS does. +_MIN_CELLS = 10 +_PRIOR_COUNT = 2.0 +# Coordinates within this relative distance of a component's largest absolute one tie +# for it; the first of them sets the sign. +_NEAR_TIE = 1e-9 +# Centred values below this multiple of the largest value are rounding errors. +_ROUNDING = 16.0 * float(np.finfo(np.float64).eps) + + +@dataclass(frozen=True) +class PrincipalComponents: + """The samples on their principal components. + + Rows of ``coordinates`` are the samples used, in the order of ``samples``; columns + are the components, in decreasing order of the variance they explain. There are + min(n_samples - 1, n_genes) components: the centred values have no more + dimensions. The arrays are read-only. + + Attributes: + coordinates: ``(n_samples, n_components)`` float64, the samples' coordinates, + the centred values times the loadings; DESeq2's ``PC1``, ``PC2``, ... up + to sign. + variance_explained: ``(n_components,)`` float64, each component's share of the + total variance of the genes used, summing to 1; DESeq2's ``percentVar``. + loadings: ``(n_genes_used, n_components)`` float64, the unit vectors of the + components over the genes used, in the order of ``genes``. + genes: ``(n_genes_used,)`` positions of the genes used, in decreasing order of + their variance: columns of the matrix, or features of the cube. + samples: One row per sample used: ``sample_id`` and, for a cube, the columns + of its sample table; for a matrix, the ids ``"1"``, ``"2"``, ... of its + rows. + settings: ``top``, and for a cube its ``group`` and ``min_cells``. + """ + + coordinates: FloatArray + variance_explained: FloatArray + loadings: FloatArray + genes: IntArray + samples: Table + settings: Mapping[str, Any] + + @property + def n_components(self) -> int: + """The number of components.""" + return int(self.coordinates.shape[1]) + + def __repr__(self) -> str: + n_samples, n_components = self.coordinates.shape + shares = ", ".join(f"{100 * v:.1f}%" for v in self.variance_explained[:3]) + return ( + f"PrincipalComponents(n_samples={n_samples}, n_genes={self.genes.size}, " + f"n_components={n_components}, variance_explained=[{shares}" + f"{', ...' if n_components > 3 else ''}])" + ) + + +def principal_components( + data: Any, + *, + group: str | None = None, + top: int = _TOP, + min_cells: int | None = None, +) -> PrincipalComponents: + """Place samples on the principal components of their most variable genes. + + As DESeq2's ``plotPCA``: the ``top`` genes of largest variance over the samples + (every gene when there are fewer, genes of equal variance in their order), each + centred on its mean over the samples, decomposed by singular values. The module's + docstring describes the values a cube gives and the sign of each component. + + Args: + data: ``(n_samples, n_genes)`` log-expression values, finite; a + :class:`~cytepool.stats.voom.VoomResult`, whose log-CPM values are used; + or a :class:`~cytepool.cube.Cube`, of which one group's samples are used. + group: The group (cell type) of a cube; it may be left out when the cube has + one group. Only cubes take it. + top: The number of genes of largest variance used, at least 1 (DESeq2's + ``ntop``). + min_cells: The fewest cells a sample of a cube needs in the group, 10 by + default as for ``limma_voom``; for a cube that records cell numbers. + + Returns: + The coordinates, shares of variance and loadings, with the genes and samples + used. + + Raises: + TypeError: If an argument has the wrong type, or is or holds a numpy masked + array. + ValueError: If the values are not finite or not 2-D, fewer than two samples + remain, the group of a cube of several groups is not given, ``group`` or + ``min_cells`` is given for data other than a cube, no gene of a cube + passes the filter, the genes used do not vary, or an option is out of + range. + KeyError: If the group is not in the cube. + """ + count = positive_integer(top, "top") + values, genes, samples, settings = _values(data, group, min_cells) + n_samples, n_genes = values.shape + if n_samples < 2: + raise ValueError( + f"principal components need at least two samples, got {n_samples}" + ) + # The genes of largest variance; equal variances keep the genes' order (found by + # probing DESeq2 1.46.0). + variances = np.var(values, axis=0, ddof=1) + chosen = np.argsort(-variances, kind="stable")[: min(count, n_genes)] + used = values[:, chosen] + centred = used - used.mean(axis=0) + left, singular, right = np.linalg.svd(centred, full_matrices=False) + n_components = min(n_samples - 1, chosen.size) + singular = singular[:n_components] + total = float(np.sum(singular**2)) + # Centred values that are rounding errors of their means, as constant genes + # leave them, do not vary. + noise = (_ROUNDING * float(np.max(np.abs(used)))) ** 2 * centred.size + if not total > noise: + raise ValueError( + "the values of the genes used do not vary between the samples, so they " + "have no principal components" + ) + coordinates = left[:, :n_components] * singular + loadings = right[:n_components].T.copy() + signs = _signs(coordinates) + coordinates *= signs + loadings *= signs + settings["top"] = count + return PrincipalComponents( + coordinates=read_only(np.ascontiguousarray(coordinates)), + variance_explained=read_only(singular**2 / total), + loadings=read_only(np.ascontiguousarray(loadings)), + genes=read_only(genes[chosen]), + samples=samples, + settings=MappingProxyType(settings), + ) + + +def _signs(coordinates: FloatArray) -> FloatArray: + """The sign of each component that makes its largest absolute coordinate positive; + of near ties, the first.""" + magnitude = np.abs(coordinates) + largest = np.argmax(magnitude >= magnitude.max(axis=0) * (1.0 - _NEAR_TIE), axis=0) + signs = np.sign(coordinates[largest, np.arange(coordinates.shape[1])]) + result: FloatArray = np.where(signs == 0.0, 1.0, signs) + return result + + +def _values( + data: Any, group: Any, min_cells: Any +) -> tuple[FloatArray, IntArray, Table, dict[str, Any]]: + """The values (samples x genes), the positions of their genes among the data's, + the samples and the settings that chose them.""" + if isinstance(data, Cube): + return _cube_values(data, group, min_cells) + if group is not None or min_cells is not None: + given = "group" if group is not None else "min_cells" + if given == "group" and not isinstance(group, str): + reject_masked(group, "group") + raise ValueError( + f"{given} applies to a cube, and data is a {type(data).__name__}" + ) + source = data.log_cpm if isinstance(data, VoomResult) else data + values = real_array(unmasked_array(source, "data"), "data", boolean=False) + if values.ndim != 2: + raise ValueError( + "data must be a Cube, a VoomResult or a 2-D array of shape " + f"(n_samples, n_genes), got shape {values.shape}" + ) + if values.shape[1] == 0: + raise ValueError("data has no genes") + if not np.isfinite(values).all(): + raise ValueError( + "data must be finite; leave out the genes with missing values first" + ) + n_samples = values.shape[0] + samples = Table({"sample_id": text_column([str(i + 1) for i in range(n_samples)])}) + return values, np.arange(values.shape[1]), samples, {} + + +def _cube_values( + cube: Cube, group: Any, min_cells: Any +) -> tuple[FloatArray, IntArray, Table, dict[str, Any]]: + """A cube group's log-CPM values, of its samples with enough cells and counts, + over the genes that filterByExpr keeps.""" + fewest = ( + _MIN_CELLS if min_cells is None else integral_number(min_cells, "min_cells") + ) + if fewest < 0: + raise ValueError(f"min_cells must be at least 0, got {fewest}") + ids = cube.groups["group_id"].tolist() + if group is None: + if cube.n_groups != 1: + raise ValueError( + f"the cube has {cube.n_groups} groups; name the one to use with group, " + f"one of {quoted(ids)}" + ) + group = ids[0] + if not isinstance(group, str): + reject_masked(group, "group") + raise TypeError(f"group must be a group id (str), got {type(group).__name__}") + counts = np.asarray(cube.layer("counts", group), dtype=np.float64) + used = counts.sum(axis=1) > 0.0 + if cube.n_cells is not None: + used &= np.asarray(cube.n_cells)[:, ids.index(group)] >= fewest + rows = np.flatnonzero(used) + if rows.size < 2: + raise ValueError( + f"{rows.size} sample(s) of group {group!r} have counts and at least " + f"{fewest} cells; principal components need two or more" + ) + kept = counts[rows] + keep = filter_by_expression(kept, group=np.zeros(rows.size)) + if not keep.any(): + raise ValueError( + f"no gene of group {group!r} passes filter_by_expression, so there is " + "nothing to decompose" + ) + kept = kept[:, keep] + values = log_cpm(kept, norm_factors=tmm_factors(kept), prior_count=_PRIOR_COUNT) + settings: dict[str, Any] = {"group": group, "min_cells": fewest} + return values, np.flatnonzero(keep), cube.samples.take(rows), settings diff --git a/cytepool/stats/_scalar_special_functions.py b/cytepool/stats/_scalar_special_functions.py index 20d5c59..8f4ee2a 100644 --- a/cytepool/stats/_scalar_special_functions.py +++ b/cytepool/stats/_scalar_special_functions.py @@ -8,8 +8,9 @@ those that :mod:`cytepool.stats.special` describes: Dekker's (1971) exact products, Loader's (2000) deviance form of Stirling's series, Maechler's (2012) log(1 - exp(-a)), the power series and continued fractions of DLMF 8.7, 8.9 and 8.17 evaluated by the -modified Lentz method (Lentz 1976; Thompson and Barnett 1986), and Temme's (1979) -uniform expansion of the incomplete gamma function. +modified Lentz method (Lentz 1976; Thompson and Barnett 1986), and Temme's uniform +expansions of the incomplete gamma function (Temme 1979) and of the incomplete beta +function (Temme 1982). ``elementwise`` runs a compiled loop over broadcast arguments: every argument reaches the loop as a read-only, contiguous, one-dimensional float64 array, so that each loop @@ -39,6 +40,10 @@ Temme NM (1979). The asymptotic expansion of the incomplete gamma functions. SIAM Journal on Mathematical Analysis 10, 757-766. doi:10.1137/0510071 + + Temme NM (1982). The uniform asymptotic expansion of a class of integrals related to + cumulative distribution functions. SIAM Journal on Mathematical Analysis 13, + 239-253. doi:10.1137/0513017 """ import math @@ -73,6 +78,36 @@ _MAX_ITERATIONS = 10_000_000 +# Arguments beyond the range of the usual forms. Deviance terms halve their arguments +# from this sum (2^1023). two_product splits factors below 1e290 only; a larger factor +# is scaled by 2^-600 first and its parts by 2^600 after. From shapes of 2^500, the +# products of the incomplete beta function's continued fraction and prefactor could +# overflow, so such shapes take forms of them ordered to stay in range. +_HALF_RANGE = 8.98846567431158e307 +_SPLIT_LIMIT = 1e290 +_SCALE_DOWN = 2.0**-600 +_SCALE_UP = 2.0**600 +_HUGE_SHAPE = 2.0**500 +# Below this shape, the continued fraction's first elements, a times numbers of order +# one, could be subnormal; it is evaluated divided by a. +_TINY_SHAPE = 1e-280 +# From this ratio of two shapes, the direct form of a deviance term subtracts +# logarithms that can be 40 times its own. +_FAR_SHAPES = 1e17 +# Below this shape, and where the exponent s of the small-shape series is below the +# second number in absolute value, the series takes s / a instead of s, whose terms +# the Stirling corrections it adds and subtracts would absorb. +_SMALL_SERIES_SHAPE = 1e-3 +_ABSORBED_EXPONENT = 1e-4 + +# Near its center, where x is within about one standard deviation of a / (a + b), the +# incomplete beta function of two shapes from 1e10 comes from Temme's uniform expansion +# instead of the continued fraction, whose number of steps there grows like the cube +# root of the shapes (about 11,600 at 1e10, found by counting them). The band is +# D = bd0(a, (a + b) x) + bd0(b, (a + b)(1 - x)) <= 1/2. +_UNIFORM_BETA_SHAPE = 1e10 +_UNIFORM_BETA_BAND = 0.5 + # log Gamma(1 + a) = -euler a + sum_{k>=2} c_k a^k with c_k = (-1)^k zeta(k) / k # (DLMF 5.7.3), k = 2..30; used for |a| <= 1/4. _LGAMMA1P_TAYLOR = np.array( @@ -126,6 +161,23 @@ ] ) +# The coefficients B_2k / (2k (2k - 1)) of Stirling's series for delta(z), k = 1..10 +# (DLMF 5.11.1), which stirling_error sums for z >= 10. +_STIRLING_SERIES = np.array( + [ + 1.0 / 12.0, + -1.0 / 360.0, + 1.0 / 1260.0, + -1.0 / 1680.0, + 1.0 / 1188.0, + -691.0 / 360360.0, + 1.0 / 156.0, + -3617.0 / 122400.0, + 43867.0 / 244188.0, + -174611.0 / 125400.0, + ] +) + # Temme's uniform expansion of the incomplete gamma function (DLMF 8.12.3, 8.12.8): # Q(a, x) = erfc(eta sqrt(a/2)) / 2 + exp(-a eta^2 / 2) / sqrt(2 pi a) sum_k c_k(eta) a^-k, # with eta^2 / 2 = lambda - 1 - log(lambda), lambda = x / a and eta of the sign of x - a. @@ -369,11 +421,16 @@ def deviance_term(x: float, m: float, d: float) -> float: """Return Loader's bd0(x, m) = x log(x / m) + m - x, given d = x - m. The caller forms d accurately; near x = m the series in v = d / (x + m) avoids the - cancellation of the direct form (Loader 2000). + cancellation of the direct form (Loader 2000). bd0 is homogeneous of degree one: + where x + m, or 2 x, would overflow, x, m and d are halved and the result doubled. """ if x == 0.0: return m + scale = 1.0 s = x + m + if not s < _HALF_RANGE: + x, m, d, scale = 0.5 * x, 0.5 * m, 0.5 * d, 2.0 + s = x + m if abs(d) < 0.5 * s: v = d / s v2 = v * v @@ -384,10 +441,10 @@ def deviance_term(x: float, m: float, d: float) -> float: term *= v2 updated = total + term / (2.0 * j + 1.0) if updated == total: - return total + return scale * total total = updated j += 1.0 - return x * math.log(x / m) - d + return scale * (x * math.log(x / m) - d) @jit_numpy_errors @@ -624,22 +681,21 @@ def log_beta_prefactor( deviance terms bd0(a, (a + b) x) and bd0(b, (a + b) y) replace the logarithms, whose differences a - (a + b) x = lam and b - (a + b) y = -lam are known accurately. With one large parameter L and a small one, the term c log L of the gamma ratio is merged with - c log v of the small parameter's power into c log(L v). + c log v of the small parameter's power into c log(L v). Where a + b overflows, a b / + (a + b) is formed as b / (1 + b / a) and the deviance terms at half scale. """ if a >= 10.0 and b >= 10.0: n = a + b - m = n * x - if abs(lam) < 0.5 * (a + m): - dev_a = deviance_term(a, m, lam) + if n < INF: + dev_a = _shape_deviance(a, b, n * x, lam, log_x) + dev_b = _shape_deviance(b, a, n * y, -lam, log_y) + half_log = 0.5 * math.log(a / n * b) else: - dev_a = a * (-math.log1p(b / a) - log_x) - lam - m = n * y - if abs(lam) < 0.5 * (b + m): - dev_b = deviance_term(b, m, -lam) - else: - dev_b = b * (-math.log1p(a / b) - log_y) + lam + dev_a = _halved_shape_deviance(a, b, lam, log_x) + dev_b = _halved_shape_deviance(b, a, -lam, log_y) + half_log = 0.5 * math.log(b / (1.0 + b / a)) return ( - 0.5 * math.log(a / n * b) + half_log - _LN_SQRT_2PI - dev_a - dev_b @@ -668,6 +724,33 @@ def log_beta_prefactor( return a * log_x + b * log_y - log_beta_positive(a, b) +@jit_numpy_errors +def _shape_deviance(a: float, b: float, m: float, lam: float, log_x: float) -> float: + """Return bd0(a, m) for m = (a + b) x, lam = a - m and log_x = log x, a and b + from 10, as log_beta_prefactor takes it. + + Near m = a, Loader's series. Further, a log(a / m) - lam, with a log(a / m) as + a (-log1p(b / a) - log x), which needs no m; but where b exceeds 1e17 a and x is + tiny the two logarithms cancel, and m, a product rounded once, gives the quotient + better. + """ + if abs(lam) < 0.5 * (a + m): + return deviance_term(a, m, lam) + if b > _FAR_SHAPES * a and m >= _SMALLEST_NORMAL and a / m < INF: + return a * math.log(a / m) - lam + return a * (-math.log1p(b / a) - log_x) - lam + + +@jit_numpy_errors +def _halved_shape_deviance(a: float, b: float, lam: float, log_x: float) -> float: + """Return bd0(a, (a + b) x) for a + b beyond the double range, from bd0(a / 2, + (a - lam) / 2) = bd0(a, a - lam) / 2, which stays in range.""" + half = 0.5 * a - 0.5 * lam # (a + b) x / 2 + if abs(lam) < 0.5 * a + half: + return 2.0 * deviance_term(0.5 * a, half, 0.5 * lam) + return a * (-math.log1p(b / a) - log_x) - lam + + @jit_numpy_errors def _beta_continued_fraction( a: float, b: float, x: float, y: float, lam: float @@ -718,6 +801,126 @@ def _beta_continued_fraction( return NAN +@jit_numpy_errors +def _log_beta_continued_fraction_in_range( + a: float, b: float, x: float, y: float, lam: float +) -> float: + """Return log K, for the K of :func:`_beta_continued_fraction`, for shapes from 2^500, + or for a below 1e-280. + + The same fraction, its products ordered so that none leaves the float64 range: + a (1 + lam) / (a + 1) as (1 + lam) (a / (a + 1)); (a + b + k - 1) x / (s - 1) as + x (a + k - 1) / (s - 1) + x b / (s - 1); lam (a - 2k / (s - 1)) / (s + 1) with the + quotient first. The modified Lentz method keeps its second sequence as the + denominators e_k = den + num / e_(k-1) instead of their reciprocals, which would be + subnormal, and so inexact, for denominators beyond 2^1022; the value is scaled by + 2^-600 when it is that large. For a tiny a, whose products with the first two + elements, both proportional to a, would be subnormal, it evaluates K / a: their + quotients by a and the later elements as they are. A step whose terms are not finite + ends it with NaN, instead of the remaining steps running on NaN. + """ + tiny = a < _TINY_SHAPE + if tiny: + f = (1.0 + lam) / (a + 1.0) + scale = math.log(a) + else: + f = (1.0 + lam) * (a / (a + 1.0)) + scale = 0.0 + c = f + if f > _HUGE_SHAPE: + f *= _SCALE_DOWN + scale += 600.0 * _LN2 + e = INF + k = 1.0 + while k < _MAX_ITERATIONS: + s = a + 2.0 * k + ratio = (a + k - 1.0) / (s - 1.0) + first = 1.0 / (s - 1.0) if tiny and k == 1.0 else ratio + num = first * (x * ratio + x * (b / (s - 1.0))) * (k * ((b - k) * x)) + den = ( + 2.0 * k * (1.0 + y) * ((a + k) / (s - 1.0)) * (s / (s + 1.0)) + + ((a - 1.0) / (s - 1.0)) * (s / (s + 1.0)) + + lam * ((a - 2.0 * k / (s - 1.0)) / (s + 1.0)) + ) + if not (abs(num) < INF and abs(den) < INF): + return NAN + e = den + num / e + if e == 0.0: + e = _TINY + c = den + num / c + if c == 0.0: + c = _TINY + delta = c / e + f *= delta + if abs(delta - 1.0) < _CF_TOLERANCE: + return math.log(f) + scale + k += 1.0 + return NAN + + +@jit_numpy_errors +def _log_beta_fraction(a: float, b: float, x: float, y: float, lam: float) -> float: + """Return log K, for the K of :func:`_beta_continued_fraction`, for any shapes.""" + if _TINY_SHAPE <= a < _HUGE_SHAPE and b < _HUGE_SHAPE: + return math.log(_beta_continued_fraction(a, b, x, y, lam)) + return _log_beta_continued_fraction_in_range(a, b, x, y, lam) + + +@jit_numpy_errors +def _log_betainc_uniform( + a: float, b: float, lam: float, deviance: float +) -> tuple[float, float]: + """Return log I_x(a, b) and log(1 - I_x(a, b)) near the center, for a, b >= 1e10. + + Temme's (1982) uniform expansion, as the incomplete gamma function's above. With + r = a + b, x0 = a / r and eta of the sign of x - x0 (so of -lam, lam = r (x0 - x)) + defined by r eta^2 / 2 = bd0(a, r x) + bd0(b, r (1 - x)) = deviance, + + I_x(a, b) = erfc(-eta sqrt(r / 2)) / 2 - G exp(-deviance) c0(eta) / sqrt(2 pi r) + + up to terms of relative order min(a, b)^(-3/2), below 1e-15 from 1e10; G is + Gamma*(r) / (Gamma*(a) Gamma*(b)), the ratio of the Stirling series' corrections, + and 1 - I has the opposite signs. The density in eta is + sqrt(r / (2 pi)) G exp(-r eta^2 / 2) f(eta) with f(eta) = sqrt(x0 (1 - x0)) eta / + (t - x0) at the t of eta, and c0(eta) = (f(eta) - 1) / eta. Its Taylor series, + derived by reverting the series of eta in t - x0 (and checked against mpmath), is + + c0(eta) = (p - q) / (3 s) + (1 - s^2) / (12 s^2) eta + + (p - q)(2 + s^2) / (135 s^3) eta^2 + O(eta^3), + + with p = x0, q = 1 - x0 and s^2 = p q. In the band, r eta^2 / 2 <= 1/2, so that the + omitted terms are below 1e-15 of c0, itself of order min(a, b)^(-1/2) of I. Every + quantity is formed from p, q and a q = a b / r, so that r itself, which may + overflow, is never needed. + """ + if a <= b: + ratio = a / b + p = ratio / (1.0 + ratio) + q = 1.0 / (1.0 + ratio) + else: + ratio = b / a + p = 1.0 / (1.0 + ratio) + q = ratio / (1.0 + ratio) + pq = p * q + pqr = a * q # a b / r + difference = (0.5 * a - 0.5 * b) / (0.5 * a + 0.5 * b) # p - q + root = math.sqrt(deviance) # |eta| sqrt(r / 2) + side = 1.0 if lam > 0.0 else -1.0 # x below the center: I is the smaller tail + # c0(eta) / sqrt(r), with eta / sqrt(r) = -side sqrt(2 deviance) / r + c0 = ( + difference / (3.0 * math.sqrt(pqr)) + - side * math.sqrt(2.0 * deviance) * (1.0 - pq) / (12.0 * pqr) + + difference * (2.0 + pq) * 2.0 * deviance / (135.0 * pqr * math.sqrt(pqr)) + ) + gamma_ratio = math.exp( + stirling_error(a + b) - stirling_error(a) - stirling_error(b) + ) + correction = gamma_ratio * math.exp(-deviance) * c0 * _INV_SQRT_2PI + lower = 0.5 * math.erfc(side * root) - correction + upper = 0.5 * math.erfc(-side * root) + correction + return math.log(lower), math.log(upper) + + @jit_numpy_errors def _beta_upper_small_shape(a: float, b: float, x: float, log_x: float) -> float: """Return 1 - I_x(a, b) for 0 < a < 1 and b x < 2, without forming 1 - I. @@ -728,18 +931,31 @@ def _beta_upper_small_shape(a: float, b: float, x: float, log_x: float) -> float point of the continued fraction, where b x < a + 1, so 1 - I would lose accuracy; here the two parts cancel by at most a factor of about ten. """ + s = _small_shape_exponent_sum(a, b, x, log_x) + return -math.expm1(s) - math.exp(s) * a * _small_shape_series(a, b, x) + + +@jit_numpy_errors +def _small_shape_exponent_sum(a: float, b: float, x: float, log_x: float) -> float: + """Return the s = a log x - log(a B(a, b)) of :func:`_beta_upper_small_shape`, as a + sum of its terms.""" if b >= 10.0: # log Gamma(a + b) - log Gamma(b) = a log b + ..., merged with a log x into # a log(b x), as in the prefactor, so that no large logarithms cancel - s = ( + return ( a * _log_product(b, x, log_x) + ((a + b - 0.5) * math.log1p(a / b) - a) + stirling_error(a + b) - stirling_error(b) - _lgamma1p(a) ) - else: - s = a * log_x - (_lgamma1p(a) - log_gamma_increment(b, a)) + return a * log_x - (_lgamma1p(a) - log_gamma_increment(b, a)) + + +@jit_numpy_errors +def _small_shape_series(a: float, b: float, x: float) -> float: + """Return sum_{n>=1} (1 - b)_n x^n / (n! (a + n)), the series of + :func:`_beta_upper_small_shape`.""" total = 0.0 term = 1.0 # (1 - b)_n x^n / n! n = 1.0 @@ -750,7 +966,94 @@ def _beta_upper_small_shape(a: float, b: float, x: float, log_x: float) -> float if abs(contribution) <= 0.25 * EPS * abs(total): break n += 1.0 - return -math.expm1(s) - math.exp(s) * a * total + return total + + +@jit_numpy_errors +def _log_beta_upper_small_shape(a: float, b: float, x: float, log_x: float) -> float: + """Return log(1 - I_x(a, b)) for 0 < a < 1 and b x < 2. + + The logarithm of :func:`_beta_upper_small_shape`, unless a is below 1e-3 and s + below 1e-4 in absolute value. The sum that forms s adds and subtracts Stirling + corrections of up to 0.0084, whose rounding leaves an error of about 1e-18 in s: + a relative error of 1e-14 or more there, all of s at shapes below about 1e-20, + and a subnormal a rounds its terms too. There the same formula is taken as + 1 - I = a [-sigma expm1(a sigma) / (a sigma) - exp(a sigma) S], with sigma = s / a + and S the series; sigma is log x - log Gamma(1 + a) / a + + (log Gamma(b + a) - log Gamma(b)) / a, with every quotient by a taken from its + series. + """ + s = _small_shape_exponent_sum(a, b, x, log_x) + sigma = NAN + if a < _SMALL_SERIES_SHAPE and abs(s) < _ABSORBED_EXPONENT: + sigma = _small_shape_exponent(a, b, x, log_x) + if not abs(sigma) < INF: # not needed, or beyond the range for subnormal shapes + return math.log(-math.expm1(s) - math.exp(s) * a * _small_shape_series(a, b, x)) + t = a * sigma + ratio = math.expm1(t) / t if t != 0.0 else 1.0 + bracket = -sigma * ratio - math.exp(t) * _small_shape_series(a, b, x) + return math.log(a) + math.log(bracket) + + +@jit_numpy_errors +def _small_shape_exponent(a: float, b: float, x: float, log_x: float) -> float: + """Return s / a = log x - log(a B(a, b)) / a for 0 < a < 1e-3 (any b > 0). + + log(a B(a, b)) = log Gamma(1 + a) - (log Gamma(b + a) - log Gamma(b)). The first + quotient comes from the Taylor series of log Gamma(1 + a); the second, for b from + 10, from log Gamma(b + a) - log Gamma(b) = a log b + (b + a - 1/2) log1p(a / b) - a + + delta(b + a) - delta(b), its log b merged with log x into log(b x) as in the + prefactor; for b below 10, b is shifted up to 10 or more first, as + :func:`log_gamma_increment` shifts it. + """ + c = _LGAMMA1P_TAYLOR + poly: float = c[c.size - 1] + for k in range(c.size - 2, -1, -1): + poly = poly * a + c[k] + lgamma1p_over = a * poly - _EULER # log Gamma(1 + a) / a + if b >= 10.0: + return _log_product(b, x, log_x) + _gamma_ratio_slope(b, a) - lgamma1p_over + total = 0.0 # sum_j log1p(a / (b + j)) / a over the shifts + z = b + while z < 10.0: + total += _log1p_ratio(a / z) / z + z += 1.0 + return log_x + ((math.log(z) + _gamma_ratio_slope(z, a)) - total) - lgamma1p_over + + +@jit_numpy_errors +def _gamma_ratio_slope(z: float, a: float) -> float: + """Return (log Gamma(z + a) - log Gamma(z)) / a - log z for z >= 10 and 0 < a < 1. + + That is (z + a - 1/2) log1p(u) / (u z) - 1 + (delta(z + a) - delta(z)) / a with + u = a / z, each quotient from a form that stays accurate as u goes to 0, where + their limits are -1/(2 z) and delta'(z). + """ + u = a / z + g = _log1p_ratio(u) + # (z + a - 1/2) g / z - 1 = (g - 1) + g (a - 1/2) / z + g_minus_one = -u * (0.5 - u / 3.0) if u < 1e-8 else g - 1.0 + # (delta(z + a) - delta(z)) / a = sum_k c_k z^-2k [(1 + u)^-(2k-1) - 1] / u + slope = 0.0 + inverse_square = 1.0 / (z * z) + power = inverse_square + for k in range(_STIRLING_SERIES.size): + m = 2.0 * k + 1.0 + if u < 1e-8: + h = -m * (1.0 - 0.5 * (m + 1.0) * u) + else: + h = math.expm1(-m * math.log1p(u)) / u + slope += _STIRLING_SERIES[k] * power * h + power *= inverse_square + return g_minus_one + g * (a - 0.5) / z + slope + + +@jit_numpy_errors +def _log1p_ratio(u: float) -> float: + """Return log1p(u) / u for u >= 0, 1 at u = 0.""" + if u < 1e-8: + return 1.0 - u * (0.5 - u / 3.0) + return math.log1p(u) / u @jit_numpy_errors @@ -767,25 +1070,30 @@ def log_betainc_pair( return -INF, 0.0 if log_y == -INF: return 0.0, -INF + if a >= _UNIFORM_BETA_SHAPE and b >= _UNIFORM_BETA_SHAPE: + # (a + b) x = a - lam and (a + b) y = b + lam + deviance = deviance_term(a, a - lam, lam) + deviance_term(b, b + lam, -lam) + if deviance <= _UNIFORM_BETA_BAND: + return _log_betainc_uniform(a, b, lam, deviance) # The continued fraction gives the tail on its own side of the branch point; the # logarithm of the other tail is taken through it, log(1 - s) = log1p(-s), unless the # computed tail exceeds 1/2 with a shape below 1, where the other tail comes from the # small-shape series instead. if lam > x - y: # x < (a + 1) / (a + b + 2) - log_i = log_beta_prefactor(a, b, x, y, log_x, log_y, lam) - math.log( - _beta_continued_fraction(a, b, x, y, lam) + log_i = log_beta_prefactor(a, b, x, y, log_x, log_y, lam) - _log_beta_fraction( + a, b, x, y, lam ) log_i = min(log_i, 0.0) if a < 1.0 and log_i > -_LN2: - log_c = math.log(_beta_upper_small_shape(a, b, x, log_x)) + log_c = _log_beta_upper_small_shape(a, b, x, log_x) return log1mexp(-log_c), log_c return log_i, log1mexp(-log_i) - log_c = log_beta_prefactor(b, a, y, x, log_y, log_x, -lam) - math.log( - _beta_continued_fraction(b, a, y, x, -lam) + log_c = log_beta_prefactor(b, a, y, x, log_y, log_x, -lam) - _log_beta_fraction( + b, a, y, x, -lam ) log_c = min(log_c, 0.0) if b < 1.0 and log_c > -_LN2: - log_i = math.log(_beta_upper_small_shape(b, a, y, log_y)) + log_i = _log_beta_upper_small_shape(b, a, y, log_y) return log_i, log1mexp(-log_i) return log1mexp(-log_c), log_c @@ -802,12 +1110,27 @@ def log_betainc_pair_at(a: float, b: float, x: float) -> tuple[float, float]: y_low = 0.0 # 1 - x is exact log_x = math.log1p(-y) log_y = math.log(y) - p, p_low = two_product(a, y) - q, q_low = two_product(b, x) + if a < _SPLIT_LIMIT and b < _SPLIT_LIMIT: + p, p_low = two_product(a, y) + q, q_low = two_product(b, x) + else: + p, p_low = _scaled_two_product(a, y) + q, q_low = _scaled_two_product(b, x) lam = (p - q) + ((p_low - q_low) + a * y_low) return log_betainc_pair(a, b, x, y, log_x, log_y, lam) +@jit_numpy_errors +def _scaled_two_product(c: float, v: float) -> tuple[float, float]: + """Return two_product(c, v) for finite c >= 0 and 0 <= v <= 1, c up to the largest + double: from 1e290, c is scaled by 2^-600 before Dekker's split, which would + overflow, and both parts are scaled back, all exactly.""" + if c < _SPLIT_LIMIT: + return two_product(c, v) + p, e = two_product(c * _SCALE_DOWN, v) + return p * _SCALE_UP, e * _SCALE_UP + + def elementwise( loop: Callable[..., Any], kind: int, **arguments: ArrayLike ) -> float | NDArray[np.float64]: diff --git a/cytepool/stats/special.py b/cytepool/stats/special.py index 9a4ae6a..02a12bc 100644 --- a/cytepool/stats/special.py +++ b/cytepool/stats/special.py @@ -43,7 +43,13 @@ t and F distributions with many degrees of freedom). For a shape below 1 the tail given by the fraction can be close to 1; the other tail then comes from a term-by-term series of its own rather than from 1 minus it (the analogue of the - incomplete gamma case). + incomplete gamma case). For two shapes from 1e10 within about one standard + deviation of the center, where the fraction's steps grow like the cube root of + the shapes, Temme's (1982) uniform expansion replaces it, as for the incomplete + gamma function. Shapes from 2^500, and shapes below 1e-280, take forms of the + fraction and its prefactor whose products and quotients stay in the double range; + a shape below 1e-3 whose series exponent would be absorbed by the Stirling + corrections takes it divided by the shape. Limits: Accuracy is checked against 50-digit mpmath values (``validation/distributions.py``) @@ -51,11 +57,13 @@ the t and F distributions with degrees of freedom from 0.5 to 1e12 and infinity: relative errors stay below 1e-12 down to probabilities of 1e-300. In deep tails they are of order eps |log p| (about 2e-13 at p = 1e-300), because a probability is - evaluated through its logarithm. Larger shapes keep working: the incomplete beta - fraction needs a few thousand steps at shapes of 1e12, and the uniform expansion - serves the incomplete gamma function near the center for every shape from 1e3. A - series or continued fraction that has not converged after 1e7 steps gives NaN rather - than an inaccurate value. + evaluated through its logarithm. Larger shapes keep working: the uniform expansions + serve the incomplete gamma function near the center for every shape from 1e3, and + the incomplete beta function for two shapes from 1e10. The incomplete beta function + is checked, too, at every combination of shapes and arguments from the smallest + subnormal to the largest double (relative errors below 1e-12 there as well), and + takes under a millisecond for any of them. A series or continued fraction that has + not converged after 1e7 steps gives NaN rather than an inaccurate value. References: Olver FWJ, Olde Daalhuis AB, Lozier DW, Schneider BI, Boisvert RF, Clark CW, @@ -82,6 +90,10 @@ Temme NM (1979). The asymptotic expansion of the incomplete gamma functions. SIAM Journal on Mathematical Analysis 10, 757-766. doi:10.1137/0510071 + Temme NM (1982). The uniform asymptotic expansion of a class of integrals related to + cumulative distribution functions. SIAM Journal on Mathematical Analysis 13, + 239-253. doi:10.1137/0513017 + Smyth GK (2004). Linear models and empirical Bayes methods for assessing differential expression in microarray experiments. Statistical Applications in Genetics and Molecular Biology 3, article 3. doi:10.2202/1544-6115.1027 diff --git a/docs/coverage_exemptions.toml b/docs/coverage_exemptions.toml index 0d3a8fc..35d213b 100644 --- a/docs/coverage_exemptions.toml +++ b/docs/coverage_exemptions.toml @@ -66,4 +66,6 @@ "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_pca" = "guides/normalization_and_quality and guides/bulk_rna_seq: samples on their principal components" +"cytepool.plotting.plot_dotplot" = "guides/pool_cells: the markers of each cell type as a dot plot, from the cube" "cytepool.plotting.plot_quasi_dispersions" = "guides/bulk_rna_seq: edgeR's quasi-dispersions" diff --git a/docs/source/reference/api/index.md b/docs/source/reference/api/index.md index ea0f072..2e58e4b 100644 --- a/docs/source/reference/api/index.md +++ b/docs/source/reference/api/index.md @@ -1,7 +1,7 @@ --- myst: html_meta: - description: The CytePool Python API for pseudobulk cubes, their conversion to and from AnnData and MuData, designs, normalization, linear models, mixed models, voom, dream, sample quality weights, negative binomial models, quasi-likelihood F-tests, DESeq2's dispersions and tests, empirical Bayes moderation, gene-set tests (including rotation tests), pathway scores, networks of regulators and pathways and their activities, multiple testing, differential expression, cell-type composition, meta-analysis, combined p-values and replication across studies, SumRank and held-out prediction, and distributions. + description: The CytePool Python API for pseudobulk cubes, their conversion to and from AnnData and MuData, designs, normalization, linear models, mixed models, voom, dream, sample quality weights, negative binomial models, quasi-likelihood F-tests, DESeq2's dispersions and tests, empirical Bayes moderation, gene-set tests (including rotation tests), pathway scores, networks of regulators and pathways and their activities, multiple testing, differential expression, principal components of samples, cell-type composition, meta-analysis, combined p-values and replication across studies, SumRank and held-out prediction, and distributions. --- # CytePool API reference @@ -46,6 +46,7 @@ multiple_testing combined_p_values differential_expression variance_partitioning +principal_components composition meta_analysis replication diff --git a/docs/source/reference/api/principal_components.md b/docs/source/reference/api/principal_components.md new file mode 100644 index 0000000..3a6c56a --- /dev/null +++ b/docs/source/reference/api/principal_components.md @@ -0,0 +1,20 @@ +--- +myst: + html_meta: + description: Principal components of samples in CytePool, as DESeq2's plotPCA computes them from the genes of largest variance, for a matrix of log-expression values or a cell type of a pseudobulk cube, with their coordinates, shares of variance and loadings, and their figure. +--- + +# Principal components + +{func}`~cytepool.principal_components.principal_components` places samples on the +principal components of their most variable genes, as DESeq2's `plotPCA` does: the +`top` genes of largest variance over the samples, each centred, decomposed by singular +values. It takes a samples x genes matrix of log-expression values or a cell type of a +cube, whose samples are prepared as edgeR prepares counts for `plotMDS`. The result, a +{class}`~cytepool.principal_components.PrincipalComponents`, holds the samples' +coordinates, each component's share of the variance and the loadings of the genes, and +{func}`~cytepool.plotting.plot_pca` draws it. + +```{eval-rst} +.. automodule:: cytepool.principal_components +``` diff --git a/docs/source/reference/bibliography.md b/docs/source/reference/bibliography.md index d94a018..a48ae2d 100644 --- a/docs/source/reference/bibliography.md +++ b/docs/source/reference/bibliography.md @@ -78,6 +78,7 @@ SumRank and the prediction of held-out studies [@nakatsuka2025; @irwin1927; @hal [@cleveland1979; @cleveland1981; @deboor1978; @hastie1992; @scott1979; @freedman1981; @abramowitz1964; @johnson1995; @olver2026; @wichura1988; @paulson1942; @temme1979; +@temme1982; @lentz1976; @thompson1986; @loader2000; @maechler2012; @dekker1971; @knuth1997; @ogita2005; @neumaier1974; @kendall1945; @golub2013; @giraud2005; @berrut2004; @brent1973; @dowell1971; @graham1969; @briggs1993; @goldberg1991; @hoare1961; @@ -87,7 +88,8 @@ SumRank and the prediction of held-out studies [@nakatsuka2025; @irwin1927; @hal Classical scaling for MDS plots, forest plots, violin plots and the bandwidth of their densities, and the check of the colour palette for colour-vision deficiency [@gower1966; -@lewis2001; @hintze1998; @silverman1986; @machado2009]. +@lewis2001; @hintze1998; @silverman1986; @machado2009]; principal components +[@pearson1901]. ## Data structures, file formats and standards diff --git a/docs/source/reference/capabilities.md b/docs/source/reference/capabilities.md index 64dd9b3..376de05 100644 --- a/docs/source/reference/capabilities.md +++ b/docs/source/reference/capabilities.md @@ -55,6 +55,8 @@ size; its sections give every check. | The proportion of true nulls, pi0, and pi1 (qvalue's `pi0est`, its smoother and bootstrap) | {api}`proportion_of_true_nulls` | qvalue 2.38.0 | exact | [`combined_p_values_and_replication.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/combined_p_values_and_replication.py) | | Replication across the studies of a collection: combined and partial conjunction p-values, sign agreement, heterogeneity, pi1 between studies, genes found in k studies, and its figures | {api}`differential_expression_replication`, {api}`plot_study_agreement`, {api}`plot_replicated_genes` | metap 1.1, qvalue 2.38.0, metafor 5.0.1 | exact; close for REML | [`combined_p_values_and_replication.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/combined_p_values_and_replication.py) | | SumRank (Nakatsuka et al. 2025): every gene's best ranks summed across studies, with Irwin-Hall p-values and permutations that re-run each study's analysis; prediction of cases in held-out studies by genes chosen from the others, by SumRank and by meta-analysis; their figures | {api}`sumrank`, {api}`differential_expression_sumrank`, {api}`held_out_prediction`, {api}`plot_sumrank`, {api}`plot_held_out_auc` | the authors' code (commit 46392f66), R 4.4.3, mpmath 1.3.0, scikit-learn 1.9.1 | exact; calibrated for p-values from random permutations | [`sumrank.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/sumrank.py) | +| Principal components of samples (DESeq2's `plotPCA`), of a matrix or a cell type of a cube, and their figure | {api}`principal_components`, {api}`plot_pca` | DESeq2 1.46.0, edgeR 4.4.0 | exact, up to the sign of each component | [`pca_and_dot_plots.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/pca_and_dot_plots.py) | +| Dot plots of the fraction of cells that express each gene and their mean expression, per cell type, from a cube (scanpy's `dotplot`) | {api}`plot_dotplot` | scanpy 1.12.4 | exact | [`pca_and_dot_plots.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/pca_and_dot_plots.py) | | Figures of limma, edgeR, DESeq2, speckle and metafor (MDS, voom's trend, plotSA, MD and volcano plots, barcode plots, dispersions, composition, forest plots) | {api}`cytepool.plotting` | limma 3.62.1, edgeR 4.4.0, DESeq2 1.46.0, speckle 1.6.0, metafor 5.0.1 | exact for what is drawn; close for plotSA's prior | [`plotting.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/plotting.py) | | Tails and quantiles of distributions, special functions | {api}`cytepool.stats.distributions`, {api}`cytepool.stats.special` | mpmath 1.3.0 | exact (1e-15 to 1e-12) | [`distributions.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/distributions.py) | | LOWESS, natural spline bases | {api}`lowess`, {api}`natural_spline_basis` | R 4.4.3 | close for LOWESS; exact for spline bases | [`smoothing.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/smoothing.py) | diff --git a/docs/source/references.bib b/docs/source/references.bib index 513c06c..35e8117 100644 --- a/docs/source/references.bib +++ b/docs/source/references.bib @@ -1467,6 +1467,17 @@ @article{paulson1942 doi = {10.1214/aoms/1177731608}, } +@article{pearson1901, + author = {Pearson, Karl}, + title = {On lines and planes of closest fit to systems of points in space}, + journal = {The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science}, + year = {1901}, + volume = {2}, + number = {11}, + pages = {559--572}, + doi = {10.1080/14786440109462720}, +} + @article{pearson1933, author = {Pearson, K.}, title = {On a method of determining whether a sample of size n supposed to have been drawn from a parent population having a known probability integral has probably been drawn at random}, @@ -1856,6 +1867,17 @@ @article{temme1979 doi = {10.1137/0510071}, } +@article{temme1982, + author = {Temme, N. M.}, + title = {The uniform asymptotic expansion of a class of integrals related to cumulative distribution functions}, + journal = {SIAM Journal on Mathematical Analysis}, + year = {1982}, + volume = {13}, + number = {2}, + pages = {239--253}, + doi = {10.1137/0513017}, +} + @article{thompson1986, author = {Thompson, I. J. and Barnett, A. R.}, title = {Coulomb and Bessel functions of complex arguments and order}, diff --git a/docs/source/statistics/cubes.md b/docs/source/statistics/cubes.md index d943190..d66caa0 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, how merging cell types and pooling samples stay exact, and how the cubes of several studies are aligned. +description: What a CytePool cube holds, why every layer is a sum over cells, how merging cell types and pooling samples stay exact, the principal components of its samples and the dot plot of its cells, and how the cubes of several studies are aligned. --- # The cube and its additive layers @@ -56,6 +56,38 @@ features. Integer layers are summed exactly, and a sum that would exceed the lar are not sums, cannot be merged this way; that is why the cube keeps sums and derives them. +## Samples and cells at a glance + +Two figures show a cube's data before any test. {api}`principal_components` places the +samples of one cell type on the principal components of their most variable genes +[@pearson1901], as DESeq2's `plotPCA` does [@love2014]. The samples are prepared as +for an MDS plot: those with fewer than 10 cells are left out, genes are kept by +{api}`filter_by_expression`, and the values are log-CPM with TMM factors and an average +count of 2 added. The `top` genes of largest variance over the samples (500 by +default), each centred on its mean, form a samples x genes matrix $X = U S V^\top$; +the samples' coordinates are $U S$, each component's share of the variance is its +squared singular value over their sum, and the result, a {api}`PrincipalComponents`, +keeps the loadings $V$ of the genes. The sign of a component is arbitrary; the +coordinate of largest absolute value is made positive. `plot_pca` draws two +components, labelled with their shares as `plotPCA` labels them. + +The dot plot shows, for chosen genes and cell types, the fraction of the cells that +express each gene and their mean expression, as scanpy draws it [@wolf2018]. A cube +holds what it needs without the cells. Because every layer is a sum, the fraction and +the mean over the cells of all the samples pooled, + +$$ +f_{gf} = \frac{\sum_j \text{detected}_{gjf}}{\sum_j n_{jg}}, \qquad +m_{gf} = \frac{\sum_j \text{log\_sum}_{gjf}}{\sum_j n_{jg}}, +$$ + +are those computed from the cells themselves, and so is the mean over the cells that +express the gene, $\sum_j \text{log\_sum}_{gjf} / \sum_j \text{detected}_{gjf}$, +since $\log(1 + 0) = 0$ for the others. {api}`Cube.fraction_detected` and {api}`Cube.mean` +give the same quantities per sample. `plot_dotplot` draws them, pooled over every sample +or some, with scanpy's scaling of the means by gene or by cell type and its groups of +genes. + ## What a cube carries - **Tables.** Samples, groups and features each have a table ({api}`Table`) with a diff --git a/docs/source/validation.md b/docs/source/validation.md index a3661aa..f5209c2 100644 --- a/docs/source/validation.md +++ b/docs/source/validation.md @@ -18,8 +18,9 @@ This page summarises the tables the scripts printed when they were run against t 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`, `gsea.py`, `cube_collection.py`, `pathway_scores.py`, -`combined_p_values_and_replication.py`, `activity.py` and `sumrank.py`, and for -`dream.py`, which printed the same table again after changes to the code it checks; on +`combined_p_values_and_replication.py`, `activity.py`, `sumrank.py`, +`pca_and_dot_plots.py` and `distributions.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`, `cube_fingerprints.py` and `anndata_conversion.py`), for `linear_models.py`, @@ -73,9 +74,10 @@ covers and the largest difference among them. Run a script to see every check. | [`combined_p_values_and_replication.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/combined_p_values_and_replication.py) | Combined p-values and replication | scipy 1.18.1, mpmath 1.3.0, R 4.4.3, metap 1.1, qvalue 2.38.0, metafor 5.0.1 | 54 | 47 pass, 7 deliberate differences | 67.4 s | 444 MiB | | [`sumrank.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/sumrank.py) | SumRank and held-out prediction | R 4.4.3, the authors' code (commit 46392f66), mpmath 1.3.0, scipy 1.18.1, scikit-learn 1.9.1 | 23 | 13 pass, 9 information, 1 deliberate difference | 44.9 s | 399 MiB | | [`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 | +| [`pca_and_dot_plots.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/pca_and_dot_plots.py) | Principal components and dot plots | R 4.4.3, DESeq2 1.46.0, edgeR 4.4.0, ggplot2 4.0.3, scanpy 1.12.4 | 13 | 11 pass, 1 deliberate difference, 1 information | 19.8 s | 671 MiB | +| [`distributions.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/distributions.py) | Distributions and special functions | mpmath 1.3.0 | 51 | 51 pass | 12 min 28 s | 174 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,379** | **2,049 pass, none failed; 185 information, 143 deliberate differences, 2 controls** | | | +| **All** | | | **2,397** | **2,065 pass, none failed; 186 information, 144 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 @@ -2248,25 +2250,75 @@ difference: a gene without residual degrees of freedom, such as one without coun has no quasi-dispersion (its deviance and its degrees of freedom are 0); edgeR draws it at 0, and CytePool leaves it out. +## Principal components and dot plots + +[`validation/pca_and_dot_plots.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/pca_and_dot_plots.py) +compares `principal_components` and `plot_pca` with DESeq2's `plotPCA`, run as a black +box with `returnData = TRUE` and drawn, on the same values, and `plot_dotplot` with +scanpy's `dotplot` on the cells a cube was pooled from. The principal components (seed +20261103) are those of the log2 counts per million of negative binomial counts of 300 +genes in 12 samples of three conditions and two batches, with the 100 genes of largest +variance and with 500, more than there are; of a matrix whose genes 19 to 22 have equal +variance, with the 20 of largest variance; and of a cube of the counts, two of whose +samples have 5 cells, which R prepares as the edgeR User's Guide does (`filterByExpr`, +`normLibSizes`, `cpm(log = TRUE, prior.count = 2)`). The cells of the dot plots (seed +20261104) are 4 cell types in 6 samples, one sample without cells of one type, with +counts of 40 genes; scanpy normalizes them (`normalize_total`, `log1p`) and +`from_cells` pools them. Reference: R 4.4.3, DESeq2 1.46.0, edgeR 4.4.0, ggplot2 +4.0.3; scanpy 1.12.4, anndata 0.13.4, matplotlib 3.11.2. + +| Check family | Reference | Largest difference | Tolerance | Outcome | +| --- | --- | --- | --- | --- | +| Principal components: the coordinates up to sign and the shares of variance, with the 100 genes of largest variance (components 1 and 2, and 3 and 4) and with 500; genes of equal variance; a cube's cell type prepared by edgeR | `plotPCA(returnData = TRUE)` | 1.1e-14 | 1e-10 | 3 pass | +| The figure: the points up to sign, the axis labels, equal scales, the samples coloured by condition and batch | `plotPCA` | 6.9e-15 | 1e-10 | 1 pass | +| Two samples: the first component; DESeq2 returns a second of rounding noise (share 9.5e-31), CytePool none | `plotPCA(returnData = TRUE)` | 8.9e-16 | n/a | 1 deliberate difference | +| Components whose sign differs from R's | `plotPCA(returnData = TRUE)` | 2 of 10 | n/a | 1 information | +| Dot plots: the areas and the fractions they encode, and the values the colours show, plain, scaled by gene or by group and over the expressing cells only, for every sample and for 3 samples and 3 cell types in a chosen order, and for genes in groups; the names of the rows and columns, and the key of the areas | `scanpy.pl.dotplot` | 1.0e-15 | 1e-10 | 6 pass | +| The brackets of genes given in groups: their columns and names | `scanpy.pl.dotplot` | identical | 0 | 1 pass | + +Differences are scaled by max(\|reference\|, 1). R's `prcomp`, whose coordinates +`plotPCA` returns, leaves the sign of each component arbitrary; CytePool makes each +component's largest coordinate positive, so the coordinates are compared after each +component takes the sign that matches R's. What DESeq2's help page leaves open was +found by probing: the genes are ranked by their variance, genes of equal variance keep +their order (the first are kept), `ntop` beyond the number of genes takes every gene, +the shares are those of all the components (`percentVar`), the axes are labelled with +them rounded to a whole percent, and the figure has equal scales. The dot plots +compare the areas scanpy draws and the values its colours show (`dot_color_df`), and +its fractions (`dot_size_df`) with those that CytePool's areas encode. + +Deliberate difference: with two samples, DESeq2's `plotPCA` returns, and draws, a +second component of rounding noise; CytePool returns one component, since the centred +values of n samples have at most n - 1 dimensions, and refuses to draw a second. + ## Distributions and special functions [`validation/distributions.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/distributions.py) checks every public function of `cytepool.stats.special` and -`cytepool.stats.distributions` on 19,493 points against mpmath evaluated at the exact +`cytepool.stats.distributions` on 20,375 points against mpmath evaluated at the exact binary values of the inputs, with 50 significant digits (70 for intermediate steps); quantile references are roots of the 70-digit tail probabilities. The grids reach probabilities of 1e-300 (and log probabilities far below), shapes from 1e-3 to 1e6 (incomplete gamma shapes to 1e18 above the center), degrees of freedom from 0.5 to 1e12 and infinity, hypergeometric populations up to 1e13, batches of 2,000 moderated t -and quasi-likelihood F statistics, and cases where textbook methods lose accuracy. -The incomplete beta reference is DLMF 8.17.22 in 70-digit arithmetic, and the -hypergeometric reference sums exact ratios in fixed-point integer arithmetic. -Reference: mpmath 1.3.0, Python 3.13.12, numpy 2.5.4. +and quasi-likelihood F statistics, and cases where textbook methods lose accuracy. For +the incomplete beta function they also hold the edges of the double range: every +combination of the shapes 5e-324 (the smallest subnormal), 1e-300, 1e-8, 0.5, +1 - 2^-53, 1, 2, 30, 1e8, 1e300 and 1.8e308 (the largest double) with the arguments 0, +5e-324, 1e-300, 1e-8, 0.5, 1 - 2^-53 and 1, and 35 points within 1.5 standard +deviations of the center of shapes from 1e10 to 1e20. The incomplete beta reference +is DLMF 8.17.22 in 70-digit arithmetic, at the edges of the double range in 1,200 to +2,400 bits, and near the center of large shapes, where the fraction converges too +slowly, the beta density integrated numerically; the hypergeometric reference sums +exact ratios in fixed-point integer arithmetic. Reference: mpmath 1.3.0, Python +3.13.12, numpy 2.5.4. | Check family | Reference | Largest difference | Tolerance | Outcome | | --- | --- | --- | --- | --- | | Log-gamma (x > 0, x < 0), digamma (x > 0, x < 0), trigamma, inverse trigamma, log-beta | mpmath | 9.99e-16 (digamma, x < 0) | 1e-15, or 2e-15 for digamma at x < 0 and log-beta | 7 pass | | Regularized incomplete gamma and beta functions, their complements and logarithms | mpmath | 2.44e-13 | 1e-12 | 8 pass | +| The regularized incomplete beta function, its complement and their logarithms at the edges of the double range and near the center of large shapes | mpmath | 2.04e-13 | 1e-12 | 4 pass | +| The time per value of the incomplete beta function at every combination of 0, 5e-324, 1e-300, 1e-8, 0.5, 1 - 2^-53, 1, 2, 30, 1e8, 1e300 and 1.8e308 as shapes and argument (1,728, those outside the domain included), the best of three calls | timed | 2.1e-4 s (median 2.5e-5 s) | 1e-2 s | 1 pass | | Normal: cdf, sf, logcdf, logsf, ppf, isf | mpmath | 6.22e-16 | 1e-12 | 6 pass | | t: cdf, sf, logsf, ppf, isf, and sf and logsf of 2,000 moderated t-statistics | mpmath | 1.95e-13 | 1e-12 | 7 pass | | F: cdf, sf, logsf, ppf, isf, and sf and logsf of 2,000 quasi-likelihood F statistics | mpmath | 2.12e-13 | 1e-12 | 7 pass | @@ -2402,8 +2454,10 @@ genes 10 times: each permutation and study took 0.21 s with `limma_voom` and 0.1 resident memory of 567 MiB for the run. `adjust_p_values` with Hommel's method took 0.006 s on the 100,000 p-values of -`multiple_testing.py`. In `distributions.py`, mpmath took 395 s to compute the -references and CytePool 3.6 s to evaluate them, including compilation. +`multiple_testing.py`. In `distributions.py`, mpmath took 746 s to compute the +references, 128 s of them at the edges of the double range, and CytePool 1.8 s to +evaluate them, including compilation; the incomplete beta function took at most +2.1e-4 s for any value at those edges, with the time of a call from Python. `cube_reads.py` read every layer of each selection of its cube file into memory with `to_memory()`. The file holds 720 MB of values in chunks of one group, every sample @@ -2472,6 +2526,7 @@ RSCRIPT=$REF/r-meta/bin/Rscript METAFOR_RSCRIPT=$REF/r-metafor/bin/Rscript \ python validation/combined_p_values_and_replication.py RSCRIPT=$REF/r/bin/Rscript METAFOR_RSCRIPT=$REF/r-metafor/bin/Rscript \ SPECKLE_RSCRIPT=$REF/r-speckle/bin/Rscript python validation/plotting.py +SCANPY_PYTHON=$REF/python-scanpy/bin/python python validation/pca_and_dot_plots.py # variancePartition, in an R environment of its own with lme4 1.1-38 (below). RSCRIPT=$REF/r-dream/bin/Rscript python validation/dream.py RSCRIPT=$REF/r-dream/bin/Rscript python validation/kenward_roger_and_f_tests.py @@ -2607,5 +2662,6 @@ The reference packages run only as separate programs that produce outputs to com against. CytePool implements the methods from their papers and documentation, and conventions that the documentation leaves open were found by running the packages on chosen inputs; the source of a package under a copyleft licence is never read, and -what CytePool adapts from permissively licensed ones, fgsea and decoupler, -`THIRD_PARTY_NOTICES.md` records. +what CytePool adapts or follows from permissively licensed ones (fgsea, decoupler, +the authors' code of SumRank, and the conventions of scanpy's dot plot that its +documentation leaves open) `THIRD_PARTY_NOTICES.md` records. diff --git a/tests/data/distributions/betainc_extreme.npz b/tests/data/distributions/betainc_extreme.npz new file mode 100644 index 0000000..af5ddc3 Binary files /dev/null and b/tests/data/distributions/betainc_extreme.npz differ diff --git a/tests/data/pca_and_dot_plots/dot_plots.json b/tests/data/pca_and_dot_plots/dot_plots.json new file mode 100644 index 0000000..be5c0b6 --- /dev/null +++ b/tests/data/pca_and_dot_plots/dot_plots.json @@ -0,0 +1,183 @@ +{ + "plain": { + "rows": [ + "B", + "T", + "NK", + "pDC" + ], + "columns": [ + "g00", + "g01", + "g02", + "g03", + "g04", + "g10", + "g20", + "g30", + "g39" + ], + "key": [ + "20", + "40", + "60", + "80", + "100" + ], + "gene_groups": null + }, + "scaled by gene": { + "rows": [ + "B", + "T", + "NK", + "pDC" + ], + "columns": [ + "g00", + "g01", + "g02", + "g03", + "g04", + "g10", + "g20", + "g30", + "g39" + ], + "key": [ + "20", + "40", + "60", + "80", + "100" + ], + "gene_groups": null + }, + "scaled by group": { + "rows": [ + "B", + "T", + "NK", + "pDC" + ], + "columns": [ + "g00", + "g01", + "g02", + "g03", + "g04", + "g10", + "g20", + "g30", + "g39" + ], + "key": [ + "20", + "40", + "60", + "80", + "100" + ], + "gene_groups": null + }, + "expressing cells": { + "rows": [ + "B", + "T", + "NK", + "pDC" + ], + "columns": [ + "g00", + "g01", + "g02", + "g03", + "g04", + "g10", + "g20", + "g30", + "g39" + ], + "key": [ + "20", + "40", + "60", + "80", + "100" + ], + "gene_groups": null + }, + "chosen samples and groups": { + "rows": [ + "NK", + "B", + "T" + ], + "columns": [ + "g00", + "g01", + "g02", + "g03", + "g04", + "g10", + "g20", + "g30", + "g39" + ], + "key": [ + "20", + "40", + "60", + "80", + "100" + ], + "gene_groups": null + }, + "genes in groups": { + "rows": [ + "B", + "T", + "NK", + "pDC" + ], + "columns": [ + "g00", + "g01", + "g02", + "g03", + "g04", + "g10", + "g20", + "g30", + "g39" + ], + "key": [ + "20", + "40", + "60", + "80", + "100" + ], + "gene_groups": [ + [ + "B cells", + "T cells", + "others" + ], + [ + [ + 0, + 2 + ], + [ + 3, + 4 + ], + [ + 5, + 8 + ] + ] + ] + } +} \ No newline at end of file diff --git a/tests/data/pca_and_dot_plots/dot_plots.npz b/tests/data/pca_and_dot_plots/dot_plots.npz new file mode 100644 index 0000000..0a941bc Binary files /dev/null and b/tests/data/pca_and_dot_plots/dot_plots.npz differ diff --git a/tests/data/pca_and_dot_plots/pca.npz b/tests/data/pca_and_dot_plots/pca.npz new file mode 100644 index 0000000..495430e Binary files /dev/null and b/tests/data/pca_and_dot_plots/pca.npz differ diff --git a/tests/test_masked_inputs.py b/tests/test_masked_inputs.py index a17e508..f710d68 100644 --- a/tests/test_masked_inputs.py +++ b/tests/test_masked_inputs.py @@ -1962,6 +1962,18 @@ def _held_out_arguments(context: Context) -> dict[str, Any]: }, {"contrasts": None}, ), + # Principal components: of an array, or of a cube's group (an object). + "cytepool.principal_components.principal_components": Entry( + lambda c: {"data": log_cpm(COUNTS), "top": 4}, + {"data": None}, + options=("top",), + ), + "cytepool.principal_components.principal_components[cube]": Entry( + lambda c: {"data": _bulk_cube(), "group": "bulk", "min_cells": 1}, + {}, + options=("group", "min_cells"), + objects=("data",), + ), # Meta-analysis. NaN marks a study that did not measure a feature. "cytepool.meta_analysis.meta_analysis": Entry( lambda c: { @@ -2094,6 +2106,53 @@ def _held_out_arguments(context: Context) -> dict[str, Any]: options=("group", "color_by", "min_cells"), objects=("data", "ax"), ), + "cytepool.plotting.plot_pca": Entry( + lambda c: { + "data": log_cpm(COUNTS), + "color_by": ["a", "b", "a", "b"], + "marker_by": ["x", "x", "y", "y"], + "labels": ["S1", "S2", "S3", "S4"], + "top": 4, + "components": (1, 2), + "ax": _axes(), + }, + { + "data": None, + "color_by": None, + "marker_by": None, + "labels": None, + "components": None, + }, + options=("top",), + objects=("ax",), + ), + "cytepool.plotting.plot_pca[cube]": Entry( + lambda c: { + "data": _bulk_cube(), + "group": "bulk", + "color_by": "condition", + "min_cells": 1, + "ax": _axes(), + }, + {}, + options=("group", "color_by", "min_cells"), + objects=("data", "ax"), + ), + # Genes, groups and samples named by id; the cube is an object. + "cytepool.plotting.plot_dotplot": Entry( + lambda c: { + "cube": c.cube, + "genes": ["feature_0", "feature_2"], + "groups": ["T", "B"], + "samples": ["S1", "S2"], + "standard_scale": "gene", + "mean_only_expressed": True, + "ax": _axes(), + }, + {"genes": None, "groups": None, "samples": None}, + options=("standard_scale", "mean_only_expressed"), + objects=("cube", "ax"), + ), "cytepool.plotting.plot_cell_numbers": Entry( lambda c: {"cube": c.cube, "min_cells": 1, "numbers": True, "ax": _axes()}, {}, @@ -2653,6 +2712,7 @@ def _held_out_arguments(context: Context) -> dict[str, Any]: "cytepool.stats.combined_p_values.CombinedPValues", "cytepool.stats.multiple_testing.TrueNullProportion", "cytepool.sumrank.SumRankFit", + "cytepool.principal_components.PrincipalComponents", "cytepool.stats.dream.ModeratedFixedEffectTests", "cytepool.stats.quasi_likelihood.QuasiLikelihoodFit", "cytepool.stats.quasi_likelihood.QuasiLikelihoodTest", diff --git a/tests/test_pca_and_dot_plots.py b/tests/test_pca_and_dot_plots.py new file mode 100644 index 0000000..6efebc1 --- /dev/null +++ b/tests/test_pca_and_dot_plots.py @@ -0,0 +1,536 @@ +"""Tests for principal components, their figure and dot plots. + +They compare with: + +- what DESeq2's ``plotPCA`` returned and drew, and what scanpy's ``dotplot`` computed + and drew from the cells, saved by ``validation/pca_and_dot_plots.py`` in + ``tests/data/pca_and_dot_plots/``; +- values derived by hand: three samples on two genes whose components are the genes, + two samples, one gene, and a gene expressed in half the cells of two samples; +- invariants that a wrong implementation would break: the order of the samples and of + the genes, the sign and the shift of the values, a group merged or samples pooled + in the cube against the same cells given one label, and scaled means that run from + 0 to 1. + +Figures are drawn on Figures that pyplot does not manage, and read back from the +artists that carry their gid. +""" + +import importlib +import json +from functools import cache +from pathlib import Path +from typing import Any + +import numpy as np +import pytest +from numpy.testing import assert_allclose, assert_array_equal + +from cytepool import plotting +from cytepool.cube import Cube, from_cells +from cytepool.principal_components import PrincipalComponents, principal_components + +FIGURE = pytest.importorskip("matplotlib.figure") + +DATA = Path(__file__).parent / "data" / "pca_and_dot_plots" +EXACT = 1e-10 +SAMPLES = [f"S{i}" for i in range(1, 13)] +CONDITIONS = ["a"] * 4 + ["b"] * 4 + ["c"] * 4 +BATCHES = ["x", "y"] * 6 +GROUPS = ["B", "T", "NK", "pDC"] +CELL_SAMPLES = [f"D{i}" for i in range(1, 7)] +MARKERS = {"B cells": ["g00", "g01", "g02"], "T cells": ["g03", "g04"]} +MARKERS["others"] = ["g10", "g20", "g30", "g39"] +GENES = [gene for genes in MARKERS.values() for gene in genes] +DOT_CASES: dict[str, dict[str, Any]] = { + "plain": {}, + "scaled by gene": {"standard_scale": "gene"}, + "scaled by group": {"standard_scale": "group"}, + "expressing cells": {"mean_only_expressed": True}, + "chosen samples and groups": { + "samples": ["D2", "D3", "D5"], + "groups": ["NK", "B", "T"], + }, + "genes in groups": {"genes": MARKERS}, +} + + +@cache +def fixture(name: str) -> dict[str, np.ndarray]: + with np.load(DATA / f"{name}.npz", allow_pickle=False) as data: + return {key: data[key] for key in data.files} + + +@cache +def dot_texts() -> dict[str, Any]: + result: dict[str, Any] = json.loads((DATA / "dot_plots.json").read_text()) + return result + + +def axes() -> Any: + return FIGURE.Figure().subplots() + + +def artists(ax: Any, gid: str) -> list[Any]: + return [artist for artist in ax.get_children() if artist.get_gid() == gid] + + +def scaled(ours: Any, reference: Any) -> float: + """Largest difference over max(|reference|, 1), as the validation script has it.""" + ours, reference = np.asarray(ours, dtype=float), np.asarray(reference, dtype=float) + assert ours.shape == reference.shape + return float(np.max(np.abs(ours - reference) / np.maximum(np.abs(reference), 1.0))) + + +def up_to_sign(ours: np.ndarray, theirs: np.ndarray) -> float: + signs = np.sign(np.sum(ours * theirs, axis=0)) + return scaled(ours * signs, theirs) + + +def pca_cube() -> Cube: + data = fixture("pca") + return Cube( + {"counts": data["counts"][np.newaxis].astype(np.uint32)}, + samples={"sample_id": SAMPLES, "condition": CONDITIONS, "batch": BATCHES}, + groups={"group_id": ["T"]}, + features={"feature_id": [f"g{j:03d}" for j in range(data["counts"].shape[1])]}, + n_cells=data["n_cells"][:, np.newaxis], + ) + + +@cache +def dot_cube() -> Cube: + data = fixture("dot_plots") + return Cube( + { + "counts": data["detected"], # counts are not drawn; any valid layer will do + "detected": data["detected"], + "log_sum": data["log_sum"], + }, + samples={"sample_id": CELL_SAMPLES}, + groups={"group_id": GROUPS}, + features={"feature_id": [f"g{j:02d}" for j in range(40)]}, + n_cells=data["n_cells"], + normalization={"method": "library_size", "target_sum": 1e4, "pseudocount": 1.0}, + ) + + +# What DESeq2 and scanpy computed and drew. + + +@pytest.mark.parametrize( + ("case", "top", "components"), + [("top100", 100, (1, 2)), ("top100_pc34", 100, (3, 4)), ("top500", 500, (1, 2))], +) +def test_components_are_those_of_plotPCA_up_to_sign( + case: str, top: int, components: tuple[int, int] +) -> None: + data = fixture("pca") + result = principal_components(data["values"], top=top) + chosen = [components[0] - 1, components[1] - 1] + + assert up_to_sign(result.coordinates[:, chosen], data[f"{case} coordinates"]) <= ( + EXACT + ) + assert scaled(result.variance_explained[chosen], data[f"{case} percent"]) <= EXACT + # 500 is more than the genes: every gene, in decreasing order of variance. + assert result.genes.size == min(top, data["values"].shape[1]) + + +def test_genes_of_equal_variance_keep_their_order_as_plotPCA_keeps_them() -> None: + data = fixture("pca") + result = principal_components(data["ties"], top=20) + + # Genes 18 to 21 (from 0) have equal variance; the 20 kept end with 18 and 19. + assert result.genes.tolist()[-2:] == [18, 19] + assert up_to_sign(result.coordinates[:, :2], data["ties coordinates"]) <= EXACT + + +def test_a_cube_group_is_prepared_as_edgeR_prepares_it_for_plotPCA() -> None: + data = fixture("pca") + result = principal_components(pca_cube()) + + # The two samples with 5 cells are left out. + assert result.samples["sample_id"].tolist() == SAMPLES[:10] + assert result.samples["condition"].tolist() == CONDITIONS[:10] + assert up_to_sign(result.coordinates[:, :2], data["cube coordinates"]) <= EXACT + assert result.settings == {"group": "T", "min_cells": 10, "top": 500} + + +def test_the_figure_draws_the_points_and_labels_of_plotPCA() -> None: + data = fixture("pca") + ax = plotting.plot_pca(data["values"], top=100, labels=True, ax=axes()) + drawn = np.array([label.xy for label in artists(ax, "sample labels")]) + + assert up_to_sign(drawn, data["figure coordinates"]) <= EXACT + assert [ax.get_xlabel(), ax.get_ylabel()] == data["figure labels"].tolist() + assert ax.get_aspect() == 1.0 # plotPCA's coord_fixed + + +@pytest.mark.parametrize("case", list(DOT_CASES)) +def test_dot_plot_areas_and_colours_are_scanpys(case: str) -> None: + options = {"genes": GENES} | DOT_CASES[case] + genes = options.pop("genes") + ax = plotting.plot_dotplot(dot_cube(), genes, ax=axes(), **options) + dots = artists(ax, "dots")[0] + reference = fixture("dot_plots") + texts = dot_texts()[case] + + assert scaled(dots.get_sizes(), reference[f"{case} sizes"]) <= EXACT + assert scaled(dots.get_array(), reference[f"{case} color"].ravel()) <= EXACT + assert [t.get_text() for t in ax.get_yticklabels()] == texts["rows"] + assert [t.get_text() for t in ax.get_xticklabels()] == texts["columns"] + legends = [a.get_legend() for a in ax.figure.axes if a.get_legend() is not None] + assert [t.get_text() for t in legends[0].get_texts()] == texts["key"] + if texts["gene_groups"] is not None: + labels, positions = texts["gene_groups"] + brackets = [line.get_xdata() for line in artists(ax, "gene groups")] + assert [[round(x[0] + 0.3), round(x[-1] - 0.3)] for x in brackets] == positions + assert [t.get_text() for t in artists(ax, "gene group names")] == labels + + +@pytest.mark.parametrize("standard_scale", [None, "gene"]) +def test_dot_plot_labels_its_colour_bar_clear_of_the_key( + standard_scale: str | None, +) -> None: + """The colour bar's label is written above the bar, as scanpy writes it, and the + key of the dots' areas below it: a label along the bar ran into the key.""" + ax = plotting.plot_dotplot( + dot_cube(), GENES, standard_scale=standard_scale, ax=axes() + ) + figure = ax.figure + canvas = importlib.import_module("matplotlib.backends.backend_agg").FigureCanvasAgg( + figure + ) + canvas.draw() + renderer = canvas.get_renderer() + bar = figure.axes[-1] + title = bar.get_title(loc="left") + written = [ + child + for child in bar.get_children() + if hasattr(child, "get_text") and child.get_text() == title + ] + label = written[0].get_window_extent(renderer) + key = bar.get_legend().get_window_extent(renderer) + frame = bar.get_window_extent(renderer) + + assert title.startswith("Mean expression") + assert not label.overlaps(key) + assert not label.overlaps(frame) and label.y0 >= frame.y1 + assert not key.overlaps(frame) and key.y1 <= frame.y0 + + +# Values derived by hand, and invariants. + + +def test_three_samples_on_two_genes_whose_components_are_the_genes() -> None: + # Centred, gene 0 is (-1, 1, 0) and gene 1 (-1/3, -1/3, 2/3): orthogonal, with + # sums of squares 2 and 2/3, so the components are the genes with shares 3/4 and + # 1/4. Of PC1's tied largest coordinates, the first sample's is made positive. + values = np.array([[0.0, 0.0], [2.0, 0.0], [1.0, 1.0]]) + result = principal_components(values) + + assert_allclose( + result.coordinates, [[1, -1 / 3], [-1, -1 / 3], [0, 2 / 3]], atol=1e-15 + ) + assert_allclose(result.variance_explained, [0.75, 0.25], rtol=1e-15) + assert_allclose(np.abs(result.loadings), np.eye(2), atol=1e-15) + assert result.genes.tolist() == [0, 1] + assert result.samples["sample_id"].tolist() == ["1", "2", "3"] + + +def test_two_samples_have_one_component_along_their_difference() -> None: + # The genes of variance 8 and 4.5 are kept, the constant one is not; centred, the + # first sample is (-2, -1.5) on them, at distance 2.5 from the centre. + values = np.array([[0.0, 0.0, 5.0], [3.0, 4.0, 5.0]]) + result = principal_components(values, top=2) + + assert result.genes.tolist() == [1, 0] + assert_allclose(result.coordinates, [[2.5], [-2.5]], rtol=1e-15) + assert_allclose(result.loadings, [[-0.8], [-0.6]], rtol=1e-15) + assert_allclose(result.variance_explained, [1.0]) + with pytest.raises(ValueError, match="need 2 components; the data have 1"): + plotting.plot_pca(result, ax=axes()) + + +def test_one_gene_gives_its_centred_values() -> None: + values = np.array([[1.0, 9.0], [4.0, 9.5], [7.0, 10.0], [0.0, 11.0]]) + result = principal_components(values, top=1) + + assert result.genes.tolist() == [0] # variance 10 against 0.729 + assert_allclose(result.coordinates[:, 0], [-2.0, 1.0, 4.0, -3.0]) + + +@cache +def random_values() -> np.ndarray: + rng = np.random.default_rng(7) + values = rng.normal(size=(9, 60)) + values[:4, :12] += 2.5 + values[::3, 30:40] -= 1.5 + return values + + +def test_the_order_of_samples_and_genes_and_the_sign_of_the_values_do_not_matter() -> ( + None +): + values = random_values() + base = principal_components(values, top=25) + samples = np.array([4, 0, 8, 2, 6, 1, 3, 7, 5]) + genes = np.random.default_rng(8).permutation(60) + + by_sample = principal_components(values[samples], top=25) + assert_allclose(by_sample.coordinates, base.coordinates[samples], atol=1e-12) + by_gene = principal_components(values[:, genes], top=25) + assert_allclose(by_gene.coordinates, base.coordinates, atol=1e-12) + assert_array_equal(genes[by_gene.genes], base.genes) + negated = principal_components(-values, top=25) + assert_allclose(negated.coordinates, base.coordinates, atol=1e-12) + assert_allclose(negated.loadings, -base.loadings, atol=1e-12) + shifted = principal_components(values * 3.0 + np.arange(60.0), top=25) + assert_allclose(shifted.coordinates, 3.0 * base.coordinates, atol=1e-11) + assert_allclose(shifted.variance_explained, base.variance_explained, rtol=1e-12) + + +def test_components_are_orthonormal_and_their_shares_sum_to_one() -> None: + result = principal_components(random_values(), top=25) + centred = random_values()[:, result.genes] - random_values()[:, result.genes].mean( + 0 + ) + + assert result.n_components == 8 # one fewer than the samples + assert_allclose(result.loadings.T @ result.loadings, np.eye(8), atol=1e-12) + assert_allclose(centred @ result.loadings, result.coordinates, atol=1e-12) + assert_allclose(result.variance_explained.sum(), 1.0, rtol=1e-15) + assert np.all(np.diff(result.variance_explained) <= 0.0) + largest = np.argmax(np.abs(result.coordinates), axis=0) + assert np.all(result.coordinates[largest, np.arange(8)] > 0.0) + + +def test_results_are_read_only_records() -> None: + result = principal_components(random_values(), top=10) + + assert isinstance(result, PrincipalComponents) + with pytest.raises(ValueError, match="read-only"): + result.coordinates[0, 0] = 1.0 + assert "n_components=8" in repr(result) + + +def test_the_figure_colours_marks_and_labels_samples_by_their_variables() -> None: + cube = pca_cube() + ax = plotting.plot_pca( + cube, color_by="condition", marker_by="batch", labels=True, ax=axes() + ) + + collections = artists(ax, "samples") + assert len(collections) == 2 # one per batch + hues = {tuple(face) for c in collections for face in c.get_facecolors().tolist()} + assert len(hues) == 3 # one per condition + assert [t.get_text() for t in artists(ax, "sample labels")] == SAMPLES[:10] + legend = [t.get_text() for t in ax.get_legend().get_texts()] + assert legend == ["condition", "a", "b", "c", "batch", "x", "y"] + # Labels given one per sample of the cube, of which the plotted ten are shown. + ax = plotting.plot_pca(cube, labels=[f"L{i}" for i in range(12)], ax=axes()) + assert [t.get_text() for t in artists(ax, "sample labels")] == [ + f"L{i}" for i in range(10) + ] + + +def test_the_figure_draws_a_result_and_chosen_components() -> None: + result = principal_components(random_values(), top=25) + ax = plotting.plot_pca(result, components=(2, 3), labels=True, ax=axes()) + drawn = np.array([label.xy for label in artists(ax, "sample labels")]) + + assert_array_equal(drawn, result.coordinates[:, 1:3]) + share = 100 * result.variance_explained + assert ax.get_xlabel() == f"PC2: {share[1]:.0f}% variance" + assert ax.get_ylabel() == f"PC3: {share[2]:.0f}% variance" + + +@pytest.mark.parametrize( + ("arguments", "error", "match"), + [ + ({"data": np.ones((1, 3))}, ValueError, "at least two samples"), + ({"data": np.ones(3)}, ValueError, "2-D array"), + ({"data": np.array([[1.0, np.nan], [2.0, 3.0]])}, ValueError, "finite"), + ({"data": np.ones((3, 4))}, ValueError, "do not vary"), + # The mean of three values of 0.1 is not 0.1 in floating point: the centred + # values are rounding errors, not variation. + ({"data": np.full((3, 4), 0.1)}, ValueError, "do not vary"), + ({"data": np.ones((3, 0))}, ValueError, "no genes"), + ({"data": np.eye(3), "top": 0}, ValueError, "top must be at least 1"), + ({"data": np.eye(3), "top": 2.5}, TypeError, "top must be an integer"), + ({"data": np.eye(3), "group": "T"}, ValueError, "group applies to a cube"), + ( + {"data": np.eye(3), "min_cells": 3}, + ValueError, + "min_cells applies to a cube", + ), + ({"data": "values"}, TypeError, "real numbers"), + ], +) +def test_principal_components_refuse_invalid_arguments( + arguments: dict[str, Any], error: type, match: str +) -> None: + with pytest.raises(error, match=match): + principal_components(**arguments) + + +def test_a_cube_needs_its_group_and_enough_samples() -> None: + cube = dot_cube() + with pytest.raises(ValueError, match="the cube has 4 groups; name the one"): + principal_components(cube) + with pytest.raises(KeyError, match="Mono"): + principal_components(cube, group="Mono") + with pytest.raises(ValueError, match="principal components need two or more"): + principal_components(cube, group="pDC", min_cells=10_000) + with pytest.raises(ValueError, match="min_cells must be at least 0"): + principal_components(cube, group="B", min_cells=-1) + + +def test_the_figure_refuses_options_a_result_has_used() -> None: + result = principal_components(random_values(), top=25) + with pytest.raises(ValueError, match="top applies to the values"): + plotting.plot_pca(result, top=10, ax=axes()) + with pytest.raises(ValueError, match="two different components"): + plotting.plot_pca(result, components=(1, 1), ax=axes()) + with pytest.raises(TypeError, match="pair of components"): + plotting.plot_pca(result, components=1, ax=axes()) + + +def hand_cube(**options: Any) -> Cube: + """Two samples of one group: gene a detected in 2 of 4 and 3 of 6 cells, gene b in + every cell.""" + return Cube( + { + "counts": np.array([[[3, 9], [5, 12]]], dtype=np.uint32), + "detected": np.array([[[2, 4], [3, 6]]], dtype=np.uint32), + "log_sum": np.array([[[2.0, 4.0], [4.0, 6.0]]]), + }, + samples={"sample_id": ["s1", "s2"]}, + groups={"group_id": ["T"]}, + features={"feature_id": ["a", "b"], "symbol": ["GA", "GB"]}, + n_cells=np.array([[4], [6]]), + normalization=options.get( + "normalization", {"method": "library_size", "pseudocount": 1.0} + ), + ) + + +def test_a_gene_expressed_in_half_the_cells_of_two_samples() -> None: + ax = plotting.plot_dotplot(hand_cube(), ["GA", "b"], ax=axes()) + dots = artists(ax, "dots")[0] + + # Fractions 5/10 and 10/10; the largest is 1, so the areas are 200 (f / 1)^1.5. + assert_allclose(dots.get_sizes(), [200.0 * 0.5**1.5, 200.0], rtol=1e-15) + # Means 6/10 and 10/10 over every cell. + assert_allclose(dots.get_array(), [0.6, 1.0], rtol=1e-15) + assert [t.get_text() for t in ax.get_xticklabels()] == ["GA", "b"] + ax = plotting.plot_dotplot(hand_cube(), ["a"], mean_only_expressed=True, ax=axes()) + # 6 over the 5 cells that express it; the largest fraction, 0.5, has the largest dot. + assert_allclose(artists(ax, "dots")[0].get_array(), [1.2], rtol=1e-15) + assert_allclose(artists(ax, "dots")[0].get_sizes(), [200.0], rtol=1e-15) + ax = plotting.plot_dotplot(hand_cube(), ["a", "b"], samples=["s1"], ax=axes()) + assert_allclose(artists(ax, "dots")[0].get_array(), [0.5, 1.0], rtol=1e-15) + + +@cache +def cells() -> dict[str, Any]: + rng = np.random.default_rng(11) + n = 600 + counts = rng.poisson(rng.uniform(0.05, 3.0, 12), size=(n, 12)).astype(np.float64) + counts[counts.sum(axis=1) == 0, 0] = 1.0 + return { + "counts": counts, + "samples": rng.choice(["s1", "s2", "s3"], n), + "groups": rng.choice(["B", "T", "NK"], n), + "features": [f"x{j}" for j in range(12)], + } + + +def test_a_merged_group_and_pooled_samples_are_the_same_cells_given_one_label() -> None: + data = cells() + cube = from_cells(**data) + genes = ["x0", "x3", "x7", "x11"] + merged = cube.regroup({"BT": ["B", "T"]}) + relabelled = from_cells( + **(data | {"groups": np.where(data["groups"] == "NK", "NK", "BT")}) + ) + one_sample = from_cells(**(data | {"samples": np.full(data["groups"].size, "all")})) + + for first, second, groups in ( + (merged, relabelled, ["BT"]), + (cube, one_sample, ["NK", "T"]), + ): + ours = artists( + plotting.plot_dotplot(first, genes, groups=groups, ax=axes()), "dots" + ) + theirs = artists( + plotting.plot_dotplot(second, genes, groups=groups, ax=axes()), "dots" + ) + assert_allclose(ours[0].get_sizes(), theirs[0].get_sizes(), rtol=1e-14) + assert_allclose(ours[0].get_array(), theirs[0].get_array(), rtol=1e-14) + + +def test_scaled_means_run_from_zero_to_one() -> None: + cube = from_cells(**cells()) + genes = ["x0", "x3", "x7", "x11"] + for scale, axis in (("gene", 0), ("group", 1)): + ax = plotting.plot_dotplot(cube, genes, standard_scale=scale, ax=axes()) + values = np.asarray(artists(ax, "dots")[0].get_array()).reshape(3, 4) + assert_allclose(values.min(axis=axis), 0.0, atol=0) + assert_allclose(values.max(axis=axis), 1.0, rtol=1e-15) + # A gene equal in every group is 0 throughout. + ax = plotting.plot_dotplot( + hand_cube(), ["a", "b"], standard_scale="gene", ax=axes() + ) + assert_array_equal(artists(ax, "dots")[0].get_array(), [0.0, 0.0]) + + +def test_dot_plots_refuse_what_they_cannot_draw() -> None: + bulk = Cube( + {"counts": np.ones((1, 2, 2), dtype=np.uint32)}, + samples={"sample_id": ["s1", "s2"]}, + groups={"group_id": ["bulk"]}, + features={"feature_id": ["a", "b"]}, + ) + with pytest.raises(ValueError, match="the cube lacks cell numbers"): + plotting.plot_dotplot(bulk, ["a"], ax=axes()) + with pytest.raises(KeyError, match="no gene 'c'"): + plotting.plot_dotplot(hand_cube(), ["c"], ax=axes()) + with pytest.raises(ValueError, match="genes is empty"): + plotting.plot_dotplot(hand_cube(), [], ax=axes()) + with pytest.raises(ValueError, match="the group of genes 'g' is empty"): + plotting.plot_dotplot(hand_cube(), {"g": []}, ax=axes()) + with pytest.raises(KeyError, match="no group 'B'"): + plotting.plot_dotplot(hand_cube(), ["a"], groups=["B"], ax=axes()) + with pytest.raises(ValueError, match="names a sample more than once"): + plotting.plot_dotplot(hand_cube(), ["a"], samples=["s1", "s1"], ax=axes()) + with pytest.raises(ValueError, match="one entry per sample"): + plotting.plot_dotplot(hand_cube(), ["a"], samples=[True], ax=axes()) + with pytest.raises(ValueError, match="standard_scale must be one of"): + plotting.plot_dotplot(hand_cube(), ["a"], standard_scale="var", ax=axes()) + without = hand_cube(normalization={"method": "library_size", "pseudocount": 0.5}) + with pytest.raises(ValueError, match="pseudocount of 1"): + plotting.plot_dotplot(without, ["a"], mean_only_expressed=True, ax=axes()) + empty = from_cells( + **(cells() | {"groups": np.where(cells()["samples"] == "s1", "A", "B")}) + ) + with pytest.raises( + ValueError, match="group 'A' has no cells in the samples pooled" + ): + plotting.plot_dotplot(empty, ["x0"], samples=["s2"], ax=axes()) + twice = Cube( + { + "counts": np.ones((1, 1, 2), dtype=np.uint32), + "detected": np.ones((1, 1, 2), dtype=np.uint32), + "log_sum": np.ones((1, 1, 2)), + }, + samples={"sample_id": ["s1"]}, + groups={"group_id": ["T"]}, + features={"feature_id": ["a", "b"], "symbol": ["G", "G"]}, + n_cells=np.array([[1]]), + normalization={"method": "library_size"}, + ) + with pytest.raises(ValueError, match="the symbol 'G' names 2 genes"): + plotting.plot_dotplot(twice, ["G"], ax=axes()) diff --git a/tests/test_special.py b/tests/test_special.py index e94658d..08762b6 100644 --- a/tests/test_special.py +++ b/tests/test_special.py @@ -1,9 +1,12 @@ +import itertools import math +import time from fractions import Fraction from pathlib import Path import numpy as np import pytest +from numpy.testing import assert_allclose, assert_array_equal from cytepool.stats import special @@ -131,6 +134,85 @@ def test_incomplete_beta_matches_high_precision_values() -> None: assert_log_close(special.log_betaincc(a, b, x), data["logC"]) +def test_incomplete_beta_matches_high_precision_values_at_extreme_arguments() -> None: + # Shapes from the smallest subnormal to the largest double, and points near the + # center of shapes from 1e10, where the uniform expansion serves them. + data = _load("betainc_extreme") + a, b, x = data["a"], data["b"], data["x"] + assert_probability_close(special.betainc(a, b, x), data["I"]) + assert_probability_close(special.betaincc(a, b, x), data["C"]) + assert_log_close(special.log_betainc(a, b, x), data["logI"]) + assert_log_close(special.log_betaincc(a, b, x), data["logC"]) + + +LARGEST = float(np.finfo(np.float64).max) +EDGES = [0.0, 5e-324, 1e-300, 1e-8, 0.5, 1.0 - 2.0**-53, 1.0, 2.0, 30.0, 1e8, 1e300] +EDGES.append(LARGEST) + + +def test_incomplete_beta_is_fast_at_every_extreme_argument() -> None: + # Some of these once ran 10^7 steps of the continued fraction (0.15 s a value), + # and with a + b beyond the double range a deviance term never returned. + a, b, x = (np.array(v) for v in zip(*itertools.product(EDGES, repeat=3))) + special.log_betainc(2.0, 3.0, 0.5) # compiled + start = time.perf_counter() + log_i = special.log_betainc(a, b, x) + log_c = special.log_betaincc(a, b, x) + elapsed = time.perf_counter() - start + + # 3,456 values take about 0.01 s; the bound leaves room for slow machines. + assert elapsed < 2.0 + domain = (a > 0.0) & (b > 0.0) & (x <= 1.0) + assert not np.isnan(log_i[domain]).any() and not np.isnan(log_c[domain]).any() + assert np.isnan(log_i[~domain]).all() and np.isnan(log_c[~domain]).all() + # One tail is the complement of the other. + both = domain & (log_i > -700.0) & (log_c > -700.0) + assert_allclose(np.exp(log_i[both]) + np.exp(log_c[both]), 1.0, rtol=4 * EPS) + + +def test_incomplete_beta_of_equal_shapes_is_one_half_at_one_half() -> None: + shapes = np.array([1e10, 3e15, 1e20, 1e300, LARGEST]) + assert_array_equal(special.betainc(shapes, shapes, 0.5), 0.5) + assert_array_equal(special.betaincc(shapes, shapes, 0.5), 0.5) + + +def test_incomplete_beta_near_the_center_of_large_shapes_is_its_reflection() -> None: + # I_x(a, b) = 1 - I_(1-x)(b, a), with x and 1 - x exact, in the band of the uniform + # expansion (about one standard deviation, 2.2e-7, of the center) and beyond it, + # where the continued fraction takes over. + a, b = 1e12, 3e12 + sd = math.sqrt(3.0 / 16.0 / (a + b + 1.0)) + x = 0.25 + np.arange(-6, 7) * 2.0**-23 + assert np.abs(x - 0.25).max() > 3.0 * sd + total = special.betainc(a, b, x) + special.betainc(b, a, 1.0 - x) + assert_allclose(total, 1.0, rtol=2 * EPS) + assert np.all(np.diff(special.betainc(a, b, x)) > 0.0) + + +def test_incomplete_gamma_returns_where_the_shape_and_the_argument_overflow() -> None: + # Loader's deviance term of a + x beyond the double range once never returned. + # mpmath, from Legendre's continued fraction at 40 digits, gives + # log Q(1e300, LARGEST) = -1.7976929347904657564e308. + assert_allclose(special.log_gammaincc(1e300, LARGEST), -1.7976929347904657564e308) + assert special.gammainc(LARGEST, LARGEST) == 0.5 # 1/2 + 1/(3 sqrt(2 pi a)) + assert special.log_gammainc(LARGEST, 1e300) == -np.inf # below -LARGEST + + +@pytest.mark.parametrize("a", [1e-8, 1e-20, 1e-300, 5e-324]) +def test_incomplete_beta_of_a_tiny_shape_against_its_closed_form(a: float) -> None: + # 1 - I_x(a, 1) = 1 - x^a = -expm1(t) with t = a log x, whose logarithm is + # log(a) + log(-log x) + log(expm1(t) / t); t is subnormal for the smallest a. + x = np.array([1e-300, 1e-8, 0.5, 0.99]) + t = a * np.log(x) + ratio = np.where(np.abs(t) > 1e-300, np.expm1(t) / np.where(t == 0.0, 1.0, t), 1.0) + expected = np.log(a) + np.log(-np.log(x)) + np.log(ratio) + assert_allclose(special.log_betaincc(a, 1.0, x), expected, rtol=4 * EPS) + # 1 - x is exact from x = 1/2. + assert_allclose( + special.log_betainc(1.0, a, 1.0 - x[2:]), expected[2:], rtol=4 * EPS + ) + + # --------------------------------------------------------------------------------------- # Closed forms # --------------------------------------------------------------------------------------- diff --git a/validation/distributions.py b/validation/distributions.py index de21f31..b51ddb4 100644 --- a/validation/distributions.py +++ b/validation/distributions.py @@ -7,8 +7,12 @@ tails (probabilities down to 1e-300, and log probabilities far below), shapes from 1e-3 to 1e6 (incomplete gamma shapes to 1e18 above the center), degrees of freedom from 0.5 to 1e12 and infinity, hypergeometric populations up to 1e13, and the hard cases where -textbook methods and scipy lose accuracy. scipy, when installed, is evaluated on the -same points for information only. +textbook methods and scipy lose accuracy. For the incomplete beta function they also +cover every combination of shapes and arguments at the edges of the double range, +from the smallest subnormal to the largest double, and points near the center of +shapes from 1e10, and time it per value at every combination of those values as +shapes and arguments, those outside the domain included. scipy, when installed, is +evaluated on the same points for information only. Requirements: Python: numpy and numba (for CytePool) and mpmath. scipy is optional. @@ -24,11 +28,15 @@ ``inputs.npz`` from the directory and writes ``reference.npz``, ``scipy.npz``, ``versions.json`` and ``notes.json`` back; it imports mpmath and scipy but not CytePool. The incomplete beta reference is DLMF 8.17.22 evaluated in 70-digit arithmetic, checked -against mpmath's own ``betainc`` where that converges quickly; the hypergeometric -reference sums exact ratios in fixed-point integer arithmetic. The script prints +against mpmath's own ``betainc`` where that converges quickly; at the edges of the +double range, it is evaluated in 1200 to 2400 bits, and near the center of large +shapes, where it converges too slowly, the beta density is integrated numerically +instead. The hypergeometric reference sums exact ratios in fixed-point integer +arithmetic. The script prints the reference versions and a markdown table, and exits with status 1 when any check exceeds its tolerance. With ``--write-fixtures`` it writes a few hundred points per -family, the hardest included, with their reference values to tests/data/distributions/. +family, the hardest included, with their reference values to tests/data/distributions/, +and every point of the extreme arguments of the incomplete beta function. """ import argparse @@ -58,9 +66,26 @@ NEAR_MACHINE = 1e-15 # Probabilities below this are checked through their logarithm. SMALLEST_CHECKED = 1e-300 +# The edges of the double range for the incomplete beta function: shapes and +# arguments of every size from the smallest subnormal to the largest double. The +# timing runs every combination of TIMED_VALUES as a, b and x, those outside the +# domain included, and every value must take at most MOST_TIME seconds (best of three). +LARGEST = 1.7976931348623157e308 +EXTREME_SHAPES = [5e-324, 1e-300, 1e-8, 0.5, 1.0 - 2.0**-53, 1.0, 2.0, 30.0, 1e8] +EXTREME_SHAPES += [1e300, LARGEST] +EXTREME_ARGUMENTS = [0.0, 5e-324, 1e-300, 1e-8, 0.5, 1.0 - 2.0**-53, 1.0] +TIMED_VALUES = [0.0, 5e-324, 1e-300, 1e-8, 0.5, 1.0 - 2.0**-53, 1.0, 2.0, 30.0, 1e8] +TIMED_VALUES += [1e300, LARGEST] +MOST_TIME = 1e-2 +# Where the reference integrates the beta density near its center: from 40 standard +# deviations below it to 40 above, split at these. +QUADRATURE_MARKS = (-40, -16, -8, -4, -2, -1, 0, 1, 2, 4, 8, 16, 40) # The tests check these quantiles through the tails, which have fixtures, and closed # forms; their families get no fixtures of their own. UNFIXTURED = {"f_quantile", "chi2_quantile"} +# Families whose fixtures keep every point (and draw no random sample, so that the +# samples of the other families stay as they are). +FIXTURE_EVERY_POINT = {"betainc_extreme"} # --------------------------------------------------------------------------------------- @@ -185,6 +210,27 @@ def build_inputs() -> tuple[dict[str, np.ndarray], dict[str, np.ndarray]]: inputs["betainc__a"], inputs["betainc__b"], inputs["betainc__x"] = beta_points.T hard["betainc"] = np.arange(5) + # Incomplete beta at the edges of the double range: every shape and argument of + # EXTREME_SHAPES and EXTREME_ARGUMENTS; then points near the center of large shapes, + # within about one standard deviation, where the uniform expansion serves them, and + # just outside, where the continued fraction does. + points = [ + (a_, b_, x) for a_ in EXTREME_SHAPES for b_ in EXTREME_SHAPES + for x in EXTREME_ARGUMENTS + ] # fmt: skip + first_band = len(points) + pairs = [(1e10, 1e10), (1e10, 3e10), (1e12, 1e12), (1e12, 3e12), (2e15, 7e15)] + pairs += [(1e20, 1e20), (1e10, 1e20)] + for a_, b_ in pairs: + center = a_ / (a_ + b_) + sd = math.sqrt(a_ / (a_ + b_) * (b_ / (a_ + b_)) / (a_ + b_ + 1.0)) + points += [(a_, b_, center + z * sd) for z in (0.0, 0.4, -0.9, 1.1, -1.5)] + extreme = np.array(points) + keys = ("betainc_extreme__a", "betainc_extreme__b", "betainc_extreme__x") + for key, column in zip(keys, extreme.T): + inputs[key] = column + hard["betainc_extreme"] = np.arange(first_band, len(points)) + # Normal distribution and quantiles. tails = np.concatenate([np.linspace(36.5, 38.6, 43), _log_grid(1, 150, 60)]) inputs["norm__x"] = np.concatenate([np.linspace(-40.0, 40.0, 801), tails, -tails]) @@ -402,11 +448,21 @@ def reference_main(directory: Path) -> None: no_convergence = mp.libmp.libhyper.NoConvergence inputs = dict(np.load(directory / "inputs.npz")) - def ibeta_fraction(A: Any, B: Any, X: Any, Y: Any) -> Any: + def ibeta_fraction( + A: Any, + B: Any, + X: Any, + Y: Any, + log_x: Any = None, + log_y: Any = None, + digits: int = 55, + ) -> Any: """I_X(A, B) for X < (A + 1) / (A + B + 2) from DLMF 8.17.22. The plain continued fraction is evaluated backward in 70-digit arithmetic, with - the depth doubled until two values agree to 55 digits. + the depth doubled until two values agree to 55 digits (or ``digits``). log_x and + log_y, when given, replace log(X) and log(Y) in the prefactor (Y = 1 - X rounds + to 1 at 70 digits when X is tiny, but its logarithm need not). """ def value(depth: int) -> Any: @@ -425,13 +481,75 @@ def value(depth: int) -> Any: while True: depth *= 2 current = value(depth) - if abs(current - previous) <= mp.mpf(10) ** -55 * abs(current): + if abs(current - previous) <= mp.mpf(10) ** -digits * abs(current): break previous = current log_beta = mp.loggamma(A) + mp.loggamma(B) - mp.loggamma(A + B) - log_pre = A * mp.log(X) + B * mp.log(Y) - mp.log(A) - log_beta + log_x = mp.log(X) if log_x is None else log_x + log_y = mp.log(Y) if log_y is None else log_y + log_pre = A * log_x + B * log_y - mp.log(A) - log_beta return mp.exp(log_pre) * current + def ibeta_by_quadrature(A: Any, B: Any, X: Any) -> tuple[Any, Any]: + """Both tails of the beta distribution near the center of large shapes. + + The density t^(A-1) (1-t)^(B-1) / B(A, B) is integrated numerically on each side + of X, out to 40 standard deviations from the center, beyond which the mass is + below 1e-340, with enough digits for A log t to keep 50 after its cancellation. + """ + digits = int(70 + mp.log10(A + B)) + with mp.workdps(digits): + log_norm = mp.loggamma(A + B) - mp.loggamma(A) - mp.loggamma(B) + center = A / (A + B) + sd = mp.sqrt(A * B / (A + B) ** 2 / (A + B + 1)) + + def density(t: Any) -> Any: + return mp.exp(log_norm + (A - 1) * mp.log(t) + (B - 1) * mp.log1p(-t)) + + marks = [center + k * sd for k in QUADRATURE_MARKS] + below = [max(mp.mpf(0), marks[0])] + [s for s in marks if s < X] + [X] + above = [X] + [s for s in marks if s > X] + [min(mp.mpf(1), marks[-1])] + lower, upper = mp.quad(density, below), mp.quad(density, above) + if abs(lower + upper - 1) > mp.mpf(10) ** -45: + raise RuntimeError(f"the quadrature lost accuracy at {A}, {B}, {X}") + return lower, upper + + def ibeta_extreme(A: Any, B: Any, X: Any) -> tuple[Any, Any]: + """I_X(A, B) and 1 - I_X(A, B) anywhere in the domain, at 50 digits. + + I is 1/2 at the center of equal shapes; near the center of shapes from 1e9 the + density is integrated; elsewhere the continued fraction converges, as in + ibeta_pair, with 1200 bits, so that 1 - X is exact down to the smallest + subnormal X. When the tail it gives is within 1e-30 of 1, its complement would + keep fewer than 25 digits, and the fraction is evaluated again with 2400 bits, + enough for A + B when one shape is the smallest subnormal and the other the + largest double, and to 400 digits, which leaves more than 50 for a tail of + 1e-324. + """ + if X == 0: + return mp.mpf(0), mp.mpf(1) + if X == 1: + return mp.mpf(1), mp.mpf(0) + if A == B and X == mp.mpf(1) / 2: + return mp.mpf(1) / 2, mp.mpf(1) / 2 + with mp.workdps(work): + center = A / (A + B) + sd = mp.sqrt(A * B / (A + B) ** 2 / (A + B + 1)) + if min(A, B) >= 1e9 and abs(X - center) <= 3 * sd: + return ibeta_by_quadrature(A, B, X) + for bits, digits in ((1200, 55), (2400, 400)): + with mp.workprec(bits): + Y, log_x, log_y = 1 - X, mp.log(X), mp.log1p(-X) + if X * (A + B + 2) < A + 1: + lower = ibeta_fraction(A, B, X, Y, log_x, log_y, digits) + upper = 1 - lower + else: + upper = ibeta_fraction(B, A, Y, X, log_y, log_x, digits) + lower = 1 - upper + if min(lower, upper) > mp.mpf(10) ** -30: + break + return lower, upper + def ibeta_pair(A: Any, B: Any, X: Any, Y: Any) -> tuple[Any, Any]: """I_X(A, B) and 1 - I_X(A, B) for Y = 1 - X. @@ -830,6 +948,16 @@ def store( rows.append((float(lower), float(upper), *logs(lower, upper))) store("betainc", ("I", "C", "logI", "logC"), rows) + rows = [] + for a_, b_, x in zip( + inputs["betainc_extreme__a"], + inputs["betainc_extreme__b"], + inputs["betainc_extreme__x"], + ): + lower, upper = ibeta_extreme(mp.mpf(a_), mp.mpf(b_), mp.mpf(x)) + rows.append((float(lower), float(upper), *logs(lower, upper))) + store("betainc_extreme", ("I", "C", "logI", "logC"), rows) + rows = [] for x in inputs["norm__x"]: with mp.workdps(work): @@ -1040,6 +1168,11 @@ def cytepool_values(inputs: dict[str, np.ndarray]) -> dict[str, np.ndarray]: hyper = [i[f"hypergeom__{c}"] for c in "kMnN"] a_g, x_g = i["gammainc__a"], i["gammainc__x"] abx = (i["betainc__a"], i["betainc__b"], i["betainc__x"]) + edges = ( + i["betainc_extreme__a"], + i["betainc_extreme__b"], + i["betainc_extreme__x"], + ) t_args = (i["t__x"], i["t__df"]) tg_args = (i["t_genes__x"], i["t_genes__df"]) f_args = (i["f__x"], i["f__dfn"], i["f__dfd"]) @@ -1064,6 +1197,10 @@ def cytepool_values(inputs: dict[str, np.ndarray]) -> dict[str, np.ndarray]: "betainc__C": s.betaincc(*abx), "betainc__logI": s.log_betainc(*abx), "betainc__logC": s.log_betaincc(*abx), + "betainc_extreme__I": s.betainc(*edges), + "betainc_extreme__C": s.betaincc(*edges), + "betainc_extreme__logI": s.log_betainc(*edges), + "betainc_extreme__logC": s.log_betaincc(*edges), "norm__cdf": d.norm_cdf(i["norm__x"]), "norm__sf": d.norm_sf(i["norm__x"]), "norm__logcdf": d.norm_logcdf(i["norm__x"]), @@ -1180,6 +1317,17 @@ def add( ]: add(name, key, "probability", CONTRACT) add(f"log_{name}", key.replace("__", "__log"), "log", CONTRACT) + for name, key in [ + ("betainc", "betainc_extreme__I"), + ("betaincc", "betainc_extreme__C"), + ]: + add(f"{name}, extreme arguments", key, "probability", CONTRACT) + add( + f"log_{name}, extreme arguments", + key.replace("__", "__log"), + "log", + CONTRACT, + ) distributions = [ ("norm", ("cdf", "sf", "logcdf", "logsf")), ("t", ("cdf", "sf", "logsf")), @@ -1329,6 +1477,25 @@ def compare( return table, worst, scipy_table, passed, errors_by_family +def time_extreme_arguments() -> np.ndarray: + """The time per value of log_betainc at every combination of TIMED_VALUES as a, b + and x, the best of three calls each, in seconds (after one call that compiles).""" + import itertools + + from cytepool.stats import special as s + + s.log_betainc(2.0, 3.0, 0.5) + times = [] + for a_, b_, x in itertools.product(TIMED_VALUES, repeat=3): + best = math.inf + for _ in range(3): + start = time.perf_counter() + s.log_betainc(a_, b_, x) + best = min(best, time.perf_counter() - start) + times.append(best) + return np.array(times) + + def _markdown(header: list[str], rows: list[list[str]]) -> str: lines = [ "| " + " | ".join(header) + " |", @@ -1352,13 +1519,16 @@ def write_fixtures( keys_in = [k for k in inputs if _family(k) == family] keys_ref = [k for k in reference if _family(k) == family] size = len(inputs[keys_in[0]]) - chosen = set(hard.get(family, np.array([], dtype=int)).tolist()) - errors = errors_by_family.get(family, np.zeros(size)) - chosen |= set(np.argsort(-errors)[:40].tolist()) - remaining = np.setdiff1d(np.arange(size), np.fromiter(chosen, dtype=int)) - extra = max(0, min(FIXTURE_POINTS - len(chosen), remaining.size)) - chosen |= set(rng.choice(remaining, extra, replace=False).tolist()) - index = np.array(sorted(chosen)) + if family in FIXTURE_EVERY_POINT: + index = np.arange(size) + else: + chosen = set(hard.get(family, np.array([], dtype=int)).tolist()) + errors = errors_by_family.get(family, np.zeros(size)) + chosen |= set(np.argsort(-errors)[:40].tolist()) + remaining = np.setdiff1d(np.arange(size), np.fromiter(chosen, dtype=int)) + extra = max(0, min(FIXTURE_POINTS - len(chosen), remaining.size)) + chosen |= set(rng.choice(remaining, extra, replace=False).tolist()) + index = np.array(sorted(chosen)) arrays = {k.split("__")[1]: inputs[k][index] for k in keys_in} arrays |= {k.split("__")[1]: reference[k][index] for k in keys_ref} np.savez_compressed(FIXTURE_DIRECTORY / f"{family}.npz", **arrays) @@ -1399,6 +1569,19 @@ def main() -> int: table, worst, scipy_table, passed, errors_by_family = compare( inputs, reference, ours, scipy_values ) + times = time_extreme_arguments() + fast = float(times.max()) <= MOST_TIME + passed &= fast + table.append( + [ + f"betainc, time per value at {times.size} extreme arguments", + "timed", + f"largest (median {np.median(times):.1e} s)", + f"{times.max():.2e} s", + f"{MOST_TIME:.0e} s", + "pass" if fast else "FAIL", + ] + ) print("\nReference versions: " + ", ".join(f"{k} {v}" for k, v in versions.items())) print( diff --git a/validation/pca_and_dot_plots.py b/validation/pca_and_dot_plots.py new file mode 100644 index 0000000..e607628 --- /dev/null +++ b/validation/pca_and_dot_plots.py @@ -0,0 +1,654 @@ +"""Check principal components and their figure against DESeq2, and dot plots against +scanpy. + +``principal_components`` and ``plot_pca`` are compared with DESeq2's +``plotPCA(returnData = TRUE)`` and the ggplot it draws, on the same values: 12 samples' +log2 counts per million (three conditions, two batches, negative binomial counts of 300 +genes, seed 20261103), with the 100 genes of largest variance (components 1 and 2, and +3 and 4) and with 500, more than there are; a matrix whose genes 19 to 22 have equal +variance, with the 20 of largest variance, which keeps the first of them; and the +samples of a cube's group, two of which have fewer than 10 cells, which R prepares as +the edgeR User's Guide does (``filterByExpr``, ``normLibSizes``, ``cpm(log = TRUE, +prior.count = 2)``). Coordinates are compared up to the sign of each component, which +R's ``prcomp`` leaves arbitrary; the shares of variance, the axis labels, the equal +scales of the axes and the samples' groups of colour must agree. With two samples +DESeq2 also returns a second component, of rounding noise, which CytePool leaves out. + +``plot_dotplot`` is compared with ``scanpy.pl.dotplot`` on the cells that the cube was +pooled from: 4 cell types in 6 samples (one sample without cells of one type, seed +20261104), 40 genes, the cells normalized by scanpy (``normalize_total`` and ``log1p``) +and the cube by ``from_cells``. The areas of the dots and the fractions they encode, the +values their colours show, plain, scaled by gene or by group, and over the cells that +express a gene only, for the cells of every sample and for three samples and three cell +types in a chosen order, and the brackets of genes given in groups and the key of the +areas, must agree. + +Tiers: exact (1e-10, relative to max(1, abs(reference))) for values both programs +compute by the same deterministic steps; texts and groupings must be identical. + +Requirements: + Python 3.13 with numpy, numba and matplotlib, run from the repository root so that + ``cytepool`` imports from the working tree. + R 4.4 or newer with DESeq2, edgeR, ggplot2 and jsonlite (``RSCRIPT``, default + ``Rscript``); R runs with 8 GB of address space at most. + Python with scanpy 1.12 (``SCANPY_PYTHON``, default this interpreter). + +Usage: + PYTHONPATH=. python validation/pca_and_dot_plots.py + PYTHONPATH=. python validation/pca_and_dot_plots.py --write-fixtures +""" + +import argparse +import json +import os +import resource +import subprocess +import sys +import tempfile +import textwrap +from pathlib import Path +from typing import Any + +import numpy as np +from matplotlib.figure import Figure + +from cytepool import plotting +from cytepool.cube import Cube, from_cells +from cytepool.normalization import log_cpm, tmm_factors +from cytepool.principal_components import principal_components + +FIXTURE_DIRECTORY = ( + Path(__file__).resolve().parents[1] / "tests" / "data" / "pca_and_dot_plots" +) +EXACT = 1e-10 +# scanpy's dot plots: the area of the largest dot and the power of the fractions. +LARGEST_DOT, SIZE_EXPONENT = 200.0, 1.5 +MEMORY = 8_000_000_000 # bytes of address space for R +PCA_SEED, DOT_SEED = 20261103, 20261104 +CONDITIONS = ["a"] * 4 + ["b"] * 4 + ["c"] * 4 +BATCHES = ["x", "y"] * 6 +SAMPLES = [f"S{i}" for i in range(1, 13)] +PCA_CASES = {"top100": (100, (1, 2)), "top100_pc34": (100, (3, 4))} +PCA_CASES["top500"] = (500, (1, 2)) +GROUPS = ["B", "T", "NK", "pDC"] +CELL_SAMPLES = [f"D{i}" for i in range(1, 7)] +MARKERS = {"B cells": ["g00", "g01", "g02"], "T cells": ["g03", "g04"]} +MARKERS["others"] = ["g10", "g20", "g30", "g39"] +GENES = [gene for genes in MARKERS.values() for gene in genes] +DOT_CASES: dict[str, dict[str, Any]] = { + "plain": {}, + "scaled by gene": {"standard_scale": "gene"}, + "scaled by group": {"standard_scale": "group"}, + "expressing cells": {"mean_only_expressed": True}, + "chosen samples and groups": { + "samples": ["D2", "D3", "D5"], + "groups": ["NK", "B", "T"], + }, + "genes in groups": {"genes": MARKERS}, +} + +R_PROGRAM = r""" +suppressPackageStartupMessages({ + library(DESeq2); library(SummarizedExperiment); library(edgeR); library(ggplot2) +}) +root <- commandArgs(trailingOnly = TRUE)[[1]] +rd <- function(name) readBin(file.path(root, name), "double", + n = file.size(file.path(root, name)) / 8, size = 8) +samples <- jsonlite::read_json(file.path(root, "samples.json"), simplifyVector = TRUE) +info <- data.frame(condition = factor(samples$condition), batch = factor(samples$batch), + row.names = samples$sample_id) +transformed <- function(values, rows = seq_len(nrow(info))) { + colnames(values) <- rownames(info)[rows] + DESeqTransform(SummarizedExperiment(list(values), colData = info[rows, , drop = FALSE])) +} +pca <- function(object, ntop, pcs = 1:2) { + d <- suppressMessages(plotPCA(object, intgroup = c("condition", "batch"), + ntop = ntop, returnData = TRUE, pcsToUse = pcs)) + list(x = d[[1]], y = d[[2]], percent = attr(d, "percentVar"), + name = as.character(d$name), group = as.character(d$group)) +} +values <- matrix(rd("values.bin"), ncol = 12) # genes x samples +out <- list() +object <- transformed(values) +out$top100 <- pca(object, 100) +out$top100_pc34 <- pca(object, 100, 3:4) +out$top500 <- pca(object, 500) +ties <- matrix(rd("ties.bin"), ncol = 12) +out$ties <- pca(transformed(ties), 20) +out$two_samples <- pca(transformed(values[, 1:2], 1:2), 100) +p <- suppressWarnings(suppressMessages( + plotPCA(object, intgroup = c("condition", "batch"), ntop = 100))) +layer <- ggplot_build(p)$data[[1]] +out$figure <- list(x = layer$x, y = layer$y, colour = layer$colour, + xlabel = p$labels$x, ylabel = p$labels$y, + ratio = p$coordinates$ratio) +counts <- matrix(rd("counts.bin"), ncol = 12) +used <- as.logical(rd("used.bin")) +y <- DGEList(counts[, used]) +y <- normLibSizes(y[suppressWarnings(filterByExpr(y)), , keep.lib.sizes = FALSE]) +out$cube <- pca(transformed(cpm(y, log = TRUE, prior.count = 2), which(used)), 500) +out$versions <- c(paste("R", getRversion()), paste("DESeq2", packageVersion("DESeq2")), + paste("edgeR", packageVersion("edgeR")), + paste("ggplot2", packageVersion("ggplot2"))) +jsonlite::write_json(out, file.path(root, "out.json"), digits = NA, auto_unbox = TRUE) +""" + +SCANPY_PROGRAM = textwrap.dedent( + ''' + """scanpy's dot plots of the cells, normalized by scanpy.""" + + import json + import sys + import warnings + from importlib.metadata import version + + import matplotlib + + matplotlib.use("Agg") + import anndata + import numpy as np + import pandas as pd + import scanpy as sc + import scipy.sparse as sp + + warnings.filterwarnings("ignore") + directory = sys.argv[1] + data = np.load(f"{directory}/cells.npz") + cases = json.loads(open(f"{directory}/cases.json").read()) + obs = pd.DataFrame( + { + "sample": pd.Categorical(data["samples"]), + "cell_type": pd.Categorical(data["groups"], categories=list(data["order"])), + }, + index=[f"cell{i}" for i in range(data["groups"].size)], + ) + cells = anndata.AnnData( + X=sp.csr_matrix(data["counts"].astype(np.float64)), + obs=obs, + var=pd.DataFrame(index=list(data["genes"])), + ) + sc.pp.normalize_total(cells, target_sum=1e4) + sc.pp.log1p(cells) + out = {} + for name, case in cases.items(): + subset = cells + if case["samples"] is not None: + subset = subset[subset.obs["sample"].isin(case["samples"])] + if case["groups"] is not None: + subset = subset[subset.obs["cell_type"].isin(case["groups"])] + subset = subset.copy() + for column in ("sample", "cell_type"): + subset.obs[column] = subset.obs[column].cat.remove_unused_categories() + plot = sc.pl.dotplot( + subset, + case["genes"], + groupby="cell_type", + categories_order=case["groups"], + standard_scale={"gene": "var", "group": "group", None: None}[ + case["standard_scale"] + ], + mean_only_expressed=case["mean_only_expressed"], + return_fig=True, + show=False, + ) + plot.make_figure() + legend = plot.ax_dict["size_legend_ax"].get_xticklabels() + groups = plot.var_groups + out[name] = { + "rows": [str(v) for v in plot.dot_color_df.index], + "columns": [str(v) for v in plot.dot_color_df.columns], + "color": plot.dot_color_df.to_numpy().tolist(), + "fraction": plot.dot_size_df.to_numpy().tolist(), + "sizes": plot.ax_dict["mainplot_ax"].collections[0].get_sizes().tolist(), + "key": [text.get_text() for text in legend], + "gene_groups": None + if groups is None + else [list(groups.labels), [list(p) for p in groups.positions]], + } + out["versions"] = { + package: version(package) for package in ("scanpy", "anndata", "matplotlib") + } + with open(f"{directory}/out.json", "w") as handle: + json.dump(out, handle) + ''' +) + + +# Data. + + +def pca_data() -> dict[str, np.ndarray]: + """Counts of 12 samples, their log-CPM values, a matrix with tied variances, and + the cell numbers of a cube of the counts.""" + rng = np.random.default_rng(PCA_SEED) + n_genes = 300 + means = np.exp(rng.normal(4.0, 1.5, n_genes)) + effects = np.zeros((12, n_genes)) + for k, level in enumerate("abc"): + changed = rng.choice(n_genes, 40, replace=False) + rows = [i for i, value in enumerate(CONDITIONS) if value == level] + effects[np.ix_(rows, changed)] = rng.normal(0.0, 1.0 + 0.3 * k, 40) + effects[1::2, rng.choice(n_genes, 30, replace=False)] += 0.7 # batch y + mu = means * np.exp(effects) * rng.uniform(0.6, 1.6, (12, 1)) + counts = rng.negative_binomial(8.0, 8.0 / (8.0 + mu)).astype(np.float64) + values = log_cpm(counts, norm_factors=tmm_factors(counts)) + # 18 genes of distinct, large variance, then 4 of variance 12 (permutations of + # one profile, whose variances are equal to the last bit), then 6 small ones. + base = np.array([0.0, 1.0, 3.0, 4.0, 7.0, 9.0] * 2) + ties = np.vstack( + [rng.normal(0.0, 10.0 + j, 12) for j in range(18)] # fmt: skip + + [base[order] for order in (np.arange(12), np.arange(12)[::-1])] + + [base[np.roll(np.arange(12), shift)] for shift in (1, 5)] + + [rng.normal(0.0, 0.1, 12) for _ in range(6)] + ).T # samples x genes + n_cells = np.where(np.arange(12) < 10, 60, 5) + return {"counts": counts, "values": values, "ties": ties, "n_cells": n_cells} + + +def dot_data() -> dict[str, np.ndarray]: + """Cells of 4 types in 6 samples with their counts of 40 genes; sample D6 has no + pDC.""" + rng = np.random.default_rng(DOT_SEED) + n_genes = 40 + sizes = rng.integers(15, 120, (len(CELL_SAMPLES), len(GROUPS))) + sizes[-1, -1] = 0 + samples, groups = [], [] + for i, sample in enumerate(CELL_SAMPLES): + for g, group in enumerate(GROUPS): + samples += [sample] * int(sizes[i, g]) + groups += [group] * int(sizes[i, g]) + n_cells = len(groups) + base = np.exp(rng.normal(-1.0, 1.2, n_genes)) + level = np.ones((len(GROUPS), n_genes)) + for g in range(len(GROUPS)): # markers: genes 0-2 in B, 3-4 in T, and so on + level[g, rng.choice(n_genes, 6, replace=False)] *= rng.uniform(3.0, 12.0, 6) + level[0, :3] *= 20.0 + level[1, 3:5] *= 20.0 + code = np.array([GROUPS.index(group) for group in groups]) + depth = rng.uniform(0.5, 2.0, n_cells)[:, np.newaxis] + mu = base * level[code] * depth + counts = rng.negative_binomial(2.0, 2.0 / (2.0 + mu)).astype(np.float64) + counts[counts.sum(axis=1) == 0, 39] = 1.0 # every cell has counts + return { + "counts": counts, + "samples": np.array(samples), + "groups": np.array(groups), + "order": np.array(GROUPS), + "genes": np.array([f"g{j:02d}" for j in range(n_genes)]), + } + + +def dot_cube(data: dict[str, np.ndarray]) -> Cube: + return from_cells( + data["counts"], + samples=data["samples"], + groups=data["groups"], + sample_table={"sample_id": CELL_SAMPLES}, + group_table={"group_id": GROUPS}, + features=data["genes"].tolist(), + layers=["counts", "detected", "log_sum"], + ) + + +def pca_cube(data: dict[str, np.ndarray]) -> Cube: + return Cube( + {"counts": data["counts"][np.newaxis].astype(np.uint32)}, + samples={"sample_id": SAMPLES, "condition": CONDITIONS, "batch": BATCHES}, + groups={"group_id": ["T"]}, + features={"feature_id": [f"g{j:03d}" for j in range(data["counts"].shape[1])]}, + n_cells=data["n_cells"][:, np.newaxis], + ) + + +# The references, each run as a separate program. + + +def run_r(data: dict[str, np.ndarray]) -> dict[str, Any]: + with tempfile.TemporaryDirectory() as temporary: + root = Path(temporary) + for name in ("values", "ties", "counts"): + (root / f"{name}.bin").write_bytes( + np.ascontiguousarray(data[name], dtype="= 10).astype(np.float64) + (root / "used.bin").write_bytes(used.astype(" dict[str, dict[str, Any]]: + """Every case with each option, its default where the case does not set it.""" + defaults = { + "genes": GENES, + "samples": None, + "groups": None, + "standard_scale": None, + "mean_only_expressed": False, + } + return {name: defaults | case for name, case in DOT_CASES.items()} + + +def run_scanpy(data: dict[str, np.ndarray]) -> dict[str, Any]: + with tempfile.TemporaryDirectory() as temporary: + root = Path(temporary) + np.savez(root / "cells.npz", **data) + (root / "cases.json").write_text(json.dumps(dot_cases())) + (root / "program.py").write_text(SCANPY_PROGRAM) + python = os.environ.get("SCANPY_PYTHON", sys.executable) + subprocess.run([python, "-I", "program.py", str(root)], check=True, cwd=root) + result: dict[str, Any] = json.loads((root / "out.json").read_text()) + return result + + +# CytePool's figures, and the comparisons. + + +def axes() -> Any: + return Figure().subplots() + + +def artists(ax: Any, gid: str) -> list[Any]: + return [a for a in ax.get_children() if a.get_gid() == gid] + + +def scaled(ours: Any, reference: Any) -> float: + """Largest |ours - reference| over max(|reference|, 1); inf for unequal shapes.""" + ours, reference = np.asarray(ours, dtype=float), np.asarray(reference, dtype=float) + if ours.shape != reference.shape: + return np.inf + if ours.size == 0: + return 0.0 + return float(np.max(np.abs(ours - reference) / np.maximum(np.abs(reference), 1.0))) + + +def up_to_sign(ours: np.ndarray, theirs: np.ndarray) -> tuple[float, int]: + """The difference of two coordinate matrices after each column of ours takes the + sign that matches theirs, and the number of columns whose sign differed.""" + signs = np.sign(np.sum(ours * theirs, axis=0)) + return scaled(ours * signs, theirs), int(np.sum(signs < 0)) + + +def figure_samples(ax: Any) -> tuple[np.ndarray, list[tuple[float, ...]]]: + """The drawn coordinates of the samples, in their order (that of their labels), + and the colour of each.""" + positions = np.array([label.xy for label in artists(ax, "sample labels")]) + colours: list[tuple[float, ...]] = [()] * len(positions) + for collection in artists(ax, "samples"): + faces = collection.get_facecolors() + for offset, face in zip(collection.get_offsets(), faces, strict=True): + i = int(np.argmin(np.abs(positions - offset).sum(axis=1))) + colours[i] = tuple(face.tolist()) + return positions, colours + + +def partition(labels: list[Any]) -> set[frozenset[int]]: + """The samples grouped by their labels, whatever the labels are.""" + groups: dict[Any, set[int]] = {} + for i, label in enumerate(labels): + groups.setdefault(label, set()).add(i) + return {frozenset(members) for members in groups.values()} + + +def row(check: str, reference: str, value: float, tolerance: float) -> tuple: + return ( + check, + reference, + "max scaled difference", + f"{value:.1e}", + f"{tolerance:.0e}", + "pass" if value <= tolerance else "FAIL", + ) + + +def pca_rows(data: dict[str, np.ndarray], r: dict[str, Any]) -> list[tuple]: + rows = [] + flipped = 0 + worst = 0.0 + for case, (top, (first, second)) in PCA_CASES.items(): + result = principal_components(data["values"], top=top) + ours = result.coordinates[:, [first - 1, second - 1]] + theirs = np.column_stack([r[case]["x"], r[case]["y"]]) + difference, flips = up_to_sign(ours, theirs) + shares = result.variance_explained[[first - 1, second - 1]] + worst = max(worst, difference, scaled(shares, r[case]["percent"])) + flipped += flips + rows.append( + row( + "Log-CPM values of 12 samples, the 100 genes of largest variance " + "(components 1 and 2, 3 and 4) and 500, more than the 300 genes: " + "coordinates up to sign, shares of variance", + "plotPCA(returnData = TRUE)", + worst, + EXACT, + ) + ) + result = principal_components(data["ties"], top=20) + theirs = np.column_stack([r["ties"]["x"], r["ties"]["y"]]) + difference, flips = up_to_sign(result.coordinates[:, :2], theirs) + flipped += flips + keeps_first = result.genes.tolist()[-2:] == [18, 19] + rows.append( + row( + "Genes 19 to 22 of equal variance, the 20 of largest variance: the first " + "two of them kept; coordinates up to sign", + "plotPCA(returnData = TRUE)", + difference if keeps_first else np.inf, + EXACT, + ) + ) + cube = pca_cube(data) + result = principal_components(cube) + theirs = np.column_stack([r["cube"]["x"], r["cube"]["y"]]) + difference, flips = up_to_sign(result.coordinates[:, :2], theirs) + flipped += flips + same_samples = result.samples["sample_id"].tolist() == r["cube"]["name"] + rows.append( + row( + "A cube's group: the 10 samples with at least 10 cells, filterByExpr, TMM, " + "log-CPM with an average count of 2; coordinates up to sign, the samples", + "edgeR, plotPCA(returnData = TRUE)", + difference if same_samples else np.inf, + EXACT, + ) + ) + figure = r["figure"] + groups = [f"{c}:{b}" for c, b in zip(CONDITIONS, BATCHES, strict=True)] + ax = plotting.plot_pca( + data["values"], top=100, color_by=groups, labels=True, ax=axes() + ) + positions, colours = figure_samples(ax) + theirs = np.column_stack([figure["x"], figure["y"]]) + difference, _ = up_to_sign(positions, theirs) + labels = (ax.get_xlabel(), ax.get_ylabel()) == (figure["xlabel"], figure["ylabel"]) + equal_scales = ax.get_aspect() == float(figure["ratio"]) + same_groups = ( + partition(colours) + == partition(figure["colour"]) + == partition(r["top100"]["group"]) + ) + rows.append( + row( + "The figure: the points up to sign, the axis labels, equal scales, and the " + "samples' colours by condition and batch", + "plotPCA", + difference if (labels and equal_scales and same_groups) else np.inf, + EXACT, + ) + ) + two = principal_components(data["values"][:2], top=100) + pc1, _ = up_to_sign( + two.coordinates[:, :1], np.array(r["two_samples"]["x"])[:, None] + ) + noise = float(np.asarray(r["two_samples"]["percent"])[1]) + rows.append( + ( + "Two samples: PC1; DESeq2 returns a second component of rounding noise " + f"(share {noise:.1e}), CytePool one component", + "plotPCA(returnData = TRUE)", + "max scaled difference of PC1", + f"{pc1:.1e}", + "n/a", + "deliberate difference" + if two.n_components == 1 and pc1 <= EXACT + else "FAIL", + ) + ) + rows.append( + ( + "Components whose sign differs from R's (prcomp's signs are arbitrary)", + "plotPCA(returnData = TRUE)", + "count", + str(flipped), + "n/a", + "information", + ) + ) + return rows + + +def dot_rows(cube: Cube, scanpy: dict[str, Any]) -> list[tuple]: + rows = [] + for name, case in dot_cases().items(): + theirs = scanpy[name] + ax = plotting.plot_dotplot( + cube, + case["genes"], + groups=case["groups"], + samples=case["samples"], + standard_scale=case["standard_scale"], + mean_only_expressed=case["mean_only_expressed"], + ax=axes(), + ) + dots = artists(ax, "dots")[0] + sizes = dots.get_sizes() + fraction = np.array(theirs["fraction"]) + # The fractions that the areas encode, with scanpy's largest fraction shown. + f_max = np.ceil(fraction.max() * 10) / 10 + encoded = f_max * (sizes / LARGEST_DOT) ** (1.0 / SIZE_EXPONENT) + difference = max( + scaled(sizes, theirs["sizes"]), + scaled(encoded, fraction.ravel()), + scaled(np.asarray(dots.get_array()), np.ravel(theirs["color"])), + ) + names = ( + [t.get_text() for t in ax.get_yticklabels()] == theirs["rows"] + and [t.get_text() for t in ax.get_xticklabels()] == theirs["columns"] + and key_labels(ax) == theirs["key"] + ) + what = { + "plain": "fractions and means of 9 genes in 4 cell types, pooled over 6 " + "samples", + "scaled by gene": "means scaled by gene (standard_scale = 'var')", + "scaled by group": "means scaled by group (standard_scale = 'group')", + "expressing cells": "means over the cells that express the gene " + "(mean_only_expressed)", + "chosen samples and groups": "3 samples and 3 cell types, in a chosen order", + "genes in groups": "genes given in 3 groups", + }[name] + rows.append( + row( + f"Dot plot, {what}: the areas, the fractions they encode, the colour " + "values, the names of the rows and columns, the key of the areas", + "scanpy.pl.dotplot", + difference if names else np.inf, + EXACT, + ) + ) + if theirs["gene_groups"] is not None: + labels, positions = theirs["gene_groups"] + brackets = [line.get_xdata() for line in artists(ax, "gene groups")] + ours = [[round(x[0] + 0.3), round(x[-1] - 0.3)] for x in brackets] + texts = [t.get_text() for t in artists(ax, "gene group names")] + same = ours == positions and texts == labels + rows.append( + ( + "Dot plot, genes given in 3 groups: the columns of each bracket " + "and its name", + "scanpy.pl.dotplot", + "identical", + "yes" if same else "no", + "0", + "pass" if same else "FAIL", + ) + ) + return rows + + +def key_labels(ax: Any) -> list[str]: + """The fractions of the key of a dot plot's areas, as its legend writes them.""" + legends = [a.get_legend() for a in ax.figure.axes if a.get_legend() is not None] + return [text.get_text() for text in legends[0].get_texts()] + + +def write_fixtures( + pca: dict[str, np.ndarray], + r: dict[str, Any], + cube: Cube, + scanpy: dict[str, Any], +) -> None: + FIXTURE_DIRECTORY.mkdir(parents=True, exist_ok=True) + arrays = { + "values": pca["values"], + "ties": pca["ties"], + "counts": pca["counts"].astype(np.int32), + "n_cells": pca["n_cells"], + } + for case in (*PCA_CASES, "ties", "cube"): + arrays[f"{case} coordinates"] = np.column_stack([r[case]["x"], r[case]["y"]]) + arrays[f"{case} percent"] = np.array(r[case]["percent"], dtype=float) + figure = r["figure"] + arrays["figure coordinates"] = np.column_stack([figure["x"], figure["y"]]) + arrays["figure labels"] = np.array([figure["xlabel"], figure["ylabel"]]) + np.savez_compressed(FIXTURE_DIRECTORY / "pca.npz", **arrays) + arrays = {name: np.asarray(cube.layer(name)) for name in ("detected", "log_sum")} + arrays["n_cells"] = np.asarray(cube.n_cells) + for name in DOT_CASES: + for field in ("color", "fraction", "sizes"): + arrays[f"{name} {field}"] = np.array(scanpy[name][field], dtype=float) + np.savez_compressed(FIXTURE_DIRECTORY / "dot_plots.npz", **arrays) + texts = { + name: {k: scanpy[name][k] for k in ("rows", "columns", "key", "gene_groups")} + for name in DOT_CASES + } + (FIXTURE_DIRECTORY / "dot_plots.json").write_text(json.dumps(texts, indent=1)) + print(f"wrote {FIXTURE_DIRECTORY}/pca.npz, dot_plots.npz and dot_plots.json") + + +def main() -> int: + parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + parser.add_argument( + "--write-fixtures", action="store_true", help="write test fixtures" + ) + args = parser.parse_args() + pca = pca_data() + r = run_r(pca) + cells = dot_data() + cube = dot_cube(cells) + scanpy = run_scanpy(cells) + versions = list(r["versions"]) + [f"{k} {v}" for k, v in scanpy["versions"].items()] + print("Reference versions: " + ", ".join(versions) + "\n") + rows = pca_rows(pca, r) + dot_rows(cube, scanpy) + header = ["check", "reference", "metric", "value", "tolerance", "result"] + print("| " + " | ".join(header) + " |") + print("|" + "|".join(" --- " for _ in header) + "|") + for entry in rows: + print("| " + " | ".join(entry) + " |") + failed = any(entry[-1] == "FAIL" for entry in rows) + if args.write_fixtures: + write_fixtures(pca, r, cube, scanpy) + print("\nSOME CHECKS FAILED" if failed else "\nall checks passed") + return 1 if failed else 0 + + +if __name__ == "__main__": + sys.exit(main())