From d812f2ea4b9908501a4938fc61bc47fa89589a72 Mon Sep 17 00:00:00 2001 From: parashardhapola Date: Mon, 26 Oct 2026 16:47:51 +0100 Subject: [PATCH] Add gene set enrichment analysis of ranked genes, as fgsea gsea ports fgsea 1.32.2 (MIT licence; notice in THIRD_PARTY_NOTICES.md): enrichment scores and leading edges, normalized scores, and p-values by fgsea's multilevel method or by permutations of the genes, from seeds that a numpy generator draws. Given the seeds fgsea draws, it returns fgsea's table. plot_enrichment and plot_gsea_table draw fgsea's enrichment plot and table. Checked against fgsea 1.32.2 (validation/gsea.py: 24 checks, none failed). --- README.md | 2 +- THIRD_PARTY_NOTICES.md | 53 + cytepool/_options.py | 39 +- cytepool/_ranked_genes.py | 329 ++++ cytepool/gsea.py | 1685 +++++++++++++++++++ cytepool/over_representation.py | 42 +- cytepool/plotting.py | 301 +++- cytepool/stats/rotation_gene_set_tests.py | 24 +- docs/source/index.md | 2 +- docs/source/reference/api/gsea.md | 21 + docs/source/reference/api/index.md | 1 + docs/source/reference/bibliography.md | 3 + docs/source/reference/capabilities.md | 1 + docs/source/reference/citing.md | 1 + docs/source/reference/glossary.md | 5 + docs/source/reference/related_software.md | 2 +- docs/source/references.bib | 41 + docs/source/statistics/choosing_a_method.md | 1 + docs/source/statistics/gene_set_tests.md | 58 +- docs/source/validation.md | 85 +- tests/data/gsea/fgsea.json | 1 + tests/test_gsea.py | 629 +++++++ tests/test_masked_inputs.py | 119 ++ validation/gsea.py | 476 ++++++ 24 files changed, 3841 insertions(+), 80 deletions(-) create mode 100644 cytepool/_ranked_genes.py create mode 100644 cytepool/gsea.py create mode 100644 docs/source/reference/api/gsea.md create mode 100644 tests/data/gsea/fgsea.json create mode 100644 tests/test_gsea.py create mode 100644 validation/gsea.py diff --git a/README.md b/README.md index 1b28eed..2e493dd 100644 --- a/README.md +++ b/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 | diff --git a/THIRD_PARTY_NOTICES.md b/THIRD_PARTY_NOTICES.md index 225073e..1a010bc 100644 --- a/THIRD_PARTY_NOTICES.md +++ b/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 diff --git a/cytepool/_options.py b/cytepool/_options.py index 4a16064..0ce2bf8 100644 --- a/cytepool/_options.py +++ b/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) diff --git a/cytepool/_ranked_genes.py b/cytepool/_ranked_genes.py new file mode 100644 index 0000000..4d6557d --- /dev/null +++ b/cytepool/_ranked_genes.py @@ -0,0 +1,329 @@ +"""Genes ranked by a statistic, gene sets matched to them, and the running sum of gene +set enrichment analysis along them. + +Gene set enrichment analysis, its figures and over-representation analysis read their +genes the same way: from a group and contrast of a +:class:`~cytepool.differential_expression.DifferentialExpressionResult`, from a table +of results, or, for a ranking, from a mapping of gene id to statistic, and their gene +sets from a :class:`~cytepool.gene_sets.GeneSets` or a mapping of set name to genes. + +The ranking: + Genes are ranked by decreasing statistic. Genes with equal statistics keep the + order in which they are given, as R's ``sort(decreasing = TRUE)``, with which + fgsea ranks them, keeps them (found by probing R 4.4.3). Rows without a gene id + (None, NaN or an empty string) and rows without a statistic (NaN) are left out; + an infinite statistic, or a gene id on two rows, is refused. + +The running sum: + For a set of m of the N ranked genes, at ranks p_1 < ... < p_m, with weights + w_i = |r_i|^p of their statistics r_i (Subramanian et al. 2005), the running sum + after the set's i-th gene is the share of the set's weight reached, minus the + share of the N - m other genes passed: + + top_i = (w_1 + ... + w_i) / W - (p_i - i) / (N - m), W = w_1 + ... + w_m, + + and just before it, bottom_i = top_i - w_i / W; when every weight is zero, each + gene weighs 1 / m. The enrichment score is the extreme of the walk, the largest + top or the smallest bottom (fgsea's ``calcGseaStat``; Korotkevich et al. 2021). + The partial sums and W are rounded once from their exact values, as R's + ``cumsum`` and ``sum``, which accumulate in extended precision, round them + (found by probing R 4.4.3), so that ties between steps of the walk, which decide + the leading edge, fall as in fgsea. + +References: + Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich + A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP (2005). Gene set enrichment + analysis: a knowledge-based approach for interpreting genome-wide expression + profiles. Proceedings of the National Academy of Sciences 102, 15545-15550. + doi:10.1073/pnas.0506580102 + + 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 +""" + +from collections.abc import Mapping +from typing import Any, NamedTuple + +import numpy as np +import numpy.typing as npt + +from ._columns import is_data_frame +from ._gene_ids import is_missing +from ._jit import jit +from ._masked_arrays import real_array, reject_masked +from ._messages import quoted +from ._summation import two_sum +from .differential_expression import DifferentialExpressionResult +from .gene_sets import GeneSets +from .table import Table + +FloatArray = npt.NDArray[np.float64] +IntArray = npt.NDArray[np.int64] + +# The names of statistics, for the axes of figures. +_STATISTIC_LABELS = {"t": "Moderated t-statistic", "statistic": "Statistic"} + + +class RankedGenes(NamedTuple): + """Genes in decreasing order of their statistic.""" + + ids: list[str] + values: FloatArray + label: str # what the statistic is, for the axes of figures + + +def gene_set_collection(gene_sets: Any) -> GeneSets: + """``gene_sets`` as a GeneSets, refusing an empty collection.""" + if isinstance(gene_sets, GeneSets): + collection = gene_sets + elif isinstance(gene_sets, Mapping): + collection = GeneSets(gene_sets) + else: + raise TypeError( + "gene_sets must be a GeneSets or a mapping of set name to gene ids, such " + f"as {{'interferon': ['ISG15', 'IFIT1']}}; got {type(gene_sets).__name__}" + ) + if not len(collection): + raise ValueError("gene_sets is empty; give at least one gene set") + return collection + + +def group_and_contrast( + result: DifferentialExpressionResult, group: Any, contrast: Any +) -> tuple[str, str]: + """The group and contrast of a result that are meant: those given, or its only + ones.""" + chosen = [] + for value, name, options in ( + (group, "group", list(result.groups)), + (contrast, "contrast", list(result.contrasts.names)), + ): + if value is None: + if len(options) != 1: + raise ValueError( + f"the result has {len(options)} {name}s; name one with {name}, " + f"one of {quoted(options)}" + ) + value = options[0] + elif not isinstance(value, str): + reject_masked(value, name) + raise TypeError(f"{name} must be a string, got {type(value).__name__}") + chosen.append(value) + return chosen[0], chosen[1] + + +def ranked_genes( + statistics: Any, + *, + group: str | None, + contrast: str | None, + id_column: str, + statistic_column: str, +) -> RankedGenes: + """The genes of ``statistics`` in decreasing order of their statistic. + + Args: + statistics: A DifferentialExpressionResult, with ``group`` and ``contrast``; + a Table or pandas DataFrame with the columns ``id_column`` and + ``statistic_column``; a pandas Series of statistics indexed by gene id; + or a mapping of gene id to statistic. + group: The group of a DifferentialExpressionResult; None for its only one. + contrast: Its contrast; None for its only one. + id_column: The column of gene ids of a table. + statistic_column: The column of statistics of a table. + + Raises: + TypeError: If ``statistics`` is none of these, a gene id is not a string, a + statistic is not a real number, or an argument is or holds a numpy masked + array. + ValueError: If ``group`` or ``contrast`` is given without a + DifferentialExpressionResult or left out of one that has several, a + table lacks a column, a statistic is infinite, a gene id is on two rows, + or fewer than two genes have an id and a statistic. + KeyError: If the result has no table of the group and contrast. + """ + for value, name in ( + (id_column, "id_column"), + (statistic_column, "statistic_column"), + ): + reject_masked(value, name) + if not isinstance(value, str): + raise TypeError(f"{name} must be a column name (a string), got {value!r}") + label = _STATISTIC_LABELS.get(statistic_column, statistic_column) + if isinstance(statistics, DifferentialExpressionResult): + table = statistics.table(*group_and_contrast(statistics, group, contrast)) + ids, values = _table_columns(table, id_column, statistic_column) + else: + if group is not None or contrast is not None: + reject_masked(group, "group") + reject_masked(contrast, "contrast") + raise ValueError( + "group and contrast choose a table of a DifferentialExpressionResult; " + f"statistics is a {type(statistics).__name__}" + ) + if isinstance(statistics, Table) or is_data_frame(statistics): + table = statistics if isinstance(statistics, Table) else Table(statistics) + ids, values = _table_columns(table, id_column, statistic_column) + elif _is_series(statistics): + ids = statistics.index.tolist() + values = _real_values(statistics.to_numpy(), "statistics") + label = str(statistics.name) if statistics.name is not None else "Statistic" + elif isinstance(statistics, Mapping): + reject_masked(statistics, "statistics", missing="NaN") + ids = list(statistics.keys()) + values = _real_values(list(statistics.values()), "statistics") + label = "Statistic" + else: + reject_masked(statistics, "statistics") + raise TypeError( + "statistics must be a DifferentialExpressionResult, a table of results " + "(a Table or a pandas DataFrame), a pandas Series indexed by gene id or " + "a mapping of gene id to statistic, got " + f"{type(statistics).__name__}" + ) + return _ranking(ids, values, label) + + +def matched_positions(genes: tuple[str, ...], index: Mapping[str, int]) -> IntArray: + """The positions of a set's genes in the ranking, in the set's order; genes that + are not ranked are left out, as fgsea leaves them out.""" + return np.array([index[gene] for gene in genes if gene in index], dtype=np.int64) + + +@jit +def running_sum( + weights: FloatArray, positions: IntArray +) -> tuple[FloatArray, FloatArray]: + """The running sum at and just before each gene of a set, fgsea's tops and + bottoms. + + Args: + weights: The weights |r|^p of all genes, in ranked order. + positions: The 0-based positions of the set's genes, increasing; fewer than + all the genes. + + Returns: + ``(tops, bottoms)``, one value per gene of the set. + """ + n = weights.size + m = positions.size + # The partial sums of the weights, each rounded once from the exact sum. + partial = np.empty(m) + high = 0.0 + low = 0.0 + for i in range(m): + high, error = two_sum(high, weights[positions[i]]) + low += error + partial[i] = high + low + total = partial[m - 1] if m else 0.0 + tops = np.empty(m) + bottoms = np.empty(m) + for i in range(m): + reached = (i + 1) / m if total == 0.0 else partial[i] / total + tops[i] = reached - (positions[i] - i) / (n - m) + step = 1.0 / m if total == 0.0 else weights[positions[i]] / total + bottoms[i] = tops[i] - step + return tops, bottoms + + +# Reading the statistics. + + +def _table_columns( + table: Table, id_column: str, statistic_column: str +) -> tuple[list[Any], FloatArray]: + missing = [name for name in (id_column, statistic_column) if name not in table] + if missing: + raise ValueError( + f"the table has no column {quoted(missing)}; its columns are " + f"{quoted(list(table.columns))}. Name the columns of the gene ids and of " + "the statistics with id_column and statistic_column; a pandas DataFrame " + "keeps its row index, such as gene ids, as a column only when the index " + "has a name" + ) + ids = table[id_column] + if ids.dtype.kind != "O": + raise TypeError( + f"the column {id_column!r} must hold gene ids as text, got dtype " + f"{ids.dtype}; convert numeric ids, such as Entrez ids, to text first" + ) + try: + values = real_array( + table[statistic_column], + f"the column {statistic_column!r}", + boolean=False, + objects=False, + ) + except TypeError: + raise TypeError( + f"the column {statistic_column!r} must hold numbers, got dtype " + f"{table[statistic_column].dtype}" + ) from None + return ids.tolist(), values + + +def _real_values(values: Any, name: str) -> FloatArray: + try: + array = real_array(values, name, boolean=False, objects=False) + except TypeError: + raise TypeError(f"{name} must map gene ids to numbers") from None + if array.ndim != 1: + raise ValueError(f"{name} must hold one number per gene") + return array + + +def _is_series(value: object) -> bool: + """Whether ``value`` is a pandas Series, recognised without importing pandas.""" + return all( + hasattr(value, name) for name in ("index", "to_numpy", "name", "dtype") + ) and not hasattr(value, "columns") + + +def _ranking(ids: list[Any], values: FloatArray, label: str) -> RankedGenes: + """The genes with an id and a statistic, in decreasing order of the statistic.""" + if len(ids) != values.size: + raise ValueError("statistics must hold one statistic per gene id") + named = np.empty(len(ids), dtype=bool) + for i, gene in enumerate(ids): + if isinstance(gene, str): + named[i] = bool(gene) + elif is_missing(gene): + named[i] = False + else: + raise TypeError( + "gene ids must be strings, got " + f"{type(gene).__name__} {gene!r}; convert numeric ids, such as Entrez " + "ids, to text first" + ) + infinite = np.isinf(values) + if infinite.any(): + where = int(np.flatnonzero(infinite)[0]) + raise ValueError( + f"the statistic of {ids[where]!r} is {values[where]}; statistics must be " + "finite, or NaN for a gene that was not tested" + ) + kept = np.flatnonzero(named & ~np.isnan(values)) + genes = [str(ids[i]) for i in kept.tolist()] + seen: dict[str, int] = {} + for gene in genes: + seen[gene] = seen.get(gene, 0) + 1 + if seen[gene] == 2: + raise ValueError( + f"the gene id {gene!r} is on more than one row; each gene needs one " + "statistic: keep one row per gene id (for example the one with the " + "largest absolute statistic), or rank the genes by another id" + ) + if len(genes) < 2: + raise ValueError( + f"{len(genes)} gene(s) have an id and a statistic; a ranking needs at " + "least two" + ) + kept_values = values[kept] + # Decreasing, with ties in the order given, as R's sort(decreasing = TRUE). + order = np.argsort(-kept_values, kind="stable") + return RankedGenes( + [genes[i] for i in order.tolist()], + np.ascontiguousarray(kept_values[order]), + label, + ) diff --git a/cytepool/gsea.py b/cytepool/gsea.py new file mode 100644 index 0000000..55c118d --- /dev/null +++ b/cytepool/gsea.py @@ -0,0 +1,1685 @@ +"""Gene set enrichment analysis (GSEA) of genes ranked by a statistic, as fgsea computes +it. + +:func:`gsea` asks of each gene set whether its genes gather at the top or at the +bottom of a ranking of all genes by a statistic, such as the moderated t-statistic of +a contrast: preranked GSEA (Subramanian et al. 2005) as fgsea 1.32.2 computes it +(Korotkevich et al. 2021), with p-values by fgsea's default multilevel method +(``fgseaMultilevel``) or by permutations of the genes (``fgseaSimple``). Results match +fgsea's exactly when CytePool is given the random seeds that fgsea draws, and within +Monte Carlo error otherwise, except where noted under "Deliberate differences". The +routines port fgsea's (MIT licence), as ``THIRD_PARTY_NOTICES.md`` records. +:func:`~cytepool.plotting.plot_enrichment` and :func:`~cytepool.plotting.plot_gsea_table` +draw its figures. + +The genes and the sets: + The statistics come from a group and contrast of a + :class:`~cytepool.differential_expression.DifferentialExpressionResult`, from a + table of results or from a mapping of gene id to statistic. Genes are ranked by + decreasing statistic, and genes with equal statistics keep the order in which they + are given, as in fgsea: the score of a set whose genes tie with other genes + depends on that order, and a warning says how many statistics tie. Rows without a + gene id or without a statistic (NaN) are left out, where fgsea refuses them; + infinite statistics and ids on two rows are refused. Each set's genes are matched + to the ranked genes, and a set is tested when ``min_size`` to ``max_size`` of its + genes are ranked (fgsea's ``minSize`` and ``maxSize``); a set of every ranked gene + is never tested, as its score is not defined. The table keeps the order of + ``gene_sets``, where fgsea orders the sets by name. + +The enrichment score: + With weights |r|^p of the statistics r, p the ``exponent`` (fgsea's + ``gseaParam``), the walk down the ranking rises at each gene of a set by its share + of the set's weight and falls at each other gene by 1 / (N - m), for m of the N + genes in the set. The enrichment score (ES) is the walk's largest deviation from + zero: with ``score_type="std"``, its maximum where that exceeds the absolute + minimum, its minimum where the absolute minimum is larger, and 0 where they are + equal; ``"pos"`` and ``"neg"`` take the maximum and the minimum, for one-sided + tests. An exponent of 0 gives the unweighted statistic of Mootha et al. (2003), 1 + the weighted one of Subramanian et al. (2005). The leading edge is the set's genes + from the top of the ranking down to the maximum, for a positive score, or from + the minimum to the bottom, for a negative one, listed from the most extreme + statistic inwards: the genes that drive the score. + +Random gene sets: + The null distribution of a set's score is that of random sets of the same size + (gene permutation). As fgsea does, each of the ``n_permutations`` permutations + draws one random sample of as many genes as the largest set, and the scores of all + its first genes, one score for each size, are computed together in + O(K^(3/2)) time for K genes, from the upper convex hulls of the walk in blocks + of about sqrt(K) genes (fgsea's ``calcGseaStatCumulative``). + + The normalized enrichment score (NES) divides the ES by the mean of the random + scores of its sign: of those at least 0 for a positive ES, of the absolute values + of those at most 0 otherwise; it is NaN where that mean is 0 or there is no such + score. The permutation p-value is (1 + b) / (1 + n), with b the random scores at + least as extreme as the ES and n those of its sign, so that it is conditional on + the sign of the score, as GSEA's is. + +P-values by the multilevel method: + ``method="multilevel"`` (fgsea's default) starts from the permutation p-values of + ``n_permutations`` permutations. A set with fewer than 10 random scores of the sign + of its ES gets NaN for its p-value and NES. A set goes to adaptive multilevel + splitting (Korotkevich et al. 2021) when its permutation p-value would be less + accurate than the multilevel estimate: with b' of the n random scores beyond its ES + on its side, when half the larger distance from log2 of (b' + 1) / (n + 1) to the + ends of a 95% interval of the proportion (from the 2.5% and 97.5% quantiles of beta + distributions) exceeds the expected error of the multilevel estimate, as fgsea + decides; with fgsea's defaults, when b' is at most 23. For each set size, + ``sample_size`` random sets + are drawn; at each level the half with the lower positive scores (the largest + deviation of the walk above zero) is replaced by copies of the other half, each + of which a Markov chain then moves, swapping one gene at a time and keeping its + score above the level's median, about ``sample_size`` times the size moves in all. + The levels continue until the median reaches the ES, the median stops rising, or + the p-value would fall below ``eps``. The probability of each level, about 1/2, is + estimated on the log scale by the expected logarithm of an order statistic of + uniform variables, psi(k) - psi(n + 1) for the k-th of n with psi the digamma + function, and the p-value is the product over the levels, times the share of the + sampled sets beyond the ES whose signed score is positive, divided by the share of + random scores of the ES's sign among the permutations. ``log2_error`` gives the + expected standard deviation of log2 of the p-value: from the number of levels for + the multilevel method, and from the trigamma function of the counts for the + permutation p-values (fgsea's ``log2err``). It is NaN where fewer than half of the + last level's sets have a positive signed score, which makes the p-value likely + too large, and where the p-value is below ``eps``, which is then its value: when + the multilevel method has run, p-values below ``eps`` are reported as ``eps``. + +Random numbers: + ``rng`` draws the seeds of the random sets: two distinct integers from 1 to 10^9 + for ``method="multilevel"``, the first for the permutations and the second for the + multilevel method, and one for each chunk of the permutations for + ``method="permutation"``, which, as ``fgseaSimple``, splits the permutations into + chunks of max(1000, ceiling(n / 128)), run in parallel. Each seed starts the + Mersenne Twister MT19937 (Matsumoto and Nishimura 1998), initialised as C++'s + ``std::mt19937`` that fgsea uses, from which random integers are drawn by + rejection as fgsea draws them. Given fgsea's seeds, which R draws with + ``sample.int(10^9, 1)``, the results are fgsea's: ``validation/gsea.py`` replays + them. A seeded ``rng`` gives the same results every time, with any number of + threads. + +Deliberate differences from fgsea: + - Adjusted p-values are Benjamini and Hochberg's over the sets with a p-value. + When no set goes to the multilevel method, fgsea adjusts them over the sets + before it removes those with fewer than 10 random scores of their sign, as + though those sets had been tested; its adjusted p-values are then smaller or + larger than those of the sets it reports. + - fgsea tests a single set with random sets from R's ``sample.int``, and several + sets with its own sampler; CytePool uses fgsea's sampler for one set too, so + that the two agree within Monte Carlo error. + - Random samples of at least half of the genes are put in random order by the + Fisher-Yates shuffle (Durstenfeld 1964) with fgsea's random integers, where + fgsea calls C++'s ``std::shuffle``, whose algorithm the standard leaves to each + library; the results agree within Monte Carlo error. + +References: + Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich + A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP (2005). Gene set enrichment + analysis: a knowledge-based approach for interpreting genome-wide expression + profiles. Proceedings of the National Academy of Sciences 102, 15545-15550. + doi:10.1073/pnas.0506580102 + + 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 + + Mootha VK, Lindgren CM, Eriksson KF, Subramanian A, Sihag S, Lehar J, Puigserver P, + Carlsson E, Ridderstrale M, Laurila E, Houstis N, Daly MJ, Patterson N, Mesirov JP, + Golub TR, Tamayo P, Spiegelman B, Lander ES, Hirschhorn JN, Altshuler D, Groop LC + (2003). PGC-1alpha-responsive genes involved in oxidative phosphorylation are + coordinately downregulated in human diabetes. Nature Genetics 34, 267-273. + doi:10.1038/ng1180 + + Matsumoto M, Nishimura T (1998). Mersenne twister: a 623-dimensionally + equidistributed uniform pseudo-random number generator. ACM Transactions on + Modeling and Computer Simulation 8, 3-30. doi:10.1145/272991.272995 + + Durstenfeld R (1964). Algorithm 235: random permutation. Communications of the + ACM 7, 420. doi:10.1145/364520.364540 + + 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 +""" + +import math +import warnings +from collections.abc import Mapping +from typing import Any, NamedTuple + +import numpy as np +import numpy.typing as npt +from numba import prange + +from ._columns import text_column +from ._jit import by_chunks, jit, thread_limit +from ._masked_arrays import reject_masked +from ._messages import PACKAGE_DIRECTORY, quoted +from ._options import ( + checked_choice, + number_at_least, + positive_integer, + random_generator, + real_number, +) +from ._ranked_genes import ( + RankedGenes, + gene_set_collection, + matched_positions, + ranked_genes, + running_sum, +) +from ._summation import accurate_sum +from .gene_sets import GeneSets +from .stats.multiple_testing import adjust_p_values +from .stats.special import betainc, digamma_scalar, trigamma +from .table import Table + +__all__ = ["gsea"] + +FloatArray = npt.NDArray[np.float64] +IntArray = npt.NDArray[np.int64] +BoolArray = npt.NDArray[np.bool_] + +_SCORE_TYPES = ("std", "pos", "neg") +_METHODS = ("multilevel", "permutation") +_SEPARATOR = ", " # between the genes of a leading edge +_SEED_RANGE = 10**9 # fgsea's seeds: sample.int(10^9, 1), from 1 to 10^9 +# fgseaMultilevel: a set with fewer random scores of its sign gets no p-value. +_FEWEST_OF_SIGN = 10 +# fgseaSimple's chunks of permutations: at least 1000, and at most 128 chunks. +_CHUNK, _MOST_CHUNKS = 1000, 128 +# The Mersenne Twister MT19937, as C++'s std::mt19937 (Matsumoto and Nishimura 1998): +# its state of 624 words, the index of the next word kept after them. +_MT_WORDS = 624 +_LARGEST_DRAW = 4294967295 # 2^32 - 1, the largest output of MT19937 +# fgsea's constants: the tolerance of its convex hulls, of its multilevel levels, +# the share of a set's genes swapped by one round of its Markov chain, the most +# draws for one gene of a random sample, and the floor of the statistics of a +# random sample (1e-5, divided by 1024 under the smallest nonzero one). +_HULL_TOLERANCE = 1e-13 +_LEVEL_TOLERANCE = 1e-10 +_SWAP_SHARE = 0.1 +_DRAWS_PER_GENE = 100 +_STATISTIC_FLOOR, _FLOOR_DIVISOR = 1e-5, 1024.0 +_COLUMNS = ( + "gene_set", + "n_genes", + "enrichment_score", + "normalized_enrichment_score", + "p_value", + "adjusted_p_value", + "log2_error", + "leading_edge", +) + + +class _Settings(NamedTuple): + min_size: int + max_size: int | None + exponent: float + score_type: int # 0 std, 1 pos, 2 neg + multilevel: bool + n_permutations: int + sample_size: int + eps: float + + +def gsea( + statistics: Any, + gene_sets: Mapping[str, Any], + *, + group: str | None = None, + contrast: str | None = None, + id_column: str = "feature_id", + statistic_column: str = "t", + min_size: int = 1, + max_size: int | None = None, + exponent: float = 1.0, + score_type: str = "std", + method: str = "multilevel", + n_permutations: int = 1000, + sample_size: int = 101, + eps: float = 1e-50, + rng: np.random.Generator | int | None = None, + n_threads: int | None = None, +) -> Table: + """Gene set enrichment analysis of genes ranked by a statistic, as fgsea. + + Preranked GSEA (Subramanian et al. 2005) as fgsea's ``fgsea`` computes it + (Korotkevich et al. 2021): whether the genes of each set gather at the top or the + bottom of the ranking, with the enrichment score, its leading edge, the score + normalized by those of random sets of the same size, and p-values by fgsea's + multilevel method or by permutations of the genes. The module docstring gives the + method, the random numbers and the deliberate differences from fgsea. + + Args: + statistics: The statistics that rank the genes: a + :class:`~cytepool.differential_expression.DifferentialExpressionResult`, + with ``group`` and ``contrast``; a table of results, a ``Table`` or a + pandas DataFrame, with a column of gene ids and a column of statistics; + or a pandas Series of statistics indexed by gene id, or a mapping of gene + id to statistic. NaN marks a gene without a statistic. + gene_sets: A :class:`~cytepool.gene_sets.GeneSets`, such as one read by + :func:`~cytepool.gene_sets.read_gmt`, or a mapping of set name to gene ids. + group: The group (cell type) of a ``DifferentialExpressionResult``; it may be + left out when the result has one. + contrast: The contrast of a ``DifferentialExpressionResult``; it may be left + out when the result tests one. + id_column: The column of gene ids of a table, such as ``"feature_id"`` or + ``"symbol"``. + statistic_column: The column of statistics of a table, such as ``"t"`` (the + moderated t-statistic) or ``"log2_fold_change"``. + min_size: The fewest of a set's genes among the ranked genes for the set to be + tested, at least 1 (fgsea's ``minSize``). + max_size: The most of them, at least ``min_size``; None, and any value of at + least the number of ranked genes, mean one fewer than the ranked genes + (fgsea's ``maxSize``). + exponent: The power p of the weights |r|^p of the statistics, at least 0 + (fgsea's ``gseaParam``): 1 for GSEA's weighted score, 0 for the unweighted + one. + score_type: ``"std"`` for the score of the larger deviation in either + direction, ``"pos"`` for enrichment at the top only and ``"neg"`` at the + bottom only (fgsea's ``scoreType``). + method: ``"multilevel"`` (fgsea's default) or ``"permutation"`` (fgsea's + ``fgseaSimple``). + n_permutations: The random samples of genes, at least 1: for the + normalized scores and the permutation p-values, which are at least + 1 / (n_permutations + 1); fgsea's ``nPermSimple`` for the multilevel + method and ``nperm`` for permutations. + sample_size: The random sets of each level of the multilevel method, at least + 3; an even number is raised by one, as fgsea raises it (fgsea's + ``sampleSize``). + eps: The smallest p-value the multilevel method estimates, in [0, 1); 0 for no + limit (fgsea's ``eps``). + rng: A ``numpy.random.Generator`` to draw the seeds of the random sets from, + or an integer seed for one, or None for a fresh, unseeded generator. + n_threads: The most threads to use; None uses numba's current setting. The + result does not depend on it. + + Returns: + A table with one row per set tested, in the order of ``gene_sets``: + ``gene_set``; ``n_genes``, the set's genes among the ranked genes; + ``enrichment_score``; ``normalized_enrichment_score``; ``p_value``; + ``adjusted_p_value``, Benjamini-Hochberg over the sets with a p-value; + ``log2_error``, the expected standard deviation of log2 of the p-value; and + ``leading_edge``, the genes of the leading edge, from the most extreme + statistic inwards, separated by ``", "``. A set with too few random scores of + its sign has NaN normalized score, p-values and error. The table has no row + when no set has between ``min_size`` and ``max_size`` genes ranked. + + Raises: + TypeError: If an argument has the wrong type, gene ids are not strings, + statistics are not numbers, or an argument is or holds a numpy masked + array. + ValueError: If ``gene_sets`` is empty, an option is out of range, a table + lacks a column, a statistic is infinite, a gene id is on two rows, fewer + than two genes are ranked, every statistic is zero, ``group`` or + ``contrast`` is given without a ``DifferentialExpressionResult`` or left + out of one that has several, or no gene of any set is among the ranked + genes. + KeyError: If the result has no table of the group and contrast; the message + says why. + + Example: + >>> import numpy as np + >>> genes = [f"g{i}" for i in range(500)] + >>> statistics = dict(zip(genes, np.linspace(2.0, -2.0, 500))) + >>> sets = {"near the top": genes[:50:2], "spread": genes[::20]} + >>> table = gsea(statistics, sets, rng=1) + >>> table["enrichment_score"].round(3).tolist() + [0.949, 0.184] + >>> table["normalized_enrichment_score"].round(2).tolist() + [3.51, 0.68] + """ + reject_masked(gene_sets, "gene_sets") + collection = gene_set_collection(gene_sets) + settings = _settings( + min_size, + max_size, + exponent, + score_type, + method, + n_permutations, + sample_size, + eps, + ) + ranking = ranked_genes( + statistics, + group=group, + contrast=contrast, + id_column=id_column, + statistic_column=statistic_column, + ) + generator = random_generator(rng) + with thread_limit(n_threads): + return _gene_set_enrichment(ranking, collection, settings, generator) + + +# The analysis. + + +def _gene_set_enrichment( + ranking: RankedGenes, + collection: GeneSets, + settings: _Settings, + generator: np.random.Generator, +) -> Table: + n = len(ranking.ids) + if not np.any(ranking.values): + raise ValueError( + "every statistic is zero, so the statistics do not rank the genes" + ) + # fgsea's weights: abs(stats)^gseaParam, in ranked order. + weights = np.ascontiguousarray(np.abs(ranking.values) ** settings.exponent) + _warn_about_statistics(ranking.values, settings.score_type) + index = {gene: position for position, gene in enumerate(ranking.ids)} + smallest = settings.min_size + largest = n - 1 if settings.max_size is None else min(settings.max_size, n - 1) + names: list[str] = [] + members: list[IntArray] = [] + matched = 0 + for name, genes in collection.items(): + positions = matched_positions(genes, index) + matched += positions.size + if smallest <= positions.size <= largest: + names.append(name) + members.append(np.sort(positions)) + if matched == 0: + raise ValueError( + "no gene of any set is among the ranked genes; the sets and the " + "statistics must name genes the same way (for example both by symbol, " + f"or both by Ensembl id): the sets name genes such as " + f"{_example_ids(collection)}, the statistics genes such as " + f"{quoted(ranking.ids[:3])}" + ) + seeds = _seeds(generator, settings) + if not names: + empty = np.zeros(0) + return _table( + [], np.zeros(0, dtype=np.int64), empty, empty, empty, empty, empty, [] + ) + sizes = np.array([positions.size for positions in members], dtype=np.int64) + scores = np.empty(len(names)) + edges = [] + for s, positions in enumerate(members): + scores[s], edge = _enrichment_score(weights, positions, settings.score_type) + edges.append(_SEPARATOR.join(ranking.ids[i] for i in edge.tolist())) + counts = _random_counts(weights, scores, sizes, settings, seeds) + normalized, p_value, beyond, of_sign = _permutation_results( + scores, counts, settings.score_type + ) + n_perm = settings.n_permutations + log2_error = np.full(len(names), np.nan) + if not settings.multilevel: + missing = int(np.count_nonzero(np.isnan(p_value))) + if missing: + _warn( + f"{missing} gene set(s) have no random score of the sign of their " + "enrichment score, or only zeros, so that their normalized scores and " + "p-values are NaN; the statistics are mostly of one sign" + ) + tested = ~np.isnan(p_value) + log2_error[tested] = _permutation_error(beyond[tested], n_perm) + else: + unbalanced = of_sign < _FEWEST_OF_SIGN + if unbalanced.any(): + _warn( + f"{int(np.count_nonzero(unbalanced))} gene set(s) have fewer than " + f"{_FEWEST_OF_SIGN} random scores of the sign of their enrichment " + f"score among {n_perm} permutations, as statistics mostly of one sign " + "give; their normalized scores and p-values are NaN. More permutations " + f"(n_permutations={10 * n_perm}) may give them some" + ) + normalized[unbalanced] = np.nan + p_value[unbalanced] = np.nan + tested = ~np.isnan(p_value) + multilevel = tested.copy() + multilevel[tested] = _needs_multilevel( + beyond[tested], n_perm, settings.sample_size + ) + simple = tested & ~multilevel + log2_error[simple] = _permutation_error(beyond[simple], n_perm) + if multilevel.any(): + rows = np.flatnonzero(multilevel) + denominator = (of_sign[rows] + 1.0) / (n_perm + 1.0) + estimated, good = _multilevel_p_values( + weights, scores[rows], sizes[rows], denominator, settings, seeds[1] + ) + p_value[rows] = np.minimum(1.0, estimated / denominator) + log2_error[rows] = np.where( + good, _multilevel_error(p_value[rows], settings.sample_size), np.nan + ) + if not good.all(): + _warn( + f"the p-values of {int(np.count_nonzero(~good))} gene set(s) are " + "likely too large: fewer than half of the random sets of the " + "multilevel method's last level had a positive signed score; " + "their log2_error is NaN" + ) + with np.errstate(invalid="ignore"): + floored = p_value < settings.eps + if floored.any(): + p_value[floored] = settings.eps + log2_error[floored] = np.nan + _warn( + f"the p-values of {int(np.count_nonzero(floored))} gene set(s) " + f"are below eps={settings.eps:g} and are reported as eps; eps=0 " + "estimates them" + ) + adjusted = adjust_p_values(p_value, "bh") + return _table( + names, sizes, scores, normalized, p_value, adjusted, log2_error, edges + ) + + +def _enrichment_score( + weights: FloatArray, positions: IntArray, score_type: int +) -> tuple[float, IntArray]: + """fgsea's calcGseaStat: the ES of a set and its leading edge, as positions.""" + tops, bottoms = running_sum(weights, positions) + highest = float(tops.max()) + lowest = float(bottoms.min()) + if score_type == 1 or (score_type == 0 and highest > -lowest): + score = highest + edge = positions[: int(np.argmax(tops)) + 1] + elif score_type == 2 or (score_type == 0 and highest < -lowest): + score = lowest + edge = positions[int(np.argmin(bottoms)) :][::-1] + else: # the deviations are equal: no direction + score = 0.0 + edge = positions[:0] + return score, edge + + +def _seeds(generator: np.random.Generator, settings: _Settings) -> IntArray: + """The seeds of the random sets, distinct integers from 1 to 10^9.""" + count = 2 if settings.multilevel else len(_chunks(settings.n_permutations)) + drawn = generator.choice(_SEED_RANGE, size=count, replace=False) + return np.asarray(drawn, dtype=np.int64) + 1 + + +def _chunks(n_permutations: int) -> list[int]: + """fgseaSimple's chunks: of max(1000, ceiling(n / 128)) permutations, and the + rest.""" + size = max(_CHUNK, -(-n_permutations // _MOST_CHUNKS)) + chunks = [size] * (n_permutations // size) + if n_permutations - sum(chunks) > 0: + chunks.append(n_permutations - sum(chunks)) + return chunks + + +class _Counts(NamedTuple): + """The random scores of each set: below or at its ES, above or at it, at most + and at least 0, and the sums of those at most and at least 0.""" + + below: FloatArray + above: FloatArray + at_most_zero: FloatArray + at_least_zero: FloatArray + negative_sum: FloatArray + positive_sum: FloatArray + + +def _random_counts( + weights: FloatArray, + scores: FloatArray, + sizes: IntArray, + settings: _Settings, + seeds: IntArray, +) -> _Counts: + if settings.multilevel: + chunks = np.array([settings.n_permutations], dtype=np.int64) + seeds = seeds[:1] + else: + chunks = np.array(_chunks(settings.n_permutations), dtype=np.int64) + arguments = (weights, scores, sizes, chunks, seeds, settings.score_type) + by_chunk: FloatArray = by_chunks( + _counts_by_chunk_parallel, _counts_by_chunk, arguments, chunks.size + ) + # The counts are whole numbers; the sums of the chunks are rounded once, as R's + # sum() of fgsea's chunks is. + totals = [by_chunk[:, column].sum(axis=0) for column in range(4)] + sums = [ + np.array([accurate_sum(by_chunk[:, column, s])[0] for s in range(sizes.size)]) + for column in (4, 5) + ] + return _Counts(*totals, *sums) + + +def _permutation_results( + scores: FloatArray, counts: _Counts, score_type: int +) -> tuple[FloatArray, FloatArray, FloatArray, FloatArray]: + """fgseaSimpleImpl: the NES, the permutation p-values, the random scores beyond + each ES on its side, and those of its sign (fgsea's modeFraction).""" + with np.errstate(divide="ignore", invalid="ignore"): + negative_mean = counts.negative_sum / counts.at_most_zero + positive_mean = counts.positive_sum / counts.at_least_zero + positive = scores > 0.0 + usable_positive = ~np.isnan(positive_mean) & (positive_mean != 0.0) + usable_negative = ~np.isnan(negative_mean) & (negative_mean != 0.0) + if score_type == 0: + normalizable = (positive & usable_positive) | (~positive & usable_negative) + beyond = np.where(positive, counts.above, counts.below) + of_sign = np.where(scores >= 0.0, counts.at_least_zero, counts.at_most_zero) + elif score_type == 1: + normalizable = (scores >= 0.0) & usable_positive + beyond = counts.above + of_sign = counts.at_least_zero + else: + normalizable = (scores <= 0.0) & usable_negative + beyond = counts.below + of_sign = counts.at_most_zero + normalized = np.full(scores.size, np.nan) + with np.errstate(divide="ignore", invalid="ignore"): + divisor = np.where(scores > 0.0, positive_mean, np.abs(negative_mean)) + if score_type == 1: + divisor = positive_mean + elif score_type == 2: + divisor = np.abs(negative_mean) + normalized[normalizable] = scores[normalizable] / divisor[normalizable] + p_value = np.minimum( + (1.0 + counts.below) / (1.0 + counts.at_most_zero), + (1.0 + counts.above) / (1.0 + counts.at_least_zero), + ) + p_value[~normalizable] = np.nan + return normalized, p_value, beyond, of_sign + + +def _permutation_error(beyond: FloatArray, n_permutations: int) -> FloatArray: + """fgsea's log2err of a permutation p-value.""" + variance = trigamma(beyond + 1.0) - trigamma(float(n_permutations + 1)) + error: FloatArray = np.sqrt(variance) / math.log(2.0) + return error + + +def _multilevel_error(p_value: FloatArray, sample_size: int) -> FloatArray: + """fgsea's multilevelError: the expected standard deviation of log2 p from the + number of levels.""" + level = trigamma((sample_size + 1) / 2) - trigamma(float(sample_size + 1)) + with np.errstate(divide="ignore"): + levels = np.floor(-np.log2(p_value) + 1.0) + error: FloatArray = np.sqrt(levels * level) / math.log(2.0) + return error + + +def _needs_multilevel( + beyond: FloatArray, n_permutations: int, sample_size: int +) -> BoolArray: + """Whether the multilevel method estimates log2 p more accurately than the + permutations did (fgseaMultilevel's comparison of multError and simpleError).""" + with np.errstate(divide="ignore"): + crude = np.log2((beyond + 1.0) / (n_permutations + 1.0)) + left = np.log2(_beta_quantile(0.025, beyond, n_permutations - beyond + 1.0)) + right = np.log2(_beta_quantile(0.975, beyond + 1.0, n_permutations - beyond)) + simple_error = 0.5 * np.maximum(crude - left, right - crude) + multiple_error = _multilevel_error( + (beyond + 1.0) / (n_permutations + 1.0), sample_size + ) + needed: BoolArray = multiple_error < simple_error + return needed + + +def _beta_quantile(q: float, a: FloatArray, b: FloatArray) -> FloatArray: + """The q-quantile of beta distributions with shapes a and b, by bisection of the + regularized incomplete beta function; 0 for a = 0 and 1 for b = 0, the point + masses of the limits.""" + a = np.asarray(a, dtype=np.float64) + b = np.asarray(b, dtype=np.float64) + low = np.zeros(a.shape) + high = np.ones(a.shape) + inside = (a > 0.0) & (b > 0.0) + for _ in range(1100): # far beyond the 1074 halvings that reach the smallest double + middle = 0.5 * (low + high) + settled = ~inside | (middle == low) | (middle == high) + if settled.all(): + break + below = np.asarray(betainc(a, b, middle), dtype=np.float64) < q + low = np.where(~settled & below, middle, low) + high = np.where(~settled & ~below, middle, high) + quantile: FloatArray = np.where(a == 0.0, 0.0, np.where(b == 0.0, 1.0, high)) + return quantile + + +def _multilevel_p_values( + weights: FloatArray, + scores: FloatArray, + sizes: IntArray, + denominator: FloatArray, + settings: _Settings, + seed: int, +) -> tuple[FloatArray, BoolArray]: + """fgsea's multilevelImpl: the multilevel estimates of the sets of each size, + with one seed for every size, and whether each is good (isCpGeHalf).""" + by_size = np.argsort(sizes, kind="stable") + group_sizes, starts = np.unique(sizes[by_size], return_index=True) + starts = np.append(starts, by_size.size).astype(np.int64) + ordered_scores = np.ascontiguousarray(scores[by_size]) + boundaries = np.array( + [ + settings.eps * float(denominator[by_size[starts[g] : starts[g + 1]]].min()) + for g in range(group_sizes.size) + ] + ) + arguments = ( + weights, + ordered_scores, + group_sizes.astype(np.int64), + starts, + boundaries, + settings.sample_size, + int(seed), + settings.score_type != 0, + ) + estimated, good = by_chunks( + _multilevel_groups_parallel, _multilevel_groups, arguments, group_sizes.size + ) + p_value = np.empty(scores.size) + p_value[by_size] = estimated + is_good = np.empty(scores.size, dtype=bool) + is_good[by_size] = good + return p_value, is_good + + +# The Mersenne Twister and fgsea's random integers. + + +@jit +def _mt_seeded(seed: int) -> IntArray: + """The state of std::mt19937 seeded with ``seed``, words in int64.""" + state = np.empty(_MT_WORDS + 1, dtype=np.int64) + state[0] = seed & 0xFFFFFFFF + for i in range(1, _MT_WORDS): + previous = state[i - 1] + state[i] = (1812433253 * (previous ^ (previous >> 30)) + i) & 0xFFFFFFFF + state[_MT_WORDS] = _MT_WORDS + return state + + +@jit +def _mt_next(state: IntArray) -> int: + """The next 32-bit output of MT19937.""" + index = state[_MT_WORDS] + if index >= _MT_WORDS: + for i in range(_MT_WORDS): + y = (state[i] & 0x80000000) | (state[(i + 1) % _MT_WORDS] & 0x7FFFFFFF) + word = state[(i + 397) % _MT_WORDS] ^ (y >> 1) + if y & 1: + word ^= 0x9908B0DF + state[i] = word + index = 0 + y = state[index] + state[_MT_WORDS] = index + 1 + y ^= y >> 11 + y ^= (y << 7) & 0x9D2C5680 + y ^= (y << 15) & 0xEFC60000 + y ^= y >> 18 + output: int = y + return output + + +@jit +def _uniform_integer(state: IntArray, low: int, high: int) -> int: + """An integer from low to high, as fgsea's uid_wrapper draws it: by rejecting the + outputs above the largest multiple of the range, then the remainder.""" + length = high - low + 1 + complete = _LARGEST_DRAW - _LARGEST_DRAW % length + while True: + x = _mt_next(state) + if x < complete: + drawn: int = low + x % length + return drawn + + +@jit +def _random_sample( + state: IntArray, low: int, high: int, k: int, used: BoolArray, out: IntArray +) -> None: + """k distinct integers from low to high in random order, as fgsea's + ``combination`` draws them; ``used`` is all False before and after.""" + n = high - low + 1 + count = 0 + if k < n * 1.0 / 2: + for _ in range(k): + for _ in range(_DRAWS_PER_GENE): + x = _uniform_integer(state, low, high) + if not used[x - low]: + out[count] = x + count += 1 + used[x - low] = True + break + else: + # Floyd's selection, as fgsea's, then the Fisher-Yates shuffle (Durstenfeld + # 1964) where fgsea calls std::shuffle. + for r in range(n - k, n): + x = _uniform_integer(state, 0, r) + if not used[x]: + out[count] = low + x + used[x] = True + else: + out[count] = low + r + used[r] = True + count += 1 + for i in range(count - 1, 0, -1): + j = _uniform_integer(state, 0, i) + out[i], out[j] = out[j], out[i] + if count < k: # 100 draws in a row fell on genes drawn already + raise RuntimeError("could not draw a random sample of distinct genes") + for i in range(count): + used[out[i] - low] = False + + +# The scores of random samples and of all their first genes (fgsea's +# calcGseaStatCumulative), from prefix sums kept in square-root blocks. + + +@jit +def _tree_shape(size: int) -> tuple[int, int, int]: + """The block length, its base-2 logarithm and the number of blocks of fgsea's + SegmentTree of ``size`` entries.""" + length = 1 + log_length = 0 + while length * length < size: + length <<= 1 + log_length += 1 + return length, log_length, (size - 1) // length + 1 + + +@jit +def _tree_add( + within: Any, between: Any, length: int, log_length: int, p: int, delta: Any +) -> None: + """Adds ``delta`` at entry ``p``: to the sums of its block from ``p`` on, and to + the sums before every later block.""" + mask = length - 1 + end = p - (p & mask) + mask + 1 + while p < end: + within[p] += delta + p += 1 + block = p >> log_length + while block < between.size: + between[block] += delta + block += 1 + + +@jit +def _tree_sum(within: Any, between: Any, log_length: int, r: int) -> Any: + """The sum of the entries before ``r``.""" + if r == 0: + return 0 + r -= 1 + return within[r] + between[r >> log_length] + + +@jit(error_model="numpy") +def _prefix_maxima( + weights: FloatArray, + selected: IntArray, + order: IntArray, + reverse: bool, + floor: float, + out: FloatArray, +) -> None: + """fgsea's gseaStats1: for each number i of the first genes of ``selected`` + (1-based positions, in the order drawn), the largest deviation of their walk + above zero, or below zero with ``reverse``. + + The walk at the end of the j-th selected gene in ranked order is the point + (x_j, y_j) of the genes passed that are not among the first i and of the weight + of those that are; its deviation is (y_j (N - i) / W - x_j) / (N - i). Adding a + gene adds (-1, w) to the points after it, which keeps the shape of the upper + convex hull of the points of each later block, so only the block of the gene is + rebuilt, and the point of each hull farthest along the current direction moves + one way as the genes are added. + """ + n = weights.size + k = selected.size + length, log_length, n_blocks = _tree_shape(k + 1) + x_within = np.zeros(length * length, dtype=np.int64) + x_between = np.zeros(n_blocks, dtype=np.int64) + y_within = np.zeros(length * length) + y_between = np.zeros(n_blocks) + rank = np.empty(k, dtype=np.int64) + for i in range(k): + rank[order[i]] = i + if not reverse: + previous = -1 + for i in range(k): + t = selected[order[i]] - 1 + _tree_add(x_within, x_between, length, log_length, i, t - previous) + previous = t + _tree_add(x_within, x_between, length, log_length, k, n - 1 - previous) + else: + previous = n + for i in range(k): + rank[i] = k - 1 - rank[i] + t = selected[order[k - 1 - i]] - 1 + _tree_add(x_within, x_between, length, log_length, i, previous - t) + previous = t + _tree_add(x_within, x_between, length, log_length, k, previous) + + # The hull of each block of points, as links to the previous point; a point off + # the hull has no next point (-1). + to_previous = np.empty(k + 2, dtype=np.int64) + to_next = np.empty(k + 2, dtype=np.int64) + block_length = int(math.sqrt(k + 1)) + n_hulls = (k + 1) // block_length + 1 + summit = np.empty(n_hulls, dtype=np.int64) + first = np.empty(n_hulls, dtype=np.int64) + last = np.empty(n_hulls, dtype=np.int64) + i = 0 + while i <= k + 1: + block = i // block_length + first[block] = i + last[block] = min(i + block_length - 1, k + 1) + j = 1 + while i + j <= last[block]: + to_previous[i + j] = i + j - 1 + to_next[i + j - 1] = i + j + j += 1 + to_previous[i] = i + to_next[last[block]] = last[block] + summit[block] = last[block] + i += block_length + + total = 0.0 + for i in range(k): + t = selected[i] - 1 + position = rank[i] + weight = weights[t] + if weight < floor: # max(|stat|, statEps), raised to gseaParam = 1 + weight = floor + _tree_add(x_within, x_between, length, log_length, position, -1) + _tree_add(y_within, y_between, length, log_length, position, weight) + total += weight + m = i + 1 + block = (position + 1) // block_length + start = first[block] + end = last[block] + top = max(position, start) + for c in range(position + 1, end + 1): + xc = float(_tree_sum(x_within, x_between, log_length, c)) + yc = _tree_sum(y_within, y_between, log_length, c) + b = top + xb = float(_tree_sum(x_within, x_between, log_length, b)) + yb = _tree_sum(y_within, y_between, log_length, b) + while to_previous[top] != top: + a = to_previous[top] + xa = float(_tree_sum(x_within, x_between, log_length, a)) + ya = _tree_sum(y_within, y_between, log_length, a) + turn = (xb - xa) * (yc - yb) - (yb - ya) * (xc - xb) + if yc - ya < _HULL_TOLERANCE: + turn = 0.0 + if turn <= 0.0: + break + top = a + to_next[b] = -1 + b = a + xb = xa + yb = ya + to_previous[c] = top + to_next[top] = c + top = c + if to_next[c] != -1: + break + coefficient = (n - m) / total + largest = 0.0 + summit[block] = max(top, summit[block]) + for h in range(n_hulls): + point = summit[h] + distance = _tree_sum( + y_within, y_between, log_length, point + ) * coefficient - _tree_sum(x_within, x_between, log_length, point) + while True: + after = to_previous[point] + farther = _tree_sum( + y_within, y_between, log_length, after + ) * coefficient - _tree_sum(x_within, x_between, log_length, after) + if farther <= distance: + break + distance = farther + point = after + summit[h] = point + if largest < distance: + largest = distance + out[i] = largest / (n - m) + + +@jit(error_model="numpy") +def _prefix_scores( + weights: FloatArray, selected: IntArray, score_type: int, out: FloatArray +) -> None: + """fgsea's calcGseaStatCumulative: the ES of each number of the first genes of + ``selected``.""" + k = selected.size + order = np.argsort(selected) + floor = _STATISTIC_FLOOR + for i in range(k): + x = abs(weights[selected[i] - 1]) + if x > 0.0 and x < floor: + floor = x + floor /= _FLOOR_DIVISOR + if score_type == 1: + _prefix_maxima(weights, selected, order, False, floor, out) + return + down = np.empty(k) + _prefix_maxima(weights, selected, order, True, floor, down) + if score_type == 2: + for i in range(k): + out[i] = -down[i] + return + _prefix_maxima(weights, selected, order, False, floor, out) + for i in range(k): + if out[i] == down[i]: + out[i] = 0.0 + elif out[i] < down[i]: + out[i] = -down[i] + + +@jit(error_model="numpy") +def _chunk_counts( + weights: FloatArray, + scores: FloatArray, + sizes: IntArray, + n_iterations: int, + seed: int, + score_type: int, + out: FloatArray, +) -> None: + """fgsea's calcGseaStatCumulativeBatch: ``out[column, set]`` holds the six + counts and sums of :class:`_Counts` over ``n_iterations`` random samples.""" + n = weights.size + k = sizes.max() + state = _mt_seeded(seed) + used = np.zeros(n, dtype=np.bool_) + selected = np.empty(k, dtype=np.int64) + random_scores = np.empty(k) + out[:] = 0.0 + for _ in range(n_iterations): + _random_sample(state, 1, n, k, used, selected) + _prefix_scores(weights, selected, score_type, random_scores) + for s in range(sizes.size): + r = random_scores[sizes[s] - 1] + if r <= scores[s]: + out[0, s] += 1.0 + if r >= scores[s]: + out[1, s] += 1.0 + if r <= 0.0: + out[2, s] += 1.0 + if r >= 0.0: + out[3, s] += 1.0 + out[4, s] = out[4, s] + (r if r < 0.0 else 0.0) + out[5, s] = out[5, s] + (r if r > 0.0 else 0.0) + + +@jit(parallel=True) +def _counts_by_chunk_parallel( + weights: FloatArray, + scores: FloatArray, + sizes: IntArray, + chunks: IntArray, + seeds: IntArray, + score_type: int, + n_runs: int, +) -> FloatArray: + """The counts of each chunk of permutations, ``(chunk, column, set)``, in + ``n_runs`` runs of chunks on numba's threads.""" + out = np.empty((chunks.size, 6, sizes.size)) + for run in prange(n_runs): + for c in range(run * chunks.size // n_runs, (run + 1) * chunks.size // n_runs): + _chunk_counts( + weights, scores, sizes, chunks[c], seeds[c], score_type, out[c] + ) + return out + + +@jit +def _counts_by_chunk( + weights: FloatArray, + scores: FloatArray, + sizes: IntArray, + chunks: IntArray, + seeds: IntArray, + score_type: int, +) -> FloatArray: + """The counts of :func:`_counts_by_chunk_parallel` on the calling thread.""" + out = np.empty((chunks.size, 6, sizes.size)) + for c in range(chunks.size): + _chunk_counts(weights, scores, sizes, chunks[c], seeds[c], score_type, out[c]) + return out + + +# The multilevel method (fgsea's EsRuler). + + +@jit(error_model="numpy") +def _signed_score(ranks: FloatArray, sample: IntArray) -> float: + """fgsea's calcES: the deviation of the walk of a sorted sample that is largest in + absolute value, the first of equal ones.""" + n = ranks.size + k = sample.size + total = 0.0 + for position in sample: + total += ranks[position] + result = 0.0 + current = 0.0 + q1 = 1.0 / (n - k) + q2 = 1.0 / total + last = -1 + for position in sample: + current -= q1 * (position - last - 1) + if abs(current) > abs(result): + result = current + current += q2 * ranks[position] + if abs(current) > abs(result): + result = current + last = position + return result + + +@jit(error_model="numpy") +def _positive_score(ranks: FloatArray, sample: IntArray) -> float: + """fgsea's calcPositiveES: the largest deviation of the walk above zero.""" + n = ranks.size + k = sample.size + total = 0.0 + for position in sample: + total += ranks[position] + result = 0.0 + current = 0.0 + q1 = 1.0 / (n - k) + q2 = 1.0 / total + last = -1 + for position in sample: + current += q2 * ranks[position] - q1 * (position - last - 1) + if result < current: # std::max, which keeps the first value against NaN + result = current + last = position + return result + + +@jit +def _appended(values: Any, count: int, value: Any) -> Any: + """``values`` with ``value`` at ``count``, enlarged when full.""" + if count == values.size: + larger = np.empty(2 * values.size, dtype=values.dtype) + larger[:count] = values + values = larger + values[count] = value + return values + + +@jit(error_model="numpy") +def _next_level( + ranks: FloatArray, + samples: IntArray, + levels: FloatArray, + positives: IntArray, + count: int, +) -> tuple[IntArray, FloatArray, IntArray, int]: + """fgsea's duplicateSamples: records the positive scores of the lower half of the + samples, the median included, with the number of the samples above them whose + signed score is positive, and replaces that half by copies of the upper half.""" + n_samples, k = samples.shape + positive_scores = np.empty(n_samples) + signed_positive = np.zeros(n_samples, dtype=np.bool_) + n_positive = 0 + for s in range(n_samples): + positive_scores[s] = _positive_score(ranks, samples[s]) + if _signed_score(ranks, samples[s]) > 0.0: + n_positive += 1 + signed_positive[s] = True + # Sorted by score, ties by sample, as std::sort sorts pairs. + order = np.argsort(positive_scores, kind="mergesort") + i = 0 + while 2 * i < n_samples: + levels = _appended(levels, count, positive_scores[order[i]]) + if signed_positive[order[i]]: + n_positive -= 1 + positives = _appended(positives, count, n_positive) + count += 1 + i += 1 + copies = np.empty((n_samples, k), dtype=np.int64) + row = 0 + i = 0 + while 2 * i < n_samples - 2: + for _ in range(2): + copies[row] = samples[order[n_samples - 1 - i]] + row += 1 + i += 1 + copies[row] = samples[order[n_samples >> 1]] + return copies, levels, positives, count + + +@jit +def _lower_bound(values: Any, start: int, end: int, value: Any) -> int: + """The first index from ``start`` before ``end`` whose value is not below + ``value``, as std::lower_bound finds it.""" + count = end - start + while count > 0: + step = count // 2 + middle = start + step + if values[middle] < value: + start = middle + 1 + count -= step + 1 + else: + count = step + return start + + +@jit +def _upper_bound(values: IntArray, value: int) -> int: + """The first index whose value is above ``value``, as std::upper_bound finds + it.""" + start = 0 + count = values.size + while count > 0: + step = count // 2 + middle = start + step + if not value < values[middle]: + start = middle + 1 + count -= step + 1 + else: + count = step + return start + + +@jit +def _spread_chunks( + slots: IntArray, starts: IntArray, sizes: IntArray, spare: int +) -> None: + """Moves the chunks of a sample apart so that each has ``spare`` free slots after + its genes, the last one all that remain.""" + n_chunks = sizes.size + genes = np.empty(sizes.sum(), dtype=np.int64) + count = 0 + for i in range(n_chunks): + for j in range(starts[i], starts[i] + sizes[i]): + genes[count] = slots[j] + count += 1 + count = 0 + for i in range(n_chunks): + starts[i + 1] = starts[i] + sizes[i] + spare + for j in range(sizes[i]): + slots[starts[i] + j] = genes[count] + count += 1 + starts[n_chunks] = slots.size + + +@jit(error_model="numpy") +def _perturb( + ranks: FloatArray, + slots: IntArray, + starts: IntArray, + sizes: IntArray, + chunk_sums: FloatArray, + chunk_ends: IntArray, + bound: float, + state: IntArray, +) -> int: + """fgsea's perturbate: about a tenth as many swaps of one gene of a sample for + another as it has genes, each kept when the sample's positive score stays above + ``bound``; returns the swaps kept. + + The sample's genes, sorted, are split into chunks of positions, as fgsea splits + them: chunk i holds ``sizes[i]`` genes from slot ``starts[i]`` of ``slots`` on, + with free slots before the next chunk; ``chunk_sums`` holds the sum of the ranks + of each chunk and ``chunk_ends`` the position after its last possible one. A + chunk whose whole sum cannot lift the walk above the bound is passed in one + step, and the first gene above the bound is remembered, so that later swaps that + keep it above need no walk. The walk of a proposed sample passes over the gene + that leaves and takes the gene that enters in its place, so that the sample is + changed only when the swap is kept; the sums are updated, and restored, in + fgsea's order of operations. + """ + n = ranks.size + n_chunks = chunk_sums.size + k = 0 + total = 0.0 + for i in range(n_chunks): + k += sizes[i] + for j in range(starts[i], starts[i] + sizes[i]): + total += ranks[slots[j]] + q1 = 1.0 / (n - k) + iterations = max(1, int(k * _SWAP_SHARE)) + moves = 0 + candidate = -1 + has_candidate = False + candidate_x = 0 + candidate_y = 0.0 + for _ in range(iterations): + old_index = _uniform_integer(state, 0, k - 1) + old_chunk = 0 + while sizes[old_chunk] <= old_index: + old_index -= sizes[old_chunk] + old_chunk += 1 + old_slot = starts[old_chunk] + old_index + old = slots[old_slot] + new = _uniform_integer(state, 0, n - 1) + new_chunk = _upper_bound(chunk_ends, new) + end = starts[new_chunk] + sizes[new_chunk] + target = _lower_bound(slots, starts[new_chunk], end, new) + if target < end and slots[target] == new: + if new == old: + moves += 1 + continue + total = total - ranks[old] + ranks[new] + chunk_sums[old_chunk] -= ranks[old] + chunk_sums[new_chunk] += ranks[new] + if has_candidate and old == candidate: + has_candidate = False + if has_candidate: + if old < candidate: + candidate_x += 1 + candidate_y -= ranks[old] + if new < candidate: + candidate_x -= 1 + candidate_y += ranks[new] + q2 = 1.0 / total + kept = has_candidate and -q1 * candidate_x + q2 * candidate_y > bound + if not kept: + x = 0 + y = 0.0 + last = -1 + for i in range(n_chunks): + size = sizes[i] + if i == old_chunk: + size -= 1 + if i == new_chunk: + size += 1 + if q2 * (y + chunk_sums[i]) - q1 * x < bound: + y += chunk_sums[i] + x += chunk_ends[i] - last - 1 - size + last = chunk_ends[i] - 1 + continue + j = starts[i] + stop = starts[i] + sizes[i] + pending = i == new_chunk # the new gene, which enters before target + while True: + if pending and j == target: + position = new + pending = False + elif j < stop: + position = slots[j] + j += 1 + if j - 1 == old_slot: + continue + else: + break + y += ranks[position] + x += position - last - 1 + if q2 * y - q1 * x > bound: + kept = True + has_candidate = True + candidate_x = x + candidate_y = y + candidate = position + break + last = position + if kept: + break + x += chunk_ends[i] - 1 - last + last = chunk_ends[i] - 1 + if kept: + moves += 1 + # The old gene leaves its chunk, and the new one enters its own. + stop = starts[old_chunk] + sizes[old_chunk] - 1 + for j in range(old_slot, stop): + slots[j] = slots[j + 1] + sizes[old_chunk] -= 1 + if new_chunk == old_chunk and old_slot < target: + target -= 1 + if starts[new_chunk] + sizes[new_chunk] == starts[new_chunk + 1]: + # The chunk is full: spread the chunks again. + offset = target - starts[new_chunk] + spare = (slots.size - (k - 1)) // n_chunks + _spread_chunks(slots, starts, sizes, spare) + target = starts[new_chunk] + offset + for j in range(starts[new_chunk] + sizes[new_chunk], target, -1): + slots[j] = slots[j - 1] + slots[target] = new + sizes[new_chunk] += 1 + else: + total = total - ranks[new] + ranks[old] + chunk_sums[old_chunk] += ranks[old] + chunk_sums[new_chunk] -= ranks[new] + if has_candidate and new == candidate: + has_candidate = False + if has_candidate: + if old < candidate: + candidate_x -= 1 + candidate_y += ranks[old] + if new < candidate: + candidate_x += 1 + candidate_y -= ranks[new] + return moves + + +@jit(error_model="numpy") +def _ruler( + ranks: FloatArray, sample_size: int, k: int, score: float, seed: int, eps: float +) -> tuple[FloatArray, IntArray]: + """fgsea's EsRuler::extend: the positive scores recorded at each level until + their median reaches ``score``, and the counts of positive signed scores.""" + n = ranks.size + state = _mt_seeded(seed) + used = np.zeros(n, dtype=np.bool_) + samples = np.empty((sample_size, k), dtype=np.int64) + for s in range(sample_size): + _random_sample(state, 0, n - 1, k, used, samples[s]) + samples[s].sort() + n_chunks = max(1, int(math.sqrt(k))) + chunk_ends = np.empty(n_chunks, dtype=np.int64) + chunk_ends[n_chunks - 1] = n + levels = np.empty(64) + positives = np.empty(64, dtype=np.int64) + samples, levels, positives, count = _next_level( + ranks, samples, levels, positives, 0 + ) + # Each chunk holds its genes and as many free slots as genes on average. + spare = (k + n_chunks - 1) // n_chunks + 4 + slots = np.empty((sample_size, k + spare * n_chunks), dtype=np.int64) + starts = np.empty((sample_size, n_chunks + 1), dtype=np.int64) + sizes = np.empty((sample_size, n_chunks), dtype=np.int64) + chunk_sums = np.empty((sample_size, n_chunks)) + column = np.empty(sample_size, dtype=np.int64) + half = (sample_size + 1) // 2 + while levels[count - 1] <= score - _LEVEL_TOLERANCE: + position = 0 + for i in range(n_chunks - 1): + position += (k + i) // n_chunks + for s in range(sample_size): + column[s] = samples[s, position] + chunk_ends[i] = np.sort(column)[sample_size // 2] + for s in range(sample_size): + chunk_sums[s, :] = 0.0 + sizes[s, :] = 0 + chunk = 0 + for j in range(k): + position = samples[s, j] + while chunk_ends[chunk] <= position: + chunk += 1 + sizes[s, chunk] += 1 + chunk_sums[s, chunk] += ranks[position] + starts[s, 0] = 0 + for i in range(n_chunks): + starts[s, i + 1] = starts[s, i] + sizes[s, i] + slots[s, :k] = samples[s] + _spread_chunks(slots[s], starts[s], sizes[s], spare) + moves = 0 + while moves < sample_size * k: + for s in range(sample_size): + moves += _perturb( + ranks, + slots[s], + starts[s], + sizes[s], + chunk_sums[s], + chunk_ends, + levels[count - 1], + state, + ) + for s in range(sample_size): + count_in = 0 + for i in range(n_chunks): + for j in range(starts[s, i], starts[s, i] + sizes[s, i]): + samples[s, count_in] = slots[s, j] + count_in += 1 + previous = levels[count - 1] + samples, levels, positives, count = _next_level( + ranks, samples, levels, positives, count + ) + if levels[count - 1] <= previous: + break # the median did not rise + if eps != 0.0 and count // half > -math.log2(0.5 * eps): + break + return levels[:count].copy(), positives[:count].copy() + + +@jit +def _beta_mean_log(a: float, b: float) -> float: + """fgsea's betaMeanLog: E log X for X ~ Beta(a, b + 1 - a), psi(a) - psi(b + 1).""" + return digamma_scalar(a) - digamma_scalar(b + 1.0) + + +@jit +def _ruler_p_value( + levels: FloatArray, + positives: IntArray, + sample_size: int, + score: float, + one_sided: bool, +) -> tuple[float, bool]: + """fgsea's EsRuler::getPvalue: the probability of a positive score of at least + ``score``, and of a positive signed one unless ``one_sided``; and whether the + estimate is good.""" + half = (sample_size + 1) // 2 + size = levels.size + good = True + if score >= levels[size - 1]: + index = size - 1 + if score > levels[size - 1] + _LEVEL_TOLERANCE: + good = False # beyond the last level + else: + index = _lower_bound(levels, 0, size, score) + full_levels = index // half + remainder = sample_size - index % half + log_level = _beta_mean_log(float(half), float(sample_size)) + log_p = full_levels * log_level + _beta_mean_log( + float(remainder + 1), float(sample_size) + ) + if not one_sided: + correction = _beta_mean_log(float(positives[index] + 1), float(remainder)) + log_p += correction + good = good and math.exp(correction) >= 0.5 + return max(0.0, min(1.0, math.exp(log_p))), good + + +@jit(error_model="numpy") +def _multilevel_group( + weights: FloatArray, + scores: FloatArray, + k: int, + eps: float, + sample_size: int, + seed: int, + one_sided: bool, + p_out: FloatArray, + good_out: BoolArray, +) -> None: + """fgsea's fgseaMultilevelCpp for the sets of one size.""" + largest = scores.max() + smallest = scores.min() + up_levels = np.empty(0) + up_positives = np.empty(0, dtype=np.int64) + down_levels = np.empty(0) + down_positives = np.empty(0, dtype=np.int64) + if largest >= 0.0: + up_levels, up_positives = _ruler( + weights, sample_size, k, abs(largest), seed, eps + ) + if smallest < 0.0: + down_levels, down_positives = _ruler( + weights[::-1].copy(), sample_size, k, abs(smallest), seed, eps + ) + for i in range(scores.size): + if scores[i] >= 0.0: + p_out[i], good_out[i] = _ruler_p_value( + up_levels, up_positives, sample_size, abs(scores[i]), one_sided + ) + else: + p_out[i], good_out[i] = _ruler_p_value( + down_levels, down_positives, sample_size, abs(scores[i]), one_sided + ) + + +@jit(parallel=True) +def _multilevel_groups_parallel( + weights: FloatArray, + scores: FloatArray, + group_sizes: IntArray, + starts: IntArray, + boundaries: FloatArray, + sample_size: int, + seed: int, + one_sided: bool, + n_runs: int, +) -> tuple[FloatArray, BoolArray]: + """The multilevel estimates of the groups of sets of each size, in ``n_runs`` + runs of groups on numba's threads.""" + p_value = np.empty(scores.size) + good = np.empty(scores.size, dtype=np.bool_) + n_groups = group_sizes.size + for run in prange(n_runs): + for g in range(run * n_groups // n_runs, (run + 1) * n_groups // n_runs): + _multilevel_group( + weights, + scores[starts[g] : starts[g + 1]], + group_sizes[g], + boundaries[g], + sample_size, + seed, + one_sided, + p_value[starts[g] : starts[g + 1]], + good[starts[g] : starts[g + 1]], + ) + return p_value, good + + +@jit +def _multilevel_groups( + weights: FloatArray, + scores: FloatArray, + group_sizes: IntArray, + starts: IntArray, + boundaries: FloatArray, + sample_size: int, + seed: int, + one_sided: bool, +) -> tuple[FloatArray, BoolArray]: + """The estimates of :func:`_multilevel_groups_parallel` on the calling thread.""" + p_value = np.empty(scores.size) + good = np.empty(scores.size, dtype=np.bool_) + for g in range(group_sizes.size): + _multilevel_group( + weights, + scores[starts[g] : starts[g + 1]], + group_sizes[g], + boundaries[g], + sample_size, + seed, + one_sided, + p_value[starts[g] : starts[g + 1]], + good[starts[g] : starts[g + 1]], + ) + return p_value, good + + +# Arguments, warnings and the table. + + +def _settings( + min_size: Any, + max_size: Any, + exponent: Any, + score_type: Any, + method: Any, + n_permutations: Any, + sample_size: Any, + eps: Any, +) -> _Settings: + for value, name in ( + (min_size, "min_size"), + (max_size, "max_size"), + (exponent, "exponent"), + (score_type, "score_type"), + (method, "method"), + (n_permutations, "n_permutations"), + (sample_size, "sample_size"), + (eps, "eps"), + ): + reject_masked(value, name) + smallest = positive_integer(min_size, "min_size") + largest = None + if max_size is not None: + largest = positive_integer(max_size, "max_size") + if largest < smallest: + raise ValueError( + f"max_size ({largest}) is smaller than min_size ({smallest}), so no " + "set could be tested" + ) + power = number_at_least(exponent, "exponent", 0.0) + kind = _SCORE_TYPES.index(checked_choice(score_type, "score_type", _SCORE_TYPES)) + how = checked_choice(method, "method", _METHODS) + count = positive_integer(n_permutations, "n_permutations") + samples = positive_integer(sample_size, "sample_size") + if samples < 3: + raise ValueError(f"sample_size must be at least 3, got {samples}") + boundary = real_number(eps, "eps") + if not 0.0 <= boundary < 1.0: + raise ValueError(f"eps must be in [0, 1), got {boundary}") + return _Settings( + smallest, + largest, + power, + kind, + how == "multilevel", + count, + samples + 1 if samples % 2 == 0 else samples, + boundary, + ) + + +def _warn(message: str) -> None: + warnings.warn(message, RuntimeWarning, skip_file_prefixes=(PACKAGE_DIRECTORY,)) + + +def _warn_about_statistics(values: FloatArray, score_type: int) -> None: + """fgsea's warnings about the statistics: ties, and all positive with "std".""" + nonzero = values[values != 0.0] + ties = nonzero.size - np.unique(nonzero).size + if ties: + _warn( + f"{ties} of the {values.size} statistics ({100 * ties / values.size:.2f}%) " + "tie with another one; genes with equal statistics keep the order in " + "which they are given, which decides the scores of sets that hold some " + "of them" + ) + if score_type == 0 and bool(np.all(values > 0.0)): + _warn( + "all statistics are positive and score_type is 'std'; for statistics " + "that are all positive, score_type='pos' tests for enrichment at the top" + ) + + +def _example_ids(collection: GeneSets) -> str: + examples: list[str] = [] + for genes in collection.values(): + for gene in genes: + if gene not in examples: + examples.append(gene) + if len(examples) == 3: + return quoted(examples) + return quoted(examples) + + +def _table( + names: list[str], + sizes: IntArray, + scores: FloatArray, + normalized: FloatArray, + p_value: FloatArray, + adjusted: FloatArray, + log2_error: FloatArray, + edges: list[str], +) -> Table: + values = ( + text_column(names), + sizes.astype(np.int64), + scores, + normalized, + p_value, + adjusted, + log2_error, + text_column(edges), + ) + return Table(dict(zip(_COLUMNS, values, strict=True))) diff --git a/cytepool/over_representation.py b/cytepool/over_representation.py index d853ba0..60e0272 100644 --- a/cytepool/over_representation.py +++ b/cytepool/over_representation.py @@ -114,6 +114,7 @@ from ._masked_arrays import real_array, reject_masked from ._messages import listing, quoted from ._options import real_number +from ._ranked_genes import gene_set_collection, group_and_contrast from .differential_expression import DifferentialExpressionResult from .gene_sets import GeneSets from .stats.distributions import hypergeom_sf @@ -197,7 +198,7 @@ def over_representation( says why. """ reject_masked(gene_sets, "gene_sets") - collection = _gene_sets(gene_sets) + collection = gene_set_collection(gene_sets) level = _alpha(alpha) columns = ( _column_name(id_column, "id_column"), @@ -207,7 +208,7 @@ def over_representation( else _column_name(fold_change_column, "fold_change_column"), ) if isinstance(genes, DifferentialExpressionResult): - chosen = _group_and_contrast(genes, group, contrast) + chosen = group_and_contrast(genes, group, contrast) selections, tested = _selections(genes.table(*chosen), level, *columns) else: if group is not None or contrast is not None: @@ -412,43 +413,6 @@ def _real_column(table: Table, column: str) -> np.ndarray: # Arguments. -def _gene_sets(gene_sets: Any) -> GeneSets: - if isinstance(gene_sets, GeneSets): - collection = gene_sets - elif isinstance(gene_sets, Mapping): - collection = GeneSets(gene_sets) - else: - raise TypeError( - "gene_sets must be a GeneSets or a mapping of set name to gene ids, such " - f"as {{'interferon': ['ISG15', 'IFIT1']}}; got {type(gene_sets).__name__}" - ) - if not len(collection): - raise ValueError("gene_sets is empty; give at least one gene set") - return collection - - -def _group_and_contrast( - result: DifferentialExpressionResult, group: Any, contrast: Any -) -> tuple[str, str]: - chosen = [] - for value, name, options in ( - (group, "group", list(result.groups)), - (contrast, "contrast", list(result.contrasts.names)), - ): - if value is None: - if len(options) != 1: - raise ValueError( - f"the result has {len(options)} {name}s; name one with {name}, " - f"one of {quoted(options)}" - ) - value = options[0] - elif not isinstance(value, str): - reject_masked(value, name) - raise TypeError(f"{name} must be a string, got {type(value).__name__}") - chosen.append(value) - return chosen[0], chosen[1] - - def _alpha(value: Any) -> float: reject_masked(value, "alpha") number = real_number(value, "alpha") diff --git a/cytepool/plotting.py b/cytepool/plotting.py index 76bb327..df63280 100644 --- a/cytepool/plotting.py +++ b/cytepool/plotting.py @@ -1,5 +1,5 @@ """Figures of pseudobulk and bulk analyses, drawn as limma, edgeR, DESeq2, speckle, -metafor and variancePartition draw them. +metafor, variancePartition and fgsea draw them. Each function draws one figure from the objects that CytePool's analyses return: a cube, a voom result, a fit, the table of a differential-expression or gene-set test, @@ -33,7 +33,10 @@ - Gene sets: :func:`plot_barcode`, limma's ``barcodeplot``; :func:`plot_gene_sets`, the sets of a ``camera``, ``fry`` or ``mroast`` table; :func:`plot_over_representation`, the sets of an over-representation analysis, as clusterProfiler's dot plots show - them (Yu et al. 2012). + them (Yu et al. 2012); :func:`plot_enrichment`, GSEA's enrichment plot (Subramanian + et al. 2005), and :func:`plot_gsea_table`, the top sets of a GSEA with their genes + on the ranking, as fgsea's ``plotEnrichment`` and ``plotGseaTable`` draw them + (Korotkevich et al. 2021). - Count models: :func:`plot_dispersions`, DESeq2's ``plotDispEsts`` (Love et al. 2014); :func:`plot_quasi_dispersions`, edgeR's ``plotQLDisp`` (Chen et al. 2016). - Meta-analysis: :func:`plot_forest`, the forest plot (Lewis and Clarke 2001) that @@ -112,6 +115,16 @@ Yu G, Wang LG, Han Y, He QY (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16, 284-287. doi:10.1089/omi.2011.0118 + Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich + A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP (2005). Gene set enrichment + analysis: a knowledge-based approach for interpreting genome-wide expression + profiles. Proceedings of the National Academy of Sciences 102, 15545-15550. + doi:10.1073/pnas.0506580102 + + 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 + Lewis S, Clarke M (2001). Forest plots: trying to see the wood and the trees. BMJ 322, 1479-1480. doi:10.1136/bmj.322.7300.1479 @@ -149,12 +162,20 @@ checked_choice, checked_flag, integral_number, + number_at_least, positive_integer, real_number, ) +from ._ranked_genes import ( + gene_set_collection, + matched_positions, + ranked_genes, + running_sum, +) from .composition import CompositionResult from .cube import Cube from .differential_expression import DifferentialExpressionResult +from .gene_sets import GeneSets from .meta_analysis import MetaAnalysisResult from .normalization import filter_by_expression, log_cpm, tmm_factors from .stats.deseq2 import DESeq2Dispersions, DESeq2Fit, DESeq2Results @@ -178,9 +199,11 @@ "plot_cell_numbers", "plot_composition", "plot_dispersions", + "plot_enrichment", "plot_forest", "plot_gene_sets", "plot_gene_variance_fractions", + "plot_gsea_table", "plot_md", "plot_mds", "plot_mean_variance", @@ -251,6 +274,12 @@ _LARGEST_DOT = 300.0 _LIST_LABELS = {"up": "Up", "down": "Down", "either": "Either", None: "Selected"} _TINY = np.finfo(np.float64).tiny # the smallest p-value shown on a log scale +# A GSEA table: the half-height of the bars of a row; where the headers of its columns +# stand, the middle of the ranks and the starts of the numbers, and its width, in +# numbers of genes ranked. +_TABLE_BAR = 0.45 +_TABLE_COLUMNS = (0.5, 1.04, 1.17, 1.37) +_TABLE_WIDTH = 1.55 _METHODS = { "reml": "REML", "dl": "DerSimonian-Laird", @@ -2387,6 +2416,274 @@ def plot_over_representation( return ax +def plot_enrichment( + statistics: Any, + gene_set: Any, + *, + group: str | None = None, + contrast: str | None = None, + id_column: str = "feature_id", + statistic_column: str = "t", + exponent: float = 1.0, + ax: "Axes | None" = None, +) -> "Axes": + """Plot the running enrichment score of a gene set along the ranked genes. + + GSEA's enrichment plot (Subramanian et al. 2005) as fgsea's ``plotEnrichment`` + draws it (Korotkevich et al. 2021): the genes ranked by decreasing statistic, as + :func:`~cytepool.gsea.gsea` ranks them, along the x-axis, from rank 0 to N + 1; + the walk that rises at each gene of the set by its share of the set's weight + |r|^p and falls at each other gene, drawn through its values just before and at + each gene of the set; a tick at the rank of each gene of the set, spanning a + sixteenth of the walk's range above and below zero; and dashed lines at the + walk's largest and smallest values, of which the larger in absolute value is the + enrichment score. + + Draws the walk (gid ``"running sum"``), the ticks (gid ``"set"``) and the dashed + lines (gid ``"extremes"``). + + Args: + statistics: The statistics, in any form :func:`~cytepool.gsea.gsea` takes. + gene_set: The genes of the set, as gene ids, such as ``gene_sets[name]``; + ids that are not ranked are left out. + group: The group of a ``DifferentialExpressionResult``, as for ``gsea``. + contrast: Its contrast, as for ``gsea``. + id_column: The column of gene ids of a table. + statistic_column: The column of statistics of a table. + exponent: The power p of the weights, at least 0, as for ``gsea``. + ax: The Axes to draw into, or None for a new figure. + + Returns: + The Axes. + + Raises: + TypeError: If an argument has the wrong type, or is or holds a numpy masked + array. + ValueError: If the statistics cannot be ranked, as for ``gsea``, the set holds + no ranked gene or every one, or ``exponent`` is negative or not finite. + KeyError: If a ``DifferentialExpressionResult`` has no table of the group and + contrast. + """ + _check_axes(ax) + power = number_at_least(exponent, "exponent", 0.0) + reject_masked(gene_set, "gene_set") + if isinstance(gene_set, str) or isinstance(gene_set, Mapping): + raise TypeError( + "gene_set must be a list of gene ids, such as gene_sets[name], got " + f"{type(gene_set).__name__}" + ) + genes = GeneSets({"gene_set": gene_set})["gene_set"] + ranking = ranked_genes( + statistics, + group=group, + contrast=contrast, + id_column=id_column, + statistic_column=statistic_column, + ) + n = len(ranking.ids) + index = {gene: i for i, gene in enumerate(ranking.ids)} + positions = np.sort(matched_positions(genes, index)) + if positions.size == 0: + raise ValueError("no gene of the set is among the ranked genes") + if positions.size == n: + raise ValueError( + "the set holds every ranked gene, so it has no enrichment score" + ) + weights = np.ascontiguousarray(np.abs(ranking.values) ** power) + tops, bottoms = running_sum(weights, positions) + ranks = positions + 1.0 + x = np.concatenate( + [[0.0], np.column_stack([ranks - 1.0, ranks]).ravel(), [n + 1.0]] + ) + y = np.concatenate([[0.0], np.column_stack([bottoms, tops]).ravel(), [0.0]]) + highest, lowest = float(tops.max()), float(bottoms.min()) + spread = highest - lowest + + ax = _axes(ax, (7.0, 3.6)) + ax.plot(x, y, color=_AQUA, linewidth=1.5, gid="running sum") + ax.vlines( + ranks, -spread / 16.0, spread / 16.0, color=_INK, linewidth=0.6, gid="set" + ) + for value in (highest, lowest): + ax.axhline(value, color=_UP, linewidth=1.0, linestyle="--", gid="extremes") + ax.axhline(0.0, color=_INK, linewidth=0.8) + ax.set_xlim(0.0, n + 1.0) + ax.set_xlabel("Rank") + ax.set_ylabel("Enrichment score") + return ax + + +def plot_gsea_table( + statistics: Any, + gene_sets: Mapping[str, Any], + table: Table, + *, + group: str | None = None, + contrast: str | None = None, + id_column: str = "feature_id", + statistic_column: str = "t", + top: int = 10, + exponent: float = 1.0, + ax: "Axes | None" = None, +) -> "Axes": + """Plot the most significant gene sets of a GSEA, each with its genes on the + ranking. + + The table of fgsea's ``plotGseaTable`` (Korotkevich et al. 2021): a row for each + of the ``top`` sets with a positive enrichment score and the smallest p-values + of a table of :func:`~cytepool.gsea.gsea`, most significant first, followed by + the ``top`` with a negative score, most significant last, as fgsea's guide + chooses them; the ties of p-values in the order of the table. Each row shows the + set's name, its genes as bars at their ranks among the genes ranked by + decreasing statistic, of height sign(r) |r|^p scaled by its largest absolute + value over all genes, and the set's normalized enrichment score, p-value and + adjusted p-value, written as fgsea writes them. + + Draws the bars (gid ``"gene ranks"``) and the numbers (gid ``"numbers"``). + + Args: + statistics: The statistics that ``table`` was computed from, in any form + :func:`~cytepool.gsea.gsea` takes. + gene_sets: The gene sets of ``table``: a + :class:`~cytepool.gene_sets.GeneSets` or a mapping of set name to gene + ids. + table: The table of ``gsea``, with the columns ``gene_set``, + ``enrichment_score``, ``normalized_enrichment_score``, ``p_value`` and + ``adjusted_p_value``. + group: The group of a ``DifferentialExpressionResult``, as for ``gsea``. + contrast: Its contrast, as for ``gsea``. + id_column: The column of gene ids of a table. + statistic_column: The column of statistics of a table. + top: How many sets of each sign to show, at least 1. + exponent: The power p of the heights of the bars, at least 0; fgsea suggests + 0.5 to flatten them. + ax: The Axes to draw into, or None for a new figure. + + Returns: + The Axes. + + Raises: + TypeError: If an argument has the wrong type, or is or holds a numpy masked + array. + ValueError: If the statistics cannot be ranked, the table lacks a column, has + no set with a p-value, or names a set that ``gene_sets`` does not hold, or + an option is out of range. + KeyError: If a ``DifferentialExpressionResult`` has no table of the group and + contrast. + """ + _check_axes(ax) + if not isinstance(table, Table): + raise TypeError(f"table must be a Table, got {type(table).__name__}") + count = positive_integer(top, "top") + power = number_at_least(exponent, "exponent", 0.0) + reject_masked(gene_sets, "gene_sets") + collection = gene_set_collection(gene_sets) + needed = ( + "gene_set", + "enrichment_score", + "normalized_enrichment_score", + "p_value", + "adjusted_p_value", + ) + missing = [column for column in needed if column not in table] + if missing: + raise ValueError( + f"the table lacks the column(s) {quoted(missing)}; plot_gsea_table reads " + "the tables of gsea" + ) + p = np.asarray(table["p_value"], dtype=np.float64) + scores = np.asarray(table["enrichment_score"], dtype=np.float64) + tested = np.isfinite(p) + if not tested.any(): + raise ValueError("no gene set of the table has a p-value") + up = np.flatnonzero(tested & (scores > 0.0)) + down = np.flatnonzero(tested & (scores < 0.0)) + up = up[np.argsort(p[up], kind="stable")[:count]] + down = down[np.argsort(p[down], kind="stable")[:count]][::-1] + rows = np.concatenate([up, down]).astype(np.intp) + names = table["gene_set"][rows].tolist() + absent = [name for name in names if name not in collection] + if absent: + raise ValueError( + f"gene_sets has no set {quoted(absent)}; give the gene sets that the " + "table was computed from" + ) + ranking = ranked_genes( + statistics, + group=group, + contrast=contrast, + id_column=id_column, + statistic_column=statistic_column, + ) + n = len(ranking.ids) + index = {gene: i for i, gene in enumerate(ranking.ids)} + heights = np.sign(ranking.values) * np.abs(ranking.values) ** power + largest = float(np.abs(heights).max()) + if largest > 0.0: + heights = heights / largest + + shown = len(rows) + ax = _axes(ax, (8.0, min(14.0, 1.4 + 0.4 * shown))) + xs: list[np.ndarray] = [] + bases: list[np.ndarray] = [] + ends: list[np.ndarray] = [] + for row, name in enumerate(names): + positions = matched_positions(collection[name], index) + base = float(shown - 1 - row) + xs.append(positions + 1.0) + bases.append(np.full(positions.size, base)) + ends.append(base + _TABLE_BAR * heights[positions]) + ax.vlines( + np.concatenate(xs), + np.concatenate(bases), + np.concatenate(ends), + color=_INK, + linewidth=0.5, + gid="gene ranks", + ) + levels = np.arange(shown - 1, -1, -1, dtype=np.float64) + ax.set_yticks(levels, names) + ax.tick_params(axis="y", length=0) + # The numbers stand in columns to the right of the ranks, inside the Axes. + ax.set_xlim(0.0, _TABLE_WIDTH * n) + ax.set_ylim(-0.6, shown + 0.4) + ticks = np.unique(np.round(np.linspace(0, n, 5)).astype(int)) + ax.set_xticks(ticks) + ax.spines[["top", "right"]].set_visible(False) + ax.spines["bottom"].set_bounds(0.0, float(n)) + ax.set_xlabel("Rank") + ax.xaxis.label.set_x(0.5 / _TABLE_WIDTH) + columns = ( + ("Gene ranks", None), + ("NES", np.asarray(table["normalized_enrichment_score"], dtype=np.float64)), + ("p-value", p), + ("Adjusted\np-value", np.asarray(table["adjusted_p_value"], dtype=np.float64)), + ) + for (header, values), start in zip(columns, _TABLE_COLUMNS, strict=True): + ax.text( + start * n, + shown - 0.3, + header, + ha="center" if values is None else "left", + va="bottom", + fontweight="bold", + ) + if values is None: + continue + for level, value in zip(levels.tolist(), values[rows].tolist(), strict=True): + text = f"{value:.2f}" if header == "NES" else _power_of_ten(value) + ax.text(start * n, level, text, ha="left", va="center", gid="numbers") + return ax + + +def _power_of_ten(value: float) -> str: + """A p-value as fgsea's tables write it, such as 1.2·10^-5, in mathtext.""" + if not math.isfinite(value) or value <= 0.0: + return f"{value:g}" + exponent = math.floor(math.log10(value)) + return rf"${value / 10.0**exponent:.1f} \cdot 10^{{{exponent}}}$" + + # Count models. diff --git a/cytepool/stats/rotation_gene_set_tests.py b/cytepool/stats/rotation_gene_set_tests.py index 798dd53..bc15581 100644 --- a/cytepool/stats/rotation_gene_set_tests.py +++ b/cytepool/stats/rotation_gene_set_tests.py @@ -178,7 +178,6 @@ import math import numbers -import operator from collections.abc import Mapping from typing import Any, NamedTuple @@ -188,7 +187,7 @@ from cytepool._jit import by_chunks, jit, thread_limit from cytepool._masked_arrays import reject_masked, unmasked_array -from cytepool._options import checked_flag, positive_integer +from cytepool._options import checked_flag, positive_integer, random_generator from cytepool._ranks import average_ranks_in_order from cytepool._summation import observed_means from cytepool.stats._gene_sets import ( @@ -515,7 +514,7 @@ def romer( shrink = checked_flag( shrink_contrast_effects, "shrink_contrast_effects", expected="a bool" ) - generator = _generator(rng) + generator = random_generator(rng) with thread_limit(n_threads): return _romer_sets( y, @@ -625,7 +624,7 @@ def _roast_tests( n_rotations = _rotation_count(n_rotations) with_trend = checked_flag(trend, "trend", expected="a bool") with_robust = checked_flag(robust, "robust", expected="a bool") - generator = _generator(rng) + generator = random_generator(rng) with thread_limit(n_threads): return _roast_sets( y, @@ -1604,20 +1603,3 @@ def _statistic(value: Any, options: tuple[str, ...]) -> int: def _rotation_count(value: Any) -> int: reject_masked(value, "n_rotations") return positive_integer(value, "n_rotations", expected="a positive integer") - - -def _generator(rng: Any) -> np.random.Generator: - reject_masked(rng, "rng") - 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( - "rng 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"rng must be a non-negative seed, got {seed}") - return np.random.default_rng(seed) diff --git a/docs/source/index.md b/docs/source/index.md index 0eab43a..7dbedb7 100644 --- a/docs/source/index.md +++ b/docs/source/index.md @@ -34,7 +34,7 @@ what CytePool offers for each today. | --- | --- | --- | | Which genes change in each cell type? | limma-voom, limma-trend and dream on every cell type of a cube; edgeR's quasi-likelihood tests and DESeq2 on one matrix at a time | [Differential expression per cell type](differential_expression.ipynb), [Bulk RNA-seq](guides/bulk_rna_seq.ipynb), [Choosing a method](statistics/choosing_a_method.md) | | Which cell types expand or contract? | propeller. Not yet: precision-weighted mixed models or count models of composition | [Composition tests](statistics/composition.md) | -| Which pathways and regulators respond? | camera, fry, roast, mroast and romer on any contrast; over-representation analysis of the genes that change; gene sets read from GMT files. Not yet: activity scores of pathways and transcription factors | [Gene-set tests](statistics/gene_set_tests.md) | +| Which pathways and regulators respond? | camera, fry, roast, mroast and romer on any contrast; gene set enrichment analysis (GSEA) of ranked genes; over-representation analysis of the genes that change; gene sets read from GMT files. Not yet: activity scores of pathways and transcription factors | [Gene-set tests](statistics/gene_set_tests.md) | | What drives the variation between donors, and what must be adjusted for? | Design checks for confounded factors, sample quality weights, mixed models that estimate a variance per random effect for every gene, and variance partitioning: the share of each gene's variance that each variable explains, in every cell type | [Designs, linear models and contrasts](statistics/linear_models.md), [Mixed models, REML and dream](statistics/mixed_models.md) | | Does it replicate, and how should studies be combined? | Meta-analysis of every gene across studies, with the between-study variance, Q and I². Not yet: replication measures and shrinkage across studies | [Meta-analysis](statistics/meta_analysis.md) | | Can patients be stratified from their profiles? | Not yet; cubes convert to AnnData and MuData for other tools | [Capability matrix](reference/capabilities.md) | diff --git a/docs/source/reference/api/gsea.md b/docs/source/reference/api/gsea.md new file mode 100644 index 0000000..01267d7 --- /dev/null +++ b/docs/source/reference/api/gsea.md @@ -0,0 +1,21 @@ +--- +myst: + html_meta: + description: Gene set enrichment analysis (GSEA) in CytePool, as fgsea computes it, of genes ranked by a statistic of a differential-expression result, a table of results or a mapping of gene id to statistic, with p-values by fgsea's multilevel method or by permutations, and its enrichment plot and table. +--- + +# Gene set enrichment analysis + +{func}`~cytepool.gsea.gsea` ranks genes by a statistic, of a group and contrast of a +differential-expression result, of a table of results or of a mapping of gene id to +statistic, and asks of each gene set whether its genes gather at the top or at the +bottom of the ranking, as fgsea's `fgsea` computes it: the enrichment score and its +leading edge, the score normalized by those of random sets of the same size, and +p-values by fgsea's multilevel method or by permutations of the genes. +{func}`~cytepool.plotting.plot_enrichment` draws the running enrichment score of one +set, and {func}`~cytepool.plotting.plot_gsea_table` the top sets with their genes on +the ranking. + +```{eval-rst} +.. automodule:: cytepool.gsea +``` diff --git a/docs/source/reference/api/index.md b/docs/source/reference/api/index.md index cb94aa2..2950fc4 100644 --- a/docs/source/reference/api/index.md +++ b/docs/source/reference/api/index.md @@ -38,6 +38,7 @@ gene_sets gene_set_tests rotation_gene_set_tests over_representation +gsea multiple_testing differential_expression variance_partitioning diff --git a/docs/source/reference/bibliography.md b/docs/source/reference/bibliography.md index 996d291..93d26af 100644 --- a/docs/source/reference/bibliography.md +++ b/docs/source/reference/bibliography.md @@ -47,6 +47,9 @@ Gene-set files and over-representation analysis [@gsea2026; @liberzon2011; @young2010; @fisher1922; @rivals2007; @tavazoie1999; @timmons2015; @wijesooriya2022; @yu2012]. +Gene set enrichment analysis [@subramanian2005; @mootha2003; @korotkevich2021] and the +random numbers of its tests [@matsumoto1998; @durstenfeld1964]. + ## Composition [@phipson2022; @buttner2021; @mangiola2023; @hoffman2026b; @dann2022; @yi2024] diff --git a/docs/source/reference/capabilities.md b/docs/source/reference/capabilities.md index 0835634..0657884 100644 --- a/docs/source/reference/capabilities.md +++ b/docs/source/reference/capabilities.md @@ -39,6 +39,7 @@ size; its sections give every check. | Rotation gene-set tests roast, mroast and romer | {api}`roast`, {api}`mroast`, {api}`romer` | limma 3.62.1 | calibrated; exact given limma's rotations | [`rotation_gene_set_tests.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/rotation_gene_set_tests.py) | | Gene-set files (GMT) read and written; collections matched to an analysis (limma's `ids2indices`) and filtered by size | {api}`read_gmt`, {api}`write_gmt`, {api}`GeneSets` | R 4.4.3, limma 3.62.1, decoupler 2.2.0 | exact | [`gene_set_files_and_ora.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/gene_set_files_and_ora.py) | | Over-representation analysis (limma's `kegga` and `goana`) of the genes up, down or either way, and its dot plot | {api}`over_representation`, {api}`plot_over_representation` | limma 3.62.1, R 4.4.3 | exact | [`gene_set_files_and_ora.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/gene_set_files_and_ora.py) | +| Gene set enrichment analysis of ranked genes (fgsea's `fgsea`: `fgseaMultilevel` and `fgseaSimple`), with its enrichment plot and table (`plotEnrichment`, `plotGseaTable`) | {api}`gsea`, {api}`plot_enrichment`, {api}`plot_gsea_table` | fgsea 1.32.2 | calibrated; exact given fgsea's seeds, and for enrichment scores, leading edges and figures | [`gsea.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/gsea.py) | | Adjustment for multiple testing (R's `p.adjust`) | {api}`adjust_p_values` | R 4.4.3 | exact | [`multiple_testing.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/multiple_testing.py) | | limma-voom and limma-trend on every cell type of a cube | {api}`limma_voom`, {api}`limma_trend` | limma 3.62.1, edgeR 4.4.0 | exact; close with robust moderation | [`differential_expression.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/differential_expression.py) | | Options for pseudobulk data: quality weights, weights from cell numbers, a minimum count, adjustment across cell types | {api}`limma_voom`, {api}`limma_trend` | limma 3.62.1, edgeR 4.4.0 | exact; close with quality weights | [`differential_expression_options.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/differential_expression_options.py) | diff --git a/docs/source/reference/citing.md b/docs/source/reference/citing.md index 01d7d5d..fea6044 100644 --- a/docs/source/reference/citing.md +++ b/docs/source/reference/citing.md @@ -34,6 +34,7 @@ the credit for the methods belongs to their authors. | {api}`roast`, {api}`mroast` | [@wu2010] | | {api}`romer` | [@majewski2010] | | {api}`over_representation` | limma's `kegga` and `goana` [@young2010]; the hypergeometric test [@tavazoie1999] | +| {api}`gsea` | fgsea [@korotkevich2021]; GSEA [@subramanian2005] | | {api}`read_gmt` | the collection you read, as its provider asks; for MSigDB, [@subramanian2005; @liberzon2011] | | {api}`propeller` | [@phipson2022] | | {api}`meta_analysis`, {api}`differential_expression_meta_analysis` | metafor [@viechtbauer2010]; the estimator of the between-study variance you chose: REML [@viechtbauer2005], DerSimonian-Laird [@dersimonian1986] or Paule-Mandel [@paule1982]; Knapp-Hartung [@knapp2003] | diff --git a/docs/source/reference/glossary.md b/docs/source/reference/glossary.md index ccea6d6..9938090 100644 --- a/docs/source/reference/glossary.md +++ b/docs/source/reference/glossary.md @@ -58,6 +58,11 @@ records](cube_files.md) such as MSigDB, often read from a GMT file. [Gene-set tests](../statistics/gene_set_tests.md#gene-set-files) +**Gene set enrichment analysis (GSEA).** A test of whether the genes of a set gather at +the top or the bottom of a ranking of all genes by a statistic, from the largest +deviation of a walk down the ranking, its enrichment score. [Gene-set +tests](../statistics/gene_set_tests.md#gene-set-enrichment-analysis) + **Group.** The second axis of a cube, usually a cell type: the cells of one sample are pooled per group. A bulk cube has one group. diff --git a/docs/source/reference/related_software.md b/docs/source/reference/related_software.md index 7c33ec0..ddd93f4 100644 --- a/docs/source/reference/related_software.md +++ b/docs/source/reference/related_software.md @@ -25,6 +25,7 @@ current releases below. | [speckle](https://bioconductor.org/packages/speckle) | propeller tests of cell-type proportions [@phipson2022] | GPL-3 | Bioconductor 1.12.0 | | [metafor](https://cran.r-project.org/package=metafor) | Meta-analysis [@viechtbauer2010] | GPL (>= 2) | CRAN 5.2-1 | | [lme4](https://cran.r-project.org/package=lme4) and [lmerTest](https://cran.r-project.org/package=lmerTest) | Linear mixed models; Satterthwaite's and Kenward-Roger's tests [@bates2015; @kuznetsova2017] | GPL (>= 2) | CRAN 2.1-0 and 3.2-1 | +| [fgsea](https://bioconductor.org/packages/fgsea) | Fast preranked gene set enrichment analysis, with p-values by multilevel splitting [@korotkevich2021]; CytePool's `gsea` adapts its algorithms, as `THIRD_PARTY_NOTICES.md` records | MIT | Bioconductor 1.38.0 | | [PyDESeq2](https://github.com/scverse/PyDESeq2) | A Python implementation of DESeq2 [@muzellec2023]; CytePool also adapts some of its constants, as `THIRD_PARTY_NOTICES.md` records | MIT | PyPI 0.5.4 | Because limma, edgeR, DESeq2, variancePartition, speckle, metafor, lme4 and lmerTest are @@ -58,7 +59,6 @@ neither uses nor reads them. | [crumblr](https://bioconductor.org/packages/crumblr) | Precision-weighted mixed models of cell-type composition [@hoffman2026b] | Artistic-2.0 | Bioconductor 1.4.5 | | [sccomp](https://bioconductor.org/packages/sccomp) | Bayesian tests of composition and its variability [@mangiola2023] | GPL-3 | Bioconductor 2.4.0 | | [miloR](https://bioconductor.org/packages/miloR) | Differential abundance in neighbourhoods of a k-nearest neighbour graph [@dann2022] | GPL-3 | Bioconductor 2.8.1 | -| [fgsea](https://bioconductor.org/packages/fgsea) | Fast preranked gene set enrichment analysis | MIT | Bioconductor 1.38.0 | | [glmGamPoi](https://bioconductor.org/packages/glmGamPoi) | Gamma-Poisson generalized linear models for single-cell counts | GPL-3 | Bioconductor 1.24.0 | | [pertpy](https://github.com/scverse/pertpy) | Perturbation analysis in the scverse ecosystem, with composition and differential expression tools [@heumos2026] | MIT | PyPI 1.4.0 | | [statsmodels](https://github.com/statsmodels/statsmodels) | Statistical models in Python, including generalized linear and mixed models | BSD 3-Clause | PyPI 0.15.0 | diff --git a/docs/source/references.bib b/docs/source/references.bib index 238cc71..c70473b 100644 --- a/docs/source/references.bib +++ b/docs/source/references.bib @@ -436,6 +436,17 @@ @book{dunn2018 doi = {10.1007/978-1-4419-0118-7}, } +@article{durstenfeld1964, + author = {Durstenfeld, Richard}, + title = {Algorithm 235: Random permutation}, + journal = {Communications of the ACM}, + year = {1964}, + volume = {7}, + number = {7}, + pages = {420}, + doi = {10.1145/364520.364540}, +} + @article{fai1996, author = {Fai, Alex Hrong-Tai and Cornelius, Paul L.}, title = {Approximate F-tests of multiple degree of freedom hypotheses in generalized least squares analyses of unbalanced split-plot experiments}, @@ -892,6 +903,14 @@ @book{knuth1997 year = {1997}, } +@misc{korotkevich2021, + author = {Korotkevich, Gennady and Sukhov, Vladimir and Budin, Nikolay and Shpak, Boris and Artyomov, Maxim N. and Sergushichev, Alexey}, + title = {Fast gene set enrichment analysis}, + howpublished = {bioRxiv 060012, version 3}, + year = {2021}, + doi = {10.1101/060012}, +} + @article{kuznetsova2017, author = {Kuznetsova, Alexandra and Brockhoff, Per B. and Christensen, Rune H. B.}, title = {lmerTest package: tests in linear mixed effects models}, @@ -1104,6 +1123,17 @@ @article{marquardt1963 doi = {10.1137/0111030}, } +@article{matsumoto1998, + author = {Matsumoto, Makoto and Nishimura, Takuji}, + title = {Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator}, + journal = {ACM Transactions on Modeling and Computer Simulation}, + year = {1998}, + volume = {8}, + number = {1}, + pages = {3--30}, + doi = {10.1145/272991.272995}, +} + @article{mccarthy2009, author = {McCarthy, Davis J. and Smyth, Gordon K.}, title = {Testing significance relative to a fold-change threshold is a TREAT}, @@ -1133,6 +1163,17 @@ @misc{miles2025 year = {2025}, } +@article{mootha2003, + author = {Mootha, Vamsi K. and Lindgren, Cecilia M. and Eriksson, Karl-Fredrik and Subramanian, Aravind and Sihag, Smita and Lehar, Joseph and Puigserver, Pere and Carlsson, Emma and Ridderstråle, Martin and Laurila, Esa and Houstis, Nicholas and Daly, Mark J. and Patterson, Nick and Mesirov, Jill P. and Golub, Todd R. and Tamayo, Pablo and Spiegelman, Bruce and Lander, Eric S. and Hirschhorn, Joel N. and Altshuler, David and Groop, Leif C.}, + title = {PGC-1α-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes}, + journal = {Nature Genetics}, + year = {2003}, + volume = {34}, + number = {3}, + pages = {267--273}, + doi = {10.1038/ng1180}, +} + @article{murphy2022, author = {Murphy, Alan E. and Skene, Nathan G.}, title = {A balanced measure shows superior performance of pseudobulk methods in single-cell RNA-sequencing analysis}, diff --git a/docs/source/statistics/choosing_a_method.md b/docs/source/statistics/choosing_a_method.md index d1f42d5..0d1fa5f 100644 --- a/docs/source/statistics/choosing_a_method.md +++ b/docs/source/statistics/choosing_a_method.md @@ -66,6 +66,7 @@ edgeR's and DESeq2's methods run on one matrix at a time, not on every group of | Question | Functions | Explained in | | --- | --- | --- | | Do the genes of a pathway respond together? | {api}`camera`, {api}`fry`, {api}`roast`, {api}`mroast`, {api}`romer` | [Gene-set tests](gene_set_tests.md) | +| Do the genes of a pathway gather at the top or the bottom of a ranking of all genes? | {api}`gsea`, with {api}`plot_enrichment` and {api}`plot_gsea_table` | [Gene-set tests](gene_set_tests.md) | | Are pathways over-represented among the genes that change? | {api}`over_representation`, with sets from {api}`read_gmt` | [Gene-set tests](gene_set_tests.md) | | Does the share of a cell type change? | {api}`propeller` | [Composition tests](composition.md) | | Do studies agree, and what is their common effect? | {api}`differential_expression_meta_analysis`, {api}`meta_analysis` | [Meta-analysis](meta_analysis.md) | diff --git a/docs/source/statistics/gene_set_tests.md b/docs/source/statistics/gene_set_tests.md index 1a9aab9..b5db82b 100644 --- a/docs/source/statistics/gene_set_tests.md +++ b/docs/source/statistics/gene_set_tests.md @@ -1,5 +1,5 @@ --- -description: Competitive and self-contained gene-set tests on a linear model, as limma's camera, fry, roast, mroast and romer compute them; over-representation analysis of the genes called significant, as limma's kegga computes it; gene sets read from GMT files; and how to choose between the tests. +description: Competitive and self-contained gene-set tests on a linear model, as limma's camera, fry, roast, mroast and romer compute them; gene set enrichment analysis of ranked genes, as fgsea computes it; over-representation analysis of the genes called significant, as limma's kegga computes it; gene sets read from GMT files; and how to choose between the tests. --- # Gene-set tests @@ -10,16 +10,17 @@ of the linear model fitted to every gene. Two questions differ [@goeman2007]: - **Self-contained**: are the genes of the set differentially expressed at all? {api}`fry`, {api}`roast` and {api}`mroast`. - **Competitive**: are they more differentially expressed than the other genes? - {api}`camera` and {api}`romer`, and {api}`over_representation`, which asks it of the - genes called significant. + {api}`camera` and {api}`romer`; {api}`gsea`, which asks it of a ranking of all genes; + and {api}`over_representation`, which asks it of the genes called significant. A self-contained test finds a set significant when its genes respond, however many other genes respond too: in a response that touches thousands of genes, many sets are. A competitive test asks whether a set stands out from the other genes, and so suits ranking many sets, such as a pathway collection. The five tests of limma work on the linear model, so they allow any design and contrast and do not permute samples; all but -romer take voom's precision weights. Over-representation analysis needs only the genes -that an analysis called significant. +romer take voom's precision weights. Gene set enrichment analysis needs only a statistic +of every gene, and over-representation analysis only the genes that an analysis called +significant. ## The effects of each gene @@ -144,6 +145,49 @@ and highly expressed genes are called significant more often, which favours sets such genes [@young2010]; the adjustment that `kegga` and `goana` offer for this is not implemented. +## Gene set enrichment analysis + +Gene set enrichment analysis (GSEA) ranks all genes by a statistic, such as the moderated +t-statistic of a contrast, and asks whether the genes of a set gather at the top or at +the bottom of the ranking [@subramanian2005]. A walk down the ranking of $N$ genes rises +at each of the $m$ genes of the set by its share of the set's weight and falls at each +other gene by $1/(N - m)$; after the $i$-th gene of the set, at rank $p_i$, it stands at + +$$ +\mathrm{RS}_i = \frac{\sum_{j \le i} |r_j|^q}{\sum_{j = 1}^{m} |r_j|^q} +- \frac{p_i - i}{N - m}, +$$ + +with $r_j$ the statistics of the set's genes and $q$ the exponent: 1 for the weighted +score of GSEA, 0 for the unweighted one of @mootha2003. The enrichment score is the +walk's largest deviation from zero, positive when the set's genes gather at the top and +negative when they gather at the bottom, and the leading edge is the set's genes that +the walk passes before its extreme, which drive the score. {api}`gsea` computes them as +fgsea does [@korotkevich2021], from a group and contrast of a differential-expression +result, a table of results or a mapping of gene id to statistic. Genes with equal +statistics keep the order in which they are given, as in fgsea, so that the score of a +set whose genes tie with others depends on that order; a warning says how many tie. + +The scores of random sets of the same size give the null distribution, so that the test +is competitive, like camera's, and it permutes genes, not samples. Like +over-representation analysis it treats genes as independent: the genes of a set that +respond together make its p-values too small [@goeman2007; @wu2012]. Each score is +divided by the mean of the random scores of its sign, so that sets of different sizes +compare, and the p-value is the share of the random scores of its sign that reach it. +With $n$ random sets it is at least $1/(n + 1)$; fgsea's multilevel method, the default, +reaches far smaller ones [@korotkevich2021]. It draws random sets of the size of the set, +keeps at each level the half with the higher scores and moves them by a Markov chain that +keeps their scores above the level's median, so that each level about halves the +probability, until the median reaches the observed score; `log2_error` gives the expected +standard deviation of $\log_2 p$. + +The random sets come from seeds that `rng` draws, through the Mersenne Twister +[@matsumoto1998] seeded as C++ seeds it, as fgsea draws them: given fgsea's seeds, +{api}`gsea` gives fgsea's results exactly, and a seed gives the same results every time. +{api}`plot_enrichment` draws the walk of one set, as fgsea's `plotEnrichment` does, and +{api}`plot_gsea_table` the top sets with their genes on the ranking, as its +`plotGseaTable` does. + ## Choices - **Competitive or self-contained.** Rank a collection with camera; ask about one set @@ -159,3 +203,7 @@ implemented. genes called significant, so it takes the results of any method and lists from elsewhere. camera uses every gene's statistic and allows for the correlation between genes, so with a linear model at hand it is the sounder competitive test. +- **GSEA or camera.** GSEA needs only a statistic of every gene, from any method, and + uses the whole ranking rather than a cutoff; like over-representation analysis it + takes genes as independent, which camera does not. Set `min_size` and `max_size`, as + GSEA's 15 and 500, to leave out sets too small to rely on and too large to mean much. diff --git a/docs/source/validation.md b/docs/source/validation.md index 452ecfd..cfa1633 100644 --- a/docs/source/validation.md +++ b/docs/source/validation.md @@ -15,10 +15,10 @@ they compare bit for bit), runs CytePool, prints a table with one row per check, exits with an error when any check falls outside its tolerance. This page summarises the tables the scripts printed when they were run against the -package in this revision of the repository: on **2026-10-26** for `plotting.py` and -`variance_partitioning.py` and `gene_set_files_and_ora.py` and -`kenward_roger_and_f_tests.py`, and for `dream.py`, which printed the same table again -after changes to the code it checks; on +package in this revision of the repository: on **2026-10-26** for `plotting.py`, +`variance_partitioning.py`, `gene_set_files_and_ora.py`, `kenward_roger_and_f_tests.py` +and `gsea.py`, and for `dream.py`, which printed the same table again after changes to +the code it checks; on **2026-10-25** for the six scripts on cubes and their files (`cube.py`, `pool_cells.py`, `cube_files.py`, `cube_reads.py`, `cube_fingerprints.py` and `anndata_conversion.py`), for `linear_models.py`, @@ -56,6 +56,7 @@ covers and the largest difference among them. Run a script to see every check. | [`gene_set_tests.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/gene_set_tests.py) | Gene-set tests | limma 3.62.1 | 93 | 79 pass, 8 deliberate differences, 6 information | 6.9 s | 278 MiB | | [`rotation_gene_set_tests.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/rotation_gene_set_tests.py) | Rotation gene-set tests | limma 3.62.1 | 103 | 81 pass, 20 deliberate differences, 2 information | 20.9 s | 190 MiB | | [`gene_set_files_and_ora.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/gene_set_files_and_ora.py) | Gene-set files and over-representation | R 4.4.3, limma 3.62.1, decoupler 2.2.0 | 34 | 29 pass, 5 deliberate differences | 4.7 s | 297 MiB | +| [`gsea.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/gsea.py) | Gene set enrichment analysis | R 4.4.3, fgsea 1.32.2 | 24 | 23 pass, 1 deliberate difference | 7 min 3 s | 323 MiB | | [`multiple_testing.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/multiple_testing.py) | Multiple testing | R 4.4.3 | 15 | 15 pass | 2 min 56 s | 202 MiB | | [`differential_expression.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/differential_expression.py) | Differential expression | limma 3.62.1, edgeR 4.4.0 | 32 | 28 pass, 4 deliberate differences | 1.8 s | 189 MiB | | [`differential_expression_options.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/differential_expression_options.py) | Options for pseudobulk differential expression | limma 3.62.1, edgeR 4.4.0 | 66 | 66 pass | 5.0 s | 182 MiB | @@ -68,7 +69,7 @@ covers and the largest difference among them. Run a script to see every check. | [`plotting.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/plotting.py) | Figures | limma 3.62.1, edgeR 4.4.0, DESeq2 1.46.0, metafor 5.0.1, speckle 1.6.0 | 26 | 23 pass, 2 information, 1 deliberate difference | 25.4 s | 810 MiB | | [`distributions.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/distributions.py) | Distributions and special functions | mpmath 1.3.0 | 46 | 46 pass | 6 min 39 s | 221 MiB | | [`smoothing.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/smoothing.py) | LOWESS and spline bases | R 4.4.3 | 90 | 87 pass, 1 information, 2 deliberate differences | 1.5 s | 175 MiB | -| **All** | | | **2,150** | **1,851 pass, none failed; 172 information, 126 deliberate differences, 1 control** | | | +| **All** | | | **2,174** | **1,874 pass, none failed; 172 information, 127 deliberate differences, 1 control** | | | Run time is the wall time of the whole script and peak memory the largest resident set of any of its processes, reference programs included, both measured with GNU @@ -1182,6 +1183,64 @@ Deliberate differences: genes the analysis called significant. With the p-values adjusted over the features with a symbol, CytePool gives kegga's values; the tests check this. +## Gene set enrichment analysis + +[`validation/gsea.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/gsea.py) +checks gene set enrichment analysis and its figures against fgsea. Statistics of 4,000 +genes are simulated from fixed seeds like moderated t-statistics, 6% of the genes up and +5% down, with 60 genes tied at 1.25 and 30 at zero, with 90 gene sets of 1 to 1,000 +genes, a fifth of each outside the ranking: sets enriched for genes up, down or both, or +for none, a set of every gene, which neither program tests, a set of one gene, half of +the tied genes, a set that lists genes twice and two identical sets. `gsea` and fgsea's +`fgsea` run with the same options: fgsea's defaults (its multilevel method, 1,000 +permutations, 101 sets per level, eps 1e-50); gseaParam 0 and 0.5; scoreType "neg"; +scoreType "pos" on -log10 p-values, which are all positive; sets of 15 to 500 genes, +with eps 0, 51 sets per level and 2,000 permutations; permutations alone (fgsea's +`nperm`, 2,500 in three chunks); statistics mostly positive, with 6 sets spread evenly +along the ranking, a set of a few genes near its top and others spread along it, whose +p-value is about 0.05, and 3 sets of 25 genes from its bottom, whose negative scores +fewer than 10 random sets of their size reach, so that no set goes to the multilevel +method; and one set alone. + +The R script records the seeds that fgsea draws where it hands them to its compiled code +(`trace` of `fgseaSimpleImpl` and `multilevelImpl`), and `gsea` draws the same through +`rng`; the two programs then draw the same random sets, and must agree exactly. fgsea +draws the random sets of a set tested alone with R's `sample.int`, which cannot be +replayed: that case is exact for the enrichment score and leading edge, and calibrated +for the rest, with seeds of their own (R's 1 to 20 and numpy's 1 to 20): the mean +log10 p-value and the mean normalized score of each set over 20 runs must agree within +4 standard errors of their difference, for 36 sets and for the set alone. The figures: +what `plot_enrichment` draws against `plotEnrichmentData`, for three sets with gseaParam +1 and 0.5, and the bars and numbers of `plot_gsea_table` against the ticks of +`plotEnrichmentData`, scaled as `plotGseaTable` scales them, and its numbers as it +writes them. Reference: R 4.4.3, fgsea 1.32.2. + +| Check family | Reference | Largest difference | Tolerance | Outcome | +| --- | --- | --- | --- | --- | +| The sets tested, their sizes, leading edges and enrichment scores, in all 9 cases | `fgsea` | 0 | 1e-10 | 9 pass | +| Normalized scores, p-values, log2 errors and adjusted p-values, given fgsea's seeds, in 7 cases | `fgsea` | 2.2e-13 | 1e-10 | 7 pass | +| The same of the mostly positive statistics, without the adjusted p-values | `fgsea` | 4.2e-16 | 1e-10 | 1 pass | +| Mean log10 p-values and mean normalized scores over 20 runs, seeds of their own, for 36 sets and one set alone | `fgsea` | 2.4 standard errors | 4 | 4 pass | +| `plot_enrichment`: the walk, the ticks and the extremes of 6 plots | `plotEnrichmentData` | 1.7e-14 | 1e-10 | 1 pass | +| `plot_gsea_table`: the bars and numbers of 10 sets | `plotEnrichmentData`, `plotGseaTable` | 1.1e-13 | 1e-10 | 1 pass | +| Adjusted p-values of the mostly positive statistics | `fgsea` | 0.16 (relative) | n/a | 1 deliberate difference | + +Numbers are compared relative to their absolute value, or 1e-300, and NaN must stand +where fgsea gives NA; sets, sizes and leading edges must be identical. Given fgsea's +seeds, the random scores and the counts of the permutations are the same to the last +bit; the multilevel p-values differ in the last bits, through the digamma function of +each program. + +Deliberate difference: + +- The adjusted p-values of the mostly positive statistics, where some sets have too few + random scores of their sign and no set goes to the multilevel method: fgsea adjusts + the p-values over the sets before it removes those that have some, but fewer than 10, + random scores of their sign, as though they had been tested, so that its adjusted + p-values of the sets it reports come out smaller or larger; CytePool adjusts over the + sets with a p-value. The tests check that CytePool's are Benjamini and Hochberg's of + the p-values it reports. + ## Multiple testing [`validation/multiple_testing.py`](https://github.com/NygenAnalytics/CytePool/blob/master/validation/multiple_testing.py) @@ -2024,6 +2083,11 @@ the run was 424 MiB. | Two crossed: 40 donors, 6 batches | 1.50 | 0.06 | 0.052 | 17.5 | | Two nested: 20 donors x 3 samples x 2 replicates | 4.12 | 0.06 | 0.051 | 17.0 | +On fgsea's example data, 12,000 genes and 586 sets of 15 to 500 genes, with its +defaults, fgsea's `fgsea` took 2.9 s and `gsea` 4.1 s, both on one thread, after a +first call of `gsea` that compiled it; with plain permutations (10,000), fgsea took +1.9 s and `gsea` 1.5 s, on 2 threads. + `adjust_p_values` with Hommel's method took 0.006 s on the 100,000 p-values of `multiple_testing.py`. In `distributions.py`, mpmath took 395 s to compute the references and CytePool 3.6 s to evaluate them, including compilation. @@ -2074,6 +2138,7 @@ python validation/empirical_bayes.py python validation/gene_set_tests.py python validation/rotation_gene_set_tests.py python validation/gene_set_files_and_ora.py +RSCRIPT=$REF/r-gsea/bin/Rscript python validation/gsea.py python validation/voom.py python validation/sample_quality_weights.py python validation/negative_binomial_fits.py @@ -2159,6 +2224,16 @@ conda create -p $REF/r-dream -c conda-forge -c bioconda r-base=4.4.3 \ r-minqa=1.2.8 r-nloptr=2.2.1 r-jsonlite=2.0.0 ``` +fgsea 1.32.2 (Bioconductor 3.20), which `gsea.py` runs, has an R environment of its own, +with R 4.4.3, limma 3.62.1, GSVA 2.0.0, BiocParallel 1.40.0, data.table 1.18.6.1 and +jsonlite 2.0.0, from conda-forge and bioconda: + +```bash +conda create -p $REF/r-gsea -c conda-forge -c bioconda r-base=4.4.3 \ + bioconductor-fgsea=1.32.2 bioconductor-gsva=2.0.0 bioconductor-limma=3.62.1 \ + r-data.table=1.18.6.1 r-jsonlite=2.0.0 +``` + The Python environments use Python 3.13.12. scanpy has its own environment, and `cube_files.py` ran CytePool on Python 3.14.4 for the zstd cases: diff --git a/tests/data/gsea/fgsea.json b/tests/data/gsea/fgsea.json new file mode 100644 index 0000000..84f53ff --- /dev/null +++ b/tests/data/gsea/fgsea.json @@ -0,0 +1 @@ +{"genes":["G00000","G00001","G00002","G00003","G00004","G00005","G00006","G00007","G00008","G00009","G00010","G00011","G00012","G00013","G00014","G00015","G00016","G00017","G00018","G00019","G00020","G00021","G00022","G00023","G00024","G00025","G00026","G00027","G00028","G00029","G00030","G00031","G00032","G00033","G00034","G00035","G00036","G00037","G00038","G00039","G00040","G00041","G00042","G00043","G00044","G00045","G00046","G00047","G00048","G00049","G00050","G00051","G00052","G00053","G00054","G00055","G00056","G00057","G00058","G00059","G00060","G00061","G00062","G00063","G00064","G00065","G00066","G00067","G00068","G00069","G00070","G00071","G00072","G00073","G00074","G00075","G00076","G00077","G00078","G00079","G00080","G00081","G00082","G00083","G00084","G00085","G00086","G00087","G00088","G00089","G00090","G00091","G00092","G00093","G00094","G00095","G00096","G00097","G00098","G00099","G00100","G00101","G00102","G00103","G00104","G00105","G00106","G00107","G00108","G00109","G00110","G00111","G00112","G00113","G00114","G00115","G00116","G00117","G00118","G00119","G00120","G00121","G00122","G00123","G00124","G00125","G00126","G00127","G00128","G00129","G00130","G00131","G00132","G00133","G00134","G00135","G00136","G00137","G00138","G00139","G00140","G00141","G00142","G00143","G00144","G00145","G00146","G00147","G00148","G00149","G00150","G00151","G00152","G00153","G00154","G00155","G00156","G00157","G00158","G00159","G00160","G00161","G00162","G00163","G00164","G00165","G00166","G00167","G00168","G00169","G00170","G00171","G00172","G00173","G00174","G00175","G00176","G00177","G00178","G00179","G00180","G00181","G00182","G00183","G00184","G00185","G00186","G00187","G00188","G00189","G00190","G00191","G00192","G00193","G00194","G00195","G00196","G00197","G00198","G00199","G00200","G00201","G00202","G00203","G00204","G00205","G00206","G00207","G00208","G00209","G00210","G00211","G00212","G00213","G00214","G00215","G00216","G00217","G00218","G00219","G00220","G00221","G00222","G00223","G00224","G00225","G00226","G00227","G00228","G00229","G00230","G00231","G00232","G00233","G00234","G00235","G00236","G00237","G00238","G00239","G00240","G00241","G00242","G00243","G00244","G00245","G00246","G00247","G00248","G00249","G00250","G00251","G00252","G00253","G00254","G00255","G00256","G00257","G00258","G00259","G00260","G00261","G00262","G00263","G00264","G00265","G00266","G00267","G00268","G00269","G00270","G00271","G00272","G00273","G00274","G00275","G00276","G00277","G00278","G00279","G00280","G00281","G00282","G00283","G00284","G00285","G00286","G00287","G00288","G00289","G00290","G00291","G00292","G00293","G00294","G00295","G00296","G00297","G00298","G00299"],"statistics":{"main":[-0.35369100406995135,-0.09370134210905152,0.2996351329719697,0.9316414831543933,0.23495758419997342,0.3852148162163202,0.36332605435434645,-0.8643802997981022,-0.16300590949862617,1.9990620116227478,0.20458944115376337,1.3561798124040774,1.463214900549598,-0.4584253564842988,2.858982587873056,2.3806547037343826,-1.406406778882225,0.6689414883315744,0.796771097507228,0.05559063743181377,-0.26979229136696303,0.09291183425079222,-0.8948405470984498,0.4382627691898982,0.3897849644648652,-0.6556750637804669,1.0422663039899254,1.2134338020297202,0.5227489630344816,-0.7625618260870751,-1.365939675567073,-1.1637633990841405,-0.752075172116772,-0.35464489731456383,-1.6111519426070073,-0.29867420398945693,0.16819556446247194,-2.5677891945764886,-1.5328902733214573,-0.008309182045553013,-2.152598955911945,1.1012017618722534,0.7571780760613831,-0.298254507880412,-1.4506862429861176,-1.672000347689788,-0.4715345336214637,0.46242883527058487,-0.22622048677754192,-0.3663121506080286,-0.40327668033404634,-0.10737864885336555,-0.19937092084129107,-2.575875610865868,1.25,-2.103757888057524,0.8402645073937265,-1.1361277492988506,-0.06374124221776316,-0.4972649369408953,-0.2108349300108254,0.6235863703099181,-0.1029825091293878,-0.6334396333780871,0.06081986764564956,1.5798647776673074,0.4946740338478745,0.379719916646331,-0.11272048765052439,0.5390412348909127,-0.8869434511917398,0.21496230037898867,-1.6401095084608632,0.12032302759780308,-0.11218678817617227,-0.02211990002135577,0.6639794497872558,0.18868274423784961,-0.6165155257573919,-0.9411251959094078,-0.040015399395260616,-0.30142632060630525,0.027160678730828974,-0.3231587334288999,0.26756033076393754,0.057932308179900055,-1.6212393004898111,0.9970578713214933,-0.525439490451896,0.40810121814604555,0.8533426674355258,0.8285854965204899,-2.1547404662137377,-0.4672778171820807,-3.9561279049063236,-0.36285659701293876,-0.5290177880719159,-0.496268222894669,-0.20723344167006913,-3.947435694301442,1.2898156026952259,-0.821185833456974,2.041903593675783,2.344559317171158,-0.11267051242102888,-2.0423089581785616,-0.4075687942896876,0.6473980505179693,0.6442084447611951,-0.025510957111514045,1.25,2.5853020659477925,-1.0365489026853167,-1.387475652464937,0.88807837484986,-0.1539565714085795,-0.8717974408858271,-4.516965390393576,1.1383789585439275,1.25,-0.9959957188010335,1.744463783410603,-2.2366371761993253,0.29175411922230804,-0.5810013610049998,-1.6223321847416063,-1.090172927061745,0.7765918025228651,1.214853736103794,-0.5153724137456338,-0.6991990569119073,-0.3209637358168869,-0.716684951029454,-0.34221558513601247,0.6303758498813837,-1.7818770915325346,-0.0785797275047747,-0.19273203406746525,-1.413884188009986,0.8360781789718202,-0.17920435028605555,2.9578242843642086,-0.8132829380231347,-0.6622848883285487,0.4546576924551794,-0.1445934169398969,-0.3863428143471749,-0.888961614120101,1.2388368485121817,-2.22971622550836,3.381155121703756,1.2192585033447296,0.6177325777117513,-0.6518558893598082,0.13472542200318482,-3.0072868192084927,2.569796207828408,-0.032002278662530026,1.58910735798479,3.866388280410967,-0.25221824446668956,2.2263364148135594,0.48726642944452936,-0.4196437280502446,-0.9302367549554909,-1.2069310802168967,-1.404164160268967,0.0,0.0,-0.1705002780145501,1.5373380345929262,-1.8479609800916446,0.2236945980047382,-1.3197489551253152,1.25,0.926132134700727,4.067052402503987,0.3317187043424216,-1.2339514739170616,-0.15173113446423198,-1.129826809423703,0.27549611256275824,0.45937701859144026,-1.137210015759904,0.8882021907998767,-0.7684305617820149,-4.26494454769387,-0.15589021446442322,-0.35684180910144525,2.0126574491454146,1.1769548964951804,0.08254991972983117,-0.027619327154997944,0.44413946675014865,1.3685168955380385,0.2484432133441479,1.2042443467561783,-0.0795742739398219,-0.7410837959480786,-1.6825783216159715,-1.410884481025365,-2.5611404069137587,-1.6397162022186866,1.6440640706684755,0.12956664882957922,2.194787608587694,0.10145334707241493,0.46057478072675534,-0.05415561499295025,0.5168296746571438,-1.0041736615609294,0.1696121391489852,0.9070788029218811,0.3516800611183043,-2.252601828943572,-0.20987731730516704,-1.10100201516091,-0.09888282756316191,-0.7422604744860672,-1.7503873545797248,2.4845153017702413,-2.2263850852595044,-0.4808797955993551,0.44519654058154884,2.1044356675938594,-0.06247122111356282,-0.08430784769088345,0.08642362072388582,-1.2198235269601996,-0.513636274361791,-0.2754755540257223,-1.4468988505039486,1.1586041432648564,1.0626703067496444,-2.1471968306104543,0.7846541144201231,-0.3801419023755764,0.014158273231564522,1.8091957324310926,-0.13576149227825862,-4.354556247929478,0.5795508292754066,3.433698334604742,-1.0366727536526206,-1.4581118551551666,0.012310285220118634,0.26639797122884873,1.2107017504918338,1.0068956803547726,-0.12214846466124739,-1.267450298839527,-1.688252422176669,-1.1290989258545734,-0.20158583219840787,1.268077383634115,0.353398614186582,1.0891782855984844,-0.30419929060217776,0.40057162133814456,1.2139046248984657,3.111155220162803,2.1834935432999982,2.9794237874096816,-0.8840010289614669,1.9463315709472044,-0.8572987267912882,0.6442001630442079,-0.3185081082252342,-1.0770726221987386,-1.238623060489314,-0.28303774119168795,3.673756114394913,1.6682699701089256,-0.1698974719945466,1.5506247307148124,3.9940027040662986,-1.6474029095703862,1.016365784974541,-2.793602601896253,0.16538990376654453,1.0777869175372632,0.2755716193321196,0.28515683970184963,0.4830064989332704,-0.5269314229441465,-0.10732028730065293,0.6843198401784071,-1.312697010871007,0.5652751621444687,-0.7978452414432053,1.3301931953102002,-0.7692102443631312,-0.02370353887002776,-1.8281269308399941,0.2791563314458031,-0.29451421983179354,0.5867299384204425,-2.581720529747996,0.48427792028313377,0.08209531388455367],"positive":[0.023259917296875533,0.26021922367666916,0.883225508095387,1.1372397403502832,0.18341946499977324,0.4259339470221224,1.1813385325480863,0.11883852026528731,0.1420470449838126,1.4909015281927622,0.14352128179857057,1.0434961089062615,0.009247610100192321,0.34458272353859276,1.5480607912320292,9.478373902229107,0.09028352187751375,0.8134008828178237,0.8798802332508289,0.0542946650627607,0.5528829336622734,0.6185662083624157,0.14872157254318108,0.44797202106864986,0.2570873938402919,0.053217304288642756,0.09665212614324972,0.14625867467671474,0.585258464248589,0.15510452972432345,0.3106351075063907,0.17377571210568074,0.10883289121337808,0.04935720685012223,0.24196857976058764,0.15463306694236842,0.5799599479278913,0.2807366067624958,0.508208086118584,0.13168078234738972,0.32245978397862857,0.9258311276635677,0.7458647310478936,0.18771073928532345,0.015131353315945752,0.4717469034340957,0.08052535535610945,0.4782390856857752,0.5289697830669731,0.4312212156790785,1.5107349830932464,0.01173416453009193,0.1254716847096276,0.8875046711685448,0.18390396522085495,0.005308164336382238,0.5073598847533144,0.06505806706149375,0.14716227770269508,0.8591011572096322,0.8945016315885206,0.06015163147671821,0.016148287086301363,0.40747956508896943,0.9028558860299929,0.059954510995156196,0.09361857135837007,0.7191673867099803,1.182386520693624,1.7463405280027173,0.682717292785221,0.33691509551568316,0.10484130330274678,1.1982954359292126,0.9956126160157436,1.0177985755275285,0.1791458914033884,0.6874515570674976,0.4728143616264503,0.8677288686154065,0.7298429771904984,0.1352764572975328,0.7278711810837406,0.4019369045955894,0.34826051132207086,0.7971118438754279,0.5644335012081221,1.0025871062641993,0.3461404097550055,0.6894621867064272,0.5477159591975556,0.3549891218091639,0.3311641740184261,0.177055206951258,0.10837627036228974,0.0781273549784751,1.8021614513942044,0.7088160093073906,0.9076658606919955,0.028203937171013924,0.17541736813940173,0.9101321095564607,2.378839787778623,0.6008281193500945,0.06302531935380726,0.433160994984071,0.005987118602239944,0.7822038821847161,1.1304841143355626,0.11628922705493876,0.12110031797048808,4.8432977293609065,0.20032869939892461,0.9204750076624899,0.2488658421491846,0.3870559432680166,0.2674187365147732,0.11136068659426142,2.3006311703965148,0.17817772173071328,0.45355204642174696,0.3503500865448041,0.061253128650704536,1.2090256338036531,0.23585361643400904,0.336707812704694,1.3518815694413535,0.17652989781800704,0.22104327397789225,0.125999037815667,0.10291065484518451,0.5533087497000562,0.06546287118788582,0.18372705423889213,0.07132506310770136,0.12891158286218674,0.6745732007539069,0.9789112823797651,1.283882849328625,0.05757936751611846,0.2577692170821223,2.3814091194900096,0.2478765863155909,1.1611367662112557,0.4673178135277266,0.25535188242182355,0.4160447940074956,0.022964881104563486,0.37697903379385955,0.4453587714953313,2.519546806220476,0.11473984840167752,0.22296909884851962,0.016302571681527097,0.5376989158286833,1.6349993402906537,4.961572888926007,0.15280505602280334,2.166026867736307,2.792668024588668,0.949777176497899,2.5219258314109365,0.41355581381203,0.04811992896758778,4.333652634467043,0.2705586774053712,0.17776928836458494,0.373141233677939,0.9906822820023178,0.6062670918022265,0.5577929107459128,0.3870829186028067,0.2950040627069547,1.7270641037123673,0.3790386266198119,0.050702535839359374,1.2650968546349917,0.3500446451255198,0.16816642491979822,0.17564185240028343,0.2586408348588401,0.11990850941582387,0.52659701449641,0.11784033803024968,1.679194992416098,0.6974963938093752,0.10755760929261239,0.042629122215988746,0.17301783315426553,0.1735466203263214,0.3973571208678083,0.3286434338205717,1.0947660057590038,0.12884594084698034,1.865581499793273,0.1443841994539627,0.850258758661687,0.413859927079418,0.009804545186388678,0.1616097583240465,0.4702859268510473,0.02418624084137464,0.2503099919302931,0.1538635961678379,0.35532686656989865,0.46772320363164754,0.9351032992819488,1.0919679820634058,0.6667667042686491,0.11274246179345251,0.28072810450558,0.5529346819468279,0.7252078170997833,0.04596695737965358,1.0303244568209422,0.3843727683177149,0.31759236172028704,0.06843000902538007,0.21183989315097232,0.8031140970280458,0.36103172771650505,0.22868303453495747,0.7458354015376584,0.979853653753938,0.1039149394468023,0.4228052624186362,0.10514651708356254,0.19480661802017418,0.08104223824150494,0.21949699329935374,0.0685494662486174,1.2942487154731412,0.008892653820538304,0.1095388768675113,0.5487867138481616,0.8326018008282453,0.5400056648627126,0.66407289079944,1.1196793668783573,0.18284669748293036,0.12097730503887992,0.1337737645581225,4.932095449961142,0.10278019457257667,0.11921509955917253,0.5585728792137234,0.2375688011344938,0.21053954323705684,1.7756767176052937e-05,0.05468321018571054,0.10054328660556738,0.31223205619892075,1.4355474442280343,0.44597462815978317,0.02526390585089547,0.7090210255330746,1.1095012447772457,0.4063306177296234,0.35057318382476976,0.3270919547092307,0.3544986102394903,0.4091148226687107,0.28877776225247037,0.20343146975271467,0.2350626154282019,1.03644911453718,1.6073632300225023,0.17498408673966237,0.8234250279986063,0.5222097536589081,0.11918464831064109,3.2223616133707567,0.17006622607095503,0.5020755297955516,0.496872248170146,1.630286589371635,0.08288627869017315,0.22227319779130397,0.31119061004554643,0.009744873322680836,0.06926237670911144,0.4794981709470576,0.6698080508652562,0.1680630883170706,0.29829241963164455,0.008527692636419065,0.13250921676218771,0.07387260250859777,0.027472794685164895,0.20714494027346042,0.14934742343969595,0.13572018151088883,1.5580478996215181,0.18623455996078966,0.12852535391017478,0.2852571519336471,0.14188400665712209,0.062147568137936336,0.044130446521723876,0.19751893881693947],"unbalanced":[0.35369100406995135,0.09370134210905152,0.2996351329719697,0.9316414831543933,0.23495758419997342,0.3852148162163202,0.36332605435434645,0.8643802997981022,0.16300590949862617,1.9990620116227478,0.20458944115376337,1.3561798124040774,1.463214900549598,0.4584253564842988,-2.858982587873056,2.3806547037343826,1.406406778882225,0.6689414883315744,0.796771097507228,0.05559063743181377,0.26979229136696303,0.09291183425079222,0.8948405470984498,0.4382627691898982,0.3897849644648652,0.6556750637804669,1.0422663039899254,1.2134338020297202,0.5227489630344816,0.7625618260870751,-1.365939675567073,1.1637633990841405,0.752075172116772,0.35464489731456383,1.6111519426070073,0.29867420398945693,0.16819556446247194,2.5677891945764886,1.5328902733214573,0.008309182045553013,2.152598955911945,1.1012017618722534,0.7571780760613831,0.298254507880412,1.4506862429861176,1.672000347689788,0.4715345336214637,0.46242883527058487,0.22622048677754192,0.3663121506080286,0.40327668033404634,0.10737864885336555,0.19937092084129107,2.575875610865868,1.25,2.103757888057524,0.8402645073937265,1.1361277492988506,0.06374124221776316,0.4972649369408953,0.2108349300108254,0.6235863703099181,-0.1029825091293878,0.6334396333780871,0.06081986764564956,1.5798647776673074,0.4946740338478745,0.379719916646331,0.11272048765052439,0.5390412348909127,0.8869434511917398,0.21496230037898867,1.6401095084608632,0.12032302759780308,0.11218678817617227,0.02211990002135577,0.6639794497872558,0.18868274423784961,0.6165155257573919,0.9411251959094078,0.040015399395260616,0.30142632060630525,0.027160678730828974,0.3231587334288999,0.26756033076393754,0.057932308179900055,1.6212393004898111,0.9970578713214933,0.525439490451896,0.40810121814604555,0.8533426674355258,0.8285854965204899,2.1547404662137377,0.4672778171820807,3.9561279049063236,0.36285659701293876,0.5290177880719159,0.496268222894669,0.20723344167006913,3.947435694301442,1.2898156026952259,0.821185833456974,2.041903593675783,2.344559317171158,0.11267051242102888,2.0423089581785616,0.4075687942896876,0.6473980505179693,0.6442084447611951,0.025510957111514045,1.25,2.5853020659477925,1.0365489026853167,1.387475652464937,0.88807837484986,0.1539565714085795,0.8717974408858271,4.516965390393576,1.1383789585439275,1.25,0.9959957188010335,1.744463783410603,2.2366371761993253,0.29175411922230804,0.5810013610049998,1.6223321847416063,1.090172927061745,0.7765918025228651,1.214853736103794,0.5153724137456338,0.6991990569119073,0.3209637358168869,0.716684951029454,0.34221558513601247,0.6303758498813837,1.7818770915325346,0.0785797275047747,0.19273203406746525,-1.413884188009986,0.8360781789718202,0.17920435028605555,2.9578242843642086,0.8132829380231347,0.6622848883285487,0.4546576924551794,0.1445934169398969,0.3863428143471749,0.888961614120101,1.2388368485121817,2.22971622550836,3.381155121703756,1.2192585033447296,0.6177325777117513,0.6518558893598082,0.13472542200318482,3.0072868192084927,2.569796207828408,0.032002278662530026,1.58910735798479,3.866388280410967,0.25221824446668956,2.2263364148135594,0.48726642944452936,0.4196437280502446,0.9302367549554909,1.2069310802168967,1.404164160268967,0.0,0.0,0.1705002780145501,1.5373380345929262,1.8479609800916446,0.2236945980047382,1.3197489551253152,1.25,0.926132134700727,4.067052402503987,0.3317187043424216,1.2339514739170616,0.15173113446423198,1.129826809423703,0.27549611256275824,0.45937701859144026,1.137210015759904,0.8882021907998767,0.7684305617820149,4.26494454769387,0.15589021446442322,0.35684180910144525,2.0126574491454146,1.1769548964951804,0.08254991972983117,0.027619327154997944,0.44413946675014865,1.3685168955380385,0.2484432133441479,1.2042443467561783,0.0795742739398219,0.7410837959480786,1.6825783216159715,1.410884481025365,2.5611404069137587,1.6397162022186866,1.6440640706684755,0.12956664882957922,2.194787608587694,0.10145334707241493,0.46057478072675534,0.05415561499295025,0.5168296746571438,1.0041736615609294,0.1696121391489852,0.9070788029218811,0.3516800611183043,2.252601828943572,0.20987731730516704,1.10100201516091,0.09888282756316191,0.7422604744860672,1.7503873545797248,2.4845153017702413,2.2263850852595044,0.4808797955993551,0.44519654058154884,2.1044356675938594,0.06247122111356282,0.08430784769088345,0.08642362072388582,1.2198235269601996,0.513636274361791,0.2754755540257223,1.4468988505039486,1.1586041432648564,1.0626703067496444,2.1471968306104543,0.7846541144201231,0.3801419023755764,-0.014158273231564522,1.8091957324310926,0.13576149227825862,4.354556247929478,0.5795508292754066,3.433698334604742,1.0366727536526206,1.4581118551551666,0.012310285220118634,0.26639797122884873,1.2107017504918338,1.0068956803547726,0.12214846466124739,1.267450298839527,1.688252422176669,1.1290989258545734,0.20158583219840787,1.268077383634115,0.353398614186582,-1.0891782855984844,0.30419929060217776,0.40057162133814456,1.2139046248984657,3.111155220162803,2.1834935432999982,2.9794237874096816,0.8840010289614669,1.9463315709472044,0.8572987267912882,0.6442001630442079,0.3185081082252342,1.0770726221987386,1.238623060489314,0.28303774119168795,3.673756114394913,1.6682699701089256,0.1698974719945466,1.5506247307148124,3.9940027040662986,1.6474029095703862,1.016365784974541,2.793602601896253,0.16538990376654453,1.0777869175372632,0.2755716193321196,0.28515683970184963,0.4830064989332704,0.5269314229441465,0.10732028730065293,-0.6843198401784071,1.312697010871007,0.5652751621444687,0.7978452414432053,1.3301931953102002,0.7692102443631312,0.02370353887002776,1.8281269308399941,0.2791563314458031,0.29451421983179354,0.5867299384204425,2.581720529747996,0.48427792028313377,0.08209531388455367]},"sets":{"main":{"set_00_up":["G00164","G00159","G00194","G00118","G00251","G00285","G00202","G00110","OUTSIDE3976","OUTSIDE5003"],"set_01_none":["G00277","G00236","G00074","G00029","G00043","G00294","G00056","G00254","G00021","G00250","G00090","G00078","G00048","G00131","G00101","G00004","G00098","G00135","G00083","G00145","G00050","G00149","G00268","G00211","OUTSIDE8290","OUTSIDE1333","OUTSIDE5502","OUTSIDE2945","OUTSIDE3298","OUTSIDE2153","OUTSIDE2692","OUTSIDE8884","OUTSIDE1129","OUTSIDE4722"],"set_02_up":["G00164","G00111","G00242","G00021","G00263","G00288","OUTSIDE6951","OUTSIDE7107"],"set_03_none":["G00158","G00155","G00249","OUTSIDE9328"],"set_04_none":["G00283","G00006","G00257","G00206","G00248","G00040","G00133","G00102","G00056","G00246","G00232","OUTSIDE1781","OUTSIDE6814","OUTSIDE1638","OUTSIDE3591","OUTSIDE2409"],"set_05_up":["G00102","G00271","G00176","G00150","G00014","G00161","G00275","G00164","G00159","G00242","G00192","G00043","G00214","G00210","G00005","G00282","G00092","G00071","G00126","G00083","OUTSIDE1654","OUTSIDE7505","OUTSIDE7239","OUTSIDE1550","OUTSIDE3093"],"set_06_none":["G00060","G00038","G00180","G00018","G00216","G00141","G00089","G00097","G00045","G00268","G00196","G00041","G00057","G00020","G00048","G00295","G00190","G00199","G00226","G00278","OUTSIDE1043","OUTSIDE4285","OUTSIDE8484","OUTSIDE691","OUTSIDE5181","OUTSIDE9357","OUTSIDE464","OUTSIDE5332","OUTSIDE5225"],"set_07_none":["G00136","G00288","G00222","G00212","G00239","OUTSIDE1402","OUTSIDE1664"],"one_gene":["G00156"],"tied_block":["G00174","G00110"],"listed_twice":["G00156","G00161","G00159","G00111","G00242","G00275","G00014","G00194","G00158","G00176","G00118","G00141","G00015","G00164","G00102","G00156","G00161","G00159","G00111","G00242","G00275","G00014","G00194","G00158","G00176","G00118","G00141","G00015","G00164","G00102"],"identical_a":["G00199","G00071","G00094","G00240","G00099","G00155","G00297","G00251","G00201","G00040","G00221","G00037","G00186","G00117","G00038"],"identical_b":["G00199","G00071","G00094","G00240","G00099","G00155","G00297","G00251","G00201","G00040","G00221","G00037","G00186","G00117","G00038"],"every_gene":["G00000","G00001","G00002","G00003","G00004","G00005","G00006","G00007","G00008","G00009","G00010","G00011","G00012","G00013","G00014","G00015","G00016","G00017","G00018","G00019","G00020","G00021","G00022","G00023","G00024","G00025","G00026","G00027","G00028","G00029","G00030","G00031","G00032","G00033","G00034","G00035","G00036","G00037","G00038","G00039","G00040","G00041","G00042","G00043","G00044","G00045","G00046","G00047","G00048","G00049","G00050","G00051","G00052","G00053","G00054","G00055","G00056","G00057","G00058","G00059","G00060","G00061","G00062","G00063","G00064","G00065","G00066","G00067","G00068","G00069","G00070","G00071","G00072","G00073","G00074","G00075","G00076","G00077","G00078","G00079","G00080","G00081","G00082","G00083","G00084","G00085","G00086","G00087","G00088","G00089","G00090","G00091","G00092","G00093","G00094","G00095","G00096","G00097","G00098","G00099","G00100","G00101","G00102","G00103","G00104","G00105","G00106","G00107","G00108","G00109","G00110","G00111","G00112","G00113","G00114","G00115","G00116","G00117","G00118","G00119","G00120","G00121","G00122","G00123","G00124","G00125","G00126","G00127","G00128","G00129","G00130","G00131","G00132","G00133","G00134","G00135","G00136","G00137","G00138","G00139","G00140","G00141","G00142","G00143","G00144","G00145","G00146","G00147","G00148","G00149","G00150","G00151","G00152","G00153","G00154","G00155","G00156","G00157","G00158","G00159","G00160","G00161","G00162","G00163","G00164","G00165","G00166","G00167","G00168","G00169","G00170","G00171","G00172","G00173","G00174","G00175","G00176","G00177","G00178","G00179","G00180","G00181","G00182","G00183","G00184","G00185","G00186","G00187","G00188","G00189","G00190","G00191","G00192","G00193","G00194","G00195","G00196","G00197","G00198","G00199","G00200","G00201","G00202","G00203","G00204","G00205","G00206","G00207","G00208","G00209","G00210","G00211","G00212","G00213","G00214","G00215","G00216","G00217","G00218","G00219","G00220","G00221","G00222","G00223","G00224","G00225","G00226","G00227","G00228","G00229","G00230","G00231","G00232","G00233","G00234","G00235","G00236","G00237","G00238","G00239","G00240","G00241","G00242","G00243","G00244","G00245","G00246","G00247","G00248","G00249","G00250","G00251","G00252","G00253","G00254","G00255","G00256","G00257","G00258","G00259","G00260","G00261","G00262","G00263","G00264","G00265","G00266","G00267","G00268","G00269","G00270","G00271","G00272","G00273","G00274","G00275","G00276","G00277","G00278","G00279","G00280","G00281","G00282","G00283","G00284","G00285","G00286","G00287","G00288","G00289","G00290","G00291","G00292","G00293","G00294","G00295","G00296","G00297","G00298","G00299"]},"null":{"spread_0":["G00117","G00155","G00103","G00055","G00251","G00065","G00194","G00148","G00031","G00233","G00175","G00056","G00042","G00108","G00096","G00283","G00163","G00095","G00081","G00020","G00010","G00187","G00051","G00058","G00075"],"spread_1":["G00240","G00262","G00214","G00105","G00199","G00274","G00011","G00269","G00232","G00026","G00212","G00139","G00032","G00266","G00284","G00222","G00089","G00188","G00002","G00084","G00253","G00115","G00285","G00225","G00245"],"spread_2":["G00186","G00141","G00122","G00102","G00045","G00170","G00290","G00178","G00118","G00243","G00022","G00091","G00218","G00063","G00088","G00046","G00106","G00033","G00035","G00246","G00052","G00179","G00206","G00064","G00039"],"spread_3":["G00176","G00278","G00149","G00189","G00272","G00038","G00173","G00228","G00183","G00112","G00147","G00101","G00198","G00134","G00028","G00093","G00050","G00000","G00043","G00160","G00137","G00145","G00217","G00085","G00167"],"spread_4":["G00275","G00111","G00221","G00009","G00276","G00012","G00287","G00151","G00057","G00277","G00184","G00142","G00132","G00061","G00209","G00047","G00258","G00255","G00295","G00195","G00077","G00239","G00001","G00019","G00168"],"spread_5":["G00094","G00297","G00161","G00264","G00203","G00244","G00100","G00128","G00180","G00248","G00114","G00289","G00130","G00152","G00129","G00207","G00024","G00213","G00123","G00004","G00140","G00154","G00021","G00208","G00237"],"tilted":["G00117","G00176","G00099","G00242","G00155","G00278","G00053","G00161","G00009","G00272","G00170","G00011","G00148","G00190","G00280","G00079","G00116","G00235","G00017","G00152","G00209","G00093","G00106","G00188","G00081","G00230","G00215","G00036","G00073","G00226"],"bottom_0":["G00281","G00020","G00160","G00048","G00060","G00010","G00137","G00169","G00036","G00187","G00145","G00204","G00068","G00051","G00217","G00227","G00299","G00058","G00085","G00080","G00082","G00075","G00167","G00062","G00030"],"bottom_1":["G00181","G00084","G00195","G00172","G00215","G00253","G00077","G00273","G00279","G00115","G00239","G00249","G00104","G00285","G00001","G00226","G00197","G00225","G00019","G00157","G00109","G00245","G00168","G00286","G00138"],"bottom_2":["G00230","G00246","G00004","G00071","G00098","G00052","G00140","G00211","G00008","G00179","G00154","G00073","G00074","G00206","G00021","G00191","G00136","G00064","G00208","G00192","G00292","G00039","G00237","G00256","G00014"]},"single":{"set_05_up":["G00102","G00271","G00176","G00150","G00014","G00161","G00275","G00164","G00159","G00242","G00192","G00043","G00214","G00210","G00005","G00282","G00092","G00071","G00126","G00083","OUTSIDE1654","OUTSIDE7505","OUTSIDE7239","OUTSIDE1550","OUTSIDE3093"]}},"cases":[{"name":"defaults","statistics":"main","sets":"main","options":{},"seed":11},{"name":"gseaParam 0","statistics":"main","sets":"main","options":{"gseaParam":0},"seed":12},{"name":"gseaParam 0.5","statistics":"main","sets":"main","options":{"gseaParam":0.5},"seed":13},{"name":"scoreType neg","statistics":"main","sets":"main","options":{"scoreType":"neg"},"seed":14},{"name":"scoreType pos, all positive","statistics":"positive","sets":"main","options":{"scoreType":"pos"},"seed":15},{"name":"sizes 15 to 500, eps 0, sampleSize 51, nPermSimple 2000","statistics":"main","sets":"main","options":{"minSize":15,"maxSize":500,"eps":0,"sampleSize":51,"nPermSimple":2000},"seed":16},{"name":"permutations, nperm 2500","statistics":"main","sets":"main","options":{"nperm":2500},"seed":17},{"name":"mostly positive statistics","statistics":"unbalanced","sets":"null","options":{},"seed":18},{"name":"one set","statistics":"main","sets":"single","options":{},"seed":19}],"figures":[{"set":["set_05_up"],"gseaParam":[1],"rank":[0,0,1,1,2,2,3,3,4,4,5,5,6,9,10,15,16,19,20,101,102,109,110,118,119,144,145,183,184,189,190,236,237,239,240,243,244,284,285,288,289,301],"ES":[0,0,0.10561147900528468,0.1056114790052847,0.20932603465444433,0.20932603465444433,0.3097267534776145,0.30972675347761447,0.40512528252215735,0.4051252825221574,0.4942900937137521,0.4942900937137521,0.5820904852035632,0.5713761994892775,0.6456170377360952,0.6277598948789523,0.6855724476228071,0.6748581619085215,0.7278814419581062,0.4385957276723918,0.4485988214240993,0.4235988214240993,0.43100365239476585,0.4024322238233372,0.4080142729898112,0.31872855870409694,0.31944576559038496,0.1837314798760993,0.1914764254976139,0.17361928264047102,0.18201093054987294,0.017725216264158555,0.04188120684235708,0.034738349699500085,0.060814302276434695,0.05010001656214888,0.0784091613195721,-0.06444798153757075,-0.008494601795214107,-0.019208887509499936,0.03928571428571426,0],"ticks":[1,2,3,4,5,6,10,16,20,102,110,119,145,184,190,237,240,244,285,289],"posES":[0.7278814419581062],"negES":[-0.06444798153757075]},{"set":["set_05_up"],"gseaParam":[0.5],"rank":[0,0,1,1,2,2,3,3,4,4,5,5,6,9,10,15,16,19,20,101,102,109,110,118,119,144,145,183,184,189,190,236,237,239,240,243,244,284,285,288,289,301],"ES":[0,0,0.08040893216938863,0.0804089321693886,0.16009246633950164,0.16009246633950167,0.23849266145160025,0.23849266145160025,0.3149148677559093,0.3149148677559094,0.3887980307355278,0.3887980307355278,0.4621137262007146,0.451399440486429,0.5188165558173587,0.5009594129602158,0.5604515823153909,0.5497372966011053,0.6067119940161583,0.3174262797304439,0.3421729004417166,0.31717290044171664,0.33846440280761075,0.309892974236182,0.32837908543590566,0.2390933711501914,0.24571966906783393,0.11000538335354824,0.13178037107777785,0.11392322822063508,0.13658909643316475,-0.027696617852549614,0.01075912084853059,0.0036162637056734517,0.04357105106553827,0.032856765351252586,0.07448731136503561,-0.06836983149210724,-0.009842072620423581,-0.02055635833470934,0.03928571428571426,0],"ticks":[1,2,3,4,5,6,10,16,20,102,110,119,145,184,190,237,240,244,285,289],"posES":[0.6067119940161583],"negES":[-0.06836983149210724]},{"set":["identical_a"],"gseaParam":[1],"rank":[0,118,119,266,267,274,275,275,276,283,284,285,286,289,290,290,291,292,293,294,295,295,296,296,297,297,298,298,299,299,300,301],"ES":[0,-0.41403508771929826,-0.4088245921082636,-0.9246140657924742,-0.8874581643750745,-0.9120195678838464,-0.8712353608645445,-0.8712353608645445,-0.8303136186773736,-0.8548750221861456,-0.8026979316343492,-0.8062067035641738,-0.7522411026436157,-0.7627674184330894,-0.7006876452545363,-0.7006876452545362,-0.638446711356758,-0.6419554832865827,-0.5793768661734732,-0.5828856381032977,-0.5099916713706093,-0.5099916713706092,-0.4143093292867339,-0.414309329286734,-0.3184162957228779,-0.3184162957228779,-0.2150378204837975,-0.21503782048379744,-0.1094872370659723,-0.10948723706597242,0,0],"ticks":[119,267,275,276,284,286,290,291,293,295,296,297,298,299,300],"posES":[0],"negES":[-0.9246140657924742]},{"set":["identical_a"],"gseaParam":[0.5],"rank":[0,118,119,266,267,274,275,275,276,283,284,285,286,289,290,290,291,292,293,294,295,295,296,296,297,297,298,298,299,299,300,301],"ES":[0,-0.41403508771929826,-0.39476032841253444,-0.9105498020967451,-0.8590787142422667,-0.8836401177510385,-0.8297144591870793,-0.8297144591870794,-0.7756979513358984,-0.8002593548446704,-0.739265037405314,-0.7427738093351384,-0.680742927523105,-0.6912692433125788,-0.6247381734804988,-0.6247381734804988,-0.5581208014217887,-0.5616295733516131,-0.49483173212655196,-0.49834050405637664,-0.4262471487706624,-0.42624714877066233,-0.34364994501750856,-0.3436499450175086,-0.2609618521974105,-0.26096185219741047,-0.17510706670557763,-0.17510706670557763,-0.08835501287496295,-0.08835501287496288,0,0],"ticks":[119,267,275,276,284,286,290,291,293,295,296,297,298,299,300],"posES":[0],"negES":[-0.9105498020967451]},{"set":["tied_block"],"gseaParam":[1],"rank":[0,38,39,40,41,301],"ES":[0,-0.12751677852348997,0.37248322147651003,0.3691275167785235,0.8691275167785235,0],"ticks":[39,41],"posES":[0.8691275167785235],"negES":[-0.12751677852348997]},{"set":["tied_block"],"gseaParam":[0.5],"rank":[0,38,39,40,41,301],"ES":[0,-0.12751677852348997,0.37248322147651003,0.3691275167785235,0.8691275167785235,0],"ticks":[39,41],"posES":[0.8691275167785235],"negES":[-0.12751677852348997]}],"multilevel":{"n1000_size101":[true,true,true,true,true,true,true,true,true,true,true,true,true,true,true,true,true,true,true,true,true,true,true,true,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false],"n2000_size51":[true,true,true,true,true,true,true,true,true,true,true,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false,false]},"fgsea":{"defaults":{"pathway":["identical_a","identical_b","listed_twice","one_gene","set_00_up","set_01_none","set_02_up","set_03_none","set_04_none","set_05_up","set_06_none","set_07_none","tied_block"],"pval":[2.6058323008722383e-10,2.6058323008722383e-10,8.670495738328725e-08,0.07045009784735812,0.4134419551934827,0.7272727272727273,0.11020408163265306,0.49317738791423,0.6975308641975309,7.796950337557118e-05,0.3603238866396761,0.7697095435684648,0.12667946257197696],"padj":[1.693790995566955e-09,1.693790995566955e-09,3.7572148199424474e-07,0.18317025440313112,0.597193935279475,0.7697095435684648,0.23526185906224292,0.641130604288499,0.7697095435684648,0.0002534008859706063,0.5855263157894737,0.7697095435684648,0.23526185906224292],"ES":[-0.9246140657924742,-0.9246140657924742,0.8753990238136659,0.9632107023411371,0.47490772993820674,-0.27555106303994625,0.6782205213750363,-0.6204988058019103,0.3425516420205607,0.7278814419581062,-0.38929215870991873,0.40129712026824776,0.8691275167785235],"NES":[-2.4070166917337503,-2.4070166917337503,2.308373709780702,1.2895522388059706,1.045250210588498,-0.8155992046661167,1.3868932081666554,-1.0161028998999986,0.8282403189427807,2.0764162899390652,-1.0916663917813307,0.765827913826721,1.3180050842119715],"size":[15,15,15,1,8,24,6,3,11,20,20,5,2],"log2err":[0.8140358378479079,0.8140358378479079,0.7049757151672383,0.2377938344236881,0.09054289484694737,0.05986031177644726,0.19189224038483796,0.07850290473486872,0.06378453903129748,0.5384340963099162,0.09821233512174378,0.05947602869573201,0.1723243449009459],"leadingEdge":[["G00117","G00240","G00186","G00094","G00099","G00155","G00297","G00037","G00201","G00221","G00040","G00251","G00199","G00038"],["G00117","G00240","G00186","G00094","G00099","G00155","G00297","G00037","G00201","G00221","G00040","G00251","G00199","G00038"],["G00176","G00275","G00159","G00242","G00141","G00014","G00111","G00156","G00015","G00161","G00102","G00158","G00194"],["G00156"],["G00159","G00194","G00110","G00118"],["G00149","G00135","G00250","G00268","G00101","G00029","G00078","G00050","G00236","G00083","G00131","G00043","G00048","G00098"],["G00242","G00111"],["G00155"],["G00102","G00232","G00248","G00056"],["G00176","G00275","G00159","G00271","G00242","G00150","G00014","G00161","G00102"],["G00278","G00199","G00045","G00038","G00057","G00180","G00216","G00268"],["G00212","G00288"],["G00110","G00174"]],"simple":[116981794],"multilevel":[277869681]},"gseaParam 0":{"pathway":["identical_a","identical_b","listed_twice","one_gene","set_00_up","set_01_none","set_02_up","set_03_none","set_04_none","set_05_up","set_06_none","set_07_none","tied_block"],"pval":[2.975356485891835e-12,2.975356485891835e-12,3.934278926186976e-10,0.07889125799573561,0.24806201550387597,0.18962075848303392,0.5613026819923371,0.8270833333333333,0.16135084427767354,0.003819968320159054,0.14314115308151093,0.5542857142857143,0.04610226723383552],"padj":[1.9339817158296926e-11,1.9339817158296926e-11,1.7048542013476894e-09,0.17093105899076047,0.32248062015503876,0.2738966511421601,0.6080779054916985,0.8270833333333333,0.2621951219512195,0.012414897040516926,0.2621951219512195,0.6080779054916985,0.11986589480797236],"ES":[-0.8631578947368421,-0.8631578947368421,0.8035087719298246,0.9632107023411371,0.339041095890411,-0.21739130434782605,0.30272108843537415,-0.31649831649831645,0.3120478137779176,0.4107142857142857,-0.2535714285714285,0.3220338983050848,0.8691275167785235],"NES":[-3.9999999999999973,-3.9999999999999973,3.8073681438072065,1.2752767527675297,1.1989981892152481,-1.2542068187044046,0.922931470313091,-0.7120086024225928,1.282812667113924,2.2386572654813,-1.3558793319891975,0.9149221636126368,1.6497432189868357],"size":[15,15,15,1,8,24,6,3,11,20,20,5,2],"log2err":[0.8986712320616942,0.8986712320616942,0.8140358378479079,0.2343926472946861,0.11934844059101546,0.14122511964253093,0.0709609456938284,0.05629407497009682,0.14920754243532175,0.4317076958033464,0.16440575583821013,0.07130529873400379,0.3217759180753611],"leadingEdge":[["G00117","G00240","G00186","G00094","G00099","G00155","G00297","G00037","G00201","G00221","G00040","G00251","G00199","G00038"],["G00117","G00240","G00186","G00094","G00099","G00155","G00297","G00037","G00201","G00221","G00040","G00251","G00199","G00038"],["G00176","G00275","G00159","G00242","G00141","G00014","G00111","G00156","G00015","G00161","G00102","G00158","G00194","G00118"],["G00156"],["G00159","G00194","G00110","G00118"],["G00149","G00135","G00250","G00268","G00101","G00029","G00078","G00050","G00236","G00083","G00131","G00043","G00048","G00098","G00145","G00074","G00021","G00211","G00004","G00294"],["G00242","G00111"],["G00155"],["G00102","G00232","G00248","G00056","G00283","G00006","G00246","G00206"],["G00176","G00275","G00159","G00271","G00242","G00150","G00014","G00161","G00102"],["G00278","G00199","G00045","G00038","G00057","G00180","G00216","G00268","G00097","G00295","G00020","G00048","G00060"],["G00212","G00288","G00136","G00239","G00222"],["G00110","G00174"]],"simple":[297914714],"multilevel":[767861829]},"gseaParam 0.5":{"pathway":["identical_a","identical_b","listed_twice","one_gene","set_00_up","set_01_none","set_02_up","set_03_none","set_04_none","set_05_up","set_06_none","set_07_none","tied_block"],"pval":[3.789796712999177e-13,3.789796712999177e-13,3.6722661636110628e-09,0.09343936381709742,0.2742616033755274,0.5436507936507936,0.1740041928721174,0.5416666666666666,0.5010309278350515,9.917864452516481e-05,0.201171875,0.7847826086956522,0.02983323454575092],"padj":[2.463367863449465e-12,2.463367863449465e-12,1.591315337564794e-08,0.20245195493704438,0.39615564932020625,0.5889550264550264,0.3231506439053609,0.5889550264550264,0.5889550264550264,0.0003223305947067856,0.326904296875,0.7847826086956522,0.0775664098189524],"ES":[-0.9105498020967451,-0.9105498020967451,0.8429442632596492,0.9632107023411371,0.41786168846742555,-0.23435079368852113,0.5082795962761097,-0.5017128860479484,0.3073313846545718,0.6067119940161583,-0.3277077976544984,0.3220338983050848,0.8691275167785235],"NES":[-3.1386770866159415,-3.1386770866159415,2.9150027987951477,1.2807937632884472,1.1681098062343673,-0.9316627567753417,1.2778673589895624,-0.9460269958677834,0.9525485897040548,2.2887497223878923,-1.2271345841791457,0.7445657241931048,1.4652233596594295],"size":[15,15,15,1,8,24,6,3,11,20,20,5,2],"log2err":[0.932595211336209,0.932595211336209,0.7749390301364363,0.20658792269689108,0.11828752612408502,0.07436253844225256,0.15214492236756472,0.0745500777420554,0.08063884959506612,0.5384340963099162,0.13500203382863138,0.060771945185353275,0.35248785758361917],"leadingEdge":[["G00117","G00240","G00186","G00094","G00099","G00155","G00297","G00037","G00201","G00221","G00040","G00251","G00199","G00038"],["G00117","G00240","G00186","G00094","G00099","G00155","G00297","G00037","G00201","G00221","G00040","G00251","G00199","G00038"],["G00176","G00275","G00159","G00242","G00141","G00014","G00111","G00156","G00015","G00161","G00102","G00158","G00194"],["G00156"],["G00159","G00194","G00110","G00118"],["G00149","G00135","G00250","G00268","G00101","G00029","G00078","G00050","G00236","G00083","G00131","G00043","G00048","G00098"],["G00242","G00111"],["G00155"],["G00102","G00232","G00248","G00056","G00283","G00006","G00246"],["G00176","G00275","G00159","G00271","G00242","G00150","G00014","G00161","G00102"],["G00278","G00199","G00045","G00038","G00057","G00180","G00216","G00268"],["G00212","G00288","G00136","G00239","G00222"],["G00110","G00174"]],"simple":[903298819],"multilevel":[319144880]},"scoreType neg":{"pathway":["identical_a","identical_b","listed_twice","one_gene","set_00_up","set_01_none","set_02_up","set_03_none","set_04_none","set_05_up","set_06_none","set_07_none","tied_block"],"pval":[1.9103823915367012e-10,1.9103823915367012e-10,1,0.951048951048951,0.6303696303696303,0.3966033966033966,0.9530469530469531,0.22677322677322678,0.6363636363636364,0.8981018981018981,0.16883116883116883,0.6083916083916084,0.7902097902097902],"padj":[1.2417485544988558e-09,1.2417485544988558e-09,1,1,1,1,1,0.737012987012987,1,1,0.7316017316017316,1,1],"ES":[-0.9246140657924742,-0.9246140657924742,0,-0.03678929765886285,-0.1816991646199555,-0.27555106303994625,-0.013605442176870763,-0.6204988058019103,-0.18220768347827013,-0.06444798153757075,-0.38929215870991873,-0.20677966101694914,-0.12751677852348997],"NES":[-3.624684744597828,-3.624684744597828,0,-0.075889284433038,-0.6390262716869867,-1.1536652542753254,-0.04449937710943639,-1.6732426603372568,-0.6750420461838605,-0.26066839065203773,-1.5745436565024984,-0.6389265754505813,-0.31036546900217993],"size":[15,15,15,1,8,24,6,3,11,20,20,5,2],"log2err":[0.826657301178302,0.826657301178302,0,0.010350458278250074,0.03494010588484777,0.05629407497009682,0.010126386883447325,0.08431443206426534,0.03449192150940737,0.015367632878669303,0.10135074379709887,0.03660822213366891,0.023508484501328895],"leadingEdge":[["G00117","G00240","G00186","G00094","G00099","G00155","G00297","G00037","G00201","G00221","G00040","G00251","G00199","G00038"],["G00117","G00240","G00186","G00094","G00099","G00155","G00297","G00037","G00201","G00221","G00040","G00251","G00199","G00038"],["G00164","G00118","G00194","G00158","G00102","G00161","G00015","G00156","G00111","G00014","G00141","G00242","G00159","G00275","G00176"],["G00156"],["G00251","G00202"],["G00149","G00135","G00250","G00268","G00101","G00029","G00078","G00050","G00236","G00083","G00131","G00043","G00048","G00098"],["G00164","G00263","G00021","G00288","G00111","G00242"],["G00155"],["G00040"],["G00214","G00092"],["G00278","G00199","G00045","G00038","G00057","G00180","G00216","G00268"],["G00222","G00239","G00136","G00288","G00212"],["G00174","G00110"]],"simple":[17343305],"multilevel":[225704875]},"scoreType pos, all positive":{"pathway":["identical_a","identical_b","listed_twice","one_gene","set_00_up","set_01_none","set_02_up","set_03_none","set_04_none","set_05_up","set_06_none","set_07_none","tied_block"],"pval":[0.932067932067932,0.932067932067932,5.7515607162759694e-09,0.0038296666652980786,0.0007647650242599387,0.5594405594405595,0.001150516686858116,0.02598388586974631,0.4325674325674326,0.00036227604108490086,0.5074925074925075,0.3106893106893107,0.7632367632367633],"padj":[0.932067932067932,0.932067932067932,7.47702893115876e-08,0.009957133329775004,0.0033139817717930675,0.7272727272727274,0.003739179232288877,0.056298419384450335,0.702922077922078,0.002354794267051856,0.7272727272727274,0.5769944341372913,0.9020070838252657],"ES":[0.3128204134275111,0.3128204134275111,0.9473684210526316,0.9966555183946488,0.9079989658004017,0.46670337221195757,0.9364088610948951,0.9285785435564251,0.5406384180008217,0.8100753309815649,0.4859029388385642,0.6395413974027904,0.30820222076696735],"NES":[0.6301744124744575,0.6301744124744575,1.90846668729961,2.0449055775142764,1.7512752846941362,0.9712488018769926,1.7782008273921588,1.7341179506021216,1.0653430008544071,1.6662094861334633,0.9994330836515356,1.2000654276626843,0.5860052897795324],"size":[15,15,15,1,8,24,6,3,11,20,20,5,2],"log2err":[0.012316745573057225,0.012316745573057225,0.7614608014455848,0.4317076958033464,0.47727081536286214,0.04049348298595027,0.45505986738722953,0.35248785758361917,0.052269331749930564,0.49849310876659025,0.04495430676400623,0.06799225915440577,0.02541180838756501],"leadingEdge":[["G00155"],["G00155"],["G00015","G00156","G00242","G00111","G00164","G00159","G00161","G00141","G00102","G00118","G00158","G00194","G00275","G00014","G00176"],["G00156"],["G00164","G00159","G00118","G00194"],["G00050","G00074","G00101","G00098","G00268","G00021","G00131","G00211","G00090","G00236","G00048","G00056","G00078","G00149"],["G00242","G00111","G00164"],["G00158","G00155"],["G00102","G00006","G00206"],["G00242","G00164","G00271","G00159","G00161","G00150","G00102","G00275","G00014","G00126","G00176","G00192","G00214"],["G00141","G00041","G00060","G00018","G00196","G00268","G00097","G00089","G00020","G00048","G00038","G00045"],["G00222","G00212","G00136"],["G00174"]],"simple":[438579695],"multilevel":[804181209]},"sizes 15 to 500, eps 0, sampleSize 51, nPermSimple 2000":{"pathway":["identical_a","identical_b","listed_twice","set_01_none","set_05_up","set_06_none"],"pval":[2.0334204201121108e-10,2.0334204201121108e-10,4.498445953979573e-08,0.7407776669990029,4.4236481738906336e-05,0.30784708249496984],"padj":[6.100261260336332e-10,6.100261260336332e-10,8.996891907959146e-08,0.7407776669990029,6.635472260835951e-05,0.3694164989939638],"ES":[-0.9246140657924742,-0.9246140657924742,0.8753990238136659,-0.27555106303994625,0.7278814419581062,-0.38929215870991873],"NES":[-2.417884273807865,-2.417884273807865,2.304809219934198,-0.8197642939371417,2.0613426431719932,-1.113575437481634],"size":[15,15,15,24,20,20],"log2err":[1.1659933427501616,1.1659933427501616,1.0148669398294836,0.04198530409973985,0.786112551315215,0.07597737366197019],"leadingEdge":[["G00117","G00240","G00186","G00094","G00099","G00155","G00297","G00037","G00201","G00221","G00040","G00251","G00199","G00038"],["G00117","G00240","G00186","G00094","G00099","G00155","G00297","G00037","G00201","G00221","G00040","G00251","G00199","G00038"],["G00176","G00275","G00159","G00242","G00141","G00014","G00111","G00156","G00015","G00161","G00102","G00158","G00194"],["G00149","G00135","G00250","G00268","G00101","G00029","G00078","G00050","G00236","G00083","G00131","G00043","G00048","G00098"],["G00176","G00275","G00159","G00271","G00242","G00150","G00014","G00161","G00102"],["G00278","G00199","G00045","G00038","G00057","G00180","G00216","G00268"]],"simple":[786447999],"multilevel":[487477163]},"permutations, nperm 2500":{"pathway":["set_00_up","set_01_none","set_02_up","set_03_none","set_04_none","set_05_up","set_06_none","set_07_none","one_gene","tied_block","listed_twice","identical_a","identical_b"],"pval":[0.4021121039805037,0.7547619047619047,0.11864406779661017,0.49715678310316813,0.7046379170056957,0.0008110300081103001,0.3033884948778566,0.7438148443735035,0.08730800323362975,0.1357202331390508,0.0008090614886731392,0.0007898894154818325,0.0007898894154818325],"padj":[0.5808285946385053,0.7547619047619047,0.2520518615439515,0.6463038180341186,0.7547619047619047,0.0026358475263584752,0.4930063041765169,0.7547619047619047,0.22700080840743736,0.2520518615439515,0.0026358475263584752,0.0026358475263584752,0.0026358475263584752],"ES":[0.47490772993820674,-0.27555106303994625,0.6782205213750363,-0.6204988058019103,0.3425516420205607,0.7278814419581062,-0.38929215870991873,0.40129712026824776,0.9632107023411371,0.8691275167785235,0.8753990238136659,-0.9246140657924742,-0.9246140657924742],"NES":[1.0485844229075767,-0.8127782326470845,1.3839040886292506,-1.0228001441032455,0.8301079263447403,2.0717834975082803,-1.1042657668091564,0.7764616259469956,1.2861230742549936,1.3002522861544077,2.309465659922996,-2.4234565218856385,-2.4234565218856385],"size":[8,24,6,3,11,20,20,5,1,2,15,15,15],"log2err":null,"leadingEdge":[["G00159","G00194","G00110","G00118"],["G00149","G00135","G00250","G00268","G00101","G00029","G00078","G00050","G00236","G00083","G00131","G00043","G00048","G00098"],["G00242","G00111"],["G00155"],["G00102","G00232","G00248","G00056"],["G00176","G00275","G00159","G00271","G00242","G00150","G00014","G00161","G00102"],["G00278","G00199","G00045","G00038","G00057","G00180","G00216","G00268"],["G00212","G00288"],["G00156"],["G00110","G00174"],["G00176","G00275","G00159","G00242","G00141","G00014","G00111","G00156","G00015","G00161","G00102","G00158","G00194"],["G00117","G00240","G00186","G00094","G00099","G00155","G00297","G00037","G00201","G00221","G00040","G00251","G00199","G00038"],["G00117","G00240","G00186","G00094","G00099","G00155","G00297","G00037","G00201","G00221","G00040","G00251","G00199","G00038"]],"simple":[665974760,937477854,678136303],"multilevel":null},"mostly positive statistics":{"pathway":["bottom_0","bottom_1","bottom_2","spread_0","spread_1","spread_2","spread_3","spread_4","spread_5","tilted"],"pval":[null,null,null,0.5548841893252769,0.5780463242698892,0.594159113796576,0.6163141993957704,0.648539778449144,0.6576032225579054,0.055165496489468405],"padj":[null,null,null,0.6576032225579054,0.6576032225579054,0.6576032225579054,0.6576032225579054,0.6576032225579054,0.6576032225579054,0.2777777777777778],"ES":[-0.8181818181818182,-0.8218181818181818,-0.8254545454545454,0.40732241391814505,0.402008246402615,0.39862251091955453,0.3931310290319043,0.3880079176637499,0.385133747991382,0.5408753102385072],"NES":[null,null,null,0.9745700539742133,0.9618552404862453,0.9537544429368824,0.9406153825346694,0.9283576948848524,0.9214808828140746,1.3115783291276588],"size":[25,25,25,25,25,25,25,25,25,30],"log2err":[null,null,null,0.041237609054006544,0.039356195952129094,0.0380842053196415,0.036378305551827785,0.03397236136349927,0.0333092817548776,0.19002330527910832],"leadingEdge":[["G00030","G00062","G00167","G00075","G00082","G00080","G00085","G00058","G00299","G00227","G00217","G00051","G00068","G00204","G00145","G00187","G00036","G00169","G00137","G00010","G00060","G00048","G00160","G00020","G00281"],["G00138","G00286","G00168","G00245","G00109","G00157","G00019","G00225","G00197","G00226","G00001","G00285","G00104","G00249","G00239","G00115","G00279","G00273","G00077","G00253","G00215","G00172","G00195","G00084","G00181"],["G00014","G00256","G00237","G00039","G00292","G00192","G00208","G00064","G00136","G00191","G00021","G00206","G00074","G00073","G00154","G00179","G00008","G00211","G00140","G00052","G00098","G00071","G00004","G00246","G00230"],["G00117","G00155","G00103","G00055","G00251","G00065","G00194","G00148","G00031","G00233"],["G00240","G00262","G00214","G00105","G00199","G00274","G00011","G00269","G00232","G00026"],["G00186","G00141","G00122","G00102","G00045","G00170","G00290","G00178","G00118","G00243"],["G00176","G00278","G00149","G00189","G00272","G00038","G00173","G00228","G00183","G00112"],["G00275","G00111","G00221","G00009","G00276","G00012","G00287","G00151","G00057","G00277"],["G00094","G00297","G00161","G00264","G00203","G00244","G00100","G00128","G00180","G00248"],["G00117","G00176","G00099","G00242","G00155","G00278","G00053","G00161","G00009","G00272"]],"simple":[313112016],"multilevel":null},"one set":{"pathway":["set_05_up"],"pval":[0.00010242157971763445],"padj":[0.00010242157971763445],"ES":[0.7278814419581062],"NES":[2.049416352699836],"size":[20],"log2err":[0.5384340963099162],"leadingEdge":[["G00176","G00275","G00159","G00271","G00242","G00150","G00014","G00161","G00102"]],"simple":[503086058],"multilevel":[293756293]}}} \ No newline at end of file diff --git a/tests/test_gsea.py b/tests/test_gsea.py new file mode 100644 index 0000000..7d56328 --- /dev/null +++ b/tests/test_gsea.py @@ -0,0 +1,629 @@ +"""Tests for gene set enrichment analysis (GSEA) and its figures. + +Reference values in ``tests/data/gsea`` come from fgsea 1.32.2, which +``validation/gsea.py`` runs on 300 genes and 14 sets, recording the seeds it draws: +given them through ``rng``, ``gsea`` must give fgsea's tables, its choice of the +multilevel method for every count, and the data of ``plotEnrichmentData``. The other +tests check enrichment scores and leading edges derived by hand, the Mersenne +Twister against the C++ standard and numpy, the scores of all the first genes of a +random sample against a walk of each, invariants that a wrong implementation would +break (negated statistics, tied genes outside a set, the order of the sets and of +the genes, seeds and threads), sizes, the forms the statistics take, warnings, +messages and the figures. +""" + +import json +import warnings +from functools import cache +from pathlib import Path +from typing import Any + +import numpy as np +import pytest +from numpy.testing import assert_allclose, assert_array_equal + +from cytepool import gsea as module +from cytepool._ranked_genes import running_sum +from cytepool.cube import Cube +from cytepool.design import Design +from cytepool.differential_expression import limma_voom +from cytepool.gsea import gsea +from cytepool.plotting import plot_enrichment, plot_gsea_table +from cytepool.stats.multiple_testing import adjust_p_values +from cytepool.table import Table + +DATA = Path(__file__).parent / "data" / "gsea" +EXACT = 1e-10 +# fgsea's arguments and gsea's. +OPTIONS = { + "gseaParam": "exponent", + "scoreType": "score_type", + "minSize": "min_size", + "maxSize": "max_size", + "eps": "eps", + "sampleSize": "sample_size", + "nPermSimple": "n_permutations", + "nperm": "n_permutations", +} +COLUMNS = ( + "gene_set", + "n_genes", + "enrichment_score", + "normalized_enrichment_score", + "p_value", + "adjusted_p_value", + "log2_error", + "leading_edge", +) +# A ranking derived by hand: weights |r| of 3, 2, 1, 1, 2, 3. +SIX = {"a": 3.0, "b": 2.0, "c": 1.0, "d": -1.0, "e": -2.0, "f": -3.0} + + +class Replay(np.random.Generator): + """A generator that draws the given seeds, fgsea's, as gsea draws its seeds.""" + + def __init__(self, seeds: list[int]) -> None: + super().__init__(np.random.PCG64()) + self.seeds = np.asarray(seeds, dtype=np.int64) + + def choice(self, a: Any, size: Any = None, *args: Any, **kwargs: Any) -> Any: + return self.seeds[:size] - 1 # gsea's seeds run from 1, numpy's from 0 + + +@cache +def fixture() -> dict[str, Any]: + data: dict[str, Any] = json.loads((DATA / "fgsea.json").read_text()) + return data + + +def case_named(name: str) -> dict[str, Any]: + return next(case for case in fixture()["cases"] if case["name"] == name) + + +def run_case(case: dict[str, Any], rng: Any, **extra: Any) -> Table: + data = fixture() + values = data["statistics"][case["statistics"]] + options = {OPTIONS[key]: value for key, value in case["options"].items()} + if "nperm" in case["options"]: + options["method"] = "permutation" + with warnings.catch_warnings(): + warnings.simplefilter("ignore", RuntimeWarning) + return gsea( + dict(zip(data["genes"], values, strict=True)), + data["sets"][case["sets"]], + rng=rng, + **options, + **extra, + ) + + +def fgsea_seeds(name: str) -> list[int]: + reference = fixture()["fgsea"][name] + seeds = list(reference["simple"]) + if "nperm" not in case_named(name)["options"]: + seeds.append((reference["multilevel"] or [1])[0]) + return seeds + + +def new_axes() -> Any: + """The Axes of a figure that pyplot does not keep; skipped without matplotlib.""" + return pytest.importorskip("matplotlib.figure").Figure().subplots() + + +def numbers(values: list[Any]) -> np.ndarray: + return np.array([np.nan if v is None else v for v in values], dtype=float) + + +def scores(statistics: dict[str, float], sets: dict[str, Any], **options: Any) -> Table: + """gsea with few permutations, for its deterministic columns.""" + options = {"method": "permutation", "n_permutations": 10, "rng": 1} | options + with warnings.catch_warnings(): + warnings.simplefilter("ignore", RuntimeWarning) + return gsea(statistics, sets, **options) + + +# fgsea, replayed. + + +@pytest.mark.parametrize( + "name", + [ + "defaults", + "gseaParam 0", + "gseaParam 0.5", + "scoreType neg", + "scoreType pos, all positive", + "sizes 15 to 500, eps 0, sampleSize 51, nPermSimple 2000", + "permutations, nperm 2500", + "mostly positive statistics", + ], +) +def test_gsea_given_fgsea_seeds_gives_fgsea_table(name: str) -> None: + reference = fixture()["fgsea"][name] + table = run_case(case_named(name), Replay(fgsea_seeds(name))) + row = {gene_set: i for i, gene_set in enumerate(table["gene_set"].tolist())} + assert sorted(row) == sorted(reference["pathway"]) + at = [row[gene_set] for gene_set in reference["pathway"]] + assert table["n_genes"][at].tolist() == reference["size"] + edges = [table["leading_edge"][i] for i in at] + assert edges == [", ".join(genes) for genes in reference["leadingEdge"]] + assert_allclose(table["enrichment_score"][at], reference["ES"], rtol=EXACT) + pairs = [("normalized_enrichment_score", "NES"), ("p_value", "pval")] + if reference["log2err"] is not None: + pairs.append(("log2_error", "log2err")) + for column, key in pairs: + assert_allclose(table[column][at], numbers(reference[key]), rtol=EXACT) + if name != "mostly positive statistics": # the deliberate difference + assert_allclose( + table["adjusted_p_value"][at], numbers(reference["padj"]), rtol=EXACT + ) + + +def test_a_single_set_has_fgsea_enrichment_score_and_leading_edge() -> None: + # fgsea draws the random sets of a single set with R's sampler. + reference = fixture()["fgsea"]["one set"] + table = run_case(case_named("one set"), 1) + assert table["gene_set"].tolist() == reference["pathway"] + assert_allclose(table["enrichment_score"], reference["ES"], rtol=EXACT) + assert table["leading_edge"][0] == ", ".join(reference["leadingEdge"][0]) + # Within Monte Carlo error: 4 standard deviations of the difference of the log2 + # p-values, as each program estimates its own; and of the NES, whose means of + # about 500 random scores each vary by about 1.5%. + errors = np.hypot(table["log2_error"][0], reference["log2err"][0]) + difference = abs(np.log2(table["p_value"][0] / reference["pval"][0])) + assert difference < 4.0 * errors + assert_allclose(table["normalized_enrichment_score"], reference["NES"], rtol=0.1) + + +def test_adjusted_p_values_are_over_the_sets_with_a_p_value() -> None: + # fgsea adjusts over the sets before it removes those with too few random scores + # of their sign; CytePool adjusts the p-values it reports. + name = "mostly positive statistics" + table = run_case(case_named(name), Replay(fgsea_seeds(name))) + p = table["p_value"] + assert np.isnan(p).sum() == 3 + assert_array_equal(table["adjusted_p_value"], adjust_p_values(p, "bh")) + assert np.all(np.isnan(table["normalized_enrichment_score"][np.isnan(p)])) + assert np.all(np.isnan(table["log2_error"][np.isnan(p)])) + + +def test_the_choice_of_the_multilevel_method_is_fgsea() -> None: + choices = fixture()["multilevel"] + for key, n, size in (("n1000_size101", 1000, 101), ("n2000_size51", 2000, 51)): + counts = np.arange(n + 1, dtype=np.float64) + assert module._needs_multilevel(counts, n, size).tolist() == choices[key] + + +# Enrichment scores and leading edges derived by hand. + + +def test_enrichment_score_and_leading_edge_of_a_set_at_the_top() -> None: + # Weights 3 and 1 of 4, at ranks 1 and 3 of 6: the walk reaches 3/4 at a and + # 1 - 1/4 at c, and never falls below 0 before them. + weighted = scores(SIX, {"ac": ["a", "c"]}) + assert weighted["enrichment_score"][0] == 0.75 + assert weighted["leading_edge"][0] == "a" + # Unweighted: 1/2 at a and 1 - 1/4 at c, the maximum. + unweighted = scores(SIX, {"ac": ["a", "c"]}, exponent=0.0) + assert unweighted["enrichment_score"][0] == 0.75 + assert unweighted["leading_edge"][0] == "a, c" + + +def test_enrichment_score_and_leading_edge_of_a_set_at_the_bottom() -> None: + # Weights 1 and 3 at ranks 4 and 6: just before d the walk is -3/4, just before f + # 1/4 - 4/4 = -3/4 too; the first minimum, at d, starts the leading edge, which + # lists the genes from the bottom. + table = scores(SIX, {"df": ["d", "f"]}) + assert table["enrichment_score"][0] == -0.75 + assert table["leading_edge"][0] == "f, d" + + +def test_equal_deviations_give_no_score_unless_one_sided() -> None: + # Weights 3 and 3 at ranks 1 and 6: the walk reaches 1/2 and falls to -1/2. + sets = {"af": ["a", "f"]} + both = scores(SIX, sets) + assert both["enrichment_score"][0] == 0.0 and both["leading_edge"][0] == "" + up = scores(SIX, sets, score_type="pos") + assert up["enrichment_score"][0] == 0.5 and up["leading_edge"][0] == "a" + down = scores(SIX, sets, score_type="neg") + assert down["enrichment_score"][0] == -0.5 and down["leading_edge"][0] == "f" + + +def test_tied_genes_keep_the_order_given() -> None: + # b ties with a and c: second as given, it is passed after one miss. + statistics = {"a": 1.0, "b": 1.0, "c": 1.0, "d": 0.5, "e": -1.0} + assert scores(statistics, {"b": ["b"]})["enrichment_score"][0] == 0.75 + first = {"b": 1.0, "a": 1.0, "c": 1.0, "d": 0.5, "e": -1.0} + assert scores(first, {"b": ["b"]})["enrichment_score"][0] == 1.0 + + +# Invariants. + + +@cache +def ranking() -> tuple[dict[str, float], dict[str, list[str]]]: + """500 genes with distinct statistics, and sets up, down and spread.""" + rng = np.random.default_rng(3) + genes = [f"g{i}" for i in range(500)] + values = rng.normal(size=500) + values[:40] += 0.8 + values[40:70] -= 0.8 + sets = { + "up": genes[:30] + genes[200:210], + "down": genes[40:65], + "spread": genes[::17], + "small": genes[100:103], + } + return dict(zip(genes, values.tolist(), strict=True)), sets + + +def test_negated_statistics_negate_the_scores() -> None: + statistics, sets = ranking() + table = scores(statistics, sets) + negated = scores({gene: -value for gene, value in statistics.items()}, sets) + assert_allclose(negated["enrichment_score"], -table["enrichment_score"], rtol=1e-12) + # The leading edge holds the same genes, from the most extreme statistic. + assert negated["leading_edge"].tolist() == table["leading_edge"].tolist() + + +def test_tied_genes_outside_a_set_leave_its_results_unchanged() -> None: + statistics, sets = ranking() + tied = dict(statistics) + block = [f"g{i}" for i in range(300, 340)] + for gene in block: + tied[gene] = 0.25 + reordered = {gene: tied[gene] for gene in reversed(block)} + reordered |= {gene: value for gene, value in tied.items() if gene not in reordered} + outside = { + name: genes for name, genes in sets.items() if not set(genes) & set(block) + } + with pytest.warns(RuntimeWarning, match="39 of the 500 statistics"): + first = gsea(tied, outside, rng=2, n_permutations=200) + with pytest.warns(RuntimeWarning, match="tie with another one"): + second = gsea(reordered, outside, rng=2, n_permutations=200) + assert list(outside) == ["up", "down", "small"] + assert first.equals(second) + + +def test_the_order_of_genes_with_distinct_statistics_does_not_matter() -> None: + statistics, sets = ranking() + shuffled = dict(sorted(statistics.items(), key=lambda item: item[0][::-1])) + assert gsea(statistics, sets, rng=4).equals(gsea(shuffled, sets, rng=4)) + + +def test_each_set_has_the_same_results_in_any_order_of_the_sets() -> None: + # The random sets of each size, and the multilevel method of each size, do not + # depend on the other sets. + statistics, sets = ranking() + forward = gsea(statistics, sets, rng=5) + backward = gsea(statistics, dict(reversed(list(sets.items()))), rng=5) + order = [backward["gene_set"].tolist().index(name) for name in sets] + assert backward.take(order).equals(forward) + + +def test_a_seed_reproduces_results_and_threads_do_not_change_them() -> None: + statistics, sets = ranking() + first = gsea(statistics, sets, rng=7, n_threads=1) + assert first.equals(gsea(statistics, sets, rng=np.random.default_rng(7))) + assert first.equals(gsea(statistics, sets, rng=7, n_threads=2)) + other = gsea(statistics, sets, rng=8) + assert other["enrichment_score"].tolist() == first["enrichment_score"].tolist() + assert other["p_value"].tolist() != first["p_value"].tolist() + # Permutations in three chunks of their own seeds. + chunks = gsea(statistics, sets, rng=9, method="permutation", n_permutations=2500) + assert chunks.equals( + gsea( + statistics, + sets, + rng=9, + method="permutation", + n_permutations=2500, + n_threads=2, + ) + ) + + +def test_sets_are_sized_and_matched_by_the_ranked_genes() -> None: + statistics, sets = ranking() + genes = list(statistics) + padded = { + "up": [*sets["up"], "NOT_RANKED_1", *sets["up"][:5], "NOT_RANKED_2"], + "every": genes, + "one": [genes[0]], + } + table = gsea(statistics, padded, rng=1, n_permutations=100) + # A set of every ranked gene is never tested; genes listed twice count once. + assert table["gene_set"].tolist() == ["up", "one"] + assert table["n_genes"].tolist() == [40, 1] + plain = gsea( + statistics, {"up": sets["up"], "one": [genes[0]]}, rng=1, n_permutations=100 + ) + assert table.equals(plain) + assert gsea(statistics, padded, rng=1, n_permutations=100, max_size=10**6).equals( + plain + ) + limited = gsea( + statistics, padded, rng=1, n_permutations=100, min_size=2, max_size=40 + ) + assert limited["gene_set"].tolist() == ["up"] + empty = gsea(statistics, padded, rng=1, min_size=41, max_size=50) + assert empty.columns == COLUMNS and len(empty["gene_set"]) == 0 + + +# Random numbers and the compiled parts. + + +def test_mersenne_twister_is_cxx_std_mt19937() -> None: + # The C++ standard requires the 10000th output of a default std::mt19937, seeded + # with 5489, to be 4123659995. + state = module._mt_seeded(5489) + outputs = [module._mt_next(state) for _ in range(10_000)] + assert outputs[-1] == 4123659995 + # numpy's Mersenne Twister, seeded as the reference generator seeds it. + for seed in (1, 123_456_789, 10**9): + state = module._mt_seeded(seed) + reference = np.random.MT19937(0) + reference._legacy_seeding(seed) + expected = reference.random_raw(2000).astype(np.int64).tolist() + assert [module._mt_next(state) for _ in range(2000)] == expected + + +def test_scores_of_all_first_genes_are_the_scores_of_each() -> None: + # The convex hulls give, for each number of the first genes drawn, the score of + # those genes, as a walk over them gives it. + rng = np.random.default_rng(11) + weights = np.ascontiguousarray(np.sort(np.abs(rng.normal(size=400)))[::-1]) + weights[150:170] = 0.0 # genes of no weight + selected = rng.choice(400, 120, replace=False).astype(np.int64) + 1 + # fgsea raises the weights of the sample to a floor below its smallest nonzero one. + drawn = weights[selected - 1] + floored = np.maximum(weights, min(1e-5, drawn[drawn > 0].min()) / 1024) + for score_type in (0, 1, 2): + prefix = np.empty(selected.size) + module._prefix_scores(weights, selected, score_type, prefix) + expected = [] + for i in range(1, selected.size + 1): + positions = np.sort(selected[:i] - 1) + tops, bottoms = running_sum(floored, positions) + up, down = tops.max(), bottoms.min() + expected.append( + [0.0 if up == -down else up if up > -down else down, up, down][ + score_type + ] + ) + assert_allclose(prefix, expected, rtol=1e-12, atol=1e-15) + + +def test_random_samples_hold_distinct_genes_in_range() -> None: + state = module._mt_seeded(3) + used = np.zeros(50, dtype=np.bool_) + for k in (1, 10, 24, 25, 40, 49): # below half the genes, and at least half + out = np.empty(k, dtype=np.int64) + module._random_sample(state, 1, 50, k, used, out) + assert len(set(out.tolist())) == k and 1 <= out.min() and out.max() <= 50 + assert not used.any() + + +# The forms of the statistics. + + +@cache +def de_result() -> Any: + """limma-voom on simulated counts of 200 genes in 8 samples and two groups.""" + rng = np.random.default_rng(5) + mean = np.exp(rng.uniform(2.0, 7.0, 200)) + counts = rng.poisson(np.tile(mean, (2, 8, 1))) + counts[:, 4:, :30] = rng.poisson(4.0 * np.tile(mean[:30], (2, 4, 1))) + counts[:, 4:, 30:50] = rng.poisson(0.25 * np.tile(mean[30:50], (2, 4, 1))) + samples = { + "sample_id": [f"S{i}" for i in range(8)], + "condition": ["ctrl"] * 4 + ["stim"] * 4, + } + cube = Cube( + {"counts": counts.astype(np.uint32)}, + samples=samples, + groups={"group_id": ["B", "T"]}, + features={"feature_id": [f"g{i}" for i in range(200)]}, + ) + design = Design(samples, "~ condition", contrasts={"stim": "conditionstim"}) + return limma_voom(cube, design) + + +def test_a_result_its_table_and_a_mapping_rank_genes_alike() -> None: + result = de_result() + sets = { + "up": [f"g{i}" for i in range(0, 30)], + "down": [f"g{i}" for i in range(30, 50)], + } + quick = {"method": "permutation", "n_permutations": 100, "rng": 1} + tested = gsea(result, sets, group="T", contrast="stim", **quick) + table = result.table("T", "stim") + assert tested.equals(gsea(table, sets, **quick)) + mapping = dict(zip(table["feature_id"].tolist(), table["t"].tolist(), strict=True)) + assert tested.equals(gsea(mapping, sets, **quick)) + assert tested["enrichment_score"][0] > 0.5 > -0.5 > tested["enrichment_score"][1] + by_change = gsea(table, sets, statistic_column="log2_fold_change", **quick) + assert by_change["enrichment_score"][0] > 0.5 + with pytest.raises( + ValueError, match="the result has 2 groups; name one with group" + ): + gsea(result, sets) + + +def test_a_pandas_series_ranks_by_its_index() -> None: + pandas = pytest.importorskip("pandas") + statistics, sets = ranking() + series = pandas.Series(statistics, name="t") + assert gsea(series, sets, rng=1).equals(gsea(statistics, sets, rng=1)) + frame = pandas.DataFrame({"feature_id": list(statistics), "t": series.to_numpy()}) + assert gsea(frame, sets, rng=1).equals(gsea(statistics, sets, rng=1)) + + +def test_genes_without_an_id_or_a_statistic_are_left_out() -> None: + statistics, sets = ranking() + table = Table( + { + "feature_id": [*statistics, None, "", "extra"], + "t": [*statistics.values(), 5.0, 6.0, np.nan], + } + ) + assert gsea(table, sets, rng=1).equals(gsea(statistics, sets, rng=1)) + + +# Warnings. + + +@pytest.mark.filterwarnings("ignore:the p-values of .* likely too large:RuntimeWarning") +@pytest.mark.filterwarnings("ignore:3 of the 300 statistics:RuntimeWarning") +def test_warnings_about_ties_positive_statistics_and_p_values() -> None: + statistics, sets = ranking() + tied = dict(statistics) | {"g300": 0.5, "g301": 0.5, "g302": 0.0, "g303": 0.0} + with pytest.warns(RuntimeWarning, match=r"1 of the 500 statistics \(0.20%\) tie"): + gsea(tied, sets, rng=1, n_permutations=50, method="permutation") + positive = {gene: abs(value) + 0.01 for gene, value in statistics.items()} + with pytest.warns(RuntimeWarning, match="all statistics are positive"): + gsea(positive, sets, rng=1, n_permutations=50, method="permutation") + data, case = fixture(), case_named("mostly positive statistics") + unbalanced = dict( + zip(data["genes"], data["statistics"][case["statistics"]], strict=True) + ) + with pytest.warns(RuntimeWarning, match="3 gene set.s. have fewer than 10 random"): + gsea(unbalanced, data["sets"][case["sets"]], rng=1) + with pytest.warns(RuntimeWarning, match="are below eps=0.001 and are reported"): + table = gsea(statistics, sets, rng=1, eps=1e-3) + floored = table["p_value"] == 1e-3 + assert floored.any() and np.all(np.isnan(table["log2_error"][floored])) + + +# Messages. + + +@pytest.mark.parametrize( + ("options", "error", "match"), + [ + ({"min_size": 0}, ValueError, "min_size must be at least 1"), + ({"min_size": 1.5}, TypeError, "min_size must be an integer"), + ({"min_size": 5, "max_size": 4}, ValueError, r"max_size \(4\) is smaller"), + ({"exponent": -1.0}, ValueError, "exponent must be finite and at least 0"), + ({"score_type": "up"}, ValueError, "score_type must be one of 'std'"), + ({"method": "simple"}, ValueError, "method must be one of 'multilevel'"), + ({"n_permutations": 0}, ValueError, "n_permutations must be at least 1"), + ({"sample_size": 2}, ValueError, "sample_size must be at least 3"), + ({"eps": 1.0}, ValueError, r"eps must be in \[0, 1\)"), + ({"rng": -1}, ValueError, "rng must be a non-negative seed"), + ({"rng": "seed"}, TypeError, "rng must be a numpy.random.Generator"), + ({"group": "T"}, ValueError, "group and contrast choose a table"), + ({"n_threads": 0}, ValueError, "n_threads must be at least 1"), + ], +) +def test_options_are_checked(options: dict[str, Any], error: type, match: str) -> None: + statistics, sets = ranking() + with pytest.raises(error, match=match): + gsea(statistics, sets, **options) + + +@pytest.mark.filterwarnings("ignore:all statistics are positive:RuntimeWarning") +@pytest.mark.parametrize( + ("statistics", "sets", "error", "match"), + [ + ( + {"a": 1.0, "b": np.inf}, + {"s": ["a"]}, + ValueError, + "the statistic of 'b' is inf", + ), + ({"a": 1.0}, {"s": ["a"]}, ValueError, "a ranking needs at least two"), + ({"a": 0.0, "b": -0.0}, {"s": ["a"]}, ValueError, "every statistic is zero"), + ({"a": 1.0, "b": 2.0}, {"s": ["x"]}, ValueError, "no gene of any set is among"), + ({"a": 1.0, "b": 2.0}, {}, ValueError, "gene_sets is empty"), + ({"a": 1.0, "b": 2.0}, ["a"], TypeError, "gene_sets must be a GeneSets"), + ({1: 1.0, "b": 2.0}, {"s": ["b"]}, TypeError, "gene ids must be strings"), + ({"a": "x", "b": 2.0}, {"s": ["b"]}, TypeError, "statistics must map gene ids"), + ([1.0, 2.0], {"s": ["b"]}, TypeError, "statistics must be a DifferentialExpr"), + ], +) +def test_statistics_and_sets_are_checked( + statistics: Any, sets: Any, error: type, match: str +) -> None: + with pytest.raises(error, match=match): + gsea(statistics, sets) + + +def test_tables_need_their_columns_and_one_row_per_gene() -> None: + table = Table({"feature_id": ["a", "b", "a"], "t": [1.0, 2.0, 3.0]}) + with pytest.raises(ValueError, match="the gene id 'a' is on more than one row"): + gsea(table, {"s": ["a"]}) + with pytest.raises(ValueError, match="the table has no column 'stat'"): + gsea(table, {"s": ["a"]}, statistic_column="stat") + numeric = Table({"feature_id": [1, 2], "t": [1.0, 2.0]}) + with pytest.raises(TypeError, match="must hold gene ids as text"): + gsea(numeric, {"s": ["1"]}) + + +# The figures. + + +def test_plot_enrichment_draws_fgsea_enrichment_data() -> None: + data = fixture() + statistics = dict(zip(data["genes"], data["statistics"]["main"], strict=True)) + for figure in data["figures"]: + ax = plot_enrichment( + statistics, + data["sets"]["main"][figure["set"][0]], + exponent=figure["gseaParam"][0], + ax=new_axes(), + ) + walk = next(line for line in ax.lines if line.get_gid() == "running sum") + assert_array_equal(walk.get_xdata(), figure["rank"]) + assert_allclose(walk.get_ydata(), figure["ES"], rtol=EXACT, atol=1e-15) + ticks = next(c for c in ax.collections if c.get_gid() == "set") + assert [s[0, 0] for s in ticks.get_segments()] == figure["ticks"] + extremes = [ln.get_ydata()[0] for ln in ax.lines if ln.get_gid() == "extremes"] + assert_allclose(extremes, figure["posES"] + figure["negES"], rtol=EXACT) + + +def test_plot_enrichment_refuses_sets_without_ranked_genes() -> None: + with pytest.raises(ValueError, match="no gene of the set is among the ranked"): + plot_enrichment(SIX, ["x"]) + with pytest.raises(ValueError, match="holds every ranked gene"): + plot_enrichment(SIX, list(SIX)) + with pytest.raises(TypeError, match="gene_set must be a list of gene ids"): + plot_enrichment(SIX, "a") + + +def test_plot_gsea_table_shows_the_top_sets_of_each_sign_with_their_numbers() -> None: + statistics, sets = ranking() + table = Table( + { + "gene_set": ["up", "down", "spread", "small"], + "enrichment_score": [0.6, -0.5, 0.2, -0.3], + "normalized_enrichment_score": [2.345, -1.5, 0.71, -1.0], + "p_value": [1.23e-5, 0.0004, 0.5, 0.02], + "adjusted_p_value": [4.9e-5, 0.0008, 0.5, 0.027], + } + ) + ax = plot_gsea_table(statistics, sets, table, top=1, ax=new_axes()) + # The most significant set up, then the most significant down, last. + assert [label.get_text() for label in ax.get_yticklabels()] == ["up", "down"] + numbers = [text.get_text() for text in ax.texts if text.get_gid() == "numbers"] + assert numbers == [ + "2.35", + "-1.50", + r"$1.2 \cdot 10^{-5}$", + r"$4.0 \cdot 10^{-4}$", + r"$4.9 \cdot 10^{-5}$", + r"$8.0 \cdot 10^{-4}$", + ] + bars = next(c for c in ax.collections if c.get_gid() == "gene ranks") + ranked = sorted(statistics, key=lambda gene: -statistics[gene]) + top_row = sorted(s[0, 0] for s in bars.get_segments() if s[0, 1] == 1.0) + assert top_row == sorted(ranked.index(gene) + 1.0 for gene in sets["up"]) + ax = plot_gsea_table(statistics, sets, table, ax=new_axes()) + names = [label.get_text() for label in ax.get_yticklabels()] + assert names == ["up", "spread", "small", "down"] + with pytest.raises(ValueError, match="gene_sets has no set 'small'"): + plot_gsea_table( + statistics, {k: v for k, v in sets.items() if k != "small"}, table + ) + with pytest.raises(ValueError, match="the table lacks the column"): + plot_gsea_table(statistics, sets, table.select(["gene_set", "p_value"])) diff --git a/tests/test_masked_inputs.py b/tests/test_masked_inputs.py index 2d11ba6..8bab538 100644 --- a/tests/test_masked_inputs.py +++ b/tests/test_masked_inputs.py @@ -50,6 +50,7 @@ from cytepool.design import Design, make_contrasts from cytepool.differential_expression import limma_voom from cytepool.gene_sets import GeneSets +from cytepool.gsea import gsea from cytepool.io.anndata import to_anndata, to_mudata from cytepool.io.zarr import open_array, write_array, write_group from cytepool.meta_analysis import differential_expression_meta_analysis @@ -135,6 +136,10 @@ def message(subject: str, found: str = "array", missing: str | None = None) -> s } LEVELS = {"condition": ["ctrl", "stim"]} COLUMNS = ("(Intercept)", "conditionstim") +# Statistics of 8 genes by gene id, for gene set enrichment analysis. +RANKED = dict( + zip([f"g{j}" for j in range(8)], [0.5, -1.0, 2.0, 0.1, 3.0, -0.2, 1.1, 0.4]) +) # Cells for from_cells: 4 cells x 3 features, two samples and two groups. CELL_COUNTS = np.array([[1, 0, 2], [3, 1, 0], [0, 2, 5], [4, 0, 1]]) CELL_SAMPLES = ["S1", "S1", "S2", "S2"] @@ -1433,6 +1438,55 @@ def _forest_arguments(context: Context) -> dict[str, Any]: options=("group", "contrast"), objects=("genes", "gene_sets"), ), + # Statistics by gene id, ranked; the options that read a table are checked too. + "cytepool.gsea.gsea": Entry( + lambda c: { + "statistics": RANKED, + "gene_sets": {"set": ["g1", "g2", "g5"]}, + "id_column": "feature_id", + "statistic_column": "t", + "min_size": 1, + "max_size": 5, + "exponent": 1.0, + "score_type": "std", + "method": "permutation", + "n_permutations": 20, + "sample_size": 11, + "eps": 1e-10, + "rng": 1, + "n_threads": 1, + }, + {"statistics": "NaN", "gene_sets": None}, + options=( + "id_column", + "statistic_column", + "min_size", + "max_size", + "exponent", + "score_type", + "method", + "n_permutations", + "sample_size", + "eps", + "rng", + "n_threads", + ), + keys={"statistics": "g0", "gene_sets": "set"}, + ), + # A differential-expression result, of which a group and a contrast are ranked. + "cytepool.gsea.gsea[result]": Entry( + lambda c: { + "statistics": _tested(c), + "gene_sets": {"set": ["feature_1", "feature_2"]}, + "group": "bulk", + "contrast": "stim", + "method": "permutation", + "n_permutations": 200, + }, + {}, + options=("group", "contrast"), + objects=("statistics", "gene_sets"), + ), # n and axis go through operator.index, which reads a 0-d masked integer. "cytepool.stats.multiple_testing.adjust_p_values": Entry( lambda c: {"p": [0.01, 0.2, 0.5], "n": 5, "axis": 0}, @@ -1802,6 +1856,71 @@ def _forest_arguments(context: Context) -> dict[str, Any]: options=("top", "alpha"), objects=("table", "ax"), ), + "cytepool.plotting.plot_enrichment": Entry( + lambda c: { + "statistics": RANKED, + "gene_set": ["g1", "g2"], + "id_column": "feature_id", + "statistic_column": "t", + "exponent": 1.0, + "ax": _axes(), + }, + {"statistics": "NaN", "gene_set": None}, + options=("id_column", "statistic_column", "exponent"), + objects=("ax",), + keys={"statistics": "g0"}, + ), + "cytepool.plotting.plot_enrichment[result]": Entry( + lambda c: { + "statistics": _tested(c), + "gene_set": ["feature_1", "feature_2"], + "group": "bulk", + "contrast": "stim", + }, + {}, + options=("group", "contrast"), + objects=("statistics", "gene_set"), + ), + "cytepool.plotting.plot_gsea_table": Entry( + lambda c: { + "statistics": RANKED, + "gene_sets": {"set": ["g1", "g2", "g5"], "other": ["g4"]}, + "table": gsea( + RANKED, + {"set": ["g1", "g2", "g5"], "other": ["g4"]}, + method="permutation", + n_permutations=20, + rng=1, + ), + "id_column": "feature_id", + "statistic_column": "t", + "top": 2, + "exponent": 0.5, + "ax": _axes(), + }, + {"statistics": "NaN", "gene_sets": None}, + options=("id_column", "statistic_column", "top", "exponent"), + objects=("table", "ax"), + keys={"statistics": "g0", "gene_sets": "set"}, + ), + "cytepool.plotting.plot_gsea_table[result]": Entry( + lambda c: { + "statistics": _tested(c), + "gene_sets": {"set": ["feature_1", "feature_2"]}, + "table": gsea( + _tested(c), + {"set": ["feature_1", "feature_2"]}, + method="permutation", + n_permutations=20, + rng=1, + ), + "group": "bulk", + "contrast": "stim", + }, + {}, + options=("group", "contrast"), + objects=("statistics", "gene_sets"), + ), "cytepool.plotting.plot_dispersions": Entry( lambda c: { "dispersions": deseq2_dispersions( diff --git a/validation/gsea.py b/validation/gsea.py new file mode 100644 index 0000000..8788de1 --- /dev/null +++ b/validation/gsea.py @@ -0,0 +1,476 @@ +"""Check gene set enrichment analysis (GSEA) and its figures against fgsea. + +Inputs, simulated from fixed seeds: statistics of 4,000 genes like moderated +t-statistics (6% up, 5% down, a block of 60 tied at 1.25 and 30 at zero), and 90 gene +sets of 1 to 1,000 genes, a fifth of each outside the ranking, enriched for genes up, +down or both or for none, with a set of every gene (which no method tests), a set of +one gene, half of the tied block, a set that lists genes twice and two identical sets. +Each case runs ``fgsea`` (``fgseaMultilevel``, or ``fgseaSimple`` given ``nperm``) and +``gsea`` with the same options: fgsea's defaults; gseaParam 0 and 0.5; scoreType +"neg"; scoreType "pos" on statistics that are all positive; minSize 15, maxSize 500, +eps 0, sampleSize 51 and nPermSimple 2,000; permutations (nperm 2,500, in three +chunks); statistics mostly positive, so that some sets have too few random scores of +their sign and none needs the multilevel method; and one set alone. + +- Exact, given fgsea's seeds: R records the seeds that fgsea draws with + ``sample.int`` where it hands them to its C++ code (``trace`` of + ``fgseaSimpleImpl`` and ``multilevelImpl``), and ``gsea`` draws the same through + ``rng``. The sets tested, their sizes and leading edges must be identical, and the + enrichment scores, normalized scores, p-values, adjusted p-values and log2 errors + must agree to 1e-10 (relative), NaN where fgsea gives NA. fgsea tests one set alone + with R's sampler, which cannot be replayed: that case is exact for the enrichment + score and leading edge only. +- Calibrated, with seeds of their own (R's 1 to 20 and numpy's 1 to 20): the mean + log10 p-value and the mean NES of each set over 20 runs must agree within 4 + standard errors of their difference, for 36 sets of 2,000 genes and for one set + alone. +- Figures: what ``plot_enrichment`` draws (the walk, the ticks, the extremes) against + ``plotEnrichmentData``, for three sets and gseaParam 1 and 0.5; and the bars and + numbers of ``plot_gsea_table`` (top 5, gseaParam 0.5) against the ticks of + ``plotEnrichmentData`` scaled as ``plotGseaTable`` scales them, and the numbers as + ``plotGseaTable`` writes them. + +Deliberate difference: in the mostly positive case, fgsea adjusts the p-values over +the sets before it leaves those with too few random scores of their sign without a +p-value; CytePool adjusts over the sets with a p-value. + +With ``--write-fixtures`` the exact cases run on 300 genes and 14 sets, without the +calibrated runs, and the tests replay them. + +Requirements: + Python 3.13 with numpy, numba and matplotlib, run from the repository root so that + ``cytepool`` imports from the working tree. + R 4.4 or newer with fgsea, data.table, BiocParallel and jsonlite (``RSCRIPT``, by + default ``Rscript``). + +Usage: + PYTHONPATH=. python validation/gsea.py + PYTHONPATH=. python validation/gsea.py --write-fixtures +""" + +import argparse +import json +import os +import shutil +import subprocess +import sys +import tempfile +import warnings +from pathlib import Path + +import matplotlib.figure +import numpy as np + +from cytepool import plotting +from cytepool.gsea import gsea + +FIXTURES = Path(__file__).resolve().parents[1] / "tests/data/gsea" +EXACT = 1e-10 +CALIBRATED = 4.0 # standard errors of the difference of two means of 20 runs +RUNS = 20 +BAR = 0.45 # the half-height of a row of plot_gsea_table +# fgsea's arguments and gsea's. +OPTIONS = { + "gseaParam": "exponent", "scoreType": "score_type", "minSize": "min_size", + "maxSize": "max_size", "eps": "eps", "sampleSize": "sample_size", + "nPermSimple": "n_permutations", "nperm": "n_permutations", +} # fmt: skip + +R_SCRIPT = r""" +suppressPackageStartupMessages({ + library(fgsea); library(data.table); library(jsonlite); library(BiocParallel) +}) +register(SerialParam()) +root <- commandArgs(trailingOnly = TRUE)[[1]] +input <- fromJSON(file.path(root, "input.json"), simplifyVector = FALSE) +# The seeds fgsea draws, recorded where it hands them to its C++ code. +seen <- new.env() +invisible(trace("fgseaSimpleImpl", where = asNamespace("fgsea"), print = FALSE, + tracer = quote(assign("simple", seeds, envir = seen)))) +invisible(trace("multilevelImpl", where = asNamespace("fgsea"), print = FALSE, + tracer = quote(assign("multilevel", seed, envir = seen)))) +statistics <- function(name) setNames(unlist(input$statistics[[name]]), + unlist(input$genes)) +sets <- function(name) lapply(input$sets[[name]], unlist) +run <- function(case, seed) { + set.seed(seed) + suppressWarnings(do.call(fgsea, c(list(pathways = sets(case$sets), + stats = statistics(case$statistics)), case$options))) +} +out <- list(cases = list(), calibrated = list(), figures = list(), bars = list()) +for (case in input$cases) { + rm(list = ls(seen), envir = seen) + res <- run(case, case$seed) + out$cases[[case$name]] <- list( + pathway = res$pathway, pval = res$pval, padj = res$padj, ES = res$ES, + NES = res$NES, size = res$size, log2err = res$log2err, + leadingEdge = lapply(res$leadingEdge, as.character), + simple = seen$simple, multilevel = seen$multilevel) +} +for (case in input$calibrated) { + runs <- lapply(seq_len(input$runs), function(s) run(case, s)) + out$calibrated[[case$name]] <- list(pathway = runs[[1]]$pathway, + pval = lapply(runs, function(r) r$pval), NES = lapply(runs, function(r) r$NES)) +} +stats <- statistics("main") +for (figure in input$figures) { + data <- plotEnrichmentData(sets("main")[[figure$set]], stats, + gseaParam = figure$gseaParam) + out$figures[[length(out$figures) + 1]] <- list( + set = figure$set, gseaParam = figure$gseaParam, rank = data$curve$rank, + ES = data$curve$ES, ticks = data$ticks$rank, posES = data$posES, + negES = data$negES) +} +# plotGseaTable's bars, plotEnrichmentData's ticks over the largest absolute +# statistic, and its numbers as its labels write them: "%.2f" and x.y 10^a. +res <- out$cases[[input$table_case]] +for (name in res$pathway) { + data <- plotEnrichmentData(sets("main")[[name]], stats, gseaParam = 0.5) + out$bars[[name]] <- list(rank = data$ticks$rank, + height = data$ticks$stat / data$maxAbsStat) +} +power <- function(v) { a <- floor(log10(v)); c(sprintf("%.1f", v / 10^a), a) } +# fgseaMultilevel's choice, for each count k of random scores beyond the ES, of the +# multilevel method: when its error is below that of the permutation p-value. +multilevel <- function(n, size) { + k <- 0:n; crude <- log2((k + 1) / (n + 1)) + left <- log2(qbeta(0.025, k, n - k + 1)); right <- log2(qbeta(0.975, k + 1, n - k)) + multilevelError((k + 1) / (n + 1), size) < 0.5 * pmax(crude - left, right - crude) +} +out$multilevel <- list(n1000_size101 = multilevel(1000, 101), + n2000_size51 = multilevel(2000, 51)) +out$numbers <- list(pathway = res$pathway, NES = sprintf("%.2f", res$NES), + pval = lapply(res$pval, power), padj = lapply(res$padj, power)) +out$versions <- c(paste("R", getRversion()), paste("fgsea", packageVersion("fgsea"))) +write_json(out, file.path(root, "output.json"), digits = I(17), na = "null", + null = "null") +""" + + +class Replay(np.random.Generator): + """A generator that draws the given seeds, fgsea's, as gsea draws its seeds.""" + + def __init__(self, seeds): + super().__init__(np.random.PCG64()) + self.seeds = np.asarray(seeds, dtype=np.int64) + + def choice(self, a, size=None, replace=True, p=None, axis=0, shuffle=True): + return self.seeds[:size] - 1 # gsea's seeds run from 1, numpy's from 0 + + +# Inputs. + + +def statistics(rng, n_genes): + """Moderated t-like statistics with genes up and down and two blocks of ties, + -log10 p-values (all positive) and mostly positive statistics; the genes up, down + and tied.""" + t = rng.normal(size=n_genes) + order = rng.permutation(n_genes) + up, down = order[: n_genes * 6 // 100], order[n_genes * 95 // 100 :] + t[up] += rng.normal(2.5, 0.7, up.size) + t[down] -= rng.normal(2.5, 0.7, down.size) + tied = order[n_genes // 3 : n_genes // 3 + n_genes * 3 // 200] + t[tied] = 1.25 + t[order[n_genes // 2 : n_genes // 2 + n_genes * 3 // 400]] = 0.0 + p = rng.uniform(size=n_genes) + p[up] = rng.beta(0.3, 20.0, up.size) + mostly_positive = np.abs(t) + mostly_positive[rng.choice(n_genes, n_genes // 40, replace=False)] *= -1.0 + values = {"main": t, "positive": -np.log10(p), "unbalanced": mostly_positive} + return {name: v.tolist() for name, v in values.items()}, (up, down, tied) + + +def gene_sets(rng, genes, groups, n_sets, largest): + """Sets enriched for genes up, down, both or none, a fifth outside the ranking, + and the special sets.""" + up, down, tied = groups + everyone = np.arange(len(genes)) + pools = {"none": [everyone], "up": [up, everyone], "down": [down, everyone], + "both": [up, down]} # fmt: skip + sets = {} + for s, size in enumerate(np.exp(rng.uniform(0.7, np.log(largest), n_sets))): + kind = str(rng.choice(list(pools), p=[0.55, 0.2, 0.15, 0.1])) + drawn = [rng.choice(pool, int(size) // 2 + 1) for pool in pools[kind]] + members = [genes[i] for i in np.concatenate(drawn)] + members += [f"OUTSIDE{i}" for i in rng.choice(9999, int(size) // 5 + 1)] + sets[f"set_{s:02d}_{kind}"] = list(dict.fromkeys(members)) + sets["one_gene"] = [genes[up[0]]] + sets["tied_block"] = [genes[i] for i in tied[::2]] + sets["listed_twice"] = [genes[i] for i in up[:15]] * 2 + sets["identical_a"] = sets["identical_b"] = [genes[i] for i in down[:25]] + return sets + + +def inputs(small): + rng = np.random.default_rng(20261026) + n_genes = 300 if small else 4000 + genes = [f"G{i:05d}" for i in range(n_genes)] + values, groups = statistics(rng, n_genes) + sets = gene_sets(rng, genes, groups, 8 if small else 84, 60 if small else 1000) + sets["every_gene"] = list(genes) + # For the mostly positive statistics: sets of genes evenly spread along their + # ranking, whose permutation p-values are large; a set of a few genes near its top + # and others spread along it, whose p-value, about 0.05, is moderate but above + # those that go to the multilevel method; and sets of 25 genes from the bottom of + # the ranking, whose negative scores only a few random sets of their size reach. + ranked = np.argsort(-np.asarray(values["unbalanced"]), kind="stable") + null = {f"spread_{s}": [genes[i] for i in ranked[s :: n_genes // 25]] + for s in range(6)} # fmt: skip + near_top, step = 6 if small else 5, n_genes // 100 + tilted = [*ranked[: near_top * step : step], + *ranked[near_top * step :: (n_genes - near_top * step) // (30 - near_top)]] # fmt: skip + null["tilted"] = [genes[i] for i in tilted[:30]] + bottom = ranked[-75:] + null |= {f"bottom_{s}": [genes[i] for i in bottom[s::3]] for s in range(3)} + strongest = max( + (name for name in sets if name.endswith("_up")), key=lambda s: len(sets[s]) + ) + cases = [ + ("defaults", "main", "main", {}), + ("gseaParam 0", "main", "main", {"gseaParam": 0}), + ("gseaParam 0.5", "main", "main", {"gseaParam": 0.5}), + ("scoreType neg", "main", "main", {"scoreType": "neg"}), + ("scoreType pos, all positive", "positive", "main", {"scoreType": "pos"}), + ("sizes 15 to 500, eps 0, sampleSize 51, nPermSimple 2000", "main", "main", + {"minSize": 15, "maxSize": 500, "eps": 0, "sampleSize": 51, + "nPermSimple": 2000}), + ("permutations, nperm 2500", "main", "main", {"nperm": 2500}), + ("mostly positive statistics", "unbalanced", "null", {}), + ("one set", "main", "single", {}), + ] # fmt: skip + data = { + "genes": genes, "statistics": values, "runs": RUNS, "table_case": "defaults", + "sets": {"main": sets, "null": null, "single": {strongest: sets[strongest]}}, + "cases": [{"name": name, "statistics": s, "sets": g, "options": o, + "seed": 11 + i} for i, (name, s, g, o) in enumerate(cases)], + "figures": [{"set": s, "gseaParam": g} for s in + (strongest, "identical_a", "tied_block") for g in (1, 0.5)], + "calibrated": [], + } # fmt: skip + if not small: # 36 sets of up to 300 genes, and one set alone, seeds of their own + data["sets"]["calibrated"] = gene_sets(rng, genes, groups, 31, 300) + data["calibrated"] = [ + {"name": "36 sets", "statistics": "main", "sets": "calibrated", + "options": {}}, + {"name": "one set alone", "statistics": "main", "sets": "single", + "options": {}}, + ] # fmt: skip + return data + + +def run_r(root, data): + (root / "input.json").write_text(json.dumps(data)) + (root / "compare.R").write_text(R_SCRIPT) + rscript = os.environ.get("RSCRIPT", "Rscript") + subprocess.run( + [rscript, "--vanilla", str(root / "compare.R"), str(root)], check=True + ) + return json.loads((root / "output.json").read_text()) + + +# Comparisons. + + +def numbers(values): + return np.array([np.nan if v is None else v for v in values], dtype=float) + + +def relative(ours, reference): + """Largest difference relative to max(|reference|, 1e-300); NaN must match.""" + ours, reference = np.asarray(ours, float), np.asarray(reference, float) + if ours.shape != reference.shape or not np.array_equal( + np.isnan(ours), np.isnan(reference) + ): + return np.inf + keep = ~np.isnan(reference) + scale = np.maximum(np.abs(reference[keep]), 1e-300) + return float(np.max(np.abs(ours[keep] - reference[keep]) / scale, initial=0.0)) + + +def run_gsea(data, case, rng): + values = data["statistics"][case["statistics"]] + stats = dict(zip(data["genes"], values, strict=True)) + options = {OPTIONS[key]: value for key, value in case["options"].items()} + if "nperm" in case["options"]: + options["method"] = "permutation" + with warnings.catch_warnings(): + warnings.simplefilter("ignore", RuntimeWarning) + return gsea(stats, data["sets"][case["sets"]], rng=rng, **options) + + +def replayed(data, case, reference): + """gsea's table, given the seeds that fgsea drew.""" + seeds = list(reference["simple"]) + if "nperm" not in case["options"]: # the permutations' seed and the multilevel's + seeds.append((reference["multilevel"] or [1])[0]) + return run_gsea(data, case, Replay(seeds)) + + +def check_case(data, case, reference, table): + ours = replayed(data, case, reference) + row = {name: i for i, name in enumerate(ours["gene_set"].tolist())} + names = reference["pathway"] + at = [row.get(name, 0) for name in names] + edges = [e.split(", ") if e else [] for e in ours["leading_edge"].tolist()] + same = ( + sorted(row) == sorted(names) + and ours["n_genes"][at].tolist() == reference["size"] + and all(edges[i] == reference["leadingEdge"][j] for j, i in enumerate(at)) + ) + label = f"{case['name']}: {len(names)} set{'s' if len(names) > 1 else ''}" + value = relative(ours["enrichment_score"][at], reference["ES"]) + table.append((f"{label}: sets tested, sizes, leading edges and enrichment scores", + "fgsea", "max relative difference", value if same else np.inf, + EXACT)) # fmt: skip + if case["sets"] == "single": # fgsea's sampler for one set is R's + return + columns = [("normalized_enrichment_score", "NES"), ("p_value", "pval")] + if "nperm" not in case["options"]: + columns.append(("log2_error", "log2err")) + value = max(relative(ours[c][at], numbers(reference[r])) for c, r in columns) + adjusted = relative(ours["adjusted_p_value"][at], numbers(reference["padj"])) + what = "NES, p-values" + (", log2 errors" if len(columns) == 3 else "") + if case["sets"] == "null": + table.append((f"{label}: {what}, given fgsea's seeds", "fgsea", + "max relative difference", value, EXACT)) # fmt: skip + table.append((f"{label}: adjusted p-values, over the sets with a p-value", + "fgsea", "max relative difference", adjusted, None)) # fmt: skip + else: + table.append((f"{label}: {what} and adjusted p-values, given fgsea's seeds", + "fgsea", "max relative difference", max(value, adjusted), + EXACT)) # fmt: skip + + +def check_calibrated(data, case, reference, table): + names = reference["pathway"] + runs = [run_gsea(data, case, seed) for seed in range(1, RUNS + 1)] + for column, key, label in ( + ("p_value", "pval", "mean log10 p-value"), + ("normalized_enrichment_score", "NES", "mean NES"), + ): + theirs = np.array([numbers(run) for run in reference[key]]) + mine = np.array( + [ + [run[column][run["gene_set"].tolist().index(s)] for s in names] + for run in runs + ] + ) + if key == "pval": + theirs, mine = np.log10(theirs), np.log10(mine) + difference = np.abs(mine.mean(axis=0) - theirs.mean(axis=0)) + error = np.sqrt((mine.var(0, ddof=1) + theirs.var(0, ddof=1)) / RUNS) + z = np.where(difference == 0.0, 0.0, difference / np.maximum(error, 1e-300)) + table.append((f"{case['name']}: {label} of {RUNS} runs, seeds of their own", + "fgsea", "largest |difference| / SE", float(z.max()), + CALIBRATED)) # fmt: skip + + +def check_figures(data, output, table): + stats = dict(zip(data["genes"], data["statistics"]["main"], strict=True)) + worst, differ = 0.0, 0 + for figure in output["figures"]: + ax = plotting.plot_enrichment( + stats, data["sets"]["main"][figure["set"][0]], + exponent=figure["gseaParam"][0], ax=matplotlib.figure.Figure().subplots(), + ) # fmt: skip + walk = next(line for line in ax.lines if line.get_gid() == "running sum") + ticks = next(c for c in ax.collections if c.get_gid() == "set") + extremes = [ln.get_ydata()[0] for ln in ax.lines if ln.get_gid() == "extremes"] + differ += not np.array_equal(walk.get_xdata(), figure["rank"]) + differ += [s[0, 0] for s in ticks.get_segments()] != figure["ticks"] + worst = max(worst, relative(walk.get_ydata(), figure["ES"]), + relative(extremes, figure["posES"] + figure["negES"])) # fmt: skip + check = f"plot_enrichment: the walk, ticks and extremes of {len(output['figures'])} plots" + table.append((check, "plotEnrichmentData", "max relative difference", + worst if differ == 0 else np.inf, EXACT)) # fmt: skip + + case = next(c for c in data["cases"] if c["name"] == data["table_case"]) + result = replayed(data, case, output["cases"][case["name"]]) + ax = plotting.plot_gsea_table( + stats, data["sets"]["main"], result, top=5, exponent=0.5, + ax=matplotlib.figure.Figure().subplots(), + ) # fmt: skip + names = [label.get_text() for label in ax.get_yticklabels()] + segments = next(c for c in ax.collections if c.get_gid() == "gene ranks") + segments = segments.get_segments() + worst, differ = 0.0, 0 + for row, name in enumerate(names): + level = len(names) - 1 - row + mine = sorted( + (s[0, 0], (s[1, 1] - level) / BAR) for s in segments if s[0, 1] == level + ) + bars = output["bars"][name] + differ += [x for x, _ in mine] != bars["rank"] + worst = max(worst, relative([h for _, h in mine], bars["height"])) + texts = output["numbers"] + at = [texts["pathway"].index(name) for name in names] + written = [texts["NES"][i] for i in at] + for key in ("pval", "padj"): + written += [ + rf"${texts[key][i][0]} \cdot 10^{{{texts[key][i][1]}}}$" for i in at + ] + differ += [t.get_text() for t in ax.texts if t.get_gid() == "numbers"] != written + check = f"plot_gsea_table: the bars and numbers of {len(names)} sets" + table.append((check, "plotEnrichmentData, plotGseaTable", "max relative difference", + worst if differ == 0 else np.inf, EXACT)) # fmt: skip + + +# Fixtures. + + +def write_fixtures(data, output): + if FIXTURES.exists(): + shutil.rmtree(FIXTURES) + FIXTURES.mkdir(parents=True) + keep = ("pathway", "pval", "padj", "ES", "NES", "size", "log2err", "leadingEdge", + "simple", "multilevel") # fmt: skip + fixture = { + "genes": data["genes"], "statistics": data["statistics"], "sets": data["sets"], + "cases": data["cases"], "figures": output["figures"], + "multilevel": output["multilevel"], + "fgsea": {name: {key: reference.get(key) for key in keep} + for name, reference in output["cases"].items()}, + } # fmt: skip + (FIXTURES / "fgsea.json").write_text(json.dumps(fixture, separators=(",", ":"))) + print(f"wrote the fixtures to {FIXTURES}") + + +def main(): + parser = argparse.ArgumentParser(description=__doc__.split("\n")[0]) + parser.add_argument( + "--write-fixtures", action="store_true", help="write test fixtures" + ) + small = parser.parse_args().write_fixtures + data = inputs(small) + table = [] + with tempfile.TemporaryDirectory() as temporary: + output = run_r(Path(temporary), data) + print("Reference: " + ", ".join(output["versions"])) + for case in data["cases"]: + check_case(data, case, output["cases"][case["name"]], table) + for case in data["calibrated"]: + check_calibrated(data, case, output["calibrated"][case["name"]], table) + check_figures(data, output, table) + failures = 0 + print("| check | reference | metric | value | tolerance | result |") + print("| --- | --- | --- | --- | --- | --- |") + for check, source, metric, value, tolerance in table: + if tolerance is None: + result, shown = "deliberate difference", "n/a" + else: + result = "pass" if value <= tolerance else "FAIL" + shown = f"{tolerance:g}" + failures += result == "FAIL" + print(f"| {check} | {source} | {metric} | {value:.1e} | {shown} | {result} |") + print(f"{len(table)} checks, {failures} failed") + if small: + if failures: + print("not writing fixtures: some checks failed") + else: + write_fixtures(data, output) + return 1 if failures else 0 + + +if __name__ == "__main__": + sys.exit(main())