Add gene set enrichment analysis of ranked genes, as fgsea

README.md

@@ -27,7 +27,7 @@ what is planned.
 | Quasi-likelihood F-tests | `fit_quasi_likelihood` and `quasi_likelihood_f_test`: edgeR 4's `glmQLFit` and `glmQLFTest` with their defaults: the common negative binomial dispersion of the most highly expressed genes by Cox-Reid adjusted profile likelihood (`common_dispersion`), deviances adjusted for small counts with exact moments of the unit deviance (`negative_binomial_unit_deviance_moments`), the average quasi-dispersion, limma's prior for unequal degrees of freedom with a trend in the average log CPM (`average_log_cpm`) or robust, prior weights and offsets, and F-tests of coefficients and named contrasts with edgeR's columns (logFC, logCPM, F, PValue, FDR) | edgeR, limma, mpmath |
 | DESeq2 | `fit_deseq2`, `deseq2_wald_test`, `deseq2_likelihood_ratio_test` and `deseq2_dispersions`: DESeq2's analysis with its defaults, from gene-wise dispersions (Cox-Reid adjusted profile likelihood), the parametric trend and the maximum a posteriori dispersions to its fits of every gene, Wald tests of coefficients and contrasts, and likelihood ratio tests of a reduced design (all the levels of a factor at once, a time course), with the columns of `results()`, its independent filtering and its Cook's distance flags | DESeq2, PyDESeq2 |
 | Empirical Bayes moderation | Moderated t and F statistics, B-statistics, TREAT, with constant, trend and robust priors, estimated as limma estimates them by default: by the method of moments, or, when the residual degrees of freedom differ between genes, by weighted moments and profile likelihood (limma's `fitFDistUnequalDF1`); either can be chosen | limma |
-| Gene-set tests | `camera` (competitive, allowing for the correlation between the genes of a set, given or estimated per set; parametric or rank-based) and `fry` (self-contained, directional and mixed p-values), for one contrast of a linear model, with voom or other precision weights and trend or robust moderation; the rotation tests `roast` and `mroast` (self-contained, with the mean, floor-mean, mean-50 and mean-square set statistics, contribution weights per gene or per set, and the proportions of active genes) and `romer` (competitive, by the ranks of the genes), with rotations from a seeded generator; gene sets named by gene id or position. Gene sets read from and written to GMT files (MSigDB, Reactome, WikiPathways, the Gene Ontology), with the source, version and licence of the collection, matched to an analysis and filtered by size; over-representation analysis of the genes up, down or either way in a differential-expression result, a table or a list of genes (limma's `kegga` and `goana`), with its dot plot | limma, R, decoupler |
+| Gene-set tests | `camera` (competitive, allowing for the correlation between the genes of a set, given or estimated per set; parametric or rank-based) and `fry` (self-contained, directional and mixed p-values), for one contrast of a linear model, with voom or other precision weights and trend or robust moderation; the rotation tests `roast` and `mroast` (self-contained, with the mean, floor-mean, mean-50 and mean-square set statistics, contribution weights per gene or per set, and the proportions of active genes) and `romer` (competitive, by the ranks of the genes), with rotations from a seeded generator; gene sets named by gene id or position. Gene sets read from and written to GMT files (MSigDB, Reactome, WikiPathways, the Gene Ontology), with the source, version and licence of the collection, matched to an analysis and filtered by size; over-representation analysis of the genes up, down or either way in a differential-expression result, a table or a list of genes (limma's `kegga` and `goana`), with its dot plot; gene set enrichment analysis (`gsea`) of genes ranked by a statistic, as fgsea computes it: enrichment scores, leading edges, normalized scores, and p-values by fgsea's multilevel method or by permutations, from seeds that a numpy generator draws, with its enrichment plot and table | limma, R, decoupler, fgsea |
 | Multiple testing | Bonferroni, Holm, Hochberg, Hommel, Benjamini-Hochberg, Benjamini-Yekutieli | R |
 | Differential expression | `limma_voom` and `limma_trend`: limma-voom and limma-trend on every group (cell type) of a cube, and `dream_voom`: dream, with the design's random intercepts, such as a donor's, in every group's fits, dream's choice of Kenward-Roger's or Satterthwaite's degrees of freedom by the number of samples in each group, and F-tests of several contrasts, with the design matched to the cube by sample id, samples with too few cells left out, rank-deficient groups fitted on their column space with the contrasts they can estimate, and a provenance record for every result. Options for pseudobulk data, off by default: sample quality weights, prior weights from the number of cells, a minimum total count per sample, and p-values adjusted across cell types | limma, edgeR, variancePartition; on simulated data and on published single-cell data |
 | 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 |

THIRD_PARTY_NOTICES.md

@@ -47,6 +47,59 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 SOFTWARE.
 ```
 
+### fgsea
+
+`cytepool/gsea.py` ports the algorithms of fgsea 1.32.2
+(https://bioconductor.org/packages/fgsea, https://github.com/ctlab/fgsea), keeping their
+structure, constants and order of operations, so that given fgsea's seeds it gives
+fgsea's results: from `src/fastGSEA.cpp`, the block prefix sums (`SegmentTree`), the
+scores of all the first genes of a random sample from the upper convex hulls of their
+walk (`gseaStats1`, `calcGseaStatCumulative`, with the hull tolerance 1e-13 and the floor
+of the statistics, the smaller of 1e-5 and the smallest nonzero one, divided by 1024) and
+the counts of `calcGseaStatCumulativeBatch`; from `src/fgseaMultilevelSupplement.cpp`,
+`src/fgseaMultilevel.cpp`, `src/esCalculation.cpp` and `src/util.cpp`, the multilevel
+method (`EsRuler`: `extend`, `duplicateSamples`, `perturbate` with its chunks, its
+remembered gene and its swaps of a tenth of a sample, `getPvalue`,
+`calcLogCorrection`), `calcES`, `calcPositiveES`, the sampler of distinct genes
+(`combination`, with its 100 draws per gene), the random integers of `uid_wrapper`,
+`betaMeanLog` and `fgseaMultilevelCpp`; and from `R/fgsea.R`, `R/fgseaMultilevel.R` and
+`R/util.R`, the enrichment score and leading edge of `calcGseaStat`, the normalized
+scores and p-values of `fgseaSimpleImpl`, the chunks of permutations of `fgseaSimple`,
+and the rules of `fgseaMultilevel` (sets with fewer than 10 random scores of their sign,
+the choice of the multilevel method by its expected error, `multilevelError`, the
+denominator of its p-values and `eps`). `cytepool/_ranked_genes.py` computes the walk of
+`calcGseaStat`, and `cytepool/plotting.py` draws `plotEnrichment` from the data of
+`plotEnrichmentData` and the rows of `plotGseaTable` with its numbers written as its
+`valueToExpExpression` writes them. `validation/gsea.py` states fgseaMultilevel's choice
+of the multilevel method in R, to compare CytePool's with it. fgsea is described in
+Korotkevich G, Sukhov V, Budin N, Shpak B, Artyomov MN, Sergushichev A (2021). Fast gene
+set enrichment analysis. bioRxiv 060012, version 3. doi:10.1101/060012. It is
+distributed under the MIT License (its `DESCRIPTION` says "MIT + file LICENCE"), as its
+`LICENCE` file states:
+
+```text
+YEAR: 2016-2019
+COPYRIGHT HOLDER: Alexey Sergushichev
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+```
+
 ### numba
 
 `cytepool/_jit.py` places the cache of its compiled routines with subclasses of numba's

cytepool/_options.py

@@ -1,10 +1,11 @@
-"""Checks of the options of public functions: flags, choices among names, numbers.
+"""Checks of the options of public functions: flags, choices among names, numbers,
+random generators.
 
 Each check takes the option as the caller gave it and its name, and returns it as a
-bool, str, float or int, or raises TypeError for a value of the wrong type and
-ValueError for one out of range. The messages name the option, say what it must be
-and what it was; where the methods word a message differently, the check takes the
-phrase as an argument, so that every message stays as it was.
+bool, str, float, int or numpy random generator, or raises TypeError for a value of
+the wrong type and ValueError for one out of range. The messages name the option, say
+what it must be and what it was; where the methods word a message differently, the
+check takes the phrase as an argument, so that every message stays as it was.
 """
 
 import math
@@ -139,3 +140,31 @@ def index_integer(value: Any, name: str, *, expected: str) -> int:
         raise TypeError(
             f"{name} must be {expected}, got {type(value).__name__}"
         ) from None
+
+
+def random_generator(rng: Any, name: str = "rng") -> np.random.Generator:
+    """A numpy random generator: ``rng`` itself, one seeded by it, or a fresh one.
+
+    Args:
+        rng: A ``numpy.random.Generator``, which is used as it is; a non-negative
+            integer, which seeds a new one; or None for a fresh, unseeded generator.
+        name: The option's name, for the message.
+
+    Raises:
+        TypeError: If ``rng`` is none of these, or is a numpy masked array.
+        ValueError: If the seed is negative.
+    """
+    reject_masked(rng, name)
+    if rng is None:
+        return np.random.default_rng()
+    if isinstance(rng, np.random.Generator):
+        return rng
+    if isinstance(rng, bool | np.bool_) or not isinstance(rng, numbers.Integral):
+        raise TypeError(
+            f"{name} must be a numpy.random.Generator, an integer seed or None, got "
+            f"{type(rng).__name__}"
+        )
+    seed = operator.index(rng)
+    if seed < 0:
+        raise ValueError(f"{name} must be a non-negative seed, got {seed}")
+    return np.random.default_rng(seed)