Skip to content

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

Merged
parashardhapola merged 1 commit into
masterfrom
50-sumrank
Oct 26, 2026
Merged

parashardhapola merged 1 commit into
masterfrom
50-sumrank

Conversation

@parashardhapola

Copy link
Copy Markdown
Member

Summary

Most genes that one single-cell study finds changed are not found by the next, and a meta-analysis that weighs studies by precision can pass a gene on one or two studies alone. SumRank (Nakatsuka et al. 2025) asks instead that a gene rank high in most studies: each study ranks its genes, and a gene's best normalized ranks are summed across the studies. This adds SumRank on plain arrays and on the studies of a CubeCollection, with Irwin-Hall p-values and p-values from permutations that re-run each study's analysis; the paper's measure of replication, the AUC of cases against controls in each held-out study scored by genes chosen from the others (by SumRank and, for comparison, by the meta-analysis already built); and the two figures. It is checked against the authors' own R code, run as they publish it, mpmath, a complete enumeration of permutations and scikit-learn.

Changes

  • cytepool.sumrank.sumrank(statistics, *, proportion_of_studies, null_statistics): studies x genes signed statistics (NaN where a study did not test a gene), and optionally the statistics of permuted studies; returns SumRankFit: the statistics up and down, Irwin-Hall p-values, permutation p-values, the normalized ranks, and the numbers of studies, ranks summed and null statistics.
  • cytepool.sumrank.differential_expression_sumrank(collection, results, group, contrast, *, proportion_of_studies, n_permutations, permute, rng): SumRank of one cell type and contrast from each study's limma_voom, limma_trend or dream_voom result, read through align_study_results; with n_permutations, permutes the design variable permute among the samples of each study's fit and re-runs its analysis; returns SumRankResult (a table of the genes, a table of the studies, settings, provenance). Warns when fewer than all studies are summed without permutations.
  • cytepool.held_out_prediction.held_out_prediction(collection, results, contrast, variable, case, control, *, groups, n_genes, proportion_of_studies): leaves each study out, chooses n_genes genes up and down per cell type from the others by SumRank and by a REML meta-analysis with Knapp-Hartung, scores every case and control by UCell's rank score, and returns HeldOutPrediction with the tables auc (by study and method), scores and genes.
  • cytepool.plotting.plot_sumrank: the authors' Manhattan plot, the -log10 p-value of every gene in alphabetical order, the genes that pass coloured and named, with the threshold where they draw it; cytepool.plotting.plot_held_out_auc: the AUC of every held-out study by method, with the mean of each.
  • cytepool/_study_statistics.py (private): each study's signed statistics, log2 fold changes and standard errors as arrays, and the number of studies summed.
  • Documentation: sections "SumRank: genes that rank high in most studies" and "Prediction in held-out studies" in How the statistics work: Meta-analysis; API pages sumrank and held_out_prediction; a row of the capability matrix; five references (Irwin 1927, Hall 1927, Andreatta and Carmona 2021, Hanley and McNeil 1982, Robin et al. 2011; DOIs checked against Crossref) cited on the Bibliography page; the Validation page (row, totals, section, timing, how to run); the README's status table; the Overview.

Design notes

The authors' code (https://github.com/nathan-nakatsuka/scRNA_Reproducibility) is under the Apache-2.0 licence, so it was read, and its conventions settle what the paper leaves open; no code was copied. THIRD_PARTY_NOTICES.md needs an entry for it, proposed at the end.

