Add SumRank across studies, with permutations and held-out prediction

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
 

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

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)))