[parashardhapola]

Summary

Gene set enrichment analysis (GSEA) is the pathway method most users know: it ranks
every gene by its change and asks whether a set's genes gather at the top or at the
bottom of the ranking (Subramanian et al. 2005). This adds preranked GSEA as fgsea
1.32.2 computes it (Korotkevich et al. 2021): enrichment scores and leading edges,
normalized scores, and p-values by fgsea's default multilevel method or by
permutations, from a differential-expression result, any table of results or a mapping
of gene id to statistic, with the two figures fgsea users expect. fgsea is MIT
licensed, so its algorithms are ported with their order of operations: given the seeds
fgsea draws, gsea returns fgsea's table, and with its own seeds it agrees within
Monte Carlo error.

Changes

• cytepool.gsea.gsea(statistics, gene_sets, *, group, contrast, id_column, statistic_column, min_size, max_size, exponent, score_type, method, n_permutations, sample_size, eps, rng, n_threads):
  one row per set tested, in the order of gene_sets, with gene_set, n_genes (the
  set's genes among the ranked genes), enrichment_score,
  normalized_enrichment_score, p_value, adjusted_p_value (Benjamini-Hochberg),
  log2_error (fgsea's log2err) and leading_edge (genes separated by ", ", from
  the most extreme statistic inwards). Statistics: a DifferentialExpressionResult with
  a group and contrast, a Table or pandas DataFrame (id_column, statistic_column,
  default "t"), a pandas Series or a mapping. Sets: a GeneSets or a mapping. fgsea's
  options under their own names: exponent (gseaParam), score_type ("std",
  "pos", "neg"), min_size/max_size, method ("multilevel", fgseaMultilevel;
  "permutation", fgseaSimple), n_permutations (nPermSimple or nperm),
  sample_size, eps.
• cytepool.plotting.plot_enrichment(statistics, gene_set, *, group, contrast, id_column, statistic_column, exponent, ax):
  fgsea's plotEnrichment, the running sum with a tick at each gene of the set and
  dashed lines at its extremes.
• cytepool.plotting.plot_gsea_table(statistics, gene_sets, table, *, group, contrast, id_column, statistic_column, top, exponent, ax):
  fgsea's plotGseaTable for the top sets up and down (as fgsea's guide chooses them),
  each with its genes as bars at their ranks, its NES, p-value and adjusted p-value.
• cytepool._ranked_genes (private): reading statistics into a ranking, matching sets
  to it, and the running sum; it also now holds the group-and-contrast and collection
  checks that over_representation had, with the same messages.
• cytepool._options.random_generator (private): the check of rng, moved from the
  rotation tests, which use it unchanged.
• Documentation: API page, capability matrix row, a GSEA section in "Gene-set tests",
  rows in "Choosing a method", the overview, citing and glossary, fgsea moved among the
  tools CytePool is validated against, four new references (Korotkevich 2021, Mootha
  2003, Matsumoto 1998, Durstenfeld 1964; DOIs checked against Crossref, the preprint
  version against bioRxiv) on the Bibliography page, the Validation page (section,
  summary row, totals, the R environment) and the README status row.
• THIRD_PARTY_NOTICES.md: what is adapted from fgsea, file by file, with its licence.

Design notes

• A port, not a reimplementation. The enrichment score and leading edge follow
  calcGseaStat; random sets follow calcGseaStatCumulative (the scores of every
  prefix of one random sample of the largest set size, from upper convex hulls in
  square-root blocks, O(K^1.5) per permutation, including fgsea's floor of 1e-5/1024
  under the smallest nonzero statistic); the multilevel method follows EsRuler
  (duplication of the upper half, the Markov chain with chunked walks and a remembered
  gene, getPvalue with digamma-based level probabilities and the sign correction); the
  R-level rules follow fgseaSimpleImpl and fgseaMultilevel (NES, the (1 + b)/(1 + n)
  p-values, NaN for fewer than 10 random scores of a set's sign, the choice of the
  multilevel method by comparing a beta-quantile interval of the permutation estimate
  with multilevelError, eps). The order of floating-point operations is kept, so
  that decisions in the Markov chain fall the same way: the p-values agree to 1e-13.
• Random numbers. rng (a numpy Generator, a seed, or None for a fresh one) draws
  distinct seeds from 1 to 10^9: two for the multilevel method, as fgsea's two
  sample.int(10^9, 1), or one per chunk of max(1000, ceiling(n / 128)) permutations,
  as fgseaSimple. Each starts MT19937 seeded as C++'s std::mt19937 (checked against
  the C++ standard's 10,000th output and against numpy's MT19937), with fgsea's
  rejection sampler for integers. Permutation chunks and multilevel set sizes run in
  parallel with seeds of their own, so results do not depend on the number of threads.
• The running sum of the observed sets rounds its partial sums once, as R's
  extended-precision cumsum and sum do (found by probing), so that ties between steps
  of the walk, which decide the leading edge, fall as in fgsea.
• Ties. Genes with equal statistics keep the order given, as R's
  sort(decreasing = TRUE) keeps them (found by probing); a warning gives the share of
  tied nonzero statistics, as fgsea's does. fgsea's other warnings are kept: all
  statistics positive with "std", sets with too few random scores of their sign,
  p-values likely overestimated, p-values below eps.
• Inputs. Rows without a gene id or without a statistic (NaN) are left out (fgsea
  refuses them, so users drop them first); infinite statistics, ids on two rows, fewer
  than two ranked genes and all-zero statistics are refused. A set of every ranked gene
  is never tested, as in fgsea. No set in the size range gives an empty table, as in
  fgsea; no gene of any set among the ranked genes is an error naming examples of both.
• Deliberate differences. (1) Adjusted p-values are BH over the sets with a p-value;
  when no set goes to the multilevel method, fgsea adjusts over the sets before it
  removes those with 1 to 9 random scores of their sign, so its values come out
  smaller or larger. (2) fgsea tests a set alone with R's sample.int; CytePool uses
  fgsea's sampler for one set as for many (agreement within Monte Carlo error).
  (3) Random samples of at least half the genes are shuffled by Fisher-Yates with
  fgsea's integers, where fgsea calls std::shuffle, whose algorithm is left to each
  C++ library. The table keeps the order of gene_sets; fgsea sorts by name.
• Performance. On fgsea's example (12,000 genes, 586 sets of 15 to 500 genes) with
  the defaults: fgsea 2.9 s, gsea 4.1 s on one thread after compilation; 10,000
  permutations: fgsea 1.9 s, gsea 1.5 s on 2 threads.
• Not included: collapsePathways, fora, label permutation (fgseaLabel) and
  GESECA.

Validation

validation/gsea.py against fgsea 1.32.2 (R 4.4.3), on 4,000 simulated genes (6% up,
5% down, 60 tied at 1.25, 30 at zero) and 90 sets of 1 to 1,000 genes (a fifth outside
the ranking; enriched up, down, both or none; every gene, one gene, the tied block, a
set listing genes twice, two identical sets). The R script records the seeds fgsea
draws where it hands them to its compiled code (trace of fgsea's MIT-licensed
fgseaSimpleImpl and multilevelImpl), and gsea draws the same through rng.

