SumRank

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. differential_expression_sumrank() does it for one cell type and contrast of the studies of a 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. plot_sumrank() draws the p-value of every gene, the genes in alphabetical order, as the authors of SumRank draw it.

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. sumrank() computes the statistic and its p-values from a signed statistic of every gene in every study; differential_expression_sumrank() does it for the results of the studies of a 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); 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 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 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.

sumrank() takes the statistics of the permutations as an array; 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

class cytepool.sumrank.SumRankFit(up_statistic, down_statistic, up_irwin_hall_p_value, down_irwin_hall_p_value, up_permutation_p_value, down_permutation_p_value, normalized_ranks, ranked, n_studies, n_studies_summed, proportion_of_studies, n_permutations, n_null_values)#

SumRank of every gene across studies, from 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.

Variables:

  • 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.

Explained in: Meta-analysis

class cytepool.sumrank.SumRankResult(table, studies, provenance)#

SumRank of the differential expression of the studies of a collection.

Built by 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.

Explained in: Meta-analysis

property table#

The statistics and p-values of every gene ranked.

Explained in: Meta-analysis

property studies#

Every study of the collection, whether it analysed the group and contrast, and why not.

Explained in: Meta-analysis

property settings#

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.

property provenance#

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.

cytepool.sumrank.differential_expression_sumrank(collection, results, group, contrast, *, proportion_of_studies=1.0, n_permutations=0, permute=None, rng=None)#

SumRank of one group and contrast across the studies of a collection.

Each study’s result is the analysis of its aligned cube (aligned_cube()), such as limma_voom() with the study’s own design, as for 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 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).

Parameters:

  • 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 limma_voom(), limma_trend() or 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.

Used in: Replication across studies

Explained in: Meta-analysis

cytepool.sumrank.sumrank(statistics, *, proportion_of_studies=1.0, null_statistics=None)#

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.

Parameters:

  • 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]

Explained in: Meta-analysis