What the paper and the authors' code specify, and CytePool follows:

  • Each study ranks the genes by signed -log10 p (negative where the log2 fold change is negative), ties by their average rank (R's rank); u = (r - 1) / (G - 1); the m smallest of a gene's k normalized ranks are summed, m = floor(p k) for the proportion p of the studies; a sum above m / 2 is set to m / 2. Genes down are tested by testing controls against cases: here their normalized ranks are 1 - u, the same.
  • Only genes that every study tested are ranked (the authors' GetCommonGenes), and a study without the cell type is left out, k counting the others.
  • The permutation null pools the statistics of every gene in every permutation, each permutation re-running every study's analysis.
  • Held-out prediction: SumRank on the other studies, the same number of genes per cell type and direction for every method, UCell scores with maxRank above the number of genes, each cell type's scores scaled to [0, 1] within the held-out study, a missing cell type at the mean, genes up minus genes down summed over cell types, and the AUC of cases against controls.

What it leaves open, and how it is settled here:

  • The Irwin-Hall p-value is P(S <= s), F_m(s), below the cap and 1 at it. It is exact for m = k under independent uniform ranks, and too small for m < k: in the validation's six null studies, summing 3 of 6, 37% of the Irwin-Hall p-values fell below 0.05 (4.8% summing all six), and the pooled permutation p-values were at 4.4%. So the Irwin-Hall p-value is anti-conservative over a subset of the studies, as expected, but not conservative over all of them: the paper's two-sided form, which doubles it, would be. The paper's final p-values are the permutation ones, and so are CytePool's when permutations run. F_m is computed as a sum of positive B-spline values (de Boor's recurrence), which keeps 1e-14 relative accuracy for 200 studies where the alternating sum cancels.
  • Ranks are formed exactly: twice a rank minus one is an integer, so sums are compared in integers and ties between the data and the permutations are exact.
  • Genes a study did not measure or did not test are left out, as the intersection the paper uses: a rank is relative to the genes ranked with it, and a gene of fewer studies would have a null of its own.
  • Permutations shuffle the values of one design variable among the samples of each study's fit of the cell type, drawn from a numpy.random.Generator permutation by permutation and study by study, and re-run limma_voom or limma_trend from the design and settings each result records, on the study's aligned cube; the result's cube fingerprint and settings are checked against the re-run. A permutation's statistics use the genes its analyses tested. Designs with random effects (dream) are not permuted: their samples are not exchangeable.
  • Cost: 0.21 s per study and permutation with limma_voom for 40 samples and 10,000 genes (0.14 s with limma_trend), so 100 permutations of 20 studies take about 7 minutes on a laptop; memory is that of one analysis and a few arrays of the genes. Pooling over genes, 100 permutations of 10,000 genes give p-values down to 1e-6. Held-out prediction of 10 studies, 2 cell types and 10,000 genes takes 1.4 s.
  • Genes are chosen by the SumRank statistic, which orders them as their permutation p-values would (every gene has the same k and m), so held-out prediction needs no permutations. The meta-analysis chooses by p-value among positive (or negative) estimates.

Deliberate differences

  • The Irwin-Hall p-value is the lower tail. The paper calls its p-values two-sided; the authors' code evaluates dirwin.hall, the Irwin-Hall density, at the sum, which is not a probability and is constant for m = 1. Both are monotone in the sum for m >= 2, so the permutation p-values are unaffected.
  • Permutation p-values are (1 + null statistics at most s) / (1 + N), never 0, and 1 at the cap. The authors' code takes the share strictly below s, 0 for a gene beyond every null statistic, and its floating-point sums can split equal statistics; the largest difference in the validation is 0.14.
  • Labels are permuted among the samples of each study's fit, keeping its numbers of cases and controls; the paper permutes all individuals of a dataset and redraws until every cell type has two of each.
  • The held-out AUC has a fixed direction, cases above controls; the paper's logistic regression with pROC's default direction = "auto" picks the direction per study, which lifts a reversed score above 0.5.
  • floor(p k) is taken of the product rounded to 9 decimals (0.57 of 100 studies is 57, where the authors' code gives 56), and m is at least 1 (their code stops for m = 1).

Validation

validation/sumrank.py, 23 checks: 13 pass, 1 deliberate difference, 9 information, none failed (44.9 s, 399 MiB peak). It downloads SumRank.R, CalibratePValueswithPermutations.R and MakingManhattanPlots.R at commit 46392f66 (SHA-256 checked) and runs them as published, with library stubbed (they load unifed only for dirwin.hall) and dirwin.hall replaced by the identity, so that their tables hold the capped sums; for the figure, plot, abline and text record what MakeManhattanPlot draws. Reference: R 4.4.3, the authors' code, mpmath 1.3.0, scipy 1.18.1, scikit-learn 1.9.1.

Check family Reference Largest difference Tolerance Outcome
Irwin-Hall distribution function, 1 to 200 studies, on a grid and at sumrank's statistics mpmath, 200 digits 6.2e-15 1e-12 2 pass
SumRank statistics of 5 limma-voom studies (2 cell types, one study without one, one without 20 genes), up and down, proportions 1, 0.65, 0.5; 4 studies with tied p-values the authors' SumRank 4.2e-15 1e-12 2 pass
Permutation p-values of 40 permutations, against the authors' null statistics of the same permutations re-analysed independently; the authors' calibration between its two counts of ties SumRank, CalibratePValueswithPermutations 0 0 2 pass
The authors' calibrated p-values CalibratePValueswithPermutations 0.14 n/a 1 deliberate difference
1,000 random permutations against all 8,000 label assignments of three 6-sample studies, 3 and 2 of 3 summed exact enumeration 0.49 standard errors 4 2 pass
Held-out scores and AUC (5 studies, 2 methods); the genes SumRank chose for each held-out study scipy's rankdata, scikit-learn's roc_auc_score, the authors' SumRank 0; 0 of 18 lists 1e-12; 0 3 pass
plot_sumrank: every gene's -log10 p-value in alphabetical order, the threshold line (theirs 1e-5 higher) and the genes named the authors' MakeManhattanPlot on the same p-values 0 1e-12 2 pass
Null calibration (Irwin-Hall and permutation, 6 of 6 and 3 of 6 summed), genes found (SumRank 35 of 60 replicated genes and none of 30 single-study ones; REML z-test 35 and none; Knapp-Hartung none), held-out AUC (SumRank 1.000, meta-analysis 0.977) simulation n/a 9 information
REF=/home/parashar/.cache/cytepool-reference
RSCRIPT=$REF/r/bin/Rscript PYTHONPATH=. $REF/python/bin/python validation/sumrank.py
RSCRIPT=$REF/r/bin/Rscript PYTHONPATH=. $REF/python/bin/python validation/sumrank.py --benchmark

Tests

tests/test_sumrank.py (43) and tests/test_held_out_prediction.py (19): 62 tests, 2.5 s warm, 73 s cold with an empty compile cache. They compare with the authors' code (fixtures in tests/data/sumrank/authors.npz, 41 KB: tied p-values; 4 studies with 20 permutations, whose permutation p-values must equal the count over the authors' null exactly); with values derived by hand (two studies with ties and a gene one study lacks, one of two studies summed, the four permutations of two three-gene studies, a permutation that tested fewer genes, infinite statistics, one study); with closed forms of the Irwin-Hall distribution (x^n / n!, symmetry, F_n(n/2) = 1/2); and with invariants: the order of studies and genes, negated statistics swapping up and down, a monotone transform of a study, explicit permutations analysed by hand against the driver, seeds reproducing permutations and threads not changing them, the held-out study's own results choosing none of its genes, and swapping cases and controls giving 1 - AUC. Held-out scores and AUC are recomputed with scipy (rankdata, mannwhitneyu), and the chosen genes with sumrank and meta_analysis on the other studies. Messages, the refusal of dream for permutations, the check of a result from another cube, the figures' data, and the masked-input inventory (38 tests) are covered.

Checks: the full suite on master, 4,816 passed, 18 skipped, also with bounds checks and with the rounding changed (OPENBLAS_CORETYPE=Sandybridge, 4 threads, NUMBA_CPU_NAME=generic, fresh cache); Python 3.14: 4,837 passed; Python 3.15: 4,731 passed, 107 skipped; ruff check and ruff format --check on cytepool, tests, validation and docs clean; mypy cytepool clean; make -C docs api and make -C docs html (Node.js 24) build; tests/test_documentation.py 74 passed.

THIRD_PARTY_NOTICES.md gains a SumRank entry under "Adapted material": the conventions
taken from the authors' Apache-2.0 code, file by file, and that no code is copied.

SumRank ranks each gene's differential expression within every study of a
CubeCollection and sums its best ranks, with Irwin-Hall p-values and a null
pooled over permutations that re-run each study's limma analysis. Held-out
prediction reports each study's AUC from genes chosen in the other studies,
by SumRank and by meta-analysis. plot_sumrank and plot_held_out_auc draw
the authors' Manhattan plot and the AUC by study.

Checked against the authors' R code (Apache-2.0), mpmath, scipy and
scikit-learn (validation/sumrank.py: 23 checks, none failed).
@parashardhapola
parashardhapola merged commit baaadaf into master Oct 26, 2026
6 checks passed
@parashardhapola
parashardhapola deleted the 50-sumrank branch October 26, 2026 21:16
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant