Over-representation analysis

Over-representation analysis#

over_representation() tests each gene set for over-representation among the genes called significant, by the one-sided hypergeometric test of limma’s kegga and goana: for the genes that go up, those that go down and those that change either way in a group and contrast of a differential-expression result or in a table of results, or for a list of genes against a universe. Its table gives each set’s genes in the universe, the overlap, the expected overlap, the p-value, the adjusted p-value and the genes of the overlap. plot_over_representation() draws the top sets.

Over-representation analysis: gene sets among the genes called significant.

over_representation() asks of each gene set whether the genes called significant in an analysis hold more of its genes than chance would give, by the one-sided hypergeometric test against a universe of genes, as limma’s kegga and goana compute it with their defaults (Young et al. 2010). From a differential-expression result it tests the genes that go up and the genes that go down separately, as kegga and goana test a fit of limma, and the genes that change in either direction together. Results match limma 3.62.1 except where noted under “Deliberate differences”.

The test:

The universe holds the N genes that could have been called significant. Of a gene set, K genes are in the universe; n genes of the universe are selected, and k of them are in the set. Were the n genes drawn at random from the universe, without regard to the set, k would follow the hypergeometric distribution of n draws without replacement from N genes of which K are in the set, whose mean K n / N is the expected overlap. The p-value is P(X >= k), the one-sided p-value of Fisher’s (1922) exact test of the 2 x 2 table of the genes of the universe, in the set or not and selected or not (Rivals et al. 2007); the hypergeometric test of functional categories goes back to Tavazoie et al. (1999). It is computed by hypergeom_sf(), which keeps its relative accuracy far into the tail. The adjusted p-values are Benjamini and Hochberg’s (1995), over the sets of each list of selected genes; limma’s kegga and goana report p-values only.

The genes and the universe:

  • From a DifferentialExpressionResult (one group and contrast) or any table of results, the genes whose adjusted p-value is below alpha are selected, as kegga selects the genes of a fit whose adjusted p-value is below its FDR (strictly, found by probing): up where the log2 fold change is positive, down where it is negative, and in either direction for both lists together. The universe is the genes tested: those with a p-value. A gene id on several rows, such as a symbol shared by two features, is one gene, in the universe when any of its rows was tested and selected up or down when any of its rows is, so that it can be in both lists, as kegga counts it; rows without an id are left out, as kegga leaves them out (both found by probing).

  • From a list of gene ids: the list is the selection, and the universe must be given.

  • A universe given replaces the default. Selected genes and genes of a set outside it are left out, and genes of the universe that are in no set stay in it, as in kegga by default (restrict.universe = FALSE).

The universe matters: the genes of the genome, or of the sets, include genes that could not have been selected, such as genes not expressed or filtered out before the tests, and make p-values too small (Timmons et al. 2015; Wijesooriya et al. 2022). The default, the genes tested, is the universe that the test assumes.

What the test assumes:

The test takes genes as independent draws. The genes of a set are often correlated, as co-regulated genes are, which makes its p-values too small (Goeman and Buehlmann 2007); camera() allows for the correlation. The test also depends on the cutoff alpha and ignores the size of each change. In RNA-seq data, long and highly expressed genes are called significant more often, which biases the test towards sets of such genes (Young et al. 2010); the adjustment of kegga and goana for such a trend (trend, Wallenius’ noncentral hypergeometric distribution) is not implemented.

Deliberate differences from limma:

  • A set without genes in the universe has a row with no p-value (NaN), which the adjustment leaves out; kegga leaves the set out of its table.

  • kegga’s method for a fit of limma keeps the genes without a p-value in the universe (found by probing); here a gene without a p-value was not tested and is left out, as the test assumes.

  • kegga’s method for a fit adjusts the p-values of the rows with an id only, leaving the others out of the Benjamini-Hochberg adjustment, so that the genes it selects depend on which rows have an id, such as a symbol (found by probing). Here the adjusted p-values of the results are used as they are, so that the genes selected are those the analysis called significant.

References

Young MD, Wakefield MJ, Smyth GK, Oshlack A (2010). Gene ontology analysis for RNA-seq: accounting for selection bias. Genome Biology 11, R14. doi:10.1186/gb-2010-11-2-r14

Fisher RA (1922). On the interpretation of chi-squared from contingency tables, and the calculation of P. Journal of the Royal Statistical Society 85, 87-94. doi:10.2307/2340521

Rivals I, Personnaz L, Taing L, Potier MC (2007). Enrichment or depletion of a GO category within a class of genes: which test? Bioinformatics 23, 401-407. doi:10.1093/bioinformatics/btl633

Tavazoie S, Hughes JD, Campbell MJ, Cho RJ, Church GM (1999). Systematic determination of genetic network architecture. Nature Genetics 22, 281-285. doi:10.1038/10343

Timmons JA, Szkop KJ, Gallagher IJ (2015). Multiple sources of bias confound functional enrichment analysis of global -omics data. Genome Biology 16, 186. doi:10.1186/s13059-015-0761-7

Wijesooriya K, Jadaan SA, Perera KL, Kaur T, Ziemann M (2022). Urgent need for consistent standards in functional enrichment analysis. PLoS Computational Biology 18, e1009935. doi:10.1371/journal.pcbi.1009935

Goeman JJ, Buehlmann P (2007). Analyzing gene expression data in terms of gene sets: methodological issues. Bioinformatics 23, 980-987. doi:10.1093/bioinformatics/btm051

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.over_representation.over_representation(genes, gene_sets, *, group=None, contrast=None, universe=None, alpha=0.05, id_column='feature_id', p_value_column='adjusted_p_value', fold_change_column='log2_fold_change')#

Test each gene set for over-representation among the genes called significant.

The one-sided hypergeometric test of limma’s kegga and goana (Young et al. 2010): whether the selected genes hold more of a set’s genes than random draws from the universe would. The module docstring gives the method, how the genes and the universe are read, and the deliberate differences from limma.

Parameters:

  • genes – The genes to test: a DifferentialExpressionResult, with group and contrast; a table of results, a Table or a pandas DataFrame, with a column of gene ids, a column of adjusted p-values and, for directions, a column of fold changes; or a list of the ids of the selected genes, with universe.

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

  • universe – The ids of the genes that could have been selected. By default, for a result or a table, the genes tested; a list of genes needs one.

  • alpha – The adjusted p-value below which a gene is selected, in (0, 1). Not used for a list of genes.

  • id_column – The column of gene ids, such as "feature_id" or "symbol".

  • p_value_column – The column of the p-values that select genes, such as "adjusted_p_value" or "adjusted_p_value_across_groups".

  • fold_change_column – The column whose sign gives the direction of a selected gene, or None to test the selected genes without directions.

Returns:

A table with a row for each gene set and each list of selected genes. For a result or a table with fold changes, the rows hold the sets in the order of gene_sets for the genes up, then for those down, then for those in either direction; for a list of genes, or without fold changes, one row per set. Columns: gene_set; direction, "up", "down" or "either", or None for genes selected without a direction; n_genes, the set’s genes in the universe (K); overlap, how many of them are selected (k); expected, the expected overlap K n / N; p_value, P(X >= k); adjusted_p_value, Benjamini-Hochberg over the sets of the same direction; n_selected, the selected genes in the universe (n); n_universe, the genes of the universe (N); and overlap_genes, the selected genes of the set, in the set’s order, separated by ", ". A set without genes in the universe has NaN p-values.

Raises:

  • TypeError – If an argument has the wrong type, gene ids are not strings, a column has the wrong type, or an argument is or holds a numpy masked array.

  • ValueError – If gene_sets is empty, alpha is not in (0, 1), a table lacks a column, has p-values outside [0, 1] or misses the fold change of a gene with a p-value, a list of genes is given without a universe, group or contrast is given without a DifferentialExpressionResult or left out of one that has several, or no gene of any set, or no selected gene, is in the universe.

  • KeyError – If the result has no table of the group and contrast; the message says why.

Explained in: Choosing a method, Gene-set tests