Gene set enrichment analysis

Contents

Gene set enrichment analysis#

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. plot_enrichment() draws the running enrichment score of one set, and plot_gsea_table() the top sets with their genes on the ranking.

Gene set enrichment analysis (GSEA) of genes ranked by a statistic, as fgsea computes it.

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. plot_enrichment() and plot_gsea_table() draw its figures.

The genes and the sets:

The statistics come from a group and contrast of a 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

cytepool.gsea.gsea(statistics, gene_sets, *, group=None, contrast=None, id_column='feature_id', statistic_column='t', min_size=1, max_size=None, exponent=1.0, score_type='std', method='multilevel', n_permutations=1000, sample_size=101, eps=1e-50, rng=None, n_threads=None)#

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.

Parameters:

  • statistics – The statistics that rank the genes: a 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 GeneSets, such as one read by 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:

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.

Return type:

A table with one row per set tested, in the order of gene_sets

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]

Used in: Gene sets and pathways

Explained in: Choosing a method, Gene-set tests

Contents