WTFgenes: What's The Function of these genes? Static sites for model-based gene set analysis

Introduction

Term Enrichment Analysis (TEA) is a common technique for finding functional patterns, specifically overrepresented ontology terms, in a set of experimentally identified genes¹. The most common approach, which we refer to as Frequentist TEA, is a one-tailed Fisher’s Exact Test (based on the hypergeometric distribution, which models the number of term-associations if the gene set was chosen by chance), with a suitable correction for multiple hypothesis testing. Frequentist TEA has been implemented many times on various platforms^1–8.

A model-based alternative to Frequentist TEA, which more directly addresses some of the multiple testing issues (for example, by modeling the ways that an observed gene list can be broken down into complementary gene sets), is Bayesian TEA. In contrast to Frequentist TEA, which just rejects a null hypothesis that genes are chosen by chance, the Bayesian TEA explicitly models the alternative hypothesis that the gene set was generated from a few random ontology terms. This approach was introduced by 9 and further developed by 10, who implemented model-based testing in Java and R¹¹. However, the model-based approach remains significantly less well-explored than frequentist approaches.

The graphical model underpinning Bayesian TEA is sketched in Figure 1. For each of the m terms there is a boolean random variable T_j (“term j is activated”). For each of the n genes there is a directly-observed boolean random variable O_i (“gene i is observed in the gene set”), and one deterministic boolean variable H_i (“gene i is activated”) defined by H_i = 1 − ${\Pi }_{j\in G_i}$ (1 − T_j), where G_i is the set of terms associated with gene i (including directly annotated terms, as well as ancestral terms implied by transitive closure of the directly annotated terms). The probability parameters are π (term activation), α (false positive) and β (false negative), and the respective hyperparameters are p = (p₀, p₁), a = (a₀, a₁) and b = (b₀, b₁).

Figure 1. Model-based explanation of observed genes (O_i) using ontology terms (T_j), following¹⁰.

Other variables and hyperparameters are defined in the text. Circular nodes indicate continuous-valued variables or hyperparameters; square nodes indicate discrete-valued (boolean) variables. Dashed lines indicate deterministic relationships; shaded nodes indicate observations. Plates (rounded rectangles) indicate replicated subgraph structures.

The model is

with π ∼ Beta(p), α ∼ Beta(a) and β ∼ Beta(b). The model of 10 is similar, but uses an ad hoc discretized prior for π, α and β .

Most Bayesian and Frequentist TEA implementations are designed for desktop use. Several Frequentist TEA implementations are designed for the web, such as DAVID-WS⁶ and Enrichr^8,12,13, which has a rich dynamic web front-end. However, web-facing Frequentist TEA implementations generally require a server-hosted back end that executes code. Further, there are no JavaScript-based Bayesian TEA implementations, and no web-facing implementations other than the Java-based Ontologizer which can be loaded via Java Web Start.

In order to further explore the model-based TEA and compare it to Frequentist TEA, and to make these investigations accessible to researchers in a way that would be easily embeddable in static websites, we developed WTFgenes, a JavaScript implementation of both approaches with (for time-sensitive applications) a parallel C++ implementation that is numerically identical.

We note in passing that Fisher’s Exact Test—which we call Frequentist TEA—was originally motivated by a blind tea-tasting challenge¹⁴.

Methods

Model

In developing our Bayesian TEA sampler, we introduce a collapsed version of the model in Figure 1 by integrating out the probability parameters. Let c_p = ${\sum }_j^m$ T_j count the number of activated terms, c_g = ${\sum }_i^n$H_i the activated genes, c_a = ${\sum }_i^n$ O_i (1 – H_i) the false positives and c_b = ${\sum }_i^n$ O_i H_i the false negatives.

Then

                                                                                P(T, O|a, b, p) = Z(c_p;m, p)Z(c_a; n – c_g, a)Z(c_b; c_g, b)

where

$$Z(k;N,\mathit{\text{A}})=\frac{B(N-k+A_0,k+A_1)}{B(A_0,A_1)}$$

is the beta-Bernoulli distribution for k ordered successes in N trials with hyperparameters A= (A₀, A₁), using the beta function

$$B(x,y\text{)}={\displaystyle {\int }_0^1{t}^{x-1}{(1-t)}^{y-1}dt=}\frac{\Gamma (x)\Gamma (y)}{\Gamma (x+y)}$$

Integrating out probability parameters improves sampling efficiency and allows for higher-dimensional models where, for example, we observe multiple gene sets and give each term its own probability π_j or each gene its own error rates (α_i, β_i). Our implementation by default uses uninformative priors with hyperparameters a = b = p = (1, 1), but this can be overridden by the user.

The MCMC sampler uses a Metropolis-Hastings kernel¹⁵. Each proposed move perturbs some subset of the term variables. The moves include flip, where a single term is toggled; step, where any activated term and any one of its unactivated ancestors or descendants are toggled; jump, where any activated term and any unactivated term are toggled; and randomize, where all term variables are uniformly randomized. The relative rates of these moves can be set by the user.

The sampler of 10 implemented only the flip move. To test the relative efficacy of the newly-introduced moves we measured the autocorrelation of the term variables for a dataset of 17 S.cerevisiae genes involved in mating (The gene IDs used in this evaluation, for purposes of reproduction, were: STE2, STE3, STE5, GPA1, SST2, STE11, STE50, STE20, STE4, STE18, FUS3, KSS1, PTP2, MSG5, DIG1, DIG2, STE12. Other representative gene sets for yeast may be obtained from the Gene Ontology website at http://geneontology.org/experimental/enrichment-genesets/yeast/ and several of these are bundled with the example dataset in the WTFgenes repository). The results, shown in Figure 2, led us to set the MCMC defaults, such that the flip, step, and jump moves are equiprobable, while randomize is disabled.

Figure 2. Autocorrelation of term variables, as a function of the number of MCMC samples, for several MCMC kernels on a set of 17 S.cerevisiae mating genes.

A rapidly-decaying curve indicates an efficiently-mixing kernel. The kernel incorporating flip, step and jump moves (defined in the text) mixes most efficiently.

Results

When compiled using clang, the C++ version of WTFgenes is about twice as fast as the JavaScript version: a benchmark of Bayesian TEA on a late-2014 iMac (4GHz Intel Core i7), using the above mentioned 17 yeast mating genes and the relevant subset of 518 GO terms, run for 1,000 samples per term, took 37.6 seconds of user time for the C++ implementation and 79.8 seconds in JavaScript.

By contrast, the Frequentist TEA approach is almost instant. However, its weaker statistical power is apparent from Figure 3, which compares the recall vs specificity of Bayesian and Frequentist methods on simulated datasets (The full workflow for this simulation is available at http://doi.org/10.5281/zenodo.400608¹⁶). For values of N from 1 to 4, we sampled N terms from the S.cerevisiae subset of the Gene Ontology, and generated a corresponding set of yeast genes with false positive rate 0.1% and false negative rate 1%. The MCMC sampler was run for 100 iterations per term, and this experiment was repeated 100 times. The model-based approach has vastly superior recall to the Fisher exact test, and the difference grows with the number of terms.

Figure 3. ROC curves for Frequentist and Bayesian TEA.

The axes are scaled per term. There are 5,919 ontology terms annotated to S.cerevisiae genes, so (for example) a false discovery rate of 0.001 corresponds to about 6 falsely reported terms.

Discussion

JavaScript genome browsers, such as JBrowse¹⁷, represent a broader web trend of producing static sites where possible, for reasons of security and performance. We have implemented such a static site generator for ontological term enrichment analysis of gene sets that offers both Bayesian and frequentist tests. In contrast with existing web services for Frequentist TEA, such as DAVID-WS or Enrichr, it requires no server resources and allows comparison of Bayesian and Frequentist approaches.

Model-based TEA is versatile: it can readily be extended to allow for datasets that are structured temporally¹⁸, spatially¹⁹, or by genomic region²⁰; to use domain-specific biological knowledge²¹; or to incorporate additional lines of evidence such as quantitative data²². We hope our development of a collapsed likelihood, and evaluation of different MCMC kernels, will assist these efforts.

Software and data availability

Latest source code: https://github.com/evoldoers/wtfgenes

Archived source code as at time of publication: http://doi.org/10.5281/zenodo.400606²³

Software license: BSD3

A demonstration for the Gene Ontology is usable at https://evoldoers.github.io/wtfgo.

A Makefile-driven simulation study underpinning results reported in this paper is available at http://doi.org/10.5281/zenodo.400608¹⁶.