Reference: R 4.4.3, fgsea 1.32.2
| check | reference | metric | value | tolerance | result |
| --- | --- | --- | --- | --- | --- |
| defaults: 89 sets: sets tested, sizes, leading edges and enrichment scores | fgsea | max relative difference | 0.0e+00 | 1e-10 | pass |
| defaults: 89 sets: NES, p-values, log2 errors and adjusted p-values, given fgsea's seeds | fgsea | max relative difference | 7.1e-15 | 1e-10 | pass |
| gseaParam 0: 89 sets: sets tested, sizes, leading edges and enrichment scores | fgsea | max relative difference | 0.0e+00 | 1e-10 | pass |
| gseaParam 0: 89 sets: NES, p-values, log2 errors and adjusted p-values, given fgsea's seeds | fgsea | max relative difference | 7.1e-15 | 1e-10 | pass |
| gseaParam 0.5: 89 sets: sets tested, sizes, leading edges and enrichment scores | fgsea | max relative difference | 0.0e+00 | 1e-10 | pass |
| gseaParam 0.5: 89 sets: NES, p-values, log2 errors and adjusted p-values, given fgsea's seeds | fgsea | max relative difference | 7.2e-15 | 1e-10 | pass |
| scoreType neg: 89 sets: sets tested, sizes, leading edges and enrichment scores | fgsea | max relative difference | 0.0e+00 | 1e-10 | pass |
| scoreType neg: 89 sets: NES, p-values, log2 errors and adjusted p-values, given fgsea's seeds | fgsea | max relative difference | 2.2e-13 | 1e-10 | pass |
| scoreType pos, all positive: 89 sets: sets tested, sizes, leading edges and enrichment scores | fgsea | max relative difference | 0.0e+00 | 1e-10 | pass |
| scoreType pos, all positive: 89 sets: NES, p-values, log2 errors and adjusted p-values, given fgsea's seeds | fgsea | max relative difference | 1.6e-13 | 1e-10 | pass |
| sizes 15 to 500, eps 0, sampleSize 51, nPermSimple 2000: 53 sets: sets tested, sizes, leading edges and enrichment scores | fgsea | max relative difference | 0.0e+00 | 1e-10 | pass |
| sizes 15 to 500, eps 0, sampleSize 51, nPermSimple 2000: 53 sets: NES, p-values, log2 errors and adjusted p-values, given fgsea's seeds | fgsea | max relative difference | 7.1e-14 | 1e-10 | pass |
| permutations, nperm 2500: 89 sets: sets tested, sizes, leading edges and enrichment scores | fgsea | max relative difference | 0.0e+00 | 1e-10 | pass |
| permutations, nperm 2500: 89 sets: NES, p-values and adjusted p-values, given fgsea's seeds | fgsea | max relative difference | 0.0e+00 | 1e-10 | pass |
| mostly positive statistics: 10 sets: sets tested, sizes, leading edges and enrichment scores | fgsea | max relative difference | 0.0e+00 | 1e-10 | pass |
| mostly positive statistics: 10 sets: NES, p-values, log2 errors, given fgsea's seeds | fgsea | max relative difference | 4.2e-16 | 1e-10 | pass |
| mostly positive statistics: 10 sets: adjusted p-values, over the sets with a p-value | fgsea | max relative difference | 1.6e-01 | n/a | deliberate difference |
| one set: 1 set: sets tested, sizes, leading edges and enrichment scores | fgsea | max relative difference | 0.0e+00 | 1e-10 | pass |
| 36 sets: mean log10 p-value of 20 runs, seeds of their own | fgsea | largest |difference| / SE | 2.3e+00 | 4 | pass |
| 36 sets: mean NES of 20 runs, seeds of their own | fgsea | largest |difference| / SE | 2.4e+00 | 4 | pass |
| one set alone: mean log10 p-value of 20 runs, seeds of their own | fgsea | largest |difference| / SE | 1.3e+00 | 4 | pass |
| one set alone: mean NES of 20 runs, seeds of their own | fgsea | largest |difference| / SE | 1.8e+00 | 4 | pass |
| plot_enrichment: the walk, ticks and extremes of 6 plots | plotEnrichmentData | max relative difference | 1.7e-14 | 1e-10 | pass |
| plot_gsea_table: the bars and numbers of 10 sets | plotEnrichmentData, plotGseaTable | max relative difference | 1.1e-13 | 1e-10 | pass |
24 checks, 0 failed

Tiers: exact (1e-10) for scores, leading edges, figures and, given fgsea's seeds,
everything; calibrated (4 standard errors of the difference of the means of 20 runs)
with seeds of their own. 7 min 3 s, 323 MiB peak, with 2 threads.

RSCRIPT=$REF/r-gsea/bin/Rscript PYTHONPATH=. python validation/gsea.py
RSCRIPT=$REF/r-gsea/bin/Rscript PYTHONPATH=. python validation/gsea.py --write-fixtures

Tests

tests/test_gsea.py, 54 tests, against fixtures that the validation writes from
fgsea (300 genes, 14 sets): every exact case replayed with fgsea's seeds (all columns
to 1e-10), fgsea's choice of the multilevel method for every count (1,000 and 2,000
permutations), the single set within Monte Carlo error, and plot_enrichment against
plotEnrichmentData. Independent of the code: enrichment scores and leading edges
derived by hand (up, down, equal deviations with each score type, the tie rule); MT19937
against the C++ standard and numpy; the hull scores of every prefix against a walk of
each prefix; invariants (negated statistics negate the scores; tied genes outside a set
and the order of genes or of sets leave results unchanged; seeds reproduce results and
threads do not change them); sizes after matching, duplicates and genes outside the
ranking; a result, its table, a mapping and a pandas Series rank alike; adjusted p-values
are BH of the reported p-values; warnings and messages; plot_gsea_table's selection,
bars and number format. The masked-input inventory gains gsea, plot_enrichment and
plot_gsea_table (each with a variant for a differential-expression result). The
module runs in about 2 to 5 s warm and 73 s cold (compiling the routines and limma-voom
for one test). Full suite, on top of the branches before it: 4,244 passed, 18 skipped;
also with bounds checks, and with NUMBA_CPU_NAME=generic,
OPENBLAS_CORETYPE=Sandybridge and 4 threads.