diff --git a/README.md b/README.md index 19fd6a9..38c7d18 100644 --- a/README.md +++ b/README.md @@ -35,6 +35,7 @@ what is planned. | Variance partitioning | `variance_fractions`: the fraction of every gene's variance that each variable of a design explains, as variancePartition's `fitExtractVarPartModel`: a mixed model of every gene by maximum likelihood (or REML), with the categorical variables as random intercepts, the continuous ones as fixed effects and precision weights such as voom's, or a linear model without random effects; `variance_partitioning`: the same on every cell type of a cube, with dream's precision weights, the samples each fit used and a provenance record | variancePartition, lme4 | | Meta-analysis | `meta_analysis`: for every gene, the fixed-effect (inverse-variance) or random-effects estimate of a common effect from the effects and standard errors of several studies, with the between-study variance by DerSimonian-Laird, REML or Paule-Mandel, Cochran's Q, I² and H², z-tests or the Knapp-Hartung adjustment, and genes measured in some studies only; `differential_expression_meta_analysis`: the same for the `limma_voom` or `limma_trend` results of several cubes, matched by cell type, contrast and gene id, with p-values adjusted over the genes and a provenance record | metafor | | Combined p-values and replication | `combine_p_values`: every gene's p-values combined over the studies that tested it, by Fisher's, Stouffer's (weighted), Tippett's, Wilkinson's, Pearson's, Mudholkar and George's or Simes's method, one-sided or two-sided with the directions of the effects; `partial_conjunction`: the p-value that at least u studies have an effect, in one direction (Benjamini and Heller); `proportion_of_true_nulls`: pi0 and pi1 as qvalue's `pi0est`, with its smoother and bootstrap; `differential_expression_replication`: for the studies of a `CubeCollection`, every gene's combined and replication p-values, the studies that find it, its sign agreement and heterogeneity, pi1 and sign agreement of every pair of studies, and the genes found in at least k studies, with their figures | scipy, metap, qvalue, metafor, mpmath | +| 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 | | 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 | @@ -57,8 +58,7 @@ These are not implemented yet: - Dose-response fits, benchmark doses and synergy - Meta-regression (meta-analysis with study-level covariates) and multivariate meta-analysis of correlated effects -- Cell types aligned up the Cell Ontology, and held-out prediction of case status from - the genes of several studies +- Cell types aligned up the Cell Ontology ## Known limits diff --git a/THIRD_PARTY_NOTICES.md b/THIRD_PARTY_NOTICES.md index cc5679a..7e36268 100644 --- a/THIRD_PARTY_NOTICES.md +++ b/THIRD_PARTY_NOTICES.md @@ -167,6 +167,26 @@ 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. ``` +### SumRank + +`cytepool/sumrank.py` and `cytepool/held_out_prediction.py` implement SumRank and its +held-out prediction from Nakatsuka et al. (2025), following conventions of the authors' +code (https://github.com/nathan-nakatsuka/scRNA_Reproducibility, commit 46392f66; +`SumRank.R`, `CalibratePValueswithPermutations.R`, `PermutingCaseControlStatus.R`, +`UCell_AUC_code.R`): ties ranked by their average rank, floor(p k) of the studies summed +and the cap at half their number, a dataset without the cell type left out, the null +pooled over genes and permutations, and, for held-out scores, each cell type's scores +scaled to [0, 1], a missing cell type at their mean and the scores of genes up minus +those of genes down; `cytepool/plotting.py` draws `plot_sumrank` as their +`MakeManhattanPlot` (`MakingManhattanPlots.R`) draws it. No code is copied. +`validation/sumrank.py` downloads `SumRank.R`, `CalibratePValueswithPermutations.R` and +`MakingManhattanPlots.R` of that commit and runs them. Nakatsuka et al. (2025) is: +Nakatsuka N, Adler D, Jiang L, Hartman A, Cheng E, Klann E, Satija R. 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. The code is distributed under the Apache License, +Version 2.0 (Copyright 2024 Nathan Nakatsuka; http://www.apache.org/licenses/LICENSE-2.0). + ### numba `cytepool/_jit.py` places the cache of its compiled routines with subclasses of numba's diff --git a/cytepool/_study_statistics.py b/cytepool/_study_statistics.py new file mode 100644 index 0000000..641b881 --- /dev/null +++ b/cytepool/_study_statistics.py @@ -0,0 +1,161 @@ +"""The results of the studies of a collection for one group and contrast, as arrays, +and the number of studies whose ranks SumRank sums. + +SumRank and held-out prediction read, from every study's differential expression of +one group and contrast, a signed statistic of every gene of the collection, with the +log2 fold changes and standard errors of a meta-analysis: one row per study and one +column per gene, NaN where the study did not test the gene. They are read from +:func:`~cytepool.meta_analysis.align_study_results`, which checks that each result +analysed its study's aligned cube and says why a study has no result for a gene. The +signed statistic is -log10 of the p-value, negative where the log2 fold change is +negative, as Nakatsuka et al. (2025) rank DESeq2's results, and a proportion p of k +studies sums floor(p k) ranks, as their code does. + +References: + Nakatsuka N, Adler D, Jiang L, Hartman A, Cheng E, Klann E, Satija R (2025). + Improving reproducibility of differentially expressed genes in single-cell + transcriptomic studies of neurodegenerative diseases through meta-analysis. + Nature Communications 16, 7436. doi:10.1038/s41467-025-62579-z +""" + +import math +from collections.abc import Mapping +from dataclasses import dataclass +from typing import Any + +import numpy as np +import numpy.typing as npt + +from ._masked_arrays import reject_masked +from ._options import real_number +from .cube import CubeCollection +from .differential_expression import DifferentialExpressionResult +from .meta_analysis import align_study_results + +FloatArray = npt.NDArray[np.float64] +BoolArray = npt.NDArray[np.bool_] +IntArray = npt.NDArray[np.int64] + + +@dataclass(frozen=True) +class StudyStatistics: + """Every study's results of one group and contrast, ``(n_studies, n_features)`` + in the collection's order. + + Attributes: + statistics: The signed -log10 p-value of every gene, NaN where the study did + not test it or its test has no p-value. + log2_fold_change: The log2 fold changes, NaN likewise. + standard_error: The standard errors, NaN likewise, or where the results have + none (dream's). + measured: Whether each study measured each gene. + analysed: Whether each study's results test the group and contrast. + reasons: Why each study does not, as its result says; empty when it does. + n_tested: The number of genes each study tested. + """ + + statistics: FloatArray + log2_fold_change: FloatArray + standard_error: FloatArray + measured: BoolArray + analysed: BoolArray + reasons: list[str] + n_tested: IntArray + + +def signed_statistics( + p_value: npt.ArrayLike, log2_fold_change: npt.ArrayLike +) -> FloatArray: + """-log10 of the p-values, negative where the log2 fold change is negative, as + the authors of SumRank sign DESeq2's p-values; infinite for a p-value of 0, and + NaN where either value is NaN.""" + p = np.asarray(p_value, dtype=np.float64) + change = np.asarray(log2_fold_change, dtype=np.float64) + with np.errstate(divide="ignore", invalid="ignore"): + magnitude = -np.log10(p) + signed: FloatArray = np.where(change < 0.0, -magnitude, magnitude) + signed[np.isnan(change)] = np.nan + return signed + + +def study_statistics( + collection: CubeCollection, + results: Mapping[str, DifferentialExpressionResult], + group: str, + contrast: str, +) -> StudyStatistics: + """Every study's signed statistics, log2 fold changes and standard errors of one + group and contrast, with which studies analysed them and why the others did not. + + Raises: + TypeError, ValueError, KeyError: As + :func:`~cytepool.meta_analysis.align_study_results` raises them, and + ValueError if ``contrast`` is an F-test, which has no direction. + """ + table = align_study_results(collection, results, group, contrast) + shape = (collection.n_studies, collection.n_features) + if "log2_fold_change" not in table: + raise ValueError( + f"{contrast!r} is an F-test, which tests several contrasts at once and " + "has no direction; SumRank ranks the genes of one contrast, up and down" + ) + status = table["status"].reshape(shape) + analysed = ((status == "tested") | (status == "not_tested")).any(axis=1) + reasons = [] + for row, flag in enumerate(analysed.tolist()): + if flag: + reasons.append("") + continue + found = table["reason"].reshape(shape)[row] + # A study that measured some gene says why it has no table; one that measured + # none of the collection's genes cannot be in a collection. + reasons.append( + next( + ( + text + for text, kind in zip( + found.tolist(), status[row].tolist(), strict=True + ) + if kind == "not_analysed" + ), + found[0], + ) + ) + tested = status == "tested" + + def column(name: str) -> FloatArray: + if name not in table: + return np.full(shape, np.nan) + values: FloatArray = np.array(table[name], dtype=np.float64).reshape(shape) + values[~tested] = np.nan + return values + + change = column("log2_fold_change") + statistics = signed_statistics(column("p_value"), change) + return StudyStatistics( + statistics=statistics, + log2_fold_change=change, + standard_error=column("standard_error"), + measured=np.array(collection.measured, dtype=bool), + analysed=analysed, + reasons=reasons, + n_tested=tested.sum(axis=1).astype(np.int64), + ) + + +def studies_summed(proportion: Any, n_studies: int) -> int: + """The number of ranks SumRank sums for each gene, m = floor(p k), at least 1; the + product is rounded to 9 decimals first, so that 0.57 of 100 studies is 57. + + Raises: + TypeError: If ``proportion`` is not a real number, or is a numpy masked array. + ValueError: If it is not in (0, 1]. + """ + reject_masked(proportion, "proportion_of_studies") + p = real_number(proportion, "proportion_of_studies", finite=True) + if not 0.0 < p <= 1.0: + raise ValueError( + f"proportion_of_studies must be in (0, 1], such as 1 to sum every study, " + f"got {p}" + ) + return max(1, math.floor(round(p * n_studies, 9))) diff --git a/cytepool/held_out_prediction.py b/cytepool/held_out_prediction.py new file mode 100644 index 0000000..6d874d3 --- /dev/null +++ b/cytepool/held_out_prediction.py @@ -0,0 +1,702 @@ +"""Prediction of case status in held-out studies, from genes chosen in the others. + +Genes found by combining studies should separate cases from controls in a study that +had no part in finding them. :func:`held_out_prediction` measures this as the paper +that introduced SumRank does (Nakatsuka et al. 2025): each study of a collection is +left out in turn; genes up and genes down in each cell type are chosen from the other +studies, by SumRank (:func:`~cytepool.sumrank.sumrank`) and, for comparison, by a +random-effects meta-analysis (:func:`~cytepool.meta_analysis.meta_analysis`); every +sample of the held-out study is scored by how high its profiles rank those genes; and +the area under the ROC curve (AUC) of the score says how well it separates the +study's cases from its controls. :func:`~cytepool.plotting.plot_held_out_auc` draws +the AUC of every held-out study by method. + +Choosing the genes: + In every group (cell type), from the studies other than the held-out one whose + results test the group and contrast, and among the genes that every one of them + tested, ``n_genes`` genes up and ``n_genes`` down, the same number for each + method, as the paper compares methods by gene sets of one size: + + - ``"sumrank"``: the genes with the smallest ``up_statistic``, and those with the + smallest ``down_statistic``, below the cap, with the given proportion of the + studies summed; ties in the collection's order of the genes. Every gene has + the same numbers of studies and ranks, so this is also the order of the + permutation p-values, which therefore need no permutations here. + - ``"meta_analysis"``: the genes with a positive estimate and the smallest + p-values, and those with a negative one, from the random-effects + meta-analysis of the studies' log2 fold changes and standard errors with REML + and the Knapp-Hartung adjustment (``method="reml"``, + ``test="knapp_hartung"``), as the paper ran metagen (``method.tau = "REML"``, + ``hakn = TRUE``); ties in the collection's order. + +Scoring a sample: + The genes of a sample's profile in a group, its counts summed over the cells, + are ranked from the highest count, tied genes sharing their average rank, among + all the genes the held-out study measured. A set of n of these genes, those of a + chosen set that the study measured, has UCell's score (Andreatta and Carmona + 2021), 1 - U / (n G), with U the sum of the set's ranks minus n (n + 1) / 2 and + G the number of genes ranked: the maximum rank is the number of genes, as the + paper set it to rank every gene (``maxRank = 16000``). Then, as the authors' + code does: + + 1. In each group, the scores of the genes up, and those of the genes down, are + scaled over the held-out study's samples to [0, 1], minus their smallest + value and divided by their range (0 for all when they are equal), which + makes the score's own scale irrelevant. + 2. A sample without a profile in the group, that is not among the samples of + the held-out study's fit of the group, takes the mean of the other samples' + scaled scores, and every sample takes 0 when none has one, as for a set + without a gene the study measured. + 3. A sample's score is the sum over the groups of its scaled scores of the + genes up, minus the sum of its scaled scores of the genes down. + +The AUC: + The probability that a case scores higher than a control, ties counting one + half: the Mann-Whitney statistic divided by the number of pairs of a case and a + control (Hanley and McNeil 1982). The samples are those of the held-out study's + design whose value of ``variable`` is ``case`` or ``control``. 0.5 says the + score does not separate them, and below 0.5 the genes predict the wrong + direction; the AUC is NaN, with the reason, for a study without a case or + without a control. + +Deliberate differences from the paper: + - The paper fits a logistic regression of case status on the score in the other + studies and takes the AUC of its predictions in the held-out study from + pROC's ``auc``, whose default ``direction = "auto"`` compares the medians of + the held-out cases and controls and takes the direction in which the cases + score higher (Robin et al. 2011). A logistic regression on one score is + monotone in it, so that AUC is the score's in the direction pROC chooses, + which puts a score that ranks controls above cases above 0.5 too. Here the + cases are expected to score higher, as genes up minus genes down should make + them, so the logistic regression is not needed. + - The profiles are the summed counts, which rank the genes of a sample as their + mean over its cells does; the paper ranked Seurat's average of the cells' + normalized expression, which weighs each cell's counts by its depth. + +Not done: + The proportion of studies is given, not chosen: the paper chose it as the one + whose genes predict best, which a loop over proportions gives. The paper's + prediction of disease severity in cases (a linear regression on the score) is + not done. + +References: + Nakatsuka N, Adler D, Jiang L, Hartman A, Cheng E, Klann E, Satija R (2025). + Improving reproducibility of differentially expressed genes in single-cell + transcriptomic studies of neurodegenerative diseases through meta-analysis. + Nature Communications 16, 7436. doi:10.1038/s41467-025-62579-z + + Andreatta M, Carmona SJ (2021). UCell: robust and scalable single-cell gene + signature scoring. Computational and Structural Biotechnology Journal 19, + 3796-3798. doi:10.1016/j.csbj.2021.06.043 + + Hanley JA, McNeil BJ (1982). The meaning and use of the area under a receiver + operating characteristic (ROC) curve. Radiology 143, 29-36. + doi:10.1148/radiology.143.1.7063747 + + Robin X, Turck N, Hainard A, Tiberti N, Lisacek F, Sanchez JC, Mueller M (2011). + pROC: an open-source package for R and S+ to analyze and compare ROC curves. BMC + Bioinformatics 12, 77. doi:10.1186/1471-2105-12-77 +""" + +import json +from collections.abc import Mapping, Sequence +from dataclasses import dataclass +from typing import Any + +import numpy as np +import numpy.typing as npt + +from ._analysis_records import software_versions +from ._columns import text_column +from ._masked_arrays import reject_masked +from ._messages import listing, quoted +from ._options import index_integer +from ._ranks import average_ranks +from ._study_statistics import StudyStatistics, study_statistics, studies_summed +from .cube import CubeCollection +from .differential_expression import DifferentialExpressionResult +from .meta_analysis import meta_analysis +from .sumrank import sumrank +from .table import Table + +__all__ = ["HeldOutPrediction", "held_out_prediction"] + +FloatArray = npt.NDArray[np.float64] +IntArray = npt.NDArray[np.intp] + +METHODS = ("sumrank", "meta_analysis") +"""The methods that choose the genes, in the order the tables list them.""" + +_DIRECTIONS = ("up", "down") +# The meta-analysis that chooses genes for comparison: the paper's metagen settings. +_META_ANALYSIS = {"method": "reml", "test": "knapp_hartung"} + + +class HeldOutPrediction: + """The prediction of case status in every held-out study, by method. + + Built by :func:`held_out_prediction`. Three tables, the studies in the + collection's order and the methods in the order ``"sumrank"``, + ``"meta_analysis"``: + + - ``auc``: one row per held-out study and method: ``study``, ``method``, + ``auc`` (NaN where it cannot be computed), ``n_cases``, ``n_controls``, + ``n_groups`` (the groups whose genes scored the study), ``n_genes_up`` and + ``n_genes_down`` (the genes chosen, summed over the groups, that the study + measured) and ``reason`` (why the AUC is NaN; empty otherwise); + - ``scores``: one row per held-out study, method and sample of its design that + is a case or a control: ``study``, ``method``, ``sample_id``, ``case`` + (whether it is a case) and ``score``; + - ``genes``: the genes chosen for each held-out study from the other studies: + ``study`` (held out), ``method``, ``group_id``, ``direction`` (``"up"`` or + ``"down"``), ``rank`` (from 1) and ``feature_id``. + + Immutable. + + Raises: + TypeError: If a table is not a ``Table``, or ``provenance`` not a mapping + that JSON can write. + ValueError: If ``provenance`` holds a NaN or an infinity. + """ + + __slots__ = ("_auc", "_scores", "_genes", "_record") + + _auc: Table + _scores: Table + _genes: Table + _record: str + + def __init__( + self, auc: Table, scores: Table, genes: Table, provenance: Mapping[str, Any] + ) -> None: + for name, given in (("auc", auc), ("scores", scores), ("genes", genes)): + if not isinstance(given, Table): + raise TypeError(f"{name} must be a Table, got {type(given).__name__}") + if not isinstance(provenance, Mapping): + raise TypeError( + f"provenance must be a mapping, got {type(provenance).__name__}" + ) + state = { + "_auc": auc, + "_scores": scores, + "_genes": genes, + # JSON holds no arrays: json.dumps raises TypeError for one. + "_record": json.dumps(provenance, allow_nan=False), + } + for name, value in state.items(): + object.__setattr__(self, name, value) + + def __setattr__(self, name: str, value: object) -> None: + raise AttributeError("HeldOutPrediction is immutable") + + def __delattr__(self, name: str) -> None: + raise AttributeError("HeldOutPrediction is immutable") + + def __reduce__(self) -> tuple[Any, tuple[Any, ...]]: + return ( + HeldOutPrediction, + (self._auc, self._scores, self._genes, json.loads(self._record)), + ) + + @property + def auc(self) -> Table: + """The AUC of every held-out study and method.""" + return self._auc + + @property + def scores(self) -> Table: + """The score of every case and control of every held-out study, by method.""" + return self._scores + + @property + def genes(self) -> Table: + """The genes chosen from the other studies for every held-out study.""" + return self._genes + + @property + def settings(self) -> dict[str, Any]: + """The settings, as a new dict (also in the provenance).""" + settings: dict[str, Any] = json.loads(self._record)["settings"] + return settings + + @property + def provenance(self) -> dict[str, Any]: + """A JSON-serialisable record of what produced every number, as a new dict. + + It holds the CytePool, numpy, numba and Python versions, the settings, and + the provenance of every study's differential-expression result, by study + name, which identifies its cube by its fingerprint and gives its design and + the samples of every group's fit. + """ + record: dict[str, Any] = json.loads(self._record) + return record + + def __repr__(self) -> str: + table = self._auc + settings = self.settings + lines = [ + f"HeldOutPrediction of {settings['case']!r} against {settings['control']!r}" + f" ({settings['variable']}), contrast {settings['contrast']!r}: " + f"{len(set(table['study'].tolist()))} held-out studies, " + f"{settings['n_genes']} genes up and down per group" + ] + for method in settings["methods"]: + values = table["auc"][(table["method"] == method) & ~np.isnan(table["auc"])] + mean = f"mean AUC {values.mean():.3f}" if values.size else "no AUC" + lines.append(f" {method}: {mean} over {values.size} studies") + return "\n".join(lines) + + +def held_out_prediction( + collection: CubeCollection, + results: Mapping[str, DifferentialExpressionResult], + contrast: str, + variable: str, + case: str, + control: str, + *, + groups: str | Sequence[str] | None = None, + n_genes: int = 50, + proportion_of_studies: float = 1.0, +) -> HeldOutPrediction: + """Leave each study out, choose genes from the others, and predict its cases. + + For every study of the collection in turn, genes up and down in each group are + chosen from the results of the other studies, by SumRank and by a + random-effects meta-analysis; each case and control of the held-out study is + scored by how high its profiles rank those genes, as UCell scores them; and the + AUC says how well the scores separate cases from controls. The module + description gives each step. The results are those that + :func:`~cytepool.sumrank.differential_expression_sumrank` takes: each study's + analysis of its aligned cube. + + Args: + collection: The collection whose aligned cubes the studies' results + analysed. + results: The differential-expression result of every study of the + collection, by study name, from + :func:`~cytepool.differential_expression.limma_voom` or + :func:`~cytepool.differential_expression.limma_trend`, whose tables + hold the standard errors of the meta-analysis. + contrast: The contrast that the studies tested, cases against controls. + variable: The variable of the studies' designs that tells cases from + controls, such as ``"condition"``. + case: Its value for a case, such as ``"AD"``. + control: Its value for a control. + groups: A group id or a list of them, whose scores are summed; by default + every group of the collection. A group that no other study tested gives + a held-out study no genes. + n_genes: The number of genes chosen up, and down, in each group, by each + method; at least 1. + proportion_of_studies: The proportion of the studies summed by SumRank, as + for :func:`~cytepool.sumrank.sumrank`. + + Returns: + The AUC of every held-out study and method, the score of every sample, and + the genes chosen. + + Raises: + TypeError: If an argument has the wrong type, or is or holds a numpy masked + array. + ValueError: If ``results`` does not give the result of every study of the + collection or a result did not analyse its study's aligned cube, the + contrast is an F-test, the results have no standard errors (dream's), a + study's design has no variable ``variable``, ``case`` equals + ``control``, a group is listed twice, or an option is invalid. + KeyError: If a group is not a group of the collection, or no study tested + ``contrast``. + """ + if not isinstance(collection, CubeCollection): + raise TypeError( + f"collection must be a CubeCollection, got {type(collection).__name__}" + ) + for name, value in ( + ("contrast", contrast), + ("variable", variable), + ("case", case), + ("control", control), + ): + reject_masked(value, name) + if not isinstance(value, str): + raise TypeError(f"{name} must be a string, got {type(value).__name__}") + if case == control: + raise ValueError( + f"case and control are both {case!r}; give the two values of {variable!r} " + "that the AUC compares" + ) + chosen = _groups(collection, groups) + count = _n_genes(n_genes) + studies_summed(proportion_of_studies, collection.n_studies) # checks it + # Every group's results, read once: the held-out studies take rows of them. + inputs = { + group: study_statistics(collection, results, group, contrast) + for group in chosen + } + for group, found in inputs.items(): + missing = found.analysed & np.isnan(found.standard_error).all(axis=1) + if missing.any(): + study = collection.studies[int(np.flatnonzero(missing)[0])] + raise ValueError( + f"the results of study {study!r} have no standard errors for group " + f"{group!r}, which the meta-analysis reads: dream_voom's tables have " + "none; give results of limma_voom or limma_trend" + ) + studies = list(collection.studies) + rows: list[dict[str, Any]] = [] + score_rows: list[dict[str, Any]] = [] + gene_rows: list[dict[str, Any]] = [] + feature_ids = collection.features["feature_id"] + for held, study in enumerate(studies): + samples = _held_out_samples(results[study], study, variable, case, control) + profiles = _profiles(collection, results[study], study, chosen) + for method in METHODS: + sets = {} + for group in chosen: + found = inputs[group] + training = found.analysed.copy() + training[held] = False + if not training.any(): + continue + chosen_genes = _chosen_genes( + method, found, training, count, proportion_of_studies + ) + sets[group] = chosen_genes + for direction, positions in zip(_DIRECTIONS, chosen_genes, strict=True): + for rank, position in enumerate(positions.tolist(), start=1): + gene_rows.append( + { + "study": study, + "method": method, + "group_id": group, + "direction": direction, + "rank": rank, + "feature_id": feature_ids[position], + } + ) + scores, summary = _scores(samples, profiles, sets, collection) + row = {"study": study, "method": method, **summary} + row |= _auc(scores, samples.case) + rows.append(row) + for sample, is_case, score in zip( + samples.ids, samples.case.tolist(), scores.tolist(), strict=True + ): + score_rows.append( + { + "study": study, + "method": method, + "sample_id": sample, + "case": is_case, + "score": score, + } + ) + settings = { + "contrast": contrast, + "variable": variable, + "case": case, + "control": control, + "groups": chosen, + "n_genes": count, + "proportion_of_studies": float(proportion_of_studies), + "methods": list(METHODS), + "meta_analysis": dict(_META_ANALYSIS), + "score": "UCell, maxRank = the number of genes ranked", + } + return HeldOutPrediction( + _auc_table(rows), + Table( + { + "study": text_column([r["study"] for r in score_rows]), + "method": text_column([r["method"] for r in score_rows]), + "sample_id": text_column([r["sample_id"] for r in score_rows]), + "case": np.array([r["case"] for r in score_rows], dtype=bool), + "score": np.array([r["score"] for r in score_rows], dtype=np.float64), + } + ), + Table( + { + "study": text_column([r["study"] for r in gene_rows]), + "method": text_column([r["method"] for r in gene_rows]), + "group_id": text_column([r["group_id"] for r in gene_rows]), + "direction": text_column([r["direction"] for r in gene_rows]), + "rank": np.array([r["rank"] for r in gene_rows], dtype=np.int64), + "feature_id": text_column([r["feature_id"] for r in gene_rows]), + } + ), + { + "analysis": "held_out_prediction", + "versions": software_versions(), + "settings": settings, + "studies": {name: results[name].provenance for name in studies}, + }, + ) + + +# Arguments. + + +def _groups(collection: CubeCollection, value: Any) -> list[str]: + reject_masked(value, "groups") + available: list[str] = collection.groups["group_id"].tolist() + if value is None: + return available + if isinstance(value, str): + listed: list[Any] = [value] + elif isinstance(value, Sequence | np.ndarray): + listed = list(value) + else: + raise TypeError( + f"groups must be a group id or a list of them, got {type(value).__name__}" + ) + if not listed: + raise ValueError("groups is empty; give at least one, or None for every group") + for group in listed: + if not isinstance(group, str): + raise TypeError( + f"groups must be group ids (str), got {type(group).__name__} {group!r}" + ) + repeated = [g for g in dict.fromkeys(listed) if listed.count(g) > 1] + if repeated: + raise ValueError(f"groups lists {quoted(repeated)} more than once") + unknown = [g for g in listed if g not in available] + if unknown: + raise KeyError( + f"no group {quoted(unknown)} in the collection; its groups are: " + f"{listing(available)}" + ) + return listed + + +def _n_genes(value: Any) -> int: + if isinstance(value, bool | np.bool_): + raise TypeError("n_genes must be a whole number of genes, not bool") + number = index_integer(value, "n_genes", expected="a whole number of genes") + if number < 1: + raise ValueError(f"n_genes must be at least 1, got {number}") + return number + + +# The genes chosen from the other studies. + + +def _chosen_genes( + method: str, + found: StudyStatistics, + training: npt.NDArray[np.bool_], + count: int, + proportion: float, +) -> tuple[IntArray, IntArray]: + """The positions, in the collection's features, of the genes chosen up and down + from the training studies, best first.""" + statistics = found.statistics[training] + universe = ~np.isnan(statistics).any(axis=0) + if np.count_nonzero(universe) < 2: + empty = np.empty(0, dtype=np.intp) + return empty, empty + chosen: list[IntArray] = [] + if method == "sumrank": + ranking = sumrank(statistics, proportion_of_studies=proportion) + cap = ranking.n_studies_summed / 2.0 + for statistic in (ranking.up_statistic, ranking.down_statistic): + candidates = np.flatnonzero(ranking.ranked & (statistic < cap)) + order = np.lexsort((candidates, statistic[candidates])) + chosen.append(candidates[order][:count].astype(np.intp)) + return chosen[0], chosen[1] + effects = found.log2_fold_change[training][:, universe] + errors = found.standard_error[training][:, universe] + combined = meta_analysis(effects, errors, method="reml", test="knapp_hartung") + positions = np.flatnonzero(universe) + for sign in (1.0, -1.0): + candidates = np.flatnonzero(sign * combined.estimate > 0.0) + order = np.lexsort((candidates, combined.p_value[candidates])) + chosen.append(positions[candidates[order][:count]].astype(np.intp)) + return chosen[0], chosen[1] + + +# The held-out study's samples, their profiles and their scores. + + +@dataclass(frozen=True) +class _Samples: + """The cases and controls of a held-out study, in its cube's order.""" + + ids: list[str] + case: npt.NDArray[np.bool_] + + +def _held_out_samples( + result: DifferentialExpressionResult, + study: str, + variable: str, + case: str, + control: str, +) -> _Samples: + design = result.provenance["design"] + samples = design["samples"] + if variable not in samples or variable == "sample_id": + variables = [name for name in samples if name != "sample_id"] + raise ValueError( + f"variable names {variable!r}, which is not a variable of the design of " + f"study {study!r}; its variables are: {quoted(variables) or '(none)'}" + ) + labels = { + sample: str(value) + for sample, value in zip(samples["sample_id"], samples[variable], strict=True) + } + # The design's samples in the cube's order, as the result's sample table has them. + table = result.samples + first = table["group_id"] == table["group_id"][0] if len(table) else [] + ordered = table["sample_id"][first].tolist() if len(table) else list(labels) + kept = [sample for sample in ordered if labels[sample] in (case, control)] + return _Samples( + ids=kept, case=np.array([labels[s] == case for s in kept], dtype=bool) + ) + + +@dataclass(frozen=True) +class _Profiles: + """A held-out study's profiles in one group: the counts of the samples of its fit + of the group, sorted within each sample for ranking.""" + + samples: list[str] + counts: FloatArray # samples x genes measured + sorted_counts: FloatArray + features: dict[str, int] # feature id: column + + +def _profiles( + collection: CubeCollection, + result: DifferentialExpressionResult, + study: str, + groups: list[str], +) -> dict[str, _Profiles]: + """The profiles of the held-out study in each group that its result analysed, + for the samples of its fit of the group.""" + table = result.samples + if not len(table): + return {} + cube = collection.aligned_cube(study) + sample_rows = {s: i for i, s in enumerate(cube.samples["sample_id"].tolist())} + features = {f: i for i, f in enumerate(cube.features["feature_id"].tolist())} + profiles = {} + analysed = set(result.groups) + for group in groups: + if group not in analysed: + continue + in_group = (table["group_id"] == group) & table["used"].astype(bool) + used = table["sample_id"][in_group].tolist() + if not used: + continue + rows = np.array([sample_rows[s] for s in used], dtype=np.intp) + counts = np.asarray(cube.layer("counts", group), dtype=np.float64)[rows] + profiles[group] = _Profiles( + samples=used, + counts=counts, + sorted_counts=np.sort(counts, axis=1), + features=features, + ) + return profiles + + +def _set_scores(profiles: _Profiles, genes: list[str]) -> FloatArray | None: + """UCell's score of a set of genes in every profile, with the maximum rank the + number of genes ranked; None when the study measured none of the genes.""" + columns = np.array( + [profiles.features[g] for g in genes if g in profiles.features], dtype=np.intp + ) + if columns.size == 0: + return None + n_genes = profiles.counts.shape[1] + n = columns.size + scores = np.empty(len(profiles.samples)) + for row in range(len(profiles.samples)): + ordered = profiles.sorted_counts[row] + values = profiles.counts[row, columns] + above = np.searchsorted(ordered, values, side="right") + below = np.searchsorted(ordered, values, side="left") + # The rank from the highest count: the genes above, plus the mean position + # among the equal ones. + ranks = (n_genes - above) + 0.5 * (above - below + 1) + u = ranks.sum() - 0.5 * n * (n + 1) + scores[row] = 1.0 - u / (n * n_genes) + return scores + + +def _scaled(scores: FloatArray) -> FloatArray: + """Scores minus their smallest value, divided by their range; 0 when equal.""" + low, high = scores.min(), scores.max() + if not high > low: + return np.zeros(scores.size) + scaled: FloatArray = (scores - low) / (high - low) + return scaled + + +def _scores( + samples: _Samples, + profiles: dict[str, _Profiles], + sets: dict[str, tuple[IntArray, IntArray]], + collection: CubeCollection, +) -> tuple[FloatArray, dict[str, Any]]: + """The score of every case and control: over the groups, the scaled scores of + the genes up minus those of the genes down, with the mean for a sample without + a profile in a group.""" + feature_ids = collection.features["feature_id"] + position = {sample: i for i, sample in enumerate(samples.ids)} + total = np.zeros(len(samples.ids)) + n_groups = 0 + measured = {"up": 0, "down": 0} + for group, chosen in sets.items(): + found = profiles.get(group) + contributed = False + for direction, sign, positions in zip( + _DIRECTIONS, (1.0, -1.0), chosen, strict=True + ): + genes = [feature_ids[p] for p in positions.tolist()] + scores = None if found is None else _set_scores(found, genes) + if scores is None: + continue # every sample takes 0 + assert found is not None + measured[direction] += sum(g in found.features for g in genes) + scaled = _scaled(scores) + values = np.full(len(samples.ids), scaled.mean()) + for sample, value in zip(found.samples, scaled.tolist(), strict=True): + if sample in position: + values[position[sample]] = value + total += sign * values + contributed = True + n_groups += contributed + return total, { + "n_groups": n_groups, + "n_genes_up": measured["up"], + "n_genes_down": measured["down"], + } + + +def _auc(scores: FloatArray, case: npt.NDArray[np.bool_]) -> dict[str, Any]: + """The AUC of the scores for cases against controls, and the numbers of each.""" + n_cases = int(np.count_nonzero(case)) + n_controls = int(case.size - n_cases) + row: dict[str, Any] = {"n_cases": n_cases, "n_controls": n_controls} + if n_cases == 0 or n_controls == 0: + row["auc"] = np.nan + row["reason"] = ( + f"the study has {n_cases} cases and {n_controls} controls; the AUC needs " + "at least one of each" + ) + return row + ranks = average_ranks(np.ascontiguousarray(scores)) + statistic = ranks[case].sum() - 0.5 * n_cases * (n_cases + 1) + row["auc"] = float(statistic / (n_cases * n_controls)) + row["reason"] = "" + return row + + +def _auc_table(rows: list[dict[str, Any]]) -> Table: + columns: dict[str, Any] = { + "study": text_column([r["study"] for r in rows]), + "method": text_column([r["method"] for r in rows]), + "auc": np.array([r["auc"] for r in rows], dtype=np.float64), + } + for name in ("n_cases", "n_controls", "n_groups", "n_genes_up", "n_genes_down"): + columns[name] = np.array([r[name] for r in rows], dtype=np.int64) + columns["reason"] = text_column([r["reason"] for r in rows]) + return Table(columns) diff --git a/cytepool/plotting.py b/cytepool/plotting.py index 2bd8b0a..91f8468 100644 --- a/cytepool/plotting.py +++ b/cytepool/plotting.py @@ -50,6 +50,10 @@ - Replication across studies: :func:`plot_study_agreement`, pi1 or the sign agreement of every pair of studies (Storey and Tibshirani 2003); :func:`plot_replicated_genes`, the genes that each number of studies finds (Nakatsuka et al. 2025). +- Several studies: :func:`plot_sumrank`, the SumRank p-value of every gene, the genes + in alphabetical order, as the authors of SumRank draw it, and + :func:`plot_held_out_auc`, the AUC of every held-out study by the method that + chose its genes (Nakatsuka et al. 2025). - Variance partitioning: :func:`plot_variance_fractions`, the violins (Hintze and Nelson 1998) of the fractions of variance of every variable over the genes, as variancePartition's ``plotVarPart`` draws them (Hoffman and Schadt 2016); @@ -152,6 +156,11 @@ Viechtbauer W (2010). Conducting meta-analyses in R with the metafor package. Journal of Statistical Software 36(3), 1-48. doi:10.18637/jss.v036.i03 + Nakatsuka N, Adler D, Jiang L, Hartman A, Cheng E, Klann E, Satija R (2025). + Improving reproducibility of differentially expressed genes in single-cell + transcriptomic studies of neurodegenerative diseases through meta-analysis. + Nature Communications 16, 7436. doi:10.1038/s41467-025-62579-z + Machado GM, Oliveira MM, Fernandes LAF (2009). A physiologically-based model for simulation of color vision deficiency. IEEE Transactions on Visualization and Computer Graphics 15, 1291-1298. doi:10.1109/TVCG.2009.113 @@ -207,6 +216,7 @@ from .cube import Cube from .differential_expression import DifferentialExpressionResult from .gene_sets import GeneSets +from .held_out_prediction import HeldOutPrediction from .meta_analysis import MetaAnalysisResult from .normalization import filter_by_expression, log_cpm, tmm_factors from .pathway_scores import PathwayScores @@ -217,6 +227,7 @@ from .stats.multiple_testing import adjust_p_values from .stats.quasi_likelihood import QuasiLikelihoodFit, QuasiLikelihoodTest from .stats.voom import VoomResult +from .sumrank import SumRankResult from .table import Table from .variance_partitioning import ( RESIDUALS, @@ -238,6 +249,7 @@ "plot_gene_sets", "plot_gene_variance_fractions", "plot_gsea_table", + "plot_held_out_auc", "plot_md", "plot_mds", "plot_mean_variance", @@ -251,6 +263,7 @@ "plot_sample_weights", "plot_significant_genes", "plot_study_agreement", + "plot_sumrank", "plot_top_activities", "plot_variance_fractions", "plot_volcano", @@ -4161,3 +4174,198 @@ def _draw_box(ax: "Axes", x: float, values: FloatArray) -> None: zorder=4, gid="outliers", ) + + +# SumRank and held-out prediction. + +_METHOD_LABELS = {"sumrank": "SumRank", "meta_analysis": "Meta-analysis"} + + +def plot_sumrank( + result: SumRankResult, + *, + direction: str = "up", + alpha: float = 0.05, + n_labels: int = 20, + ax: "Axes | None" = None, +) -> "Axes": + """Plot the SumRank p-value of every gene, the genes in alphabetical order. + + The Manhattan plot that the authors of SumRank draw (Nakatsuka et al. 2025, + Figure 4, and their ``MakeManhattanPlot``): a point per gene ranked, at -log10 of + its p-value (``up_p_value`` or ``down_p_value``, from permutations when they were + run), the genes in alphabetical order of their symbol, or of their id without + one. The genes whose adjusted p-value is below ``alpha`` are coloured, red for + genes up and blue for genes down, the ``n_labels`` with the smallest p-values are + named, and a dotted line lies at the largest -log10 p-value of the other genes, + where the authors draw it, when some gene passes. + + Draws the genes (gid ``"genes"``), those that pass (gid ``"significant"``), the + line (gid ``"threshold"``) and the names (gid ``"labels"``). + + Args: + result: The result of + :func:`~cytepool.sumrank.differential_expression_sumrank`. + direction: ``"up"`` or ``"down"``. + alpha: The adjusted p-value below which genes pass, in (0, 1). + n_labels: The most genes named, at least 0. + ax: The Axes to draw into, or None for a new figure. + + Returns: + The Axes. + + Raises: + TypeError: If ``result`` is not a SumRankResult or an option has the wrong + type. + ValueError: If ``direction`` is unknown or ``alpha`` is not in (0, 1). + """ + _check_axes(ax) + if not isinstance(result, SumRankResult): + raise TypeError(f"result must be a SumRankResult, got {type(result).__name__}") + side = checked_choice(direction, "direction", ("up", "down")) + level = _fraction(alpha, "alpha") + reject_masked(n_labels, "n_labels") + named = integral_number(n_labels, "n_labels") + if named < 0: + raise ValueError(f"n_labels must be at least 0, got {named}") + table = result.table + ids = [str(value) for value in table["feature_id"].tolist()] + names = ids + if "symbol" in table: + names = [ + str(symbol) if symbol is not None else identifier + for symbol, identifier in zip(table["symbol"].tolist(), ids, strict=True) + ] + order = sorted(range(len(names)), key=lambda i: (names[i], ids[i])) + p = np.asarray(table[f"{side}_p_value"], dtype=np.float64)[order] + adjusted = np.asarray(table[f"{side}_adjusted_p_value"], dtype=np.float64)[order] + heights = -np.log10(np.maximum(p, _TINY)) + passed = adjusted < level + x = np.arange(len(order), dtype=np.float64) + color = _UP if side == "up" else _DOWN + ax = _axes(ax, (8.0, 4.0)) + ax.scatter( + x[~passed], heights[~passed], s=_DENSE, color=_INK, linewidths=0, gid="genes" + ) + ax.scatter( + x[passed], + heights[passed], + s=_DENSE, + color=color, + linewidths=0, + gid="significant", + ) + if passed.any() and (~passed).any(): + ax.axhline( + float(heights[~passed].max()), + color=color, + linewidth=1.2, + linestyle=":", + gid="threshold", + ) + shown = np.flatnonzero(passed) + shown = shown[np.lexsort((shown, p[shown]))][:named] + for i in sorted(shown.tolist()): + ax.annotate( + names[order[i]], + (x[i], heights[i]), + xytext=(0.0, 3.0), + textcoords="offset points", + ha="center", + va="bottom", + fontsize=7, + color=color, + gid="labels", + ) + ax.set_xticks([]) + ax.set_xlim(-1.0, len(order)) + ax.set_ylim(0.0, max(1.0, float(heights.max(initial=0.0))) * 1.12) + settings = result.settings + source = ( + f"{settings['n_permutations']:,} permutations" + if settings["p_values"] == "permutation" + else "Irwin-Hall" + ) + ax.set_xlabel(f"Genes, in alphabetical order ({len(order):,})") + ax.set_ylabel(f"-log10 p-value, genes {side} ({source})") + _legend( + ax, + [ + _key(f"Adjusted p < {level:g} ({int(passed.sum()):,})", color=color), + _key("Other genes", color=_INK), + ], + ) + return ax + + +def plot_held_out_auc( + prediction: HeldOutPrediction, *, ax: "Axes | None" = None +) -> "Axes": + """Plot the AUC of every held-out study, by the method that chose its genes. + + A row per study of a :func:`~cytepool.held_out_prediction.held_out_prediction`, + in the collection's order from the top, with a point for each method at the AUC + of the study's cases against its controls, scored by the genes that the method + chose from the other studies; the legend gives each method's mean AUC, and a + dotted line marks 0.5, the AUC of a score that does not separate them, as the + paper compares its methods (Nakatsuka et al. 2025, Table 1). A study without an + AUC keeps its row, marked "no AUC". + + Draws the points of each method (gid its name, ``"sumrank"`` or + ``"meta_analysis"``) and the line (gid ``"chance"``). + + Args: + prediction: The result of + :func:`~cytepool.held_out_prediction.held_out_prediction`. + ax: The Axes to draw into, or None for a new figure. + + Returns: + The Axes. + + Raises: + TypeError: If ``prediction`` is not a HeldOutPrediction. + """ + _check_axes(ax) + if not isinstance(prediction, HeldOutPrediction): + raise TypeError( + f"prediction must be a HeldOutPrediction, got {type(prediction).__name__}" + ) + table = prediction.auc + studies = list(dict.fromkeys(table["study"].tolist())) + methods = prediction.settings["methods"] + rows = {study: len(studies) - 1 - i for i, study in enumerate(studies)} + ax = _axes(ax, (6.5, min(12.0, 1.6 + 0.4 * len(studies)))) + ax.axvline(0.5, color=_GREY, linewidth=1, linestyle=":", zorder=0, gid="chance") + handles = [] + offsets = np.linspace(-0.12, 0.12, len(methods)) if len(methods) > 1 else [0.0] + values = np.asarray(table["auc"], dtype=np.float64) + with_auc = set(table["study"][~np.isnan(values)].tolist()) + for method, hue, marker, offset in zip( + methods, _CATEGORICAL, _MARKERS, offsets, strict=False + ): + chosen = table["method"] == method + auc = values[chosen] + names = table["study"][chosen].tolist() + known = ~np.isnan(auc) + y = np.array([rows[name] for name in names], dtype=np.float64) + offset + ax.scatter( + auc[known], + y[known], + s=_SPARSE, + color=hue, + marker=marker, + zorder=3, + clip_on=False, # an AUC of 0 or 1 lies on the frame + gid=method, + ) + mean = f"mean {auc[known].mean():.2f}" if known.any() else "no AUC" + label = _METHOD_LABELS.get(method, method) + handles.append(_key(f"{label} ({mean})", color=hue, marker=marker)) + labels = [study if study in with_auc else f"{study} (no AUC)" for study in studies] + ax.set_yticks([rows[study] for study in studies], labels) + ax.tick_params(axis="y", length=0) + ax.set_ylim(-0.7, len(studies) - 0.3) + ax.set_xlim(0.0, 1.0) + ax.set_xlabel("AUC of the held-out study's cases against its controls") + _legend(ax, handles) + return ax diff --git a/cytepool/sumrank.py b/cytepool/sumrank.py new file mode 100644 index 0000000..80302a8 --- /dev/null +++ b/cytepool/sumrank.py @@ -0,0 +1,1050 @@ +"""SumRank: genes whose differential expression ranks high in study after study. + +Most genes that one single-cell study finds differentially expressed are not found by +the next (Nakatsuka et al. 2025). SumRank asks instead whether a gene ranks high in +every study: each study ranks its genes by their signed evidence of change, and a +gene's ranks are summed over the studies, so that a gene ranked high everywhere has a +small sum and a gene found by one study alone does not, however small its p-value +there. :func:`sumrank` computes the statistic and its p-values from a signed statistic +of every gene in every study; :func:`differential_expression_sumrank` does it for the +results of the studies of a :class:`~cytepool.cube.CubeCollection`, one group (cell +type) and contrast at a time, with p-values from permutations that re-run each study's +analysis on its samples with their labels permuted. + +The statistic, as the paper defines it and its authors' code (Apache-2.0) computes it: + 1. Each study ranks the genes by their signed statistic, from 1 for the largest + (the gene most up); tied genes share the average of their ranks, as R's + ``rank``, which the authors' code calls, gives them. Results of + differential expression give the signed -log10 p-value: -log10(p_value), + negative where the log2 fold change is negative, as the paper ranks DESeq2's + results. + 2. A gene of rank r among G genes has the normalized rank u = (r - 1) / (G - 1): + 0 for the gene most up and 1 for the gene most down. + 3. For each gene, the m smallest of its normalized ranks in the k studies are + summed: those of the studies that rank it best, which leaves out studies that + rank it low for reasons of their own. m = floor(p k), at least 1, for the + proportion p of the studies (``proportion_of_studies``); the product is rounded + to 9 decimals first, so that 0.57 of 100 studies is 57. The paper sums every + study when there are few (p = 1 for its 3 to 6 studies of Parkinson's, + Huntington's disease and schizophrenia) and chose p by prediction in held-out + studies when there are many (0.65 of 17 studies of Alzheimer's disease, 0.55 + of 16 of COVID-19); :func:`~cytepool.held_out_prediction.held_out_prediction` + gives that prediction for a chosen p. + 4. A sum above m / 2, the mean of a sum of m uniform ranks, is set to m / 2, so + that genes without consistent evidence share one value. + + The result is the statistic of genes up (``up_statistic``): small for a gene that + ranks high in its best m studies. Genes down have the same statistic of the + negated statistics (``down_statistic``), whose normalized ranks are 1 - u, as the + paper tests genes down by testing controls against cases. + +Genes that a study did not measure or did not test: + The genes ranked are those that every study tested (a finite or infinite + statistic, not NaN), as the paper ranks the genes that every dataset measured, + the intersection of their genes. A rank says where a gene stands among the genes + ranked with it, so the studies must rank the same genes; and a gene ranked in + fewer studies would have a statistic, and a distribution under the null, of its + own. The other genes have NaN. A study that did not analyse the group or the + contrast at all is left out by :func:`differential_expression_sumrank`, as the + authors' code leaves out a dataset without the cell type: k counts the studies + that did. + +Irwin-Hall p-values: + If every study ranked the genes in an order of its own, at random and + independently of the others, each normalized rank would be uniform on [0, 1], + and the sum of m of them would follow the Irwin-Hall distribution of a sum of m + independent uniform variables (Irwin 1927; Hall 1927), whose distribution + function is + + F_m(x) = (1/m!) sum_{j=0}^{floor(x)} (-1)^j C(m, j) (x - j)^m. + + ``up_irwin_hall_p_value`` is F_m(s) for the statistic s of a gene below the cap, + and 1 at the cap: the probability under that null of a statistic at most s. + It treats the ranks as continuous: a gene ranked first in every study (s = 0) + has p-value 0, where its probability is G^-m. + + With every study summed (m = k), the p-value is exact under that null, apart + from the steps of the ranks. A sum of the m best of k ranks (m < k) is smaller + than a sum of m ranks, so its Irwin-Hall p-values are too small: the paper takes + its p-values from permutations instead, as below, and so should an analysis + with m < k, as :func:`differential_expression_sumrank` warns. Genes whose ranks + are not uniform under the null, such as genes with a p-value of 1 in many + permutations, depart from it too. + + The alternating sum above loses every digit to cancellation near the middle of + the distribution for tens of studies. F_m is computed instead as a sum of + positive terms, F_m(x) = sum_{i=0}^{floor(x)} N_{m+1}(x - floor(x) + i), where + N_{m+1} is the cardinal B-spline of order m + 1, the density of the Irwin-Hall + distribution of m + 1 variables, which de Boor's recurrence evaluates with + positive terms only (de Boor 1978): N_{m+1}' (x) = N_m(x) - N_m(x - 1), so the + sum telescopes to the integral of N_m. It agrees with a high-precision + computation to about 1e-14 relative for up to 200 studies + (``validation/sumrank.py``). + +Permutation p-values: + The paper permutes the case and control labels of the samples of every study, + runs each study's analysis and SumRank again, and pools the statistics of every + gene in every permutation into one null distribution: with B permutations of G + genes, N = B G null statistics. The p-value of a gene of statistic s is + + p = (1 + the number of null statistics at most s) / (1 + N), + + as Phipson and Smyth (2010) form p-values from random permutations, so that it + is never zero, and a gene at the cap has p = 1. Pooling over the genes gives + small p-values from few permutations (100 permutations of 10,000 genes give a + p-value as small as 1e-6), for genes whose null distributions are alike. The + statistics of a permutation are computed as those of the data, over the genes + that every study's analysis of that permutation tested. + + :func:`sumrank` takes the statistics of the permutations as an array; + :func:`differential_expression_sumrank` permutes the values of one variable of + each study's design (``permute``, such as the condition) among the samples of + the study's fit of the group, with a ``numpy.random.Generator``, and re-runs the + study's analysis, ``limma_voom`` or ``limma_trend`` with the design and settings + its result records, on its aligned cube, for that group and contrast only. The + permutations are drawn in a fixed order, permutation after permutation and + study after study, so a seed gives the same permutations, and the analyses do + not depend on the number of threads. + +Cost: + A permutation runs the analysis of every study once, for one group: with + ``limma_voom``, 0.21 s for a study of 40 samples and 10,000 genes on a laptop on + 2 threads (``limma_trend`` 0.14 s; ``validation/sumrank.py --benchmark``), so 100 + permutations of 20 such studies take about 7 minutes; the SumRank statistics of a + permutation take milliseconds. The memory is that of one study's analysis and a + few arrays of the genes, whatever the number of permutations. + +Deliberate differences from the paper and its authors' code: + - The Irwin-Hall p-value is the lower tail of the distribution, as above. The + paper calls its Irwin-Hall p-values two-sided, which would double these, and + the authors' code computes the density of the Irwin-Hall distribution at the + statistic instead (the unifed package's ``dirwin.hall``), which is not a + probability. Its permutation p-values compare -log10 of that density, which + falls as the statistic rises up to the cap for m >= 2: they order genes as the + statistic does, so the permutation p-values here do not depend on this choice. + - The authors' code takes as the permutation p-value the share of null + statistics strictly below the gene's (``1 - ecdf`` of -log10 of the density), + without the one above: a gene beyond every null statistic has the p-value 0, + and a gene at the cap the share of null statistics below the cap rather than + 1. Here ties count as at least as extreme, which keeps the p-value valid. + - The labels are permuted among the samples of each study's fit of the group, + which keeps the number of cases and controls of the fit. The paper permutes + them among all the individuals of a dataset, for all its cell types at once, + and draws again until every cell type has at least two cases and two controls. + - The authors' code sums the normalized ranks in floating point, in each gene's + own order, so that genes with equal statistics can compare as unequal in its + calibration; here the sums are exact. It stops when it sums one study (m = 1), + and takes floor(p k) of the product in floating point, which for 0.57 of 100 + studies is 56. + +Not done: + Studies are not weighted, a limit the paper notes too. Permutations re-run + designs without random effects only (``limma_voom`` and ``limma_trend``), whose + samples are exchangeable under the null; a design with several samples of a + donor needs the donors' labels permuted instead. + +References: + Nakatsuka N, Adler D, Jiang L, Hartman A, Cheng E, Klann E, Satija R (2025). + Improving reproducibility of differentially expressed genes in single-cell + transcriptomic studies of neurodegenerative diseases through meta-analysis. + Nature Communications 16, 7436. doi:10.1038/s41467-025-62579-z + + Irwin JO (1927). On the frequency distribution of the means of samples from a + population having any law of frequency with finite moments, with special + reference to Pearson's Type II. Biometrika 19, 225-239. + doi:10.1093/biomet/19.3-4.225 + + Hall P (1927). The distribution of means for samples of size N drawn from a + population in which the variate takes values between 0 and 1, all such values + being equally probable. Biometrika 19, 240-244. doi:10.1093/biomet/19.3-4.240 + + de Boor C (1978). A Practical Guide to Splines. Springer, New York. + doi:10.1007/978-1-4612-6333-3 + + Phipson B, Smyth GK (2010). Permutation p-values should never be zero: calculating + exact p-values when permutations are randomly drawn. Statistical Applications in + Genetics and Molecular Biology 9, article 39. doi:10.2202/1544-6115.1585 + + Benjamini Y, Hochberg Y (1995). Controlling the false discovery rate: a practical + and powerful approach to multiple testing. Journal of the Royal Statistical + Society B 57, 289-300. doi:10.1111/j.2517-6161.1995.tb02031.x +""" + +import inspect +import json +import math +import numbers +import warnings +from collections.abc import Callable, Iterator, Mapping +from dataclasses import dataclass +from typing import Any + +import numpy as np +import numpy.typing as npt + +from ._analysis_records import software_versions +from ._columns import text_column +from ._jit import jit +from ._masked_arrays import read_only, reject_masked, unmasked_array +from ._messages import PACKAGE_DIRECTORY, listing, quoted +from ._options import index_integer, random_generator +from ._ranks import average_ranks +from ._study_statistics import signed_statistics, study_statistics, studies_summed +from .cube import Cube, CubeCollection +from .design import Design +from .differential_expression import ( + DifferentialExpressionResult, + limma_trend, + limma_voom, +) +from .stats.multiple_testing import adjust_p_values +from .table import Table + +__all__ = [ + "SumRankFit", + "SumRankResult", + "differential_expression_sumrank", + "sumrank", +] + +FloatArray = npt.NDArray[np.float64] +IntArray = npt.NDArray[np.int64] +BoolArray = npt.NDArray[np.bool_] + +_DIRECTIONS = ("up", "down") +# The analyses that permutations re-run: a design without random effects, whose +# samples are exchangeable under the null. +_RERUN: dict[str, Callable[..., DifferentialExpressionResult]] = { + "limma_voom": limma_voom, + "limma_trend": limma_trend, +} + + +@dataclass(frozen=True, eq=False, repr=False) +class SumRankFit: + """SumRank of every gene across studies, from :func:`sumrank`. + + The per-gene arrays hold one value per gene (feature), NaN for a gene that some + study did not test; they are read-only. + + Attributes: + up_statistic: The sum of the gene's ``n_studies_summed`` smallest normalized + ranks, set to ``n_studies_summed / 2`` where it is larger: small for a + gene up in the studies that rank it best. + down_statistic: The same for the negated statistics: small for a gene down. + up_irwin_hall_p_value: The Irwin-Hall distribution function of the sum of + ``n_studies_summed`` uniform ranks at ``up_statistic``, or 1 at the cap. + down_irwin_hall_p_value: The same for ``down_statistic``. + up_permutation_p_value: (1 + the number of null statistics at most + ``up_statistic``) / (1 + ``n_null_values``), over the statistics of + every gene in every permutation; None without permutations. + down_permutation_p_value: The same for ``down_statistic``. + normalized_ranks: ``(n_studies, n_features)``: the normalized rank of every + gene in every study, (r - 1) / (G - 1) for the rank r from 1 for the + largest statistic among the G genes ranked; NaN for a gene not ranked. + The ranks of the genes down are 1 minus these. + ranked: Whether each gene was ranked: tested by every study. + n_studies: The number of studies, k. + n_studies_summed: The number of ranks summed for each gene, m. + proportion_of_studies: The proportion of the studies that sets m. + n_permutations: The number of permutations, 0 without them. + n_null_values: The number of null statistics of each direction, the genes + ranked in each permutation, summed over the permutations; 0 without + them. + """ + + up_statistic: FloatArray + down_statistic: FloatArray + up_irwin_hall_p_value: FloatArray + down_irwin_hall_p_value: FloatArray + up_permutation_p_value: FloatArray | None + down_permutation_p_value: FloatArray | None + normalized_ranks: FloatArray + ranked: BoolArray + n_studies: int + n_studies_summed: int + proportion_of_studies: float + n_permutations: int + n_null_values: int + + def __repr__(self) -> str: + text = ( + f"SumRankFit(n_features={self.ranked.size}, ranked=" + f"{int(np.count_nonzero(self.ranked))}, n_studies={self.n_studies}, " + f"n_studies_summed={self.n_studies_summed}" + ) + if self.n_permutations: + text += f", n_permutations={self.n_permutations}" + return text + ")" + + +def sumrank( + statistics: npt.ArrayLike, + *, + proportion_of_studies: float = 1.0, + null_statistics: npt.ArrayLike | None = None, +) -> SumRankFit: + """Rank every gene in every study and sum its best ranks, as SumRank does. + + Each study ranks the genes that every study tested by their statistic, from the + largest; a gene's normalized ranks in the studies that rank it best are summed and + capped at their mean under the null, for genes up and, with the statistics + negated, for genes down. The module description gives the statistic, its + Irwin-Hall p-values and the permutation p-values. + + Args: + statistics: ``(n_studies, n_features)`` signed statistics, larger for a gene + more up, such as the signed -log10 p-values of each study's test; NaN + where a study did not test a gene, and infinite values are allowed (a + p-value of 0). Only their order within each study matters. + proportion_of_studies: The proportion p of the k studies whose ranks are + summed for each gene: the m = floor(p k) studies that rank it best, at + least one, in (0, 1]; 1 sums every study. + null_statistics: ``(n_permutations, n_studies, n_features)``: the statistics + of the same studies with their labels permuted, one study's analysis + re-run on each permutation, NaN where an analysis did not test a gene; + or None for no permutation p-values. + + Returns: + The statistics of every gene, up and down, their Irwin-Hall p-values and, + with ``null_statistics``, their permutation p-values. + + Raises: + TypeError: If an argument is not an array of real numbers, is a numpy masked + array, or ``proportion_of_studies`` is not a real number. + ValueError: If the shapes are wrong, a study tested no gene, fewer than two + genes were tested by every study (in the data or in a permutation), or + ``proportion_of_studies`` is not in (0, 1]. + + Example: + Two studies of four genes: the first gene is the largest in both, so its + ranks sum to 0; the others reach the cap of 1 (two ranks of mean 1/2). + + >>> fit = sumrank([[3.0, 1.0, -2.0, 0.5], [2.0, -1.0, 0.0, 1.5]]) + >>> fit.up_statistic.tolist() + [0.0, 1.0, 1.0, 1.0] + >>> fit.down_statistic.tolist() + [1.0, 0.6666666666666666, 0.3333333333333333, 1.0] + >>> fit.down_irwin_hall_p_value.tolist() # (1/3)^2 / 2 = 1/18 for the third + [1.0, 0.2222222222222222, 0.05555555555555555, 1.0] + """ + values = _statistics(statistics, "statistics", 2) + n_studies = values.shape[0] + summed = studies_summed(proportion_of_studies, n_studies) + null = None + if null_statistics is not None: + null = _statistics(null_statistics, "null_statistics", 3) + if null.shape[1:] != values.shape: + raise ValueError( + "null_statistics must be (n_permutations, n_studies, n_features) with " + f"the studies and features of statistics, {values.shape}; got " + f"{null.shape}" + ) + if null.shape[0] == 0: + raise ValueError( + "null_statistics holds no permutation; give at least one, or None for " + "no permutation p-values" + ) + observed = _capped_sums(values, summed, "statistics") + up = _full(observed.up, observed.ranked) + down = _full(observed.down, observed.ranked) + permutation: list[FloatArray | None] = [None, None] + n_null = 0 + if null is not None: + counts = _NullCounts(observed) + for b, permuted in enumerate(null): + counts.add(_capped_sums(permuted, summed, f"permutation {b}")) + permutation = [read_only(_full(p, observed.ranked)) for p in counts.p_values()] + n_null = counts.n_values + return SumRankFit( + up_statistic=read_only(up), + down_statistic=read_only(down), + up_irwin_hall_p_value=read_only(_full(observed.up_p, observed.ranked)), + down_irwin_hall_p_value=read_only(_full(observed.down_p, observed.ranked)), + up_permutation_p_value=permutation[0], + down_permutation_p_value=permutation[1], + normalized_ranks=read_only(observed.ranks_of_every_feature()), + ranked=read_only(observed.ranked), + n_studies=n_studies, + n_studies_summed=summed, + proportion_of_studies=float(proportion_of_studies), + n_permutations=0 if null is None else int(null.shape[0]), + n_null_values=n_null, + ) + + +# The statistic. + + +def _statistics(values: Any, name: str, ndim: int) -> FloatArray: + array = unmasked_array(values, name, missing="NaN") + if array.dtype.kind not in "iuf": + raise TypeError(f"{name} must hold real numbers, got dtype {array.dtype}") + if array.ndim != ndim: + shape = ( + "(n_studies, n_features)" + if ndim == 2 + else "(n_permutations, n_studies, n_features)" + ) + raise ValueError(f"{name} must be {shape}, got {array.ndim} dimensions") + converted: FloatArray = np.array(array, dtype=np.float64) + if ndim == 2 and converted.shape[0] == 0: + raise ValueError(f"{name} has no study: give one row per study") + return converted + + +@dataclass(frozen=True) +class _Sums: + """The capped sums of the genes every study tested, in units that are exact. + + ``up`` and ``down`` are the statistics of the ``ranked`` genes, in their order, + and ``up_p`` and ``down_p`` their Irwin-Hall p-values; ``ranks`` holds the + normalized ranks of these genes in every study. + """ + + ranked: BoolArray # per feature + up: FloatArray # per ranked gene + down: FloatArray + up_p: FloatArray + down_p: FloatArray + ranks: FloatArray # (n_studies, n_ranked) + + def ranks_of_every_feature(self) -> FloatArray: + full = np.full((self.ranks.shape[0], self.ranked.size), np.nan) + full[:, self.ranked] = self.ranks + return full + + +def _capped_sums(values: FloatArray, summed: int, where: str) -> _Sums: + """The SumRank statistics of the genes every study tested, up and down. + + The rank r of a gene among G (from 1 for the largest value) is an integer or, for + ties, a half-integer, so q = 2 (r - 1) is an integer and a sum of normalized ranks + is (sum of q) / (2 (G - 1)): it is formed exactly in integers, and compared with + the cap m (G - 1) there, so that equal sums are equal in every permutation. + """ + tested = ~np.isnan(values) + empty = np.flatnonzero(~tested.any(axis=1)) + if empty.size: + raise ValueError( + f"{where}: study {int(empty[0])} tested no gene (every value is NaN); " + "leave out a study that did not test the group and contrast" + ) + ranked = tested.all(axis=0) + n_ranked = int(np.count_nonzero(ranked)) + if n_ranked < 2: + genes = "1 gene was" if n_ranked == 1 else f"{n_ranked} genes were" + raise ValueError( + f"{where}: {genes} tested by every study, and SumRank ranks the genes " + "that every study tested against each other; it needs at least two" + ) + kept = values[:, ranked] + n_studies = kept.shape[0] + twice_up = np.empty(kept.shape, dtype=np.int64) # 2 (r - 1), r from the largest + for row in range(n_studies): + ascending = average_ranks(np.ascontiguousarray(kept[row])) # 1 = smallest + # r = G + 1 - ascending rank, so 2 (r - 1) = 2 (G - ascending rank). + twice_up[row] = np.rint(2.0 * (n_ranked - ascending)).astype(np.int64) + twice_down = 2 * (n_ranked - 1) - twice_up + scale = 2.0 * (n_ranked - 1) + statistics = [] + for twice in (twice_up, twice_down): + if summed < n_studies: + best = np.partition(twice, summed - 1, axis=0)[:summed] + else: + best = twice + total = best.sum(axis=0) + cap = summed * (n_ranked - 1) + capped = np.minimum(total, cap) + statistic = capped / scale + p_value = np.ones(statistic.size) + below = capped < cap + p_value[below] = _irwin_hall_cdf(statistic[below], summed) + statistics.append((statistic, p_value)) + return _Sums( + ranked=ranked, + up=statistics[0][0], + down=statistics[1][0], + up_p=statistics[0][1], + down_p=statistics[1][1], + ranks=twice_up / scale, + ) + + +def _full(values: FloatArray, ranked: BoolArray) -> FloatArray: + """Values of the ranked genes in an array of every feature, NaN elsewhere.""" + full = np.full(ranked.size, np.nan) + full[ranked] = values + return full + + +class _NullCounts: + """For every gene, the number of null statistics at most its own, per direction, + accumulated one permutation at a time, so that the null distribution is never + held whole.""" + + def __init__(self, observed: _Sums) -> None: + self._observed = (observed.up, observed.down) + self._counts = [np.zeros(observed.up.size, dtype=np.int64) for _ in range(2)] + self.n_values = 0 + + def add(self, permutation: _Sums) -> None: + for counts, statistic, null in zip( + self._counts, + self._observed, + (permutation.up, permutation.down), + strict=True, + ): + counts += np.searchsorted(np.sort(null), statistic, side="right") + self.n_values += permutation.up.size + + def p_values(self) -> list[FloatArray]: + return [ + np.divide(1.0 + counts, 1.0 + self.n_values, dtype=np.float64) + for counts in self._counts + ] + + +# The Irwin-Hall distribution function. + + +def _irwin_hall_cdf(x: FloatArray, n: int) -> FloatArray: + """F_n at each x, the distribution function of a sum of n uniform variables.""" + out = np.empty(x.size) + _irwin_hall_cdf_into(np.ascontiguousarray(x, dtype=np.float64), n, out) + return out + + +@jit +def _irwin_hall_cdf_into(x: FloatArray, n: int, out: FloatArray) -> None: + """Fill ``out`` with F_n(x), as a sum of values of the cardinal B-spline of order + n + 1 at t + i for the fractional part t of x, which de Boor's recurrence + evaluates with positive terms only; above n / 2 by the symmetry F_n(x) = + 1 - F_n(n - x).""" + b = np.empty(n + 1) + for g in range(x.size): + value = x[g] + if not value > 0.0: + out[g] = 0.0 if value <= 0.0 else math.nan + continue + if not value < n: + out[g] = 1.0 + continue + upper = value > 0.5 * n + if upper: + value = n - value + whole = int(math.floor(value)) + t = value - whole + # b[i] = N_k(t + i) for the cardinal B-spline N_k of order k, from k = 1 + # (the indicator of [0, 1)) to k = n + 1, updating from the top so that + # b[i - 1] still holds order k - 1. + b[0] = 1.0 + for i in range(1, n + 1): + b[i] = 0.0 + for k in range(2, n + 2): + for i in range(k - 1, -1, -1): + term = (t + i) * b[i] + if i > 0: + term += (k - t - i) * b[i - 1] + b[i] = term / (k - 1) + total = 0.0 + for i in range(whole + 1): + total += b[i] + out[g] = 1.0 - total if upper else total + + +# SumRank of the studies of a collection. + + +class SumRankResult: + """SumRank of the differential expression of the studies of a collection. + + Built by :func:`differential_expression_sumrank`, for one group (cell type) and + contrast. The table has one row per gene that every study tested, in the + collection's order of the features: + + - ``feature_id`` and, when the collection's features have it, ``symbol``; + - ``up_statistic``: the SumRank statistic of the gene up, the capped sum of its + best normalized ranks (small for a gene up in the studies that rank it best); + - ``up_p_value``: its p-value, from permutations when they were run + (``settings["p_values"]`` is ``"permutation"``), from the Irwin-Hall + distribution otherwise (``"irwin_hall"``); + - ``up_adjusted_p_value``: Benjamini-Hochberg over the genes of the table + (Benjamini and Hochberg 1995); + - ``up_irwin_hall_p_value``: the Irwin-Hall p-value, also with permutations; + - ``down_statistic``, ``down_p_value``, ``down_adjusted_p_value`` and + ``down_irwin_hall_p_value``: the same for the gene down. + + The table of studies has one row per study of the collection: ``study``, + ``analysed`` (whether its results test the group and contrast), ``reason`` (why + not; empty when analysed) and ``n_features_tested`` (0 when not analysed). + + Immutable. + + Raises: + TypeError: If ``table`` or ``studies`` is not a ``Table``, or ``provenance`` + not a mapping that JSON can write. + ValueError: If ``provenance`` holds a NaN or an infinity. + """ + + __slots__ = ("_table", "_studies", "_record") + + _table: Table + _studies: Table + _record: str + + def __init__( + self, table: Table, studies: Table, provenance: Mapping[str, Any] + ) -> None: + for name, given in (("table", table), ("studies", studies)): + if not isinstance(given, Table): + raise TypeError(f"{name} must be a Table, got {type(given).__name__}") + if not isinstance(provenance, Mapping): + raise TypeError( + f"provenance must be a mapping, got {type(provenance).__name__}" + ) + state = { + "_table": table, + "_studies": studies, + # JSON holds no arrays: json.dumps raises TypeError for one. + "_record": json.dumps(provenance, allow_nan=False), + } + for name, value in state.items(): + object.__setattr__(self, name, value) + + def __setattr__(self, name: str, value: object) -> None: + raise AttributeError("SumRankResult is immutable") + + def __delattr__(self, name: str) -> None: + raise AttributeError("SumRankResult is immutable") + + def __reduce__(self) -> tuple[Any, tuple[Any, ...]]: + return ( + SumRankResult, + (self._table, self._studies, json.loads(self._record)), + ) + + @property + def table(self) -> Table: + """The statistics and p-values of every gene ranked.""" + return self._table + + @property + def studies(self) -> Table: + """Every study of the collection, whether it analysed the group and + contrast, and why not.""" + return self._studies + + @property + def settings(self) -> dict[str, Any]: + """The settings, as a new dict (also in the provenance), with the group, the + contrast, the proportion of studies, the numbers of studies and of ranks + summed, the source of the p-values and the permutations.""" + settings: dict[str, Any] = json.loads(self._record)["settings"] + return settings + + @property + def provenance(self) -> dict[str, Any]: + """A JSON-serialisable record of what produced every number, as a new dict. + + It holds the CytePool, numpy, numba and Python versions; the settings; the + number of genes ranked and of genes left out, measured by a study but not + tested by every study that analysed the group and contrast; and the + provenance of every study's differential-expression result, by study name, + which identifies its cube by its fingerprint and gives its design and + settings and the samples of every group's fit. + """ + record: dict[str, Any] = json.loads(self._record) + return record + + def __repr__(self) -> str: + settings = self.settings + table = self._table + lines = [ + f"SumRankResult for group {settings['group']!r} and contrast " + f"{settings['contrast']!r}: {settings['n_studies']} studies, the best " + f"{settings['n_studies_summed']} summed; {len(table):,} genes" + ] + if settings["p_values"] == "permutation": + lines.append( + f" p-values from {settings['n_permutations']:,} permutations of " + f"{settings['permute']!r}" + ) + else: + lines.append(" p-values from the Irwin-Hall distribution") + counts = [ + int(np.count_nonzero(table[f"{d}_adjusted_p_value"] < 0.05)) + for d in _DIRECTIONS + ] + lines.append( + f" adjusted p-value below 0.05: {counts[0]:,} up, {counts[1]:,} down" + ) + return "\n".join(lines) + + +def differential_expression_sumrank( + collection: CubeCollection, + results: Mapping[str, DifferentialExpressionResult], + group: str, + contrast: str, + *, + proportion_of_studies: float = 1.0, + n_permutations: int = 0, + permute: str | None = None, + rng: np.random.Generator | int | None = None, +) -> SumRankResult: + """SumRank of one group and contrast across the studies of a collection. + + Each study's result is the analysis of its aligned cube + (:meth:`~cytepool.cube.CubeCollection.aligned_cube`), such as + :func:`~cytepool.differential_expression.limma_voom` with the study's own design, + as for :func:`~cytepool.meta_analysis.align_study_results`. For the group and + contrast, the studies whose results test them each rank the genes that every one + of them tested by their signed -log10 p-value, and :func:`sumrank` sums each + gene's best ranks; the module description gives the statistic and its p-values. + + With ``n_permutations``, the values of the design variable ``permute``, such as + the condition, are permuted among the samples of each study's fit of the group, + and each study's analysis is run again on its aligned cube with the design and + settings its result records (``limma_voom`` or ``limma_trend``), for this group + and contrast only; the SumRank statistics of every permutation make the null + distribution of the p-values. Each permutation costs one analysis of the group in + every study (see "Cost" in the module description). + + Args: + collection: The collection whose aligned cubes the studies' results + analysed. + results: The differential-expression result of every study of the + collection, by study name, from + :func:`~cytepool.differential_expression.limma_voom`, + :func:`~cytepool.differential_expression.limma_trend` or + :func:`~cytepool.differential_expression.dream_voom` (without + permutations). + group: An aligned group id of the collection. + contrast: A contrast name that the studies tested. + proportion_of_studies: The proportion of the studies that analysed the group + whose ranks are summed for each gene, those that rank it best, in + (0, 1]; 1 sums every study. + n_permutations: The number of permutations, at least 0; 0 gives Irwin-Hall + p-values only. + permute: The variable of the studies' designs whose values are permuted, + such as ``"condition"``; needed with permutations. + rng: A ``numpy.random.Generator`` to draw the permutations from, an integer + seed for one, or None for a fresh, unseeded generator. + + Returns: + The statistics and p-values of every gene that every study analysing the + group and contrast tested, a table of the studies, and a provenance record. + + Warns: + RuntimeWarning: If fewer than all the studies are summed without + permutations: the p-values are then the Irwin-Hall p-values, which are + too small for such sums. + + Raises: + TypeError: If an argument has the wrong type, or is or holds a numpy masked + array. + ValueError: If ``results`` does not give the result of every study of the + collection, a result did not analyse its study's aligned cube, no study + analysed the group and contrast, the contrast is an F-test, fewer than + two genes were tested by every study, an option is invalid, + ``permute`` is given without permutations or missing with them, or a + study cannot be permuted: its analysis has random effects, its design + has no variable ``permute``, or a permutation leaves the contrast + untestable. + KeyError: If ``group`` is not a group of the collection, or no study tested + ``contrast``. + """ + reject_masked(n_permutations, "n_permutations") + permutations = _n_permutations(n_permutations) + if permute is not None and not isinstance(permute, str): + reject_masked(permute, "permute") + raise TypeError( + f"permute must be a design variable name (str), got {type(permute).__name__}" + ) + if permutations and permute is None: + raise ValueError( + "permutations permute the values of one variable of each study's design, " + "such as the condition: give its name as permute" + ) + if not permutations and permute is not None: + raise ValueError( + "permute names the variable that permutations permute, but n_permutations " + "is 0; give the number of permutations" + ) + generator = random_generator(rng) if permutations else None + inputs = study_statistics(collection, results, group, contrast) + used = inputs.analysed + if not used.any(): + raise ValueError( + f"no study tested group {group!r} and contrast {contrast!r}: " + + "; ".join( + f"{study}: {reason}" + for study, reason in zip( + collection.studies, inputs.reasons, strict=True + ) + ) + ) + studies = [ + s for s, flag in zip(collection.studies, used.tolist(), strict=True) if flag + ] + statistics = inputs.statistics[used] + summed = studies_summed(proportion_of_studies, len(studies)) + if summed < len(studies) and not permutations: + warnings.warn( + f"the Irwin-Hall p-values of a sum of the best {summed} of " + f"{len(studies)} studies' ranks are too small, as the sum of the best " + "ranks is smaller than a sum of uniform ranks; the p_value columns hold " + "them, so take p-values from permutations (n_permutations and permute)", + RuntimeWarning, + skip_file_prefixes=(PACKAGE_DIRECTORY,), + ) + observed = _capped_sums(statistics, summed, "the studies' results") + counts = None + if permutations: + assert generator is not None and permute is not None + counts = _NullCounts(observed) + for b, permuted in enumerate( + _permuted_statistics( + collection, + results, + studies, + group, + contrast, + permute, + permutations, + generator, + ) + ): + counts.add(_capped_sums(permuted, summed, f"permutation {b}")) + ranked = observed.ranked + features = collection.features + columns: dict[str, Any] = {"feature_id": features["feature_id"][ranked]} + if "symbol" in features: + columns["symbol"] = features["symbol"][ranked] + p_values = [observed.up_p, observed.down_p] if counts is None else counts.p_values() + for direction, statistic, p_value, irwin_hall in zip( + _DIRECTIONS, + (observed.up, observed.down), + p_values, + (observed.up_p, observed.down_p), + strict=True, + ): + columns[f"{direction}_statistic"] = statistic + columns[f"{direction}_p_value"] = p_value + columns[f"{direction}_adjusted_p_value"] = adjust_p_values(p_value, "bh") + columns[f"{direction}_irwin_hall_p_value"] = irwin_hall + measured = inputs.measured[used].any(axis=0) + settings = { + "group": group, + "contrast": contrast, + "statistic": "signed -log10 p-value", + "proportion_of_studies": float(proportion_of_studies), + "n_studies": len(studies), + "n_studies_summed": summed, + "p_values": "irwin_hall" if counts is None else "permutation", + "n_permutations": permutations, + "permute": permute, + "seed": ( + int(rng) + if isinstance(rng, numbers.Integral) + and not isinstance(rng, bool | np.bool_) + else None + ), + "n_null_values": 0 if counts is None else counts.n_values, + "p_value_adjustment": "BH", + } + return SumRankResult( + Table(columns), + Table( + { + "study": text_column(list(collection.studies)), + "analysed": used.copy(), + "reason": text_column(inputs.reasons), + "n_features_tested": inputs.n_tested, + } + ), + { + "analysis": "differential_expression_sumrank", + "versions": software_versions(), + "settings": settings, + "n_features": int(np.count_nonzero(ranked)), + "n_features_left_out": int(np.count_nonzero(measured & ~ranked)), + "studies": {name: results[name].provenance for name in collection.studies}, + }, + ) + + +def _n_permutations(value: Any) -> int: + if isinstance(value, bool | np.bool_): + raise TypeError("n_permutations must be a whole number, not bool") + number = index_integer(value, "n_permutations", expected="a whole number") + if number < 0: + raise ValueError(f"n_permutations must be at least 0, got {number}") + return number + + +# Permutations: each study's analysis re-run with the labels of its samples permuted. + + +@dataclass(frozen=True) +class _Rerun: + """What it takes to run a study's analysis again with permuted labels.""" + + study: str + cube: Cube + function: Callable[..., DifferentialExpressionResult] + options: dict[str, Any] + formula: str + levels: dict[str, list[Any]] + samples: dict[str, list[Any]] + weights: list[float] # of the contrast on the design's columns + used: list[int] # the design rows of the samples of the group's fit + fingerprint: str + settings: dict[str, Any] + + +def _permuted_statistics( + collection: CubeCollection, + results: Mapping[str, DifferentialExpressionResult], + studies: list[str], + group: str, + contrast: str, + permute: str, + n_permutations: int, + generator: np.random.Generator, +) -> Iterator[FloatArray]: + """The signed statistics of every study for each permutation, ``(n_studies, + n_features)`` in the collection's order, NaN where an analysis did not test a + gene. The permutations are drawn permutation by permutation, study by study.""" + reruns = [ + _rerun(collection, results[study], study, group, contrast, permute) + for study in studies + ] + index = { + feature: position + for position, feature in enumerate(collection.features["feature_id"].tolist()) + } + checked = [False] * len(reruns) + for _ in range(n_permutations): + statistics = np.full((len(reruns), len(index)), np.nan) + for row, rerun in enumerate(reruns): + values = list(rerun.samples[permute]) + drawn = generator.permutation(len(rerun.used)) + for position, source in zip(rerun.used, drawn.tolist(), strict=True): + values[position] = rerun.samples[permute][rerun.used[source]] + design = Design( + {**rerun.samples, permute: values}, + rerun.formula, + levels=rerun.levels, + contrasts={contrast: rerun.weights}, + ) + result = rerun.function(rerun.cube, design, groups=[group], **rerun.options) + if not checked[row]: + _check_rerun(rerun, result) + checked[row] = True + table = result.tables.get((group, contrast)) + if table is None: + raise ValueError( + f"a permutation of {permute!r} in study {rerun.study!r} leaves " + f"contrast {contrast!r} untestable in group {group!r}: " + f"{_absence(result, group, contrast)}" + ) + positions = np.array( + [index[feature] for feature in table["feature_id"].tolist()], + dtype=np.intp, + ) + statistics[row, positions] = signed_statistics( + table["p_value"], table["log2_fold_change"] + ) + yield statistics + + +def _rerun( + collection: CubeCollection, + result: DifferentialExpressionResult, + study: str, + group: str, + contrast: str, + permute: str, +) -> _Rerun: + """How to run a study's analysis again: its function and options, its design, + and the samples of its fit of the group, from what its result records.""" + function = _RERUN.get(result.method) + if function is None: + raise ValueError( + f"the result of study {study!r} is from {result.method}, whose design has " + "random effects: its samples are not exchangeable under the null, so its " + "labels cannot be permuted sample by sample; permutations re-run " + "limma_voom and limma_trend" + ) + record = result.provenance + design = record["design"] + samples: dict[str, list[Any]] = design["samples"] + variables = [name for name in samples if name != "sample_id"] + if permute not in variables: + raise ValueError( + f"permute names {permute!r}, which is not a variable of the design of " + f"study {study!r}; its variables are: {quoted(variables) or '(none)'}" + ) + settings: dict[str, Any] = record["settings"] + parameters = inspect.signature(function).parameters + options = { + name: settings[name] + for name, parameter in parameters.items() + if parameter.kind is inspect.Parameter.KEYWORD_ONLY + and name not in ("groups", "contrasts") + and name in settings + } + if options.get("expression_filter", True) is None: + options["expression_filter"] = False # recorded as None when it was off + ids = samples["sample_id"] + rows = {sample: position for position, sample in enumerate(ids)} + table = result.samples + in_group = (table["group_id"] == group) & table["used"].astype(bool) + used = [rows[sample] for sample in table["sample_id"][in_group].tolist()] + return _Rerun( + study=study, + cube=collection.aligned_cube(study), + function=function, + options=options, + formula=design["formula"], + levels=design["levels"], + samples=samples, + weights=design["contrasts"][contrast], + used=used, + fingerprint=record["cube"]["fingerprint"], + settings=settings, + ) + + +def _check_rerun(rerun: _Rerun, result: DifferentialExpressionResult) -> None: + """Refuse a re-run on another cube, or with other settings, than the study's + analysis: its permutations would not be those of the study's results.""" + record = result.provenance + if record["cube"]["fingerprint"] != rerun.fingerprint: + raise ValueError( + f"the result of study {rerun.study!r} was computed from a cube whose " + f"fingerprint, {rerun.fingerprint}, is not that of the study's aligned " + f"cube, {record['cube']['fingerprint']}: analyse " + f"collection.aligned_cube({rerun.study!r})" + ) + if record["settings"] != rerun.settings: + changed = [ + name + for name in sorted(set(record["settings"]) | set(rerun.settings)) + if record["settings"].get(name) != rerun.settings.get(name) + ] + raise ValueError( # pragma: no cover - the options are those recorded + f"the analysis of study {rerun.study!r} could not be re-run with its " + f"recorded settings; these differ: {listing(changed)}" + ) + + +def _absence(result: DifferentialExpressionResult, group: str, contrast: str) -> str: + """Why a result has no table of a group and contrast: its own message.""" + try: + result.table(group, contrast) + except KeyError as error: + return str(error.args[0]) + raise AssertionError("a table that is not in result.tables") # pragma: no cover diff --git a/docs/source/index.md b/docs/source/index.md index c3eb678..6574942 100644 --- a/docs/source/index.md +++ b/docs/source/index.md @@ -36,7 +36,7 @@ what CytePool offers for each today. | Which cell types expand or contract? | propeller. Not yet: precision-weighted mixed models or count models of composition | [Composition tests](statistics/composition.md) | | Which pathways and regulators respond? | camera, fry, roast, mroast and romer on any contrast; gene set enrichment analysis (GSEA) of ranked genes; over-representation analysis of the genes that change; scores of pathways in every sample (GSVA, ssGSEA, combined z-scores, PLAGE); activities of transcription factors and pathways from the networks CollecTRI and PROGENy (decoupler's ULM, MLM and weighted aggregate); gene sets read from GMT files | [Gene-set tests](statistics/gene_set_tests.md) | | What drives the variation between donors, and what must be adjusted for? | Design checks for confounded factors, sample quality weights, mixed models that estimate a variance per random effect for every gene, and variance partitioning: the share of each gene's variance that each variable explains, in every cell type | [Designs, linear models and contrasts](statistics/linear_models.md), [Mixed models, REML and dream](statistics/mixed_models.md) | -| Does it replicate, and how should studies be combined? | Studies aligned by sample, gene and cell type, each analysed on its own and lined up gene by gene; meta-analysis of every gene across studies, with the between-study variance, Q and I²; combined p-values, the p-value that a change replicates in at least u studies, sign agreement, pi1 between studies and the genes found in k studies. Not yet: shrinkage across studies | [Meta-analysis](statistics/meta_analysis.md), [The cube and its additive layers](statistics/cubes.md) | +| Does it replicate, and how should studies be combined? | Studies aligned by sample, gene and cell type, each analysed on its own and lined up gene by gene; meta-analysis of every gene across studies, with the between-study variance, Q and I²; combined p-values, the p-value that a change replicates in at least u studies, sign agreement, pi1 between studies and the genes found in k studies; SumRank, which asks that a gene rank high in most studies, and the prediction of cases in held-out studies. Not yet: shrinkage across studies | [Meta-analysis](statistics/meta_analysis.md), [The cube and its additive layers](statistics/cubes.md) | | Can patients be stratified from their profiles? | Not yet; cubes convert to AnnData and MuData for other tools | [Capability matrix](reference/capabilities.md) | | Which cell-cell interactions change? | Not yet | [Related software](reference/related_software.md) | | Do effects depend on genotype, time, sex or context? | Interactions and F-tests in designs, likelihood ratio tests, random intercepts for repeated measures, with F-tests and small-sample (Kenward-Roger) tests in mixed models. Not yet: random slopes. QTL mapping is left to dedicated tools | [Designs, linear models and contrasts](statistics/linear_models.md), [Count models](statistics/count_models.md), [Mixed models, REML and dream](statistics/mixed_models.md) | diff --git a/docs/source/reference/api/held_out_prediction.md b/docs/source/reference/api/held_out_prediction.md new file mode 100644 index 0000000..a9b1fc7 --- /dev/null +++ b/docs/source/reference/api/held_out_prediction.md @@ -0,0 +1,20 @@ +--- +myst: + html_meta: + description: Held-out prediction in CytePool, as Nakatsuka et al. (2025) measure how well genes found in some studies replicate, each study left out in turn, genes chosen from the others by SumRank and by a random-effects meta-analysis, its samples scored by UCell's rank score, and the AUC of cases against controls reported by study and method. +--- + +# Held-out prediction + +{func}`~cytepool.held_out_prediction.held_out_prediction` leaves each study of a +{class}`~cytepool.cube.CubeCollection` out in turn, chooses genes up and down in each +cell type from the results of the other studies, by SumRank and by a random-effects +meta-analysis, scores every case and control of the held-out study by how high its +profiles rank those genes, as UCell scores them, and gives the area under the ROC curve +of cases against controls, by study and method, as Nakatsuka et al. (2025) compare +methods. {func}`~cytepool.plotting.plot_held_out_auc` draws the AUC of every held-out +study by method. + +```{eval-rst} +.. automodule:: cytepool.held_out_prediction +``` diff --git a/docs/source/reference/api/index.md b/docs/source/reference/api/index.md index deb65d4..ea0f072 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 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, cell-type composition, meta-analysis, combined p-values and replication across studies, SumRank and held-out prediction, and distributions. --- # CytePool API reference @@ -49,6 +49,8 @@ variance_partitioning composition meta_analysis replication +sumrank +held_out_prediction plotting distributions smoothing diff --git a/docs/source/reference/api/sumrank.md b/docs/source/reference/api/sumrank.md new file mode 100644 index 0000000..ccafe62 --- /dev/null +++ b/docs/source/reference/api/sumrank.md @@ -0,0 +1,23 @@ +--- +myst: + html_meta: + description: SumRank in CytePool, the meta-analysis of Nakatsuka et al. (2025) that ranks every gene's differential expression within each study and sums its best ranks across studies, with Irwin-Hall p-values and p-values from permutations that re-run each study's analysis, and its Manhattan plot. +--- + +# SumRank + +{func}`~cytepool.sumrank.sumrank` ranks the genes of every study by a signed statistic, +sums each gene's best normalized ranks across the studies and caps the sum at its mean, +for genes up and genes down, as SumRank does (Nakatsuka et al. 2025), with p-values from +the Irwin-Hall distribution and, given the statistics of permuted studies, from their +pooled null distribution. {func}`~cytepool.sumrank.differential_expression_sumrank` +does it for one cell type and contrast of the studies of a +{class}`~cytepool.cube.CubeCollection`, from each study's differential-expression +result, and runs the permutations itself: it permutes a variable of each study's design +among the samples of its fit, and re-runs the study's analysis on its aligned cube. +{func}`~cytepool.plotting.plot_sumrank` draws the p-value of every gene, the genes in +alphabetical order, as the authors of SumRank draw it. + +```{eval-rst} +.. automodule:: cytepool.sumrank +``` diff --git a/docs/source/reference/bibliography.md b/docs/source/reference/bibliography.md index 9a453d3..d94a018 100644 --- a/docs/source/reference/bibliography.md +++ b/docs/source/reference/bibliography.md @@ -71,6 +71,9 @@ Combined p-values, partial conjunctions and replication [@fisher1932; @stouffer1 @liptak1958; @tippett1931; @wilkinson1951; @pearson1933; @george1983; @zaykin2011; @owen2009; @benjamini2008; @heller2014; @storey2002; @storey2003; @storey2004]. +SumRank and the prediction of held-out studies [@nakatsuka2025; @irwin1927; @hall1927; +@andreatta2021; @hanley1982; @robin2011]. + ## Smoothing, distributions and numerical methods [@cleveland1979; @cleveland1981; @deboor1978; @hastie1992; @scott1979; @freedman1981; diff --git a/docs/source/reference/capabilities.md b/docs/source/reference/capabilities.md index 79859ad..64dd9b3 100644 --- a/docs/source/reference/capabilities.md +++ b/docs/source/reference/capabilities.md @@ -54,6 +54,7 @@ size; its sections give every check. | Combined p-values of several studies (Fisher, Stouffer with weights, Tippett, Wilkinson, Pearson, Mudholkar and George, Simes), one-sided or two-sided with the directions of the effects; partial conjunctions, an effect in at least u studies | {api}`combine_p_values`, {api}`partial_conjunction` | scipy 1.18.1, metap 1.1, mpmath 1.3.0 | exact | [`combined_p_values_and_replication.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/combined_p_values_and_replication.py) | | 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) | | 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) | @@ -70,7 +71,6 @@ have no CytePool function yet: - dose-response fits, benchmark doses and drug synergy; - factorization of samples across cell types; - shrinkage of effects across conditions (mashr) and meta-regression; -- cell types aligned up the Cell Ontology, and held-out prediction of case status - from the genes of several studies; +- cell types aligned up the Cell Ontology; - edgeR's quasi-likelihood F-tests and DESeq2's tests on every cell type of a cube, and DESeq2's shrunken fold changes. diff --git a/docs/source/references.bib b/docs/source/references.bib index 2188c74..513c06c 100644 --- a/docs/source/references.bib +++ b/docs/source/references.bib @@ -8,6 +8,15 @@ @book{abramowitz1964 year = {1964}, } +@article{andreatta2021, + author = {Andreatta, Massimo and Carmona, Santiago J.}, + title = {UCell: robust and scalable single-cell gene signature scoring}, + journal = {Computational and Structural Biotechnology Journal}, + year = {2021}, + volume = {19}, + pages = {3796--3798}, + doi = {10.1016/j.csbj.2021.06.043}, +} @article{anders2010, author = {Anders, Simon and Huber, Wolfgang}, title = {Differential expression analysis for sequence count data}, @@ -668,6 +677,26 @@ @article{halekoh2014 doi = {10.18637/jss.v059.i09}, } +@article{hall1927, + author = {Hall, Philip}, + title = {The distribution of means for samples of size N drawn from a population in which the variate takes values between 0 and 1, all such values being equally probable}, + journal = {Biometrika}, + year = {1927}, + volume = {19}, + number = {3-4}, + pages = {240--244}, + doi = {10.1093/biomet/19.3-4.240}, +} +@article{hanley1982, + author = {Hanley, James A. and McNeil, Barbara J.}, + title = {The meaning and use of the area under a receiver operating characteristic (ROC) curve}, + journal = {Radiology}, + year = {1982}, + volume = {143}, + number = {1}, + pages = {29--36}, + doi = {10.1148/radiology.143.1.7063747}, +} @article{hartung1999, author = {Hartung, Joachim}, title = {An alternative method for meta-analysis}, @@ -908,6 +937,16 @@ @techreport{ieee2019 doi = {10.1109/IEEESTD.2019.8766229}, } +@article{irwin1927, + author = {Irwin, J. O.}, + title = {On the frequency distribution of the means of samples from a population having any law of frequency with finite moments, with special reference to Pearson's Type II}, + journal = {Biometrika}, + year = {1927}, + volume = {19}, + number = {3-4}, + pages = {225--239}, + doi = {10.1093/biomet/19.3-4.225}, +} @misc{iso2019, author = {{ISO}}, title = {ISO 8601-1:2019, Date and time: representations for information interchange, part 1: basic rules}, @@ -1533,6 +1572,16 @@ @article{rivals2007 doi = {10.1093/bioinformatics/btl633}, } +@article{robin2011, + author = {Robin, Xavier and Turck, Natacha and Hainard, Alexandre and Tiberti, Natalia and Lisacek, Frédérique and Sanchez, Jean-Charles and Müller, Markus}, + title = {pROC: an open-source package for R and S+ to analyze and compare ROC curves}, + journal = {BMC Bioinformatics}, + year = {2011}, + volume = {12}, + number = {1}, + pages = {77}, + doi = {10.1186/1471-2105-12-77}, +} @article{robinson2008, author = {Robinson, M. D. and Smyth, G. K.}, title = {Small-sample estimation of negative binomial dispersion, with applications to SAGE data}, diff --git a/docs/source/statistics/meta_analysis.md b/docs/source/statistics/meta_analysis.md index c1df70a..bef7adc 100644 --- a/docs/source/statistics/meta_analysis.md +++ b/docs/source/statistics/meta_analysis.md @@ -1,5 +1,5 @@ --- -description: Fixed-effect and random-effects meta-analysis of every gene across studies, as metafor computes it, with the between-study variance, Cochran's Q, I-squared and the Knapp-Hartung adjustment; combined p-values, partial conjunction tests of replication, Storey's pi0 and pi1, and how well studies replicate each other. +description: Fixed-effect and random-effects meta-analysis of every gene across studies, as metafor computes it, with the between-study variance, Cochran's Q, I-squared and the Knapp-Hartung adjustment; combined p-values, partial conjunction tests of replication, Storey's pi0 and pi1, and how well studies replicate each other; SumRank, which asks that a gene rank high in most studies; and the prediction of cases in held-out studies. --- # Meta-analysis @@ -167,6 +167,99 @@ tables say how far that holds. {api}`plot_study_agreement` draws $\pi_1$ or the agreement of every pair of studies, and {api}`plot_replicated_genes` the genes that each number of studies finds. +## SumRank: genes that rank high in most studies + +A meta-analysis weighs each study by its precision, so a gene with a very small +p-value in one or two studies can pass without any sign of change in the others; such +genes predicted case status poorly in studies left out [@nakatsuka2025]. SumRank asks +instead that a gene rank high in most studies. {api}`sumrank` computes it from a signed +statistic of every gene in every study and returns a {api}`SumRankFit`; +{api}`differential_expression_sumrank` computes it from the results of one cell type +and contrast of the studies of a {api}`CubeCollection` and returns a +{api}`SumRankResult`, whose table ({api}`SumRankResult.table`) gives each gene's +statistics and p-values, up and down, and whose table of studies +({api}`SumRankResult.studies`) says which studies tested the cell type and contrast, +and why the others did not. + +Each study ranks the genes that every study tested by its signed $-\log_{10}$ +p-value, from $r = 1$ for the gene most up to $r = G$; tied genes share their average +rank. The normalized rank $u = (r - 1) / (G - 1)$ runs from 0 to 1, and a gene's +statistic sums its $m$ smallest normalized ranks over the $k$ studies, those of the +studies that rank it best, capped at their mean under the null: + +$$ +S_g = \min\Big(\sum_{j=1}^{m} u_{g(j)},\ \frac{m}{2}\Big), \qquad +m = \lfloor p\,k \rfloor, +$$ + +with $u_{g(1)} \le \dots \le u_{g(k)}$ and the proportion $p$ of the studies +(`proportion_of_studies`). The paper sums every study when there are few, and chose $p$ +by prediction in held-out studies when there are many (0.65 of 17 studies of Alzheimer's +disease, 0.55 of 16 of COVID-19). Genes down have the statistic of the ranks $1 - u$. +A gene that some study did not test, because it did not measure it or its expression +filter left it out, is not ranked, as the paper ranks the genes every study measured: +a rank says where a gene stands among the genes ranked with it. + +If each study ranked the genes at random, the normalized ranks would be independent +uniform variables, and $S_g$ below the cap would follow the Irwin-Hall distribution of a +sum of $m$ of them [@irwin1927; @hall1927], whose distribution function gives the +Irwin-Hall p-value. With every study summed it is exact under that null; the sum of the +best $m$ of $k$ ranks is smaller than a sum of $m$ ranks, so with $m < k$ these p-values +are much too small: in six simulated studies without changes, summing three, 37% of +them fell below 0.05, and {api}`differential_expression_sumrank` warns when it gives +them. The paper takes its p-values from permutations instead: each study's labels are +permuted and its analysis run again, and the statistics of every gene in every +permutation make one null distribution, $N$ statistics in all, against which a gene of +statistic $s$ has the p-value + +$$ +p = \frac{1 + \#\{\text{null statistics} \le s\}}{1 + N}, +$$ + +never zero [@phipson2010]. Pooling over the genes gives small p-values from few +permutations, 100 permutations of 10,000 genes reaching $10^{-6}$. With +`n_permutations`, {api}`differential_expression_sumrank` permutes a variable of each +study's design among the samples of its fit of the cell type, from a seeded +`numpy.random.Generator`, and re-runs the study's {api}`limma_voom` or {api}`limma_trend` +with the design and settings its result records: 0.21 s per study and permutation for +40 samples and 10,000 genes with limma-voom, so that 100 permutations of 20 such studies +take about 7 minutes on a laptop. {api}`plot_sumrank` draws every gene's p-value, the +genes in alphabetical order, as the authors of SumRank draw it. + +Where CytePool departs from the paper and its authors' code, which the validation runs: + +- The Irwin-Hall p-value is the lower tail of the distribution. The paper calls its + Irwin-Hall p-values two-sided, which would double them, and its code takes the + Irwin-Hall density at the statistic, which is not a probability; both order the genes + as the statistic does, so the permutation p-values do not depend on it. +- Permutation p-values count ties as at least as extreme and add one. The authors' code + takes the share of null statistics strictly below a gene's, which is 0 for a gene + beyond every one of them; its sums are formed in floating point, so that equal + statistics can compare as unequal, where CytePool forms them exactly. +- The labels are permuted among the samples of each study's fit, which keeps its + numbers of cases and controls; the paper permutes them among all the individuals of a + dataset and draws again until every cell type has two cases and two controls. + +## Prediction in held-out studies + +How well genes found in some studies replicate can be measured in a study that had no +part in finding them. {api}`held_out_prediction` leaves each study of a collection out +in turn, chooses genes up and genes down in each cell type from the other studies' +results, by SumRank and by the random-effects meta-analysis above with REML and the +Knapp-Hartung adjustment (the paper's comparison), scores every case and control of +the held-out study, and reports the area under the ROC curve of cases against controls +[@hanley1982], by study and method, in a {api}`HeldOutPrediction` +({api}`HeldOutPrediction.auc`). A sample's score is UCell's [@andreatta2021]: in each +cell type, the ranks of the chosen genes among all the genes of the sample's profile, +scaled to [0, 1] over the study's samples, the genes up counted for the sample and the +genes down against it, summed over the cell types, as the authors' code sums them. The genes chosen and the scores are in +{api}`HeldOutPrediction.genes` and {api}`HeldOutPrediction.scores`, and +{api}`plot_held_out_auc` draws the AUC of every study by method. The paper fits a +logistic regression on the score and lets pROC choose the direction of each study's +AUC [@robin2011], so that a score that ranks controls above cases also scores above 0.5; +here cases must score higher, and an AUC below 0.5 says the genes predict the wrong +direction. + ## Limits The studies are taken as independent, with known sampling variances. Meta-regression on @@ -176,5 +269,5 @@ assume independent studies whose p-values are uniform without effects; for discr tests they are conservative. The r-values of a primary study and its follow-up [@heller2014] are not implemented. metap's methods of Edgington (sums of p-values), Lancaster (inverse chi-squared), the inverse t, the mean z and vote counting are left -out. Held-out prediction of case status from the genes of several studies -[@nakatsuka2025] is not implemented yet. +out. SumRank weighs every study alike, as the paper notes, and its permutations re-run +designs without random effects only. diff --git a/docs/source/validation.md b/docs/source/validation.md index 1d3b245..a3661aa 100644 --- a/docs/source/validation.md +++ b/docs/source/validation.md @@ -18,8 +18,8 @@ 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` and `activity.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` and `sumrank.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`, @@ -71,10 +71,11 @@ covers and the largest difference among them. Run a script to see every check. | [`composition_tests.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/composition_tests.py) | Cell-type composition | speckle 1.6.0, limma 3.62.1 | 32 | 25 pass, 7 deliberate differences | 13.1 s | 860 MiB | | [`meta_analysis.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/meta_analysis.py) | Meta-analysis | metafor 5.0.1 | 75 | 52 pass, 20 information, 3 deliberate differences | 59.5 s | 256 MiB | | [`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 | | [`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,356** | **2,036 pass, none failed; 176 information, 142 deliberate differences, 2 controls** | | | +| **All** | | | **2,379** | **2,049 pass, none failed; 185 information, 143 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 @@ -2144,6 +2145,63 @@ the smoothing parameter: the estimates agree to about 1e-12. Deliberate differen only; CytePool keeps one at every lambda. For 96 values (0, 0.01, ..., 0.95) the estimates differed by 2.6e-5, relative. +## SumRank and held-out prediction + +[`validation/sumrank.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/sumrank.py) +runs the authors' R code of SumRank, under the Apache-2.0 licence, as they publish it: +`SumRank.R`, `CalibratePValueswithPermutations.R` and `MakingManhattanPlots.R` of their +repository at commit 46392f66, downloaded and checked by their SHA-256. `SumRank.R` loads the unifed package +only for `dirwin.hall`, the Irwin-Hall density, which it evaluates at each gene's capped +sum; the script replaces `dirwin.hall` by the identity, so that SumRank's tables hold +the sums, and the calibration, which compares -log10 of these values, orders the genes +as it would. The data (seed 20261030) are five studies of cases and controls, of 10 to +20 samples, with cell types T and B (one study without B), 300 genes (one study without +20 of them), 30 genes changed in every study and 8 changed 8-fold in one study each, +analysed by `limma_voom`; tables of 4 studies of 60 genes whose p-values take 6 values, +many tied; and three studies of 6 samples, 3 cases each, whose 20 assignments of the +labels each are all analysed. Reference: R 4.4.3, the authors' code, mpmath 1.3.0, scipy +1.18.1, scikit-learn 1.9.1, numpy 2.5.4 and numba 0.68.0. + +| Check family | Reference | Largest difference | Tolerance | Outcome | +| --- | --- | --- | --- | --- | +| The Irwin-Hall distribution function on a grid of 1 to 200 studies (118 values), and the Irwin-Hall p-values of `sumrank` for 3, 10 and 30 studies (1,199 genes) | mpmath, 200 digits | 6.2e-15 | 1e-12 | 2 pass | +| The SumRank statistics of 5 studies, T and B, up and down, summing every study, 0.65 and 0.5 of them; of 4 studies with tied p-values | the authors' SumRank | 4.2e-15 | 1e-12 | 2 pass | +| Permutation p-values of 40 permutations of the condition, up and down: (1 + the null statistics at most the gene's) / (1 + N) | the authors' null statistics | 0 | 0 | 1 pass | +| The authors' calibration, between the null statistics strictly below the gene's and those at most it, divided by N | `CalibratePValueswithPermutations` | 0 genes outside | 0 | 1 pass | +| The authors' calibrated p-values, against CytePool's | `CalibratePValueswithPermutations` | 0.14 | n/a | 1 deliberate difference | +| Permutation p-values of 1,000 random permutations, 3 and 2 of 3 studies summed | the exact distribution of every assignment of the labels, 8,000 combinations | 0.49 standard errors | 4 standard errors | 2 pass | +| The held-out score of every sample and the AUC of cases against controls, 5 studies, 2 methods | scipy's `rankdata` with UCell's formula; scikit-learn's `roc_auc_score` | 0 | 1e-12 | 2 pass | +| The genes SumRank chose for each held-out study, T and B, up and down | the authors' SumRank on the other studies | 0 of 18 lists differ | 0 | 1 pass | +| The Manhattan plot of `plot_sumrank`: every gene's -log10 p-value in alphabetical order, the threshold line and the genes named, for p-values of 20 permutations summing every study | the authors' `MakeManhattanPlot`, with `plot`, `abline` and `text` recording what it draws | 0 | 1e-12 | 2 pass | +| Calibration, what each method finds, held-out AUC (simulated studies) | the simulated changes | below | n/a | 9 information | + +The sums of the authors' code are floating-point sums of each gene's ranks in its own +order, written with 15 digits, so that equal sums can differ in their last digits; the +script compares them as the integers q of q / (2 (G - 1)), G the genes ranked, which +every sum of normalized ranks is. Its calibration then counts the null statistics +strictly below a gene's and some of those equal to it, as the rounding decides, divided +by N; CytePool counts every null statistic at most the gene's, plus one, divided by +1 + N, so that no p-value is 0 and a gene at the cap has 1, which is the deliberate +difference of 0.14. The authors' code stops when it sums the rank of one study (its +`rowSums` takes a matrix), so the held-out runs sum 0.7 of the other studies, 2 of 3 or +4. The standard error of a permutation p-value is taken as sqrt(p (1 - p) / B) for B +permutations, which bounds it whatever the correlation of the genes of a permutation. +`MakeManhattanPlot` draws its line 1e-5 above the highest gene that does not pass; +`plot_sumrank` draws it at that gene. + +The information rows, on six simulated studies of 16 samples and 1,200 genes: without +changes, the Irwin-Hall p-values summing every study fell at or below 0.05 and 0.01 in +4.8% and 0.79% of 1,400 tests, and the permutation p-values (40 permutations) in 4.4% +and 0.79%; summing the best 3 of 6, the Irwin-Hall p-values fell below them in 37% and +12%, and the permutation p-values in 4.4% and 0.50%. With 60 genes changed by a log2 +fold change of 0.5 in every study and 30 changed 8-fold in one study each, SumRank with +40 permutations found 35 of the 60 at an adjusted p-value below 0.05 and none of the 30, +the meta-analysis with REML and z-tests 35 and none, and with the Knapp-Hartung +adjustment none of either, as its t-test on 5 degrees of freedom gives no p-value small +enough for the adjustment over the genes. Chosen by SumRank, 30 genes up and 30 down +gave a held-out AUC of 1 in every study; chosen by the meta-analysis with Knapp-Hartung, +0.977 on average (0.906 to 1). + ## Figures [`validation/plotting.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/plotting.py) @@ -2338,6 +2396,11 @@ three runs after a first call that compiled it. The peak resident memory of the | Combined z-scores | 0.67 | 0.15 | | PLAGE | 1.95 | 0.96 | +`sumrank.py --benchmark` permutes the condition of 10 studies of 40 samples and 10,000 +genes 10 times: each permutation and study took 0.21 s with `limma_voom` and 0.14 s with +`limma_trend`, the analysis that a permutation runs again, on 2 threads, with a peak +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. @@ -2425,6 +2488,10 @@ SCANPY_PYTHON=$REF/python-scanpy/bin/python DECOUPLER_PYTHON=$REF/python/bin/pyt LIANA_PYTHON=$REF/python-liana/bin/python \ $REF/python-anndata/bin/python validation/anndata_conversion.py REFERENCE_PYTHON=$REF/python-scanpy/bin/python python validation/cube_collection.py +# Downloads the authors' SumRank code (13 kB) on its first run, into +# ~/.cache/cytepool-validation/sumrank unless --code or SUMRANK_CODE names another. +python validation/sumrank.py +python validation/sumrank.py --benchmark ``` Here `python` is `$REF/python/bin/python`, the environment below. With diff --git a/tests/data/sumrank/authors.npz b/tests/data/sumrank/authors.npz new file mode 100644 index 0000000..872bf6f Binary files /dev/null and b/tests/data/sumrank/authors.npz differ diff --git a/tests/test_held_out_prediction.py b/tests/test_held_out_prediction.py new file mode 100644 index 0000000..995108c --- /dev/null +++ b/tests/test_held_out_prediction.py @@ -0,0 +1,325 @@ +"""Tests of held-out prediction. + +The scores of the samples and their AUC are recomputed here from the counts with +scipy's rankdata and Mann-Whitney statistic; the genes chosen for a held-out study must +be those that sumrank and meta_analysis rank first over the other studies; and swapping +cases and controls must turn the AUC into one minus it. validation/sumrank.py checks +the AUC against scikit-learn and the genes chosen against the authors' SumRank. +""" + +import pickle +from functools import cache + +import numpy as np +import pytest +from numpy.testing import assert_allclose, assert_array_equal +from scipy.stats import mannwhitneyu, rankdata + +from cytepool.cube import Cube, CubeCollection +from cytepool.design import Design +from cytepool.differential_expression import DifferentialExpressionResult, limma_voom +from cytepool.held_out_prediction import HeldOutPrediction, held_out_prediction +from cytepool.meta_analysis import align_study_results, meta_analysis +from cytepool.sumrank import sumrank + +CONTRAST = "case_vs_ctrl" +LEVELS = {"condition": ["ctrl", "case"]} + + +def study_cube(rng: np.random.Generator, study: str, n: int, groups: list[str]) -> Cube: + genes = [f"g{j:03d}" for j in range(80)] + condition = np.array(["ctrl", "case"] * (n // 2)) + change = np.zeros(80) + change[:8], change[8:14] = 0.8, -0.8 + mean = np.exp(rng.normal(4.0, 1.0, size=(len(groups), 1, 80))) + mean = mean * np.where(condition[None, :, None] == "case", 2.0**change, 1.0) + n_cells = np.full((n, len(groups)), 30) + n_cells[0, -1] = 3 # a sample with too few cells for the last group's fit + return Cube( + {"counts": rng.negative_binomial(10, 10 / (10 + mean))}, + samples={ + "sample_id": [f"{study}{i}" for i in range(n)], + "condition": condition, + }, + groups={"group_id": groups}, + features={"feature_id": genes}, + n_cells=n_cells, + ) + + +def design(samples: object) -> Design: + return Design(samples, "~ condition", levels=LEVELS, contrasts={CONTRAST: [0, 1]}) + + +@cache +def studies() -> tuple[CubeCollection, dict[str, DifferentialExpressionResult]]: + """Four studies with groups T and B, the last without B.""" + rng = np.random.default_rng(20261031) + cubes = { + "first": study_cube(rng, "a", 10, ["T", "B"]), + "second": study_cube(rng, "b", 12, ["T", "B"]), + "third": study_cube(rng, "c", 8, ["T", "B"]), + "fourth": study_cube(rng, "d", 10, ["T"]), + } + collection = CubeCollection(cubes) + results = { + s: limma_voom( + collection.aligned_cube(s), design(collection.aligned_cube(s).samples) + ) + for s in collection.studies + } + return collection, results + + +@cache +def prediction() -> HeldOutPrediction: + collection, results = studies() + return held_out_prediction( + collection, results, CONTRAST, "condition", "case", "ctrl", n_genes=6, + proportion_of_studies=0.7, + ) # fmt: skip + + +def chosen(found: HeldOutPrediction, study: str, method: str, group: str, direction: str) -> list[str]: # fmt: skip + genes = found.genes + rows = ( + (genes["study"] == study) + & (genes["method"] == method) + & (genes["group_id"] == group) + & (genes["direction"] == direction) + ) + assert_array_equal(genes["rank"][rows], np.arange(1, rows.sum() + 1)) + return genes["feature_id"][rows].tolist() + + +def recomputed(found: HeldOutPrediction, study: str, method: str) -> dict[str, float]: + """Every sample's score, from the counts of the aligned cube.""" + collection, results = studies() + cube, samples = collection.aligned_cube(study), results[study].samples + ids = cube.samples["sample_id"].tolist() + columns = {f: i for i, f in enumerate(cube.features["feature_id"].tolist())} + total = np.zeros(len(ids)) + for group in cube.groups["group_id"].tolist(): + used = samples["sample_id"][(samples["group_id"] == group) & samples["used"]] + counts = np.asarray(cube.layer("counts", group), dtype=float) + for direction, sign in (("up", 1.0), ("down", -1.0)): + genes = [columns[g] for g in chosen(found, study, method, group, direction)] + n = len(genes) + if n == 0: + continue + raw = {} + for sample in used.tolist(): + ranks = rankdata(-counts[ids.index(sample)]) # 1 for the highest + raw[sample] = 1 - (ranks[genes].sum() - n * (n + 1) / 2) / ( + n * counts.shape[1] + ) + low, high = min(raw.values()), max(raw.values()) + spread = high - low + scaled = {s: (v - low) / spread if spread else 0.0 for s, v in raw.items()} + mean = np.mean(list(scaled.values())) + total += sign * np.array([scaled.get(s, mean) for s in ids]) + return dict(zip(ids, total.tolist(), strict=True)) + + +@pytest.mark.parametrize("method", ["sumrank", "meta_analysis"]) +def test_scores_and_auc_match_an_independent_computation(method: str) -> None: + found = prediction() + collection, _ = studies() + for study in collection.studies: + scores = found.scores.filter( + (found.scores["study"] == study) & (found.scores["method"] == method) + ) + expected = recomputed(found, study, method) + assert_allclose( + scores["score"], [expected[s] for s in scores["sample_id"].tolist()], + rtol=1e-12, atol=1e-14, + ) # fmt: skip + case = scores["case"] + statistic = mannwhitneyu(scores["score"][case], scores["score"][~case]) + row = (found.auc["study"] == study) & (found.auc["method"] == method) + n_pairs = case.sum() * (~case).sum() + assert found.auc["auc"][row][0] == pytest.approx( + statistic.statistic / n_pairs, rel=1e-12 + ) + assert found.auc["n_cases"][row][0] == case.sum() + assert found.auc["reason"][row][0] == "" + # A sample left out of a group's fit takes the mean there, so it is scored too. + assert len(found.scores) == 2 * sum( + collection.aligned_cube(s).n_samples for s in collection.studies + ) + + +def test_genes_are_chosen_from_the_other_studies() -> None: + found = prediction() + collection, results = studies() + for group in ("T", "B"): + table = align_study_results(collection, results, group, CONTRAST) + shape = (collection.n_studies, collection.n_features) + p = table["p_value"].reshape(shape) + change = table["log2_fold_change"].reshape(shape) + error = table["standard_error"].reshape(shape) + ids = np.array(collection.features["feature_id"].tolist()) + for held, study in enumerate(collection.studies): + others = [ + i for i in range(collection.n_studies) + if i != held and not np.isnan(p[i]).all() + ] # fmt: skip + statistics = np.where(change < 0, np.log10(p), -np.log10(p))[others] + fit = sumrank(statistics, proportion_of_studies=0.7) + cap = fit.n_studies_summed / 2 + for direction in ("up", "down"): + values = getattr(fit, f"{direction}_statistic") + order = np.lexsort((np.arange(values.size), values)) + best = [i for i in order if fit.ranked[i] and values[i] < cap][:6] + assert ( + chosen(found, study, "sumrank", group, direction) + == ids[best].tolist() + ) + universe = fit.ranked + combined = meta_analysis( + change[others][:, universe], error[others][:, universe], + method="reml", test="knapp_hartung", + ) # fmt: skip + positions = np.flatnonzero(universe) + for direction, sign in (("up", 1.0), ("down", -1.0)): + candidates = np.flatnonzero(sign * combined.estimate > 0) + order = candidates[ + np.lexsort((candidates, combined.p_value[candidates])) + ] + expected = ids[positions[order[:6]]].tolist() + assert ( + chosen(found, study, "meta_analysis", group, direction) == expected + ) + + +def test_the_held_out_studys_own_results_choose_none_of_its_genes() -> None: + collection, results = studies() + # The first study analysed with its labels permuted: its own genes stay, the + # others' change, as the first study is among their training studies. + cube = collection.aligned_cube("first") + labels = cube.samples["condition"][[1, 0, 3, 2, 5, 4, 7, 6, 9, 8]] + shuffled = dict(results) + shuffled["first"] = limma_voom( + cube, design({"sample_id": cube.samples["sample_id"], "condition": labels}) + ) + other = held_out_prediction( + collection, shuffled, CONTRAST, "condition", "case", "ctrl", n_genes=6, + proportion_of_studies=0.7, + ) # fmt: skip + found = prediction() + for method in ("sumrank", "meta_analysis"): + assert chosen(other, "first", method, "T", "up") == chosen( + found, "first", method, "T", "up" + ) + assert chosen(other, "second", "sumrank", "T", "up") != chosen( + found, "second", "sumrank", "T", "up" + ) + + +def test_swapping_cases_and_controls_gives_one_minus_the_auc() -> None: + collection, results = studies() + swapped = held_out_prediction( + collection, results, CONTRAST, "condition", "ctrl", "case", n_genes=6, + proportion_of_studies=0.7, + ) # fmt: skip + found = prediction() + assert_allclose(swapped.auc["auc"], 1 - found.auc["auc"], rtol=1e-13) + assert_array_equal(swapped.scores["case"], ~found.scores["case"]) + assert_array_equal(swapped.scores["score"], found.scores["score"]) + + +def test_a_study_without_controls_has_no_auc() -> None: + collection, results = studies() + found = held_out_prediction( + collection, results, CONTRAST, "condition", "case", "healthy", n_genes=3 + ) + assert np.isnan(found.auc["auc"]).all() + assert (found.auc["n_controls"] == 0).all() + assert "0 controls" in found.auc["reason"][0] + assert found.scores["case"].all() # only cases are scored + + +def test_a_group_tested_by_no_other_study_gives_no_genes() -> None: + collection, results = studies() + # Group B with the first study alone: held out, it has no other study. + kept = collection.select(studies=["first", "fourth"]) + found = held_out_prediction( + kept, + {s: results[s] for s in kept.studies}, + CONTRAST, + "condition", + "case", + "ctrl", + groups="B", + n_genes=4, + ) + row = (found.auc["study"] == "first") & (found.auc["method"] == "sumrank") + assert found.auc["n_groups"][row][0] == 0 + assert found.auc["n_genes_up"][row][0] == 0 + # Every score is then 0, and the AUC 0.5. + assert found.auc["auc"][row][0] == 0.5 + + +def test_the_result_is_a_record_that_cannot_change() -> None: + found = prediction() + assert found.settings["n_genes"] == 6 + assert found.settings["meta_analysis"] == { + "method": "reml", + "test": "knapp_hartung", + } + assert found.provenance["analysis"] == "held_out_prediction" + with pytest.raises(AttributeError, match="immutable"): + found.auc = found.scores # type: ignore[misc] + copy = pickle.loads(pickle.dumps(found)) + assert copy.auc.equals(found.auc) and copy.genes.equals(found.genes) + assert "sumrank: mean AUC" in repr(found) + with pytest.raises(TypeError, match="genes must be a Table"): + HeldOutPrediction(found.auc, found.scores, {}, {}) # type: ignore[arg-type] + + +@pytest.mark.parametrize( + ("arguments", "error", "match"), + [ + ({"case": "ctrl"}, ValueError, "both 'ctrl'"), + ({"variable": "donor"}, ValueError, "not a variable of the design"), + ({"variable": 1}, TypeError, "variable must be a string"), + ({"n_genes": 0}, ValueError, "at least 1"), + ({"n_genes": 2.5}, TypeError, "whole number"), + ({"groups": []}, ValueError, "groups is empty"), + ({"groups": ["T", "T"]}, ValueError, "more than once"), + ({"groups": "NK"}, KeyError, "no group"), + ({"proportion_of_studies": 2}, ValueError, r"in \(0, 1\]"), + ({"contrast": "other"}, KeyError, "no study tested the contrast"), + ], +) +def test_invalid_arguments_are_refused( + arguments: dict[str, object], error: type[Exception], match: str +) -> None: + collection, results = studies() + given = {"contrast": CONTRAST, "variable": "condition", "case": "case"} + given |= {"control": "ctrl"} | arguments + with pytest.raises(error, match=match): + held_out_prediction(collection, results, **given) # type: ignore[arg-type] + + +def test_the_figure_shows_every_study_and_method() -> None: + pytest.importorskip("matplotlib") + from matplotlib.figure import Figure + + from cytepool.plotting import plot_held_out_auc + + found = prediction() + ax = plot_held_out_auc(found, ax=Figure().subplots()) + studies_shown = [label.get_text() for label in ax.get_yticklabels()] + assert studies_shown == ["first", "second", "third", "fourth"] + for method in ("sumrank", "meta_analysis"): + (points,) = [c for c in ax.collections if c.get_gid() == method] + expected = found.auc["auc"][found.auc["method"] == method] + assert_allclose(points.get_offsets()[:, 0], expected) + legend = [text.get_text() for text in ax.get_legend().get_texts()] + mean = found.auc["auc"][found.auc["method"] == "sumrank"].mean() + assert legend[0] == f"SumRank (mean {mean:.2f})" + (line,) = [line for line in ax.lines if line.get_gid() == "chance"] + assert line.get_xdata()[0] == 0.5 + with pytest.raises(TypeError, match="HeldOutPrediction"): + plot_held_out_auc(found.auc) # type: ignore[arg-type] diff --git a/tests/test_masked_inputs.py b/tests/test_masked_inputs.py index 6d78296..a17e508 100644 --- a/tests/test_masked_inputs.py +++ b/tests/test_masked_inputs.py @@ -52,6 +52,7 @@ from cytepool.differential_expression import limma_voom from cytepool.gene_sets import GeneSets from cytepool.gsea import gsea +from cytepool.held_out_prediction import held_out_prediction from cytepool.io.anndata import to_anndata, to_mudata from cytepool.io.zarr import open_array, write_array, write_group from cytepool.meta_analysis import differential_expression_meta_analysis @@ -71,6 +72,7 @@ from cytepool.stats.quasi_likelihood import fit_quasi_likelihood from cytepool.stats.smoothing import lowess from cytepool.stats.voom import voom, voom_with_quality_weights +from cytepool.sumrank import differential_expression_sumrank from cytepool.table import Table from cytepool.variance_partitioning import variance_partitioning @@ -558,6 +560,55 @@ def _forest_arguments(context: Context) -> dict[str, Any]: return {"meta": meta, "studies": results, "feature": feature} +@functools.cache +def _studies() -> tuple[CubeCollection, dict[str, Any]]: + """Two studies, the bulk cube of COUNTS twice, each tested by limma_voom.""" + cube = from_bulk(COUNTS, samples=SAMPLE_TABLE) + collection = CubeCollection({"first": cube, "second": cube}) + design = Design(SAMPLE_TABLE, "~ condition", levels=LEVELS) + results = { + study: limma_voom( + collection.aligned_cube(study), + design, + contrasts={"stim": [0.0, 1.0]}, + expression_filter=False, + ) + for study in collection.studies + } + return collection, results + + +def _sumrank_arguments(context: Context) -> dict[str, Any]: + """SumRank of the two studies, with permutations and a value for every option.""" + collection, results = _studies() + return { + "collection": collection, + "results": results, + "group": "bulk", + "contrast": "stim", + "proportion_of_studies": 1.0, + "n_permutations": 2, + "permute": "condition", + "rng": 1, + } + + +def _held_out_arguments(context: Context) -> dict[str, Any]: + """Held-out prediction of the two studies, with a value for every option.""" + collection, results = _studies() + return { + "collection": collection, + "results": results, + "contrast": "stim", + "variable": "condition", + "case": "stim", + "control": "ctrl", + "groups": ["bulk"], + "n_genes": 2, + "proportion_of_studies": 1.0, + } + + TABLE = Table({"x": [1.0, 2.0], "s": ["a", "b"]}) NORM_FACTORS = [1.0, 1.1, 0.9, 1.0] LIBRARY_SIZES = COUNTS.sum(axis=1).astype(float) @@ -1985,6 +2036,29 @@ def _forest_arguments(context: Context) -> dict[str, Any]: objects=("collection", "results"), keys={"weights": "A", "null_proportion": "lambdas"}, ), + # SumRank: NaN marks a gene a study did not test. n_permutations and n_genes go + # through operator.index, which reads a 0-d masked integer. + "cytepool.sumrank.sumrank": Entry( + lambda c: { + "statistics": [[1.0, 2.0, np.nan, 0.5], [2.0, 1.0, 3.0, -1.0]], + "proportion_of_studies": 0.5, + "null_statistics": [[[0.5, 1.0, 2.0, 1.5], [1.0, 2.0, 0.0, np.nan]]], + }, + {"statistics": "NaN", "null_statistics": "NaN"}, + options=("proportion_of_studies",), + ), + "cytepool.sumrank.differential_expression_sumrank": Entry( + _sumrank_arguments, + {"n_permutations": None}, + options=("group", "contrast", "proportion_of_studies", "permute", "rng"), + objects=("collection", "results"), + ), + "cytepool.held_out_prediction.held_out_prediction": Entry( + _held_out_arguments, + {"groups": None, "n_genes": None}, + options=("contrast", "variable", "case", "control", "proportion_of_studies"), + objects=("collection", "results"), + ), # Figures, each drawn into an Axes that pyplot does not keep. A variable of the # samples is an array of labels, or a name; the other inputs are objects. "cytepool.plotting.plot_mds": Entry( @@ -2331,6 +2405,28 @@ def _forest_arguments(context: Context) -> dict[str, Any]: {}, objects=("result", "ax"), ), + "cytepool.plotting.plot_sumrank": Entry( + lambda c: { + "result": differential_expression_sumrank(*_studies(), "bulk", "stim"), + "direction": "down", + "alpha": 0.1, + "n_labels": 2, + "ax": _axes(), + }, + {}, + options=("direction", "alpha", "n_labels"), + objects=("result", "ax"), + ), + "cytepool.plotting.plot_held_out_auc": Entry( + lambda c: { + "prediction": held_out_prediction( + *_studies(), "stim", "condition", "stim", "ctrl", n_genes=2 + ), + "ax": _axes(), + }, + {}, + objects=("prediction", "ax"), + ), "cytepool.plotting.plot_variance_fractions": Entry( lambda c: { "data": _variance_partition(), @@ -2500,6 +2596,16 @@ def _forest_arguments(context: Context) -> dict[str, Any]: "Tables, which were checked when they were built, and a record kept as JSON, " "which cannot hold an array", ), + "cytepool.sumrank.SumRankResult": ( + ("table", "studies", "provenance"), + "Tables, which were checked when they were built, and a record kept as JSON, " + "which cannot hold an array", + ), + "cytepool.held_out_prediction.HeldOutPrediction": ( + ("auc", "scores", "genes", "provenance"), + "Tables, which were checked when they were built, and a record kept as JSON, " + "which cannot hold an array", + ), "cytepool.meta_analysis.align_study_results": ( ("collection", "results", "group", "contrast"), "a CubeCollection and results, which were checked when they were built, and a " @@ -2546,6 +2652,7 @@ def _forest_arguments(context: Context) -> dict[str, Any]: "cytepool.meta_analysis.MetaAnalysisFit", "cytepool.stats.combined_p_values.CombinedPValues", "cytepool.stats.multiple_testing.TrueNullProportion", + "cytepool.sumrank.SumRankFit", "cytepool.stats.dream.ModeratedFixedEffectTests", "cytepool.stats.quasi_likelihood.QuasiLikelihoodFit", "cytepool.stats.quasi_likelihood.QuasiLikelihoodTest", diff --git a/tests/test_sumrank.py b/tests/test_sumrank.py new file mode 100644 index 0000000..948126a --- /dev/null +++ b/tests/test_sumrank.py @@ -0,0 +1,678 @@ +"""Tests of SumRank. + +The reference values come from validation/sumrank.py --write-fixtures: the authors' R +code (SumRank and CalibratePValueswithPermutations, under the Apache-2.0 licence) run on +4 studies of 60 genes with tied p-values, and on 4 studies of 50 genes with 20 +permutations. Other values are worked out by hand below, or are closed forms and +invariants that a wrong ranking, sum or count would break. +""" + +import math +import pickle +from functools import cache +from pathlib import Path + +import numpy as np +import pytest +from numpy.testing import assert_allclose, assert_array_equal + +from cytepool._jit import thread_limit +from cytepool._study_statistics import studies_summed +from cytepool.cube import Cube, CubeCollection +from cytepool.design import Design +from cytepool.differential_expression import ( + DifferentialExpressionResult, + dream_voom, + limma_trend, + limma_voom, +) +from cytepool.meta_analysis import align_study_results +from cytepool.sumrank import ( + SumRankResult, + _irwin_hall_cdf, + differential_expression_sumrank, + sumrank, +) + +DATA = Path(__file__).parent / "data" / "sumrank" +CONTRAST = "case_vs_ctrl" + + +@cache +def authors() -> dict[str, np.ndarray]: + return dict(np.load(DATA / "authors.npz")) + + +def signed(p: np.ndarray, change: np.ndarray) -> np.ndarray: + """-log10 p, negative where the log2 fold change is negative.""" + return np.where(change < 0, np.log10(p), -np.log10(p)) + + +# Against the authors' code. + + +@pytest.mark.parametrize("proportion", [1, 0.5]) +def test_statistics_match_the_authors_code(proportion: float) -> None: + reference = authors() + fit = sumrank( + signed(reference["ties_p"], reference["ties_change"]), + proportion_of_studies=proportion, + ) + for direction in ("up", "down"): + # The authors' tables hold 15 significant digits. + assert_allclose( + getattr(fit, f"{direction}_statistic"), + reference[f"ties_{proportion}_{direction}_sums"], + rtol=1e-14, + ) + + +def test_permutation_p_values_count_the_authors_null_statistics() -> None: + reference = authors() + statistics = signed(reference["small_p"], reference["small_sign"]) + fit = sumrank( + statistics[0], proportion_of_studies=0.75, null_statistics=statistics[1:] + ) + assert fit.n_studies_summed == 3 + assert fit.n_permutations == 20 + assert fit.n_null_values == 20 * 50 + scale = 2 * (50 - 1) # sums are multiples of 1 / (2 (G - 1)) for G genes + for direction in ("up", "down"): + observed = np.rint(reference[f"small_{direction}_sums"] * scale) + null = np.sort(np.rint(reference[f"small_{direction}_null_sums"] * scale)) + at_most = np.searchsorted(null, observed, side="right") + below = np.searchsorted(null, observed, side="left") + ours = getattr(fit, f"{direction}_permutation_p_value") + assert_array_equal(ours, (1 + at_most) / (1 + null.size)) + # The authors' p-value counts the null statistics strictly below (and, by the + # rounding of their floating-point sums, perhaps some equal ones), without + # adding one: never larger than this one. + theirs = reference[f"small_{direction}_calibrated"] * null.size + assert np.all(theirs >= below - 1e-6) + assert np.all(theirs <= at_most + 1e-6) + assert np.all(ours > theirs / null.size) + + +# Worked out by hand. + + +def test_two_studies_by_hand() -> None: + # Study 1 ties genes b and c; study 2 did not test gene d, which is left out. + fit = sumrank([[3.0, 1.0, 1.0, 2.0, -2.0], [2.0, -1.0, 0.5, np.nan, 1.0]]) + # Over a, b, c and e (G = 4), study 1 ranks a 1, b and c 2.5, e 4, so u = (0, + # 1/2, 1/2, 1); study 2 ranks a 1, e 2, c 3, b 4, so u = (0, 1, 2/3, 1/3). + assert_array_equal(fit.ranked, [True, True, True, False, True]) + assert_allclose( + fit.normalized_ranks, + [[0, 0.5, 0.5, np.nan, 1], [0, 1, 2 / 3, np.nan, 1 / 3]], + rtol=1e-15, + ) + # Two ranks are summed and capped at 1; genes down rank 1 - u. + assert_allclose(fit.up_statistic, [0, 1, 1, np.nan, 1], rtol=1e-15) + assert_allclose(fit.down_statistic, [1, 0.5, 5 / 6, np.nan, 2 / 3], rtol=1e-15) + # F_2(x) = x^2 / 2 for x <= 1; 1 at the cap. + assert_allclose(fit.up_irwin_hall_p_value, [0, 1, 1, np.nan, 1], rtol=1e-15) + assert_allclose( + fit.down_irwin_hall_p_value, [1, 1 / 8, 25 / 72, np.nan, 2 / 9], rtol=1e-14 + ) + assert fit.up_permutation_p_value is None + assert (fit.n_studies, fit.n_studies_summed, fit.n_null_values) == (2, 2, 0) + + +def test_the_best_of_two_studies_by_hand() -> None: + fit = sumrank( + [[3.0, 1.0, 1.0, 2.0, -2.0], [2.0, -1.0, 0.5, np.nan, 1.0]], + proportion_of_studies=0.5, + ) + # One rank, the smaller, capped at 1/2: up min(u), down min(1 - u). + assert_allclose(fit.up_statistic, [0, 0.5, 0.5, np.nan, 1 / 3], rtol=1e-15) + assert_allclose(fit.down_statistic, [0.5, 0, 1 / 3, np.nan, 0], rtol=1e-15) + # F_1(x) = x. + assert_allclose(fit.up_irwin_hall_p_value, [0, 1, 1, np.nan, 1 / 3], rtol=1e-15) + assert_allclose(fit.down_irwin_hall_p_value, [1, 0, 1 / 3, np.nan, 0], rtol=1e-15) + + +def test_permutation_p_values_pool_every_gene_of_every_permutation() -> None: + # Each study ranks three genes as (3, 2, 1) or (1, 2, 3); the four combinations + # are the permutations. Their capped sums, up: (0, 1, 1), (1, 1, 1) twice and + # (1, 1, 0); so of the 12 null statistics two are 0 and ten are 1. + forward, backward = [3.0, 2.0, 1.0], [1.0, 2.0, 3.0] + null = [[a, b] for a in (forward, backward) for b in (forward, backward)] + fit = sumrank([forward, forward], null_statistics=null) + assert_allclose(fit.up_statistic, [0, 1, 1]) + assert_allclose(fit.up_permutation_p_value, [3 / 13, 1, 1], rtol=1e-15) + # Down, the data's statistics are (1, 1, 0): gene c ties the two null zeros. + assert_allclose(fit.down_permutation_p_value, [1, 1, 3 / 13], rtol=1e-15) + assert fit.n_null_values == 12 + + +def test_each_permutation_ranks_the_genes_its_studies_tested() -> None: + # The second permutation did not test gene c: it ranks two genes, and adds two + # null statistics. + statistics = [[3.0, 2.0, 1.0, 0.0], [2.0, 3.0, 1.0, 0.0]] + null = np.array([statistics, statistics]) + null[1, 0, 2] = np.nan + fit = sumrank(statistics, null_statistics=null) + assert fit.n_null_values == 4 + 3 + # Genes a and b rank 1 and 2 of 4 in one study each: they sum 0 + 1/3. The first + # permutation is the data; in the second, over a, b and d, they sum 0 + 1/2. So + # two of the 7 null statistics are at most 1/3: (1 + 2) / (1 + 7). + assert_allclose(fit.up_statistic[:2], [1 / 3, 1 / 3], rtol=1e-15) + assert_allclose(fit.up_permutation_p_value[:2], [3 / 8, 3 / 8], rtol=1e-15) + + +def test_irwin_hall_closed_forms() -> None: + x = np.linspace(0.0, 1.0, 11) + for n in (1, 2, 5, 10, 30): + # Below 1, F_n(x) = x^n / n!. + expected = x**n / math.factorial(n) + assert_allclose(_irwin_hall_cdf(x, n), expected, rtol=1e-13, atol=0) + # The distribution is symmetric about n / 2. + y = np.linspace(0.0, n, 37) + assert_allclose( + _irwin_hall_cdf(y, n) + _irwin_hall_cdf(n - y, n), 1.0, rtol=1e-13 + ) + for n in range(1, 61): + assert _irwin_hall_cdf(np.array([n / 2]), n)[0] == pytest.approx(0.5, abs=1e-14) + # Two and three uniforms at integers and halves. + assert_allclose(_irwin_hall_cdf(np.array([1.5]), 2), [1 - 0.5**2 / 2], rtol=1e-15) + assert_allclose( + _irwin_hall_cdf(np.array([1.0, 2.0]), 3), [1 / 6, 5 / 6], rtol=1e-15 + ) + assert_array_equal( + _irwin_hall_cdf(np.array([-1.0, 0.0, 3.0, 4.0]), 3), [0, 0, 1, 1] + ) + + +@pytest.mark.parametrize( + ("proportion", "n_studies", "summed"), + [ + (1, 7, 7), + (0.65, 17, 11), # the paper's 17 studies of Alzheimer's disease + (0.55, 16, 8), # and 16 of COVID-19, as the authors' code takes them + (0.57, 100, 57), # 57 although 0.57 * 100 is 56.99999999999999 + (0.1, 3, 1), # at least one + ], +) +def test_the_number_of_studies_summed( + proportion: float, n_studies: int, summed: int +) -> None: + assert studies_summed(proportion, n_studies) == summed + + +# Invariants. + + +def random_statistics(seed: int, shape: tuple[int, ...]) -> np.ndarray: + statistics = np.random.default_rng(seed).normal(size=shape) + statistics[..., 0, 5] = np.nan # a gene the first study did not test + statistics[..., 2, :3] = statistics[..., 2, 3, None] # ties + return statistics + + +def assert_fits_equal(first: object, second: object) -> None: + for name in ( + "up_statistic", + "down_statistic", + "up_irwin_hall_p_value", + "down_irwin_hall_p_value", + "up_permutation_p_value", + "down_permutation_p_value", + ): + assert_array_equal(getattr(first, name), getattr(second, name), err_msg=name) + + +@pytest.mark.parametrize("proportion", [1, 0.6]) +def test_the_order_of_studies_and_genes_does_not_matter(proportion: float) -> None: + statistics = random_statistics(1, (5, 40)) + null = random_statistics(2, (6, 5, 40)) + fit = sumrank(statistics, proportion_of_studies=proportion, null_statistics=null) + studies = [3, 0, 4, 1, 2] + genes = np.random.default_rng(3).permutation(40) + shuffled = sumrank( + statistics[studies][:, genes], + proportion_of_studies=proportion, + null_statistics=null[:, studies][:, :, genes], + ) + assert_fits_equal(shuffled, _taken(fit, genes)) + + +def _taken(fit: object, genes: np.ndarray) -> object: + class Taken: + pass + + taken = Taken() + for name in ( + "up_statistic", + "down_statistic", + "up_irwin_hall_p_value", + "down_irwin_hall_p_value", + "up_permutation_p_value", + "down_permutation_p_value", + ): + setattr(taken, name, getattr(fit, name)[genes]) + return taken + + +def test_negated_statistics_swap_up_and_down() -> None: + # Testing controls against cases negates every statistic. + statistics = random_statistics(4, (4, 30)) + null = random_statistics(5, (3, 4, 30)) + fit = sumrank(statistics, proportion_of_studies=0.75, null_statistics=null) + negated = sumrank(-statistics, proportion_of_studies=0.75, null_statistics=-null) + for name in ("statistic", "irwin_hall_p_value", "permutation_p_value"): + assert_array_equal(getattr(negated, f"up_{name}"), getattr(fit, f"down_{name}")) + assert_array_equal(getattr(negated, f"down_{name}"), getattr(fit, f"up_{name}")) + + +def test_a_monotone_transform_of_a_study_changes_nothing() -> None: + statistics = random_statistics(6, (3, 25)) + changed = statistics.copy() + changed[1] = np.exp(changed[1]) + changed[1] ** 3 # increasing + changed[2] = 10.0 * changed[2] - 4.0 + assert_fits_equal(sumrank(changed), sumrank(statistics)) + + +def test_infinite_statistics_rank_first_and_tie() -> None: + # p-values of 0 give infinite statistics: genes a and b tie at the top. + fit = sumrank([[np.inf, np.inf, 1.0, -np.inf], [5.0, 4.0, 3.0, 2.0]]) + assert_allclose(fit.normalized_ranks[0], [1 / 6, 1 / 6, 2 / 3, 1], rtol=1e-15) + assert_allclose(fit.up_statistic, [1 / 6, 1 / 2, 1, 1], rtol=1e-15) + + +def test_one_study_is_its_own_ranking() -> None: + fit = sumrank([[2.0, 0.5, -1.0, 3.0]]) + assert_allclose(fit.up_statistic, [1 / 3, 0.5, 0.5, 0.0], rtol=1e-15) + assert_allclose(fit.up_irwin_hall_p_value, [1 / 3, 1, 1, 0], rtol=1e-15) + + +def test_statistics_are_read_only() -> None: + fit = sumrank( + [[1.0, 2.0, 3.0], [3.0, 2.0, 1.0]], null_statistics=np.ones((1, 2, 3)) + ) + for name in ( + "up_statistic", + "normalized_ranks", + "ranked", + "up_permutation_p_value", + ): + with pytest.raises(ValueError, match="read-only"): + getattr(fit, name)[0] = 0 + assert repr(fit) == ( + "SumRankFit(n_features=3, ranked=3, n_studies=2, n_studies_summed=2, " + "n_permutations=1)" + ) + + +@pytest.mark.parametrize( + ("arguments", "error", "match"), + [ + ({"statistics": [1.0, 2.0]}, ValueError, r"\(n_studies, n_features\)"), + ({"statistics": np.empty((0, 3))}, ValueError, "has no study"), + ({"statistics": [["a", "b"], ["c", "d"]]}, TypeError, "real numbers"), + ({"statistics": [[1.0, 2.0], [np.nan, np.nan]]}, ValueError, "study 1 tested"), + ({"statistics": [[1.0, np.nan], [np.nan, 2.0]]}, ValueError, "0 genes were"), + ({"proportion_of_studies": 0}, ValueError, r"in \(0, 1\]"), + ({"proportion_of_studies": 1.5}, ValueError, r"in \(0, 1\]"), + ({"proportion_of_studies": "all"}, TypeError, "real number"), + ({"null_statistics": np.ones((2, 2, 3))}, ValueError, "studies and features"), + ({"null_statistics": np.ones((2, 3))}, ValueError, "n_permutations"), + ({"null_statistics": np.ones((0, 2, 2))}, ValueError, "no permutation"), + ( + {"null_statistics": [[[1.0, np.nan], [np.nan, 1.0]]]}, + ValueError, + "permutation 0: 0 genes", + ), + ], +) +def test_invalid_arguments_are_refused( + arguments: dict[str, object], error: type[Exception], match: str +) -> None: + given = {"statistics": [[1.0, 2.0], [3.0, 4.0]]} | arguments + with pytest.raises(error, match=match): + sumrank(**given) # type: ignore[arg-type] + + +# Several studies of a collection. + + +def simulated_cube( + rng: np.random.Generator, study: str, n: int, groups: list[str] +) -> Cube: + genes = [f"g{j:03d}" for j in range(60)] + condition = np.array(["ctrl", "case"] * (n // 2)) + change = np.zeros(60) + change[:6], change[6:10] = 1.5, -1.5 + mean = np.exp(rng.normal(4.0, 1.0, size=(len(groups), 1, 60))) + mean = mean * np.where(condition[None, :, None] == "case", 2.0**change, 1.0) + counts = rng.negative_binomial(10, 10 / (10 + mean)) + return Cube( + {"counts": counts}, + samples={ + "sample_id": [f"{study}{i}" for i in range(n)], + "condition": condition, + }, + groups={"group_id": groups}, + features={"feature_id": genes}, + n_cells=np.full((n, len(groups)), 30), + ) + + +def design(samples: object) -> Design: + return Design( + samples, + "~ condition", + levels={"condition": ["ctrl", "case"]}, + contrasts={CONTRAST: "conditioncase"}, + ) + + +@cache +def studies() -> tuple[CubeCollection, dict[str, DifferentialExpressionResult]]: + """Three studies, the third without group B; the second analysed by limma-trend.""" + rng = np.random.default_rng(20261030) + cubes = { + "first": simulated_cube(rng, "a", 8, ["T", "B"]), + "second": simulated_cube(rng, "b", 10, ["T", "B"]), + "third": simulated_cube(rng, "c", 8, ["T"]), + } + collection = CubeCollection(cubes) + results = {} + for study in collection.studies: + cube = collection.aligned_cube(study) + analysis = limma_trend if study == "second" else limma_voom + results[study] = analysis(cube, design(cube.samples)) + return collection, results + + +def test_studies_of_a_collection_are_ranked_as_sumrank_ranks_their_results() -> None: + collection, results = studies() + with pytest.warns(RuntimeWarning, match="best 2 of 3 studies' ranks are too small"): + result = differential_expression_sumrank( + collection, results, "T", CONTRAST, proportion_of_studies=0.7 + ) + # The signed -log10 p-values of the aligned results, as plain arrays. + table = align_study_results(collection, results, "T", CONTRAST) + shape = (collection.n_studies, collection.n_features) + p = table["p_value"].reshape(shape) + change = table["log2_fold_change"].reshape(shape) + fit = sumrank(signed(p, change), proportion_of_studies=0.7) + ranked = fit.ranked + out = result.table + assert_array_equal(out["feature_id"], collection.features["feature_id"][ranked]) + for direction in ("up", "down"): + assert_array_equal( + out[f"{direction}_statistic"], + getattr(fit, f"{direction}_statistic")[ranked], + ) + assert_array_equal( + out[f"{direction}_p_value"], + getattr(fit, f"{direction}_irwin_hall_p_value")[ranked], + ) + assert result.settings["n_studies_summed"] == 2 + assert result.settings["p_values"] == "irwin_hall" + # The genes that every study tested: those measured but left out are counted. + tested = ~np.isnan(p).any(axis=0) + assert len(out) == int(tested.sum()) + assert result.provenance["n_features_left_out"] == collection.n_features - len(out) + + +def test_a_study_without_the_group_is_left_out() -> None: + collection, results = studies() + result = differential_expression_sumrank(collection, results, "B", CONTRAST) + studies_table = result.studies + assert studies_table["analysed"].tolist() == [True, True, False] + assert "has no group 'B'" in studies_table["reason"][2] + assert studies_table["n_features_tested"][2] == 0 + assert result.settings["n_studies"] == 2 + + +def explicit_permutations( + collection: CubeCollection, + results: dict[str, DifferentialExpressionResult], + n_permutations: int, + seed: int, +) -> np.ndarray: + """The statistics of group T of the permutations that the seed draws, each + study's analysis run by hand on a design with the condition permuted among the + samples of its fit of T.""" + generator = np.random.default_rng(seed) + index = {f: i for i, f in enumerate(collection.features["feature_id"].tolist())} + null = np.full( + (n_permutations, collection.n_studies, collection.n_features), np.nan + ) + for b in range(n_permutations): + for row, study in enumerate(collection.studies): + cube, samples = collection.aligned_cube(study), results[study].samples + used = samples["sample_id"][(samples["group_id"] == "T") & samples["used"]] + ids = cube.samples["sample_id"].tolist() + condition = dict(zip(ids, cube.samples["condition"].tolist(), strict=True)) + new = dict(condition) + for i, source in enumerate(generator.permutation(len(used)).tolist()): + new[used[i]] = condition[used[source]] + analysis = limma_trend if study == "second" else limma_voom + table = analysis( + cube, + design({"sample_id": ids, "condition": [new[s] for s in ids]}), + groups=["T"], + ).table("T", CONTRAST) + positions = [index[f] for f in table["feature_id"].tolist()] + null[b, row, positions] = signed( + table["p_value"], table["log2_fold_change"] + ) + return null + + +def test_permutations_rerun_each_studys_analysis() -> None: + collection, results = studies() + result = differential_expression_sumrank( + collection, + results, + "T", + CONTRAST, + n_permutations=4, + permute="condition", + rng=7, + ) + table = align_study_results(collection, results, "T", CONTRAST) + shape = (collection.n_studies, collection.n_features) + observed = signed( + table["p_value"].reshape(shape), table["log2_fold_change"].reshape(shape) + ) + fit = sumrank( + observed, null_statistics=explicit_permutations(collection, results, 4, 7) + ) + for direction in ("up", "down"): + assert_array_equal( + result.table[f"{direction}_p_value"], + getattr(fit, f"{direction}_permutation_p_value")[fit.ranked], + ) + assert result.settings["p_values"] == "permutation" + assert result.settings["n_null_values"] == fit.n_null_values + assert result.settings["seed"] == 7 + + +def test_permutations_follow_the_seed_and_not_the_threads() -> None: + collection, results = studies() + + def p_values(rng: object, threads: int | None = None) -> np.ndarray: + with thread_limit(threads): + out = differential_expression_sumrank( + collection, results, "T", CONTRAST, n_permutations=3, + permute="condition", rng=rng, + ).table # fmt: skip + return np.concatenate([out["up_p_value"], out["down_p_value"]]) + + first = p_values(11) + assert_array_equal(p_values(11, threads=1), first) + assert_array_equal(p_values(np.random.default_rng(11)), first) + assert not np.array_equal(p_values(12), first) + + +def test_the_result_records_its_studies_and_is_immutable() -> None: + collection, results = studies() + result = differential_expression_sumrank( + collection, results, "T", CONTRAST, n_permutations=2, permute="condition", rng=1 + ) + record = result.provenance + assert record["analysis"] == "differential_expression_sumrank" + assert record["settings"]["permute"] == "condition" + assert set(record["studies"]) == set(collection.studies) + assert record["studies"]["first"] == results["first"].provenance + with pytest.raises(AttributeError, match="immutable"): + result.table = result.studies # type: ignore[misc] + copy = pickle.loads(pickle.dumps(result)) + assert copy.table.equals(result.table) and copy.provenance == record + text = repr(result) + assert "3 studies, the best 3 summed" in text + assert "2 permutations of 'condition'" in text + + +def test_permutations_need_a_variable_of_each_design() -> None: + collection, results = studies() + with pytest.raises(ValueError, match="give its name as permute"): + differential_expression_sumrank( + collection, results, "T", CONTRAST, n_permutations=2 + ) + with pytest.raises(ValueError, match="n_permutations is 0"): + differential_expression_sumrank( + collection, results, "T", CONTRAST, permute="condition" + ) + with pytest.raises(ValueError, match="not a variable of the design of study"): + differential_expression_sumrank( + collection, results, "T", CONTRAST, n_permutations=1, permute="donor" + ) + with pytest.raises(ValueError, match="at least 0"): + differential_expression_sumrank( + collection, results, "T", CONTRAST, n_permutations=-1 + ) + with pytest.raises(TypeError, match="not bool"): + differential_expression_sumrank( + collection, results, "T", CONTRAST, n_permutations=True + ) + + +def test_permutations_refuse_a_result_of_another_cube() -> None: + collection, results = studies() + # The first study's result from counts other than those of its aligned cube, with + # the same genes and groups: the line-up accepts it, a re-run would not be its own. + cube = collection.aligned_cube("first") + changed = Cube( + {"counts": cube.layer("counts") + 1}, + samples=cube.samples, + groups=cube.groups, + features=cube.features, + n_cells=cube.n_cells, + ) + other = dict(results) + other["first"] = limma_voom(changed, design(changed.samples)) + assert len(differential_expression_sumrank(collection, other, "T", CONTRAST).table) + with pytest.raises(ValueError, match="fingerprint"): + differential_expression_sumrank( + collection, other, "T", CONTRAST, n_permutations=1, permute="condition" + ) + + +def test_a_mixed_model_is_not_permuted_sample_by_sample() -> None: + rng = np.random.default_rng(3) + cube = simulated_cube(rng, "d", 8, ["T"]) + samples = cube.samples.with_columns({"donor": ["D1", "D1", "D2", "D2"] * 2}) + cube = Cube( + {"counts": cube.layer("counts")}, + samples=samples, + groups=cube.groups, + features=cube.features, + n_cells=cube.n_cells, + ) + collection = CubeCollection({"one": cube, "two": cube}) + mixed = Design( + samples, + "~ condition", + random="(1 | donor)", + levels={"condition": ["ctrl", "case"]}, + contrasts={CONTRAST: "conditioncase"}, + ) + results = { + study: dream_voom(collection.aligned_cube(study), mixed) + for study in collection.studies + } + # Without permutations dream's results are ranked; with them they are refused. + assert len( + differential_expression_sumrank(collection, results, "T", CONTRAST).table + ) + with pytest.raises(ValueError, match="random effects"): + differential_expression_sumrank( + collection, results, "T", CONTRAST, n_permutations=1, permute="condition" + ) + + +def test_invalid_studies_and_contrasts_are_refused() -> None: + collection, results = studies() + with pytest.raises(KeyError, match="no group 'NK'"): + differential_expression_sumrank(collection, results, "NK", CONTRAST) + with pytest.raises(KeyError, match="no study tested the contrast"): + differential_expression_sumrank(collection, results, "T", "other") + with pytest.raises(ValueError, match="every study of the collection"): + differential_expression_sumrank( + collection, {"first": results["first"]}, "T", CONTRAST + ) + with pytest.raises(TypeError, match="permute must be"): + differential_expression_sumrank( + collection, + results, + "T", + CONTRAST, + n_permutations=1, + permute=1, # type: ignore[arg-type] + ) + with pytest.raises(ValueError, match="non-negative seed"): + differential_expression_sumrank( + collection, results, "T", CONTRAST, n_permutations=1, + permute="condition", rng=-1, + ) # fmt: skip + + +def test_the_figure_shows_the_p_values_in_alphabetical_order() -> None: + pytest.importorskip("matplotlib") + from matplotlib.figure import Figure + + from cytepool.plotting import plot_sumrank + + collection, results = studies() + result = differential_expression_sumrank(collection, results, "T", CONTRAST) + ax = plot_sumrank(result, direction="down", n_labels=3, ax=Figure().subplots()) + table = result.table + ids = table["feature_id"].tolist() + order = np.argsort(ids, kind="stable") + # A p-value of 0 is drawn at the smallest positive float. + tiny = np.finfo(np.float64).tiny + heights = -np.log10(np.maximum(table["down_p_value"][order], tiny)) + adjusted = table["down_adjusted_p_value"][order] + shown = { + artist.get_gid(): artist.get_offsets() + for artist in ax.collections + if artist.get_gid() in ("genes", "significant") + } + passed = adjusted < 0.05 + assert passed.any() and (~passed).any() + assert_allclose(shown["significant"][:, 1], heights[passed]) + assert_allclose(shown["genes"][:, 1], heights[~passed]) + assert_allclose(shown["genes"][:, 0], np.flatnonzero(~passed)) + (line,) = [line for line in ax.lines if line.get_gid() == "threshold"] + assert line.get_ydata()[0] == pytest.approx(heights[~passed].max()) + labels = [text.get_text() for text in ax.texts if text.get_gid() == "labels"] + best = np.argsort(table["down_p_value"], kind="stable")[:3] + assert sorted(labels) == sorted(ids[i] for i in best) + with pytest.raises(ValueError, match="direction"): + plot_sumrank(result, direction="both") + with pytest.raises(TypeError, match="SumRankResult"): + plot_sumrank(table) # type: ignore[arg-type] + + +def test_a_result_is_built_from_tables_and_a_record() -> None: + collection, results = studies() + result = differential_expression_sumrank(collection, results, "T", CONTRAST) + with pytest.raises(TypeError, match="table must be a Table"): + SumRankResult({}, result.studies, {}) # type: ignore[arg-type] + with pytest.raises(ValueError, match="Out of range float values"): + SumRankResult(result.table, result.studies, {"x": float("nan")}) diff --git a/validation/sumrank.py b/validation/sumrank.py new file mode 100644 index 0000000..d089df6 --- /dev/null +++ b/validation/sumrank.py @@ -0,0 +1,779 @@ +"""Check SumRank and held-out prediction against the authors' code, mpmath, a complete +enumeration of permutations and scikit-learn. + +The authors of SumRank publish their R code under the Apache-2.0 licence +(https://github.com/nathan-nakatsuka/scRNA_Reproducibility, commit ``COMMIT``). This +script downloads ``SumRank.R``, ``CalibratePValueswithPermutations.R`` and +``MakingManhattanPlots.R`` of that commit into a directory outside the repository, +checks their SHA-256 and runs their functions in R as they are, on tables of each +study's results (``Gene``, ``p_val``, ``avg_log2FC``) that it writes. ``SumRank.R`` loads the unifed package for +``dirwin.hall``, the density of the Irwin-Hall distribution, which it evaluates at each +gene's capped sum; the R program replaces ``library`` by a function that does nothing +and ``dirwin.hall`` by the identity, so that the table SumRank writes holds the capped +sum itself. Their calibration compares -log10 of these values, which orders the genes +as -log10 of the density does (the density rises up to the cap), so its p-values are +those of the authors' code. + +Checks: + +1. The Irwin-Hall distribution function against mpmath at 200 digits (its alternating + sum) on a grid from 1e-6 to m - 0.5, for 1 to 200 studies, and the Irwin-Hall + p-values of ``sumrank`` for random statistics of 3, 10 and 30 studies. +2. SumRank statistics against the authors' SumRank: five studies of cases and controls + (seed ``SEED``), cell types T and B, one study without B and one that did not + measure 20 genes, each analysed by ``limma_voom``: ``differential_expression_sumrank`` + against SumRank on their tables, up and down (genes down with the log2 fold changes + negated, as testing controls against cases gives them), summing every study and the + proportions 0.65 and 0.5; and ``sumrank`` against SumRank on 4 studies of 60 genes + whose p-values take 6 values. +3. Permutation p-values: ``differential_expression_sumrank`` with 40 permutations of the + condition in cell type T, against the same permutations drawn here in the documented + order and analysed by ``limma_voom`` with designs built from the cubes' sample + tables, whose tables the authors' SumRank and ``CalibratePValueswithPermutations`` + take. Their sums are floating-point sums, so that equal sums can differ in their last + digits; they are compared as the integers q of q / (2 (G - 1)), G the genes ranked. + CytePool's p-value must be (1 + the null sums at most the gene's) / (1 + N); the + authors' counts those strictly below and some of those equal, as rounding decides, + without the one (a deliberate difference). +4. Permutations against a complete enumeration: three studies of 6 samples (3 cases) + and 150 genes; each study's 20 assignments of the labels analysed, and their 8,000 + combinations give the exact permutation distribution, pooled over the genes, which + 1,000 random permutations must estimate within 4 standard errors, sqrt(p (1 - p) / + 1000) (calibrated). + The authors' ``MakeManhattanPlot``, given CytePool's p-values of 20 permutations + summing every study, with ``plot``, ``abline`` and ``text`` recording what it + draws, must draw the points, the threshold (1e-5 higher) and the names of + ``plot_sumrank``. +5. Held-out prediction on the studies of 2, summing 0.7 of the other studies (2 of 3 + or 4: the authors' SumRank stops when it sums one study, as its rowSums takes a + matrix): every sample's score recomputed with scipy's ``rankdata`` from the counts + of the aligned cubes and the genes chosen, the AUC by scikit-learn's + ``roc_auc_score``, and the genes SumRank chose against the top genes of the + authors' SumRank on the other studies. +6. Information: six studies without changes, the share of Irwin-Hall and permutation + p-values at most 0.05 and 0.01, summing every study and half of them; with 60 genes + changed in every study and 30 changed 8-fold in one study only, those found at an + adjusted p-value below 0.05 by SumRank and by random-effects meta-analyses, and the + mean held-out AUC of SumRank and of the meta-analysis with Knapp-Hartung. + +Tolerances: 1e-12 relative for the Irwin-Hall distribution function and the statistics +(the authors' tables hold 15 significant digits); 0 for p-values formed from the same +counts and for the genes chosen; 1e-12 for the scores and the AUC. + +Requirements: + This script: numpy, numba, scipy, scikit-learn, mpmath and matplotlib, with + cytepool from the source tree, and access to raw.githubusercontent.com on the first + run (13 kB). R 4.4 with jsonlite: set RSCRIPT to its Rscript when it is not on PATH. + R runs with its address space limited to 8 GB. The authors' code is kept in --code, else + SUMRANK_CODE, else cytepool-validation/sumrank under XDG_CACHE_HOME or ~/.cache. + +Usage, from the repository root: + PYTHONPATH=. python validation/sumrank.py + PYTHONPATH=. python validation/sumrank.py --write-fixtures + PYTHONPATH=. python validation/sumrank.py --benchmark + +``--benchmark`` times 10 permutations of 10 studies of 40 samples and 10,000 genes, +analysed by limma_voom and by limma_trend, instead of the checks. +""" + +import argparse +import hashlib +import itertools +import json +import os +import resource +import subprocess +import sys +import tempfile +import time +import urllib.request +import warnings +from importlib.metadata import version +from pathlib import Path + +import mpmath +import numpy as np +from scipy.stats import rankdata +from sklearn.metrics import roc_auc_score + +from cytepool.cube import Cube, CubeCollection +from cytepool.design import Design +from cytepool.differential_expression import limma_trend, limma_voom +from cytepool.held_out_prediction import held_out_prediction +from cytepool.meta_analysis import differential_expression_meta_analysis +from cytepool.sumrank import _irwin_hall_cdf, differential_expression_sumrank, sumrank + +SEED = 20261030 +COMMIT = "46392f66edcbaea2614a17d90bc2409648e63951" +SOURCE = "https://raw.githubusercontent.com/nathan-nakatsuka/scRNA_Reproducibility/" +FILES = { + "SumRank.R": "78fa0054d8cfae971284733aa6fc66d3e185f4372dc4767c00518fa2f8dde163", + "CalibratePValueswithPermutations.R": ( + "d4ba83d13474c41622ce5fa656f692ccad287a45556271ab44f57a139c3b14a2" + ), + "MakingManhattanPlots.R": ( + "c370aaeeffe30264f43fe202dfbb9e7ba3962d62ec3aea82b9c1c26c38cf96aa" + ), +} +EXACT = 1e-12 +MEMORY = 8_000_000 * 1024 # bytes of address space for R +FIXTURES = Path(__file__).resolve().parent.parent / "tests" / "data" / "sumrank" +CONTRAST = "case_vs_ctrl" +DIRECTIONS = (("up", 1.0), ("down", -1.0)) +HELD_OUT = 0.7 # the proportion of the other studies of a held-out one + +R_PROGRAM = r""" +args <- commandArgs(trailingOnly = TRUE) +suppressPackageStartupMessages(library(jsonlite)) +manifest <- fromJSON(args[1], simplifyVector = FALSE) +runs <- manifest$runs +# SumRank.R loads unifed for dirwin.hall, the Irwin-Hall density at the capped sum: +# library() does nothing here, and dirwin.hall returns the sum it is given. +library <- function(...) invisible(NULL) +dirwin.hall <- function(x, n) x +source(file.path(args[2], "SumRank.R")) +# The calibration's file starts with a snippet that reads tables; take the function. +lines <- readLines(file.path(args[2], "CalibratePValueswithPermutations.R")) +start <- grep("^CalibratePValueswithPermutations <- function", lines) +eval(parse(text = lines[start:length(lines)])) +for (run in runs) { + setwd(run$directory) + datasets <- unlist(run$datasets) + presence <- data.frame(Dataset = datasets, G = unlist(run$presence)) + tables <- list() + for (suffix in unlist(run$suffixes)) { + SumRank(datasets, "G", suffix, unlist(run$genes), run$proportion, presence, + run$directory) + tables[[suffix]] <- read.table(paste0( + "G_CombinedSignedNegLogPValranksNormalized_DirwinHallPValues_Top_", + as.character(run$proportion), "_ofDatasets_", suffix, ".txt"), header = TRUE) + writeBin(as.numeric(tables[[suffix]]$P_Val), paste0("sums_", suffix, ".bin")) + } + if (length(run$null) > 0) { + comparison <- data.frame(NegLogPValue = unlist(lapply( + unlist(run$null), function(suffix) tables[[suffix]]$NegLogPValue))) + final <- CalibratePValueswithPermutations(unlist(run$genes), comparison, + tables[["Real"]]) + writeBin(as.numeric(final$PVal), "calibrated.bin") + } +} +# The authors' Manhattan plot of CytePool's p-values: plot, abline and text, which +# MakeManhattanPlot calls, record what they are given instead of drawing it. +source(file.path(args[2], "MakingManhattanPlots.R")) +drawn <- new.env() +plot <- function(x, y, ...) { drawn$x <- x; drawn$y <- y } +abline <- function(h, ...) drawn$h <- c(drawn$h, h) +text <- function(x, y, labels, ...) drawn$labels <- c(drawn$labels, as.character(labels)) +MakeManhattanPlot(read.table(manifest$manhattan, header = TRUE), + OtherNegLogPValueCutoff = 1000, TopValueCutoff = 9, Desiredggtitle = "") +setwd(args[3]) +writeBin(as.numeric(drawn$y), "manhattan_y.bin") +writeBin(as.numeric(drawn$h), "manhattan_h.bin") +writeLines(drawn$labels, "manhattan_labels.txt") +writeLines(paste("R", getRversion()), file.path(args[3], "versions.txt")) +""" + + +def fetch(directory): + """The authors' two files, downloaded once and checked by their SHA-256.""" + directory.mkdir(parents=True, exist_ok=True) + for name, checksum in FILES.items(): + path = directory / name + if not path.exists(): + print(f"downloading {SOURCE}{COMMIT}/{name} to {directory}", flush=True) + with urllib.request.urlopen(f"{SOURCE}{COMMIT}/{name}", timeout=60) as r: + path.write_bytes(r.read()) + if hashlib.sha256(path.read_bytes()).hexdigest() != checksum: + path.unlink() + sys.exit(f"{name} of commit {COMMIT} does not have the SHA-256 {checksum}") + return directory + + +def simulate(rng, sizes, n_genes, groups, changes, artefacts, lacking=(), missing=()): + """Studies of cases and controls with negative binomial counts, analysed by + limma_voom: log2 fold changes ``changes`` in every study and ``artefacts`` genes + changed 8-fold in one study only; a study in ``lacking`` has no last group, one in + ``missing`` no last 20 genes.""" + genes = [f"g{j:04d}" for j in range(n_genes)] + base = rng.lognormal(3.0, 1.3, size=(len(groups), n_genes)) + cubes, own = {}, {} + for s, n in enumerate(sizes): + study = f"s{s + 1}" + own[study] = rng.choice(np.flatnonzero(changes == 0), artefacts, replace=False) + change = changes.copy() + change[own[study]] = 3.0 + condition = np.array(["ctrl", "case"] * (n // 2)) + g = len(groups) - (study in lacking) + mean = base[:g, None, :] * rng.lognormal(0.0, 0.25, size=(g, 1, n_genes)) + mean = mean * np.where(condition[None, :, None] == "case", 2.0**change, 1.0) + counts = rng.negative_binomial(10, 10 / (10 + mean)) + kept = n_genes - 20 if study in missing else n_genes + cubes[study] = Cube( + {"counts": counts[:, :, :kept]}, + samples={"sample_id": [f"{study}_{i}" for i in range(n)], "condition": condition}, + groups={"group_id": list(groups[:g])}, + features={"feature_id": genes[:kept]}, + n_cells=np.full((n, g), 50), + ) # fmt: skip + collection = CubeCollection(cubes) + cubes = {study: collection.aligned_cube(study) for study in collection.studies} + results = { + s: limma_voom(cube, design_of(cube.samples)) for s, cube in cubes.items() + } + return collection, results, own + + +def design_of(samples): + levels = {"condition": ["ctrl", "case"]} + return Design(samples, "~ condition", levels=levels, contrasts={CONTRAST: [0, 1]}) + + +def tables_of(results, studies, group, sign): + """Each study's tested genes, p-values and signed log2 fold changes, None where it + has no table of the group.""" + found = {s: results[s].tables.get((group, CONTRAST)) for s in studies} + return { + s: None + if t is None + else (t["feature_id"].tolist(), t["p_value"], sign * t["log2_fold_change"]) + for s, t in found.items() + } + + +class Runs: + """Runs of the authors' SumRank, each in a directory of its own: ``tables`` maps a + suffix, ``"Real"`` first, to each dataset's table, None where it lacks the cell + type; the genes ranked are those of every dataset's table of the data.""" + + def __init__(self, root): + self.root, self.runs = root, {} + + def add(self, name, tables, proportion, null=()): + directory = self.root / name + directory.mkdir() + for suffix, by_dataset in tables.items(): + for dataset, table in by_dataset.items(): + if table is None: + continue + # Python floats print as the shortest text that reads back as them. + lines = [f"{g}\t{float(p)!r}\t{float(c)!r}" for g, p, c in zip(*table)] + path = directory / f"{dataset}_G_DifferentialExpression_{suffix}.txt" + path.write_text("\n".join(["Gene\tp_val\tavg_log2FC", *lines]) + "\n") + real = tables["Real"] + common = set.intersection(*(set(t[0]) for t in real.values() if t is not None)) + self.runs[name] = { + "directory": str(directory), + "datasets": list(real), + "presence": [int(t is not None) for t in real.values()], + "genes": sorted(common), + "proportion": proportion, + "suffixes": list(tables), + "null": list(null), + } + + def execute(self, code, manhattan): + """Run every run, and the authors' Manhattan plot of the table ``manhattan`` + (Gene, PVal and PVal_BH, the genes in alphabetical order).""" + lines = [f"{g}\t{float(p)!r}\t{float(q)!r}" for g, p, q in zip(*manhattan)] + path = self.root / "manhattan.tsv" + path.write_text("\n".join(["Gene\tPVal\tPVal_BH", *lines]) + "\n") + manifest = {"runs": list(self.runs.values()), "manhattan": str(path)} + (self.root / "runs.json").write_text(json.dumps(manifest)) + (self.root / "sumrank.R").write_text(R_PROGRAM) + rscript = os.environ.get("RSCRIPT", "Rscript") + files = [str(self.root / name) for name in ("sumrank.R", "runs.json")] + subprocess.run( + [rscript, "--vanilla", *files, str(code), str(self.root)], + check=True, + preexec_fn=lambda: resource.setrlimit(resource.RLIMIT_AS, (MEMORY, MEMORY)), + ) + return (self.root / "versions.txt").read_text().strip() + + def sums(self, name, suffix="Real"): + return np.fromfile(Path(self.runs[name]["directory"]) / f"sums_{suffix}.bin") + + def grid(self, name, suffix="Real"): + """A run's sums as the integers q of q / (2 (G - 1)), G the genes ranked.""" + scale = 2 * (len(self.runs[name]["genes"]) - 1) + return np.rint(self.sums(name, suffix) * scale).astype(np.int64) + + def calibrated(self, name): + return np.fromfile(Path(self.runs[name]["directory"]) / "calibrated.bin") + + +def relative(ours, reference): + ours, reference = np.asarray(ours, float), np.asarray(reference, float) + differences = np.abs(ours - reference) / np.maximum(np.abs(reference), 1e-300) + return float(np.max(np.where(ours == reference, 0.0, differences), initial=0.0)) + + +def at(table, genes, column): + """A column of a table at the given genes.""" + index = {g: i for i, g in enumerate(table["feature_id"].tolist())} + return np.asarray(table[column])[[index[g] for g in genes]] + + +def irwin_hall(): + mpmath.mp.dps = 200 + + def exact(x, m): + x = mpmath.mpf(float(x)) + terms = range(int(mpmath.floor(x)) + 1) + total = sum((-1) ** j * mpmath.binomial(m, j) * (x - j) ** m for j in terms) + return float(total / mpmath.factorial(m)) + + worst, count = 0.0, 0 + for m in (1, 2, 3, 5, 10, 20, 50, 100, 200): + fixed = [1e-6, 1e-3, 0.01, 0.5, 0.999999, 1.0, 1.000001, 1.5, 2.0] + scaled = [0.1 * m, 0.3 * m, 0.45 * m, 0.5 * m - 1e-9, 0.5 * m, 0.7 * m, m - 0.5] + grid = np.unique(np.clip(fixed + scaled, 1e-6, m - 1e-6)) + reference = np.array([exact(x, m) for x in grid]) + keep = reference > 1e-300 # below, both underflow to 0 + worst = max(worst, relative(_irwin_hall_cdf(grid[keep], m), reference[keep])) + count += int(keep.sum()) + rows = [("Irwin-Hall F_m on a grid, 1 to 200 studies", "mpmath, 200 digits", + f"max relative difference ({count} values)", worst, EXACT)] # fmt: skip + worst, count, rng = 0.0, 0, np.random.default_rng(SEED) + for k in (3, 10, 30): + fit = sumrank(rng.normal(size=(k, 400))) + for direction, _ in DIRECTIONS: + statistic = getattr(fit, f"{direction}_statistic") + p = getattr(fit, f"{direction}_irwin_hall_p_value") + below = statistic < k / 2 + reference = [exact(s, k) for s in statistic[below]] + worst = max(worst, relative(p[below], reference)) + count += int(below.sum()) + rows.append(("Irwin-Hall p-values of sumrank, 3, 10 and 30 studies", + "mpmath, 200 digits", f"max relative difference ({count} genes)", + worst, EXACT)) # fmt: skip + return rows + + +def permuted_tables(collection, results, n_permutations): + """Group T's tables of the permutations that differential_expression_sumrank draws + with the seed SEED: permutation by permutation and study by study, the condition + permuted among the samples of the study's fit of T, then limma_voom.""" + tables = {"Real": tables_of(results, collection.studies, "T", 1.0)} + generator = np.random.default_rng(SEED) + for b in range(n_permutations): + tables[f"Permutation_{b}"] = {} + for study in collection.studies: + cube, samples = collection.aligned_cube(study), results[study].samples + used = samples["sample_id"][(samples["group_id"] == "T") & samples["used"]] + ids = cube.samples["sample_id"].tolist() + condition = dict(zip(ids, cube.samples["condition"].tolist())) + new = dict(condition) + for i, source in enumerate(generator.permutation(len(used)).tolist()): + new[used[i]] = condition[used[source]] + design = design_of({"sample_id": ids, "condition": [new[s] for s in ids]}) + t = limma_voom(cube, design, groups=["T"]).tables[("T", CONTRAST)] + table = (t["feature_id"].tolist(), t["p_value"], t["log2_fold_change"]) + tables[f"Permutation_{b}"][study] = table + return tables + + +def enumeration(rng): + """Every assignment of the labels of three studies of 6 samples, against 1,000 + random permutations.""" + collection, results, _ = simulate(rng, [6] * 3, 150, ("T",), np.zeros(150), 0) + statistics = [] + for study in collection.studies: + cube = collection.aligned_cube(study) + ids, by_assignment = cube.samples["sample_id"].tolist(), [] + for cases in itertools.combinations(range(6), 3): + condition = ["case" if i in cases else "ctrl" for i in range(6)] + design = design_of({"sample_id": ids, "condition": condition}) + t = limma_voom(cube, design).tables[("T", CONTRAST)] + p, change = np.asarray(t["p_value"]), np.asarray(t["log2_fold_change"]) + signed = np.where(change < 0, np.log10(p), -np.log10(p)) + by_assignment.append(dict(zip(t["feature_id"].tolist(), signed))) + statistics.append(by_assignment) + genes = sorted(set.intersection(*(set(s[0]) for s in statistics))) + ranks = [] + for by_assignment in statistics: + values = np.array([[a[g] for g in genes] for a in by_assignment]) + ranks.append((rankdata(-values, axis=1) - 1) / (len(genes) - 1)) + rows = [] + for proportion, summed in ((1, 3), (0.7, 2)): + ours = differential_expression_sumrank( + collection, results, "T", CONTRAST, proportion_of_studies=proportion, + n_permutations=1000, permute="condition", rng=SEED + summed, + ).table # fmt: skip + worst = 0.0 + for direction, sign in DIRECTIONS: + u = [r if sign > 0 else 1 - r for r in ranks] + every = np.broadcast_arrays( + u[0][:, None, None], u[1][None, :, None], u[2][None, None, :] + ) + best = np.sort(np.stack(every), axis=0)[:summed].sum(axis=0) + null = np.sort(np.minimum(best, summed / 2).ravel()) + statistic = at(ours, genes, f"{direction}_statistic") + exact = np.searchsorted(null, statistic, side="right") / null.size + error = np.sqrt(exact * (1 - exact) / 1000) + 1 / null.size + mine = at(ours, genes, f"{direction}_p_value") + worst = max(worst, float(np.max(np.abs(mine - exact) / error))) + rows.append((f"Permutation p-values, 1,000 random permutations, {summed} of 3 " + "studies summed", "the exact distribution of the 8,000 assignments", + "max |difference| / standard error", worst, 4.0)) # fmt: skip + return rows + + +def recomputed_scores(collection, results, genes, held, method): + """The score of every sample of a held-out study, from the counts of its aligned + cube, scipy's ranks and UCell's formula, with the genes CytePool chose.""" + cube, samples = collection.aligned_cube(held), results[held].samples + ids = cube.samples["sample_id"].tolist() + features = {f: i for i, f in enumerate(cube.features["feature_id"].tolist())} + total = np.zeros(len(ids)) + for group in cube.groups["group_id"].tolist(): + used = samples["sample_id"][(samples["group_id"] == group) & samples["used"]] + counts = np.asarray(cube.layer("counts", group), float) + for direction, sign in DIRECTIONS: + chosen = genes["feature_id"][ + (genes["study"] == held) + & (genes["method"] == method) + & (genes["group_id"] == group) + & (genes["direction"] == direction) + ] + columns = [features[g] for g in chosen.tolist() if g in features] + if not columns: + continue + n, raw = len(columns), {} + for s in used.tolist(): + ranks = rankdata(-counts[ids.index(s)]) + u = ranks[columns].sum() - n * (n + 1) / 2 + raw[s] = 1 - u / (n * counts.shape[1]) + low, high = min(raw.values()), max(raw.values()) + scaled = { + s: (v - low) / (high - low) if high > low else 0.0 + for s, v in raw.items() + } + mean = np.mean(list(scaled.values())) + total += sign * np.array([scaled.get(s, mean) for s in ids]) + return dict(zip(ids, total)) + + +def held_out_rows(collection, results, prediction, runs): + genes, scores, auc = prediction.genes, prediction.scores, prediction.auc + worst_score = worst_auc = 0.0 + mismatched = compared = 0 + for held, method in itertools.product( + collection.studies, ("sumrank", "meta_analysis") + ): + expected = recomputed_scores(collection, results, genes, held, method) + mine = scores.filter((scores["study"] == held) & (scores["method"] == method)) + reference = [expected[s] for s in mine["sample_id"].tolist()] + worst_score = max(worst_score, relative(mine["score"], reference)) + value = auc["auc"][(auc["study"] == held) & (auc["method"] == method)][0] + reference_auc = roc_auc_score(mine["case"], mine["score"]) + worst_auc = max(worst_auc, abs(value - reference_auc)) + for name, run in runs.runs.items(): + parts = name.split("_") # a held-out study's runs are named s1_T_up and so on + if len(parts) != 3 or not parts[0][1:].isdigit(): + continue + held, group, direction = parts + q = runs.grid(name) + cap = int(HELD_OUT * sum(run["presence"])) * (len(run["genes"]) - 1) + keep = np.flatnonzero(q < cap) + top = [run["genes"][i] for i in keep[np.lexsort((keep, q[keep]))][:15]] + chosen = genes["feature_id"][ + (genes["study"] == held) + & (genes["method"] == "sumrank") + & (genes["group_id"] == group) + & (genes["direction"] == direction) + ] + mismatched += chosen.tolist() != top + compared += 1 + return [ + ("Held-out scores of every sample, 5 studies, 2 methods", + "scipy's rankdata, UCell's formula", "max relative difference", worst_score, EXACT), + ("Held-out AUC, 5 studies, 2 methods", "scikit-learn's roc_auc_score", + "max absolute difference", worst_auc, EXACT), + ("Genes SumRank chose for each held-out study, T and B, up and down", + "the authors' SumRank on the other studies", + f"lists that differ (of {compared})", float(mismatched), 0.0), + ] # fmt: skip + + +def manhattan_rows(result, root): + """plot_sumrank against the authors' MakeManhattanPlot of the same p-values.""" + from matplotlib.figure import Figure + + from cytepool.plotting import plot_sumrank + + ax = plot_sumrank(result, n_labels=len(result.table), ax=Figure().subplots()) + heights = np.empty(len(result.table)) + for artist in ax.collections: + offsets = artist.get_offsets() + heights[np.rint(offsets[:, 0]).astype(int)] = offsets[:, 1] + (line,) = [line for line in ax.lines if line.get_gid() == "threshold"] + names = sorted(t.get_text() for t in ax.texts if t.get_gid() == "labels") + theirs = np.fromfile(root / "manhattan_y.bin") + # MakeManhattanPlot draws its line 1e-5 above the highest gene that does not pass. + line_difference = abs( + line.get_ydata()[0] + 1e-5 - np.fromfile(root / "manhattan_h.bin")[0] + ) + same_names = names == sorted((root / "manhattan_labels.txt").read_text().split()) + return [ + ("Manhattan plot: every gene's -log10 p-value, in alphabetical order", + "the authors' MakeManhattanPlot", "max relative difference", + relative(heights, theirs), EXACT), + ("Manhattan plot: the threshold line (theirs 1e-5 higher) and the genes named", + "the authors' MakeManhattanPlot", "absolute difference; names that differ", + line_difference + (not same_names), EXACT), + ] # fmt: skip + + +def simulation(rng): + """Calibration without changes, and what each method finds with them.""" + rows = [] + collection, results, _ = simulate(rng, [16] * 6, 1200, ("T",), np.zeros(1200), 0) + for proportion in (1, 0.5): + out = differential_expression_sumrank( + collection, results, "T", CONTRAST, proportion_of_studies=proportion, + n_permutations=40, permute="condition", rng=SEED, + ) # fmt: skip + summed = out.settings["n_studies_summed"] + for source, column in ( + ("Irwin-Hall", "irwin_hall_p_value"), + ("permutation", "p_value"), + ): + p = np.concatenate([out.table[f"{d}_{column}"] for d in ("up", "down")]) + shares = ( + f"{np.mean(p <= 0.05):.4f} at 0.05, {np.mean(p <= 0.01):.4f} at 0.01" + ) + rows.append((f"No changes, {summed} of 6 studies summed: {source} p-values", + "the nominal level", f"share of {p.size:,} p-values", shares, + "information")) # fmt: skip + changes = np.zeros(1200) + changes[:30], changes[30:60] = 0.5, -0.5 + collection, results, own = simulate(rng, [16] * 6, 1200, ("T",), changes, 5) + replicated = {f"g{j:04d}" for j in range(60)} + single = {f"g{j:04d}" for genes in own.values() for j in genes} + table = differential_expression_sumrank( + collection, results, "T", CONTRAST, n_permutations=40, permute="condition", + rng=SEED, + ).table # fmt: skip + either = (table["up_adjusted_p_value"] < 0.05) | ( + table["down_adjusted_p_value"] < 0.05 + ) + found = {"SumRank, 40 permutations": set(table["feature_id"][either].tolist())} + for test in ("z", "knapp_hartung"): + meta = differential_expression_meta_analysis(results, test=test) + meta_table = meta.table("T", CONTRAST) + called = meta_table["feature_id"][meta_table["adjusted_p_value"] < 0.05] + found[f"meta-analysis, REML, {test} test"] = set(called.tolist()) + for method, genes in found.items(): + counts = f"{len(genes & replicated)}; {len(genes & single)}" + rows.append((f"60 genes changed in 6 studies, 30 in one: {method}", + "the simulated changes", "found of the 60; of the 30", counts, + "information")) # fmt: skip + prediction = held_out_prediction( + collection, results, CONTRAST, "condition", "case", "ctrl", n_genes=30 + ) + for method in ("sumrank", "meta_analysis"): + auc = prediction.auc["auc"][prediction.auc["method"] == method] + text = f"{auc.mean():.3f} ({auc.min():.3f} to {auc.max():.3f})" + rows.append((f"Held-out AUC of the 6 studies, genes chosen by {method}", + "cases against controls", "mean (smallest to largest)", text, + "information")) # fmt: skip + return rows + + +def benchmark(): + """The cost of permutations: 10 studies of 40 samples and 10,000 genes.""" + rng = np.random.default_rng(SEED) + sizes, n_genes = [40] * 10, 10_000 + collection, results, _ = simulate(rng, sizes, n_genes, ("T",), np.zeros(n_genes), 0) + for analysis in (limma_voom, limma_trend): + if analysis is limma_trend: + cubes = {s: collection.aligned_cube(s) for s in collection.studies} + results = { + s: limma_trend(c, design_of(c.samples)) for s, c in cubes.items() + } + options = {"permute": "condition", "rng": SEED} + differential_expression_sumrank( # compiles + collection, results, "T", CONTRAST, n_permutations=1, **options + ) + started = time.perf_counter() + differential_expression_sumrank( + collection, results, "T", CONTRAST, n_permutations=10, **options + ) + seconds = (time.perf_counter() - started) / (10 * len(sizes)) + print(f"{analysis.__name__}: {seconds:.3f} s per study and permutation") + + +def main(): + # The statistics of sums of the best studies are compared without permutations. + warnings.filterwarnings("ignore", "the Irwin-Hall p-values of a sum of the best") + parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + parser.add_argument("--write-fixtures", action="store_true") + parser.add_argument("--benchmark", action="store_true") + parser.add_argument("--code", type=Path) + arguments = parser.parse_args() + if arguments.benchmark: + return benchmark() + cache = Path(os.environ.get("XDG_CACHE_HOME") or Path.home() / ".cache") + default = cache / "cytepool-validation" / "sumrank" + code = fetch(arguments.code or Path(os.environ.get("SUMRANK_CODE", default))) + rows = irwin_hall() + rng = np.random.default_rng(SEED) + changes = np.zeros(300) + changes[:15], changes[15:30] = 1.0, -1.0 + studies = ([12, 16, 10, 14, 20], 300, ("T", "B"), changes, 8, ("s5",), ("s2",)) + collection, results, _ = simulate(rng, *studies) + # Tied p-values, with a log2 fold change of 0 where p is 1, as a t of 0 gives it. + ties_p = rng.choice([1.0, 0.5, 0.2, 0.05, 0.01, 1e-4], size=(4, 60)) + ties_change = rng.choice([-1.0, 1.0], size=(4, 60)) + ties_change[(ties_p == 1.0) & (rng.uniform(size=ties_p.shape) < 0.5)] = 0.0 + # For the tests: the data, then 20 permutations, each 4 studies of 50 genes. + small_p = rng.uniform(size=(21, 4, 50)) + small_sign = rng.choice([-1.0, 1.0], size=(21, 4, 50)) + small_p[0, :, :10] *= 0.01 + small_sign[0, :, :10] = 1.0 + prediction = held_out_prediction( + collection, results, CONTRAST, "condition", "case", "ctrl", n_genes=15, + proportion_of_studies=HELD_OUT, + ) # fmt: skip + proportions = (1, 0.65, 0.5) + null = [f"Permutation_{b}" for b in range(40)] + small_null = [f"Permutation_{b}" for b in range(20)] + with tempfile.TemporaryDirectory() as temporary: + runs = Runs(Path(temporary)) + permuted = permuted_tables(collection, results, len(null)) + for direction, sign in DIRECTIONS: + for group, proportion in itertools.product(("T", "B"), proportions): + tables = tables_of(results, collection.studies, group, sign) + runs.add( + f"{group}_{proportion}_{direction}", {"Real": tables}, proportion + ) + for proportion in (1, 0.5): + tables = {d: (sorted_genes(60), ties_p[i], sign * ties_change[i]) + for i, d in enumerate("abcd")} # fmt: skip + runs.add(f"ties_{proportion}_{direction}", {"Real": tables}, proportion) + signed = { + suffix: {d: (t[0], t[1], sign * t[2]) for d, t in by.items()} + for suffix, by in permuted.items() + } + runs.add(f"permutations_{direction}", signed, 0.65, null) + small = { + suffix: { + d: (sorted_genes(50), small_p[b, i], sign * small_sign[b, i]) + for i, d in enumerate("abcd") + } + for b, suffix in enumerate(["Real", *small_null]) + } + runs.add(f"small_{direction}", small, 0.75, small_null) + for held, group in itertools.product(collection.studies, ("T", "B")): + if group in collection.aligned_cube(held).groups["group_id"].tolist(): + others = [s for s in collection.studies if s != held] + tables = tables_of(results, others, group, sign) + runs.add(f"{held}_{group}_{direction}", {"Real": tables}, HELD_OUT) + # The figure of every study summed, whose p-values from 20 permutations pass + # the adjustment for some genes and are never 0 (MakeManhattanPlot jitters 0). + drawn = differential_expression_sumrank( + collection, results, "T", CONTRAST, n_permutations=20, + permute="condition", rng=SEED, + ) # fmt: skip + table = drawn.table + columns = ("feature_id", "up_p_value", "up_adjusted_p_value") + references = runs.execute(code, [table[column] for column in columns]) + rows += manhattan_rows(drawn, runs.root) + + worst = 0.0 + for group, proportion in itertools.product(("T", "B"), proportions): + ours = differential_expression_sumrank( + collection, results, group, CONTRAST, proportion_of_studies=proportion + ).table + for direction, _ in DIRECTIONS: + name = f"{group}_{proportion}_{direction}" + mine = at(ours, runs.runs[name]["genes"], f"{direction}_statistic") + worst = max(worst, relative(mine, runs.sums(name))) + rows.append(("SumRank statistics of 5 studies, T and B, up and down, " + "proportions 1, 0.65 and 0.5", "the authors' SumRank", + "max relative difference", worst, EXACT)) # fmt: skip + worst = 0.0 + signed = np.where(ties_change < 0, np.log10(ties_p), -np.log10(ties_p)) + for proportion in (1, 0.5): + fit = sumrank(signed, proportion_of_studies=proportion) + for direction, _ in DIRECTIONS: + mine = getattr(fit, f"{direction}_statistic") + worst = max( + worst, relative(mine, runs.sums(f"ties_{proportion}_{direction}")) + ) + rows.append(("sumrank of 4 studies of 60 genes with tied p-values, proportions " + "1 and 0.5", "the authors' SumRank", "max relative difference", + worst, EXACT)) # fmt: skip + + ours = differential_expression_sumrank( + collection, results, "T", CONTRAST, proportion_of_studies=0.65, + n_permutations=len(null), permute="condition", rng=SEED, + ).table # fmt: skip + differences = deliberate = 0.0 + outside = 0 + for direction, _ in DIRECTIONS: + name = f"permutations_{direction}" + q = runs.grid(name) + pooled = np.sort(np.concatenate([runs.grid(name, s) for s in null])) + at_most = np.searchsorted(pooled, q, side="right") + below = np.searchsorted(pooled, q, side="left") + mine = at(ours, runs.runs[name]["genes"], f"{direction}_p_value") + expected = (1 + at_most) / (1 + pooled.size) + differences = max(differences, float(np.max(np.abs(mine - expected)))) + theirs = runs.calibrated(name) + count = theirs * pooled.size # a count, but for the rounding of 1 - k / N + wrong = (count < below - 1e-6) | (count > at_most + 1e-6) + outside += int(np.count_nonzero(wrong)) + deliberate = max(deliberate, float(np.max(np.abs(mine - theirs)))) + rows += [ + ("Permutation p-values of 40 permutations, up and down", + "(1 + the authors' null sums at most the gene's) / (1 + N)", + "max absolute difference", differences, 0.0), + ("The authors' calibration: null sums below the gene's, and some equal " + "to it, divided by N", "CalibratePValueswithPermutations", + "genes outside these bounds", float(outside), 0.0), + ("The authors' calibrated p-values, against CytePool's", + "CalibratePValueswithPermutations", "max absolute difference", + deliberate, "deliberate difference"), + ] # fmt: skip + rows += held_out_rows(collection, results, prediction, runs) + fixtures = {"ties_p": ties_p, "ties_change": ties_change} + fixtures |= {"small_p": small_p, "small_sign": small_sign} + for name in runs.runs: + if name.startswith(("ties_", "small_")): + fixtures[f"{name}_sums"] = runs.sums(name) + if name.startswith("small_"): + sums = [runs.sums(name, suffix) for suffix in small_null] + fixtures[f"{name}_null_sums"] = np.concatenate(sums) + fixtures[f"{name}_calibrated"] = runs.calibrated(name) + rows += enumeration(rng) + rows += simulation(np.random.default_rng(SEED + 1)) + + print(f"References: {references}, the authors' code at commit {COMMIT[:8]}, " + f"mpmath {mpmath.__version__}, scikit-learn {version('scikit-learn')}, scipy " + f"{version('scipy')}; CytePool with numpy {np.__version__}, numba " + f"{version('numba')}") # fmt: skip + print("\n| check | reference | metric | value | tolerance | result |") + print("| --- | --- | --- | --- | --- | --- |") + failed = 0 + for check, reference, metric, value, tolerance in rows: + if isinstance(tolerance, str): # "information" or "deliberate difference" + result, shown = tolerance, "n/a" + else: + result = "pass" if value <= tolerance else "FAIL" + shown, failed = f"{tolerance:g}", failed + (result == "FAIL") + text = value if isinstance(value, str) else f"{value:.3g}" + print(f"| {check} | {reference} | {metric} | {text} | {shown} | {result} |") + if arguments.write_fixtures: + FIXTURES.mkdir(parents=True, exist_ok=True) + np.savez_compressed(FIXTURES / "authors.npz", **fixtures) + print(f"wrote {FIXTURES / 'authors.npz'}") + sys.exit(1 if failed else 0) + + +def sorted_genes(n): + return [f"g{j:04d}" for j in range(n)] + + +if __name__ == "__main__": + main()