Gene Set Enrichment Analysis

---

• Table of Contents
  • I. Goals
  • II. Background
  • III. Significance Analysis of Function and Expression

---

I. Goals

In this section we discuss the use of Gene Set Enrichment Analysis (GSEA) to identify pathways enriched in ranked gene lists, with a particular emphasis on ordering based on a measure of differential gene expression. We aim to convey how the approach works from an intuitive standpoint before dividing into a full discussion of the statistical underpinnings. By then end of this discussion you should:

1. Understand what you can get out of GSEA
2. Be aware of the advantages over previous methods
3. Be aware of the statistical basis of the approach

II. Background

High-throughput approaches for gene expression measurement can generate a tremendous volume of data but can be unwieldy and easily outstrip intuition. The noisy nature of biological processes compounds the difficulty in interpreting outputs of large-scale experiments. Clearly, there is a need for robust approaches that place data in a wider context that is more biologically meaningful to the scientific question at hand.

To this end, approaches collectively termed ‘Overrepresentation Analyses’ (ORA) were developed to take large lists of genes emerging from experimental results and determine whether there was evidence of enrichment for gene sets grouped on the basis of some shared theme (Khatri 2005, Khatri 2012). In simple terms, ORA approaches aim to distill which pathways are present within a list of genes. A popular source of sets is the Gene Ontology (GO) which groups of genes according to various biological processes and molecular functions.

The ‘SAFE’ framework

While tremendously useful for interpreting differential expression output, ORA approaches have three major limitations. First, the inclusion criteria for input gene lists are rather arbitrary and typically involves selecting genes that exceed some user-defined statistical cutoff. This risks excluding potentially important genes that for whatever reason fail to reach statistical significance. Second, ORA approaches use gene names but not any of the rich quantitative information associated with gene expression experiments. In this way, equal importance is assigned to each an every gene. Third, many of the ORA procedures uses statistical procedures that assume independence among genes: Changes in any gene do not affect or are not influenced by any others. Clearly, this is unrealistic for biological systems and has the effect of making ORA procedures more prone to erroneous discoveries or false-positives.

Gene Set Enrichment Analysis (GSEA) is a tool that belongs to a class of second-generation pathway analysis approaches referred to as significance analysis of function and expression (SAFE) (Barry 2005). These methods are distinguished from their forerunners in that they make use of entire data sets including quantitive data gene expression values or their proxies.

Methods that fall under the SAFE framework use a four-step approach to map gene lists onto pathways

1. Calculate a local (gene-level) statistic
2. Calculate a global (gene set or pathway-level) statistic
3. Determine significance of the global statistic
4. Adjust for multiple testing

Origins of GSEA

GSEA was first described by Mootha et al. (Mootha 2003) in an attempt to shed light on the mechanistic basis of Type 2 diabetes mellitus. They reasoned that alterations in gene expression associated with a disease can manifest at the level of biological pathways or co-regulated gene sets, rather than individual genes. The lack of power to detect true signals at the gene level may be a consequence of high-throughput expression measurements which involve heterogeneous samples, modest sample sizes and subtle but nevertheless meaningful alterations expression levels. Ultimately these confound attempts to derive reproducible associations with a biological state of interest.

In their study, Mootha et al. employed microarrays to compare gene expression in mice with diabetes mellitus (DM2) to controls with normal glucose tolerance (NGT). On one hand, two statistical analysis approaches failed to identify a single differentially expressed gene between these two groups of mice. On the other hand, GSEA identified oxidative phosphorylation (OXPHOS) as the top scoring gene set down-regulated in DM2; The next four top-scoring gene sets largely overlapped with OXPHOS.

Importantly, the prominence of OXPHOS genes provided the necessary clues for the authors to hypothesize that peroxisome proliferator-activated receptor $\gamma$ coactivator 1$\alpha$ (PGC-1$\alpha$ encoded by PPARGC1) might play a role in DM2. Indeed, follow-up experiments demonstrated that PPARGC1 expression was lower in DM2 and over-expression in mouse skeletal cells was sufficient to increase OXPHOS expression. Buoyed by the success of GSEA in this case, the authors went on to suggest the particular niche that the approach might occupy.

Single-gene methods are powerful only when the individual gene effect is marked and the variance is small across individuals, which may not be the case in many disease states. Methods like GSEA are complementary to single-gene approaches and provide a framework with which to examine changes operating at a higher level of biological organization. This may be needed if common, complex disorders typically result from modest variation in the expression or activity of multiple members of a pathway. As gene sets are systematically assembled using functional and genomic approaches, methods such as GSEA will be valuable in detecting coordinated variation in gene function that contributes to common human diseases.

Criticisms concerning the original methodology (Damian 2004) were considered in an updated version of GSEA described in detail by Subramanian et al. (Subramanian 2005). Below, we provide a description of the approach with particular emphasis on the protocol we recommend for analyzing gene lists ranked according to differential expression.

III. Significance Analysis of Function and Expression

SAFE Step 1. Local statistic

In this step, we describe a local or gene-level measure that is used to rank genes, in GSEA terminology, a ‘rank metric’. Previously, we described how to obtain and process RNA-seq datasets into a single list of genes ordered according to a function of each gene’s p-value calculated as part of differential expression testing. In this context, a p-value assigned to a gene can be interpreted as the probability of a difference in gene expression between groups at least as extreme as that observed given no inter-group difference. We simply declare this function of p-values the rank metric (Figure 1).

An example of a rank metric is the product of the sign of the ‘direction’ in the expression change (i.e. ‘up-regulation’ is positive and ‘down-regulation’ is negative) and the p-value ($P$).

$$\begin{equation*} s_i=s(P_i) = sign(\text{fold change gene }i)\cdot -log_{10}(P_i) \end{equation*}$$

Under this rank metric, up-regulated genes with relatively small p-values appear at the top of the list and down-regulated genes with small p-values at the bottom.

SAFE Step 2. Global statistic

The global statistic is at the very heart of GSEA and the rather simple algorithm used to calculate it belies a rather profound statistical method to judge each candidate gene set. We will set aside the technical details for a moment to see how GSEA calculates the enrichment score for a given pathway.

Let’s pick up the process following Figure 1 where we have a list of ranked genes. For illustrative purposes, suppose we wish to test the list for enrichment of a cell cycle pathway (Figure 2).

GSEA considers candidate gene sets one at a time. To calculate the enrichment score, GSEA starts at the top of the ranked gene list. If a gene is a member of the candidate gene set then it adds to a running sum, otherwise, it subtracts. This process is repeated for each gene in the ranked list and the enrichment score for that gene set is equal to the largest absolute value that the running sum achieved.

One aspect of this algorithm we side-stepped is the value that gets added or subtracted. In the original version of GSEA (Mootha 2003) the values were chosen specifically such that the sum over all genes would be zero. In Figure 2, this would be equivalent to the running sum meeting the horizontal axis at the end of the list. In this case, the enrichment score is the Kolmogorov-Smirnov (K-S) statistic that is used to determines if the enrichment score is statistically significant. The updated procedure described by Subramanian et al (Subramanian 2005) uses a ‘weighted’ version of the procedure whereby the increments to the running sum are proportional to the rank metric for that gene. The reasons for these choices are rather technical and we reserve this for those more mathematically inclined in the following section.

The enrichment score

Consider a single gene set $G_k$ indexed by $k$. The gene set consists of a list of $n_k$ genes ($g_{kj}$), that is $G_k=\{g_{kj}: j = 1, \ldots, n_k\}$. Note that each gene in the set must be a represented in the ranked list $L$ as you will see shortly.

Define the set of genes outside of the set as $\bar{G}_k = \{\bar{g}_{kj}: 1, \ldots, n-n_k\}$. We summarize the relevant notation up to this point

• Notation
  • Number of genes: $n$
  • Gene rank metric: $s$
  • Ranked gene list: $L$
  • Gene set: $G_k$ where $k=1,\ldots,K$
  • Genes in gene set: $G_k=\{g_{kj}: j=1,\ldots,n_k\}$
  • Genes not in gene set: $\bar{G}_k=\{\bar{g}_{kj}: j=1,\ldots,n-n_k\}$

Define the enrichment score for a given gene set as $S^{GSEA}_k$ which is the (weighted) Kolmogorov-Smirnov (K-S) statistic.

$$\begin{equation} \label{eq:1} S^{GSEA}_k = \sup_{1 \leq i \leq n} (F^{G_k}_i-F^{\bar{G}_k}_i) \end{equation}$$

Where $sup$ is the supremum and the indices $i$ represents the position or rank in $L$. The $S^{GSEA}_k$ is the largest difference in $F$ which are the (weighted) empirical cumulative distribution functions.

$$\begin{equation} \label{eq:2} F^{G_k}_i = \frac{\sum\limits_{t=1}^i |s_t|^\alpha \cdot \mathbb{1}_{\{gene_t \in G_k\}}}{\sum\limits_{t=1}^n |s_t|^\alpha \cdot \mathbb{1}_{\{gene_t \in G_k\}}} \end{equation}$$ $$\begin{equation} \label{eq:3} F^{\bar{G}_k}_i = \frac{\sum\limits_{t=1}^i \mathbb{1}_{\{gene_t \in \bar{G}_k\}}}{n-n_k} \end{equation}$$

Where $\mathbb{1}$ is the indicator function for membership in the specified gene set. We will note here that the enrichment score is a function of the gene set size, which will come into play when we discuss significance testing below. To get a better feel for what these equations mean, let’s see what happens when we vary the exponent $\alpha$.

Equal weights: The ‘classic’ Kolmogorov-Smirnov statistic

Consider the simple case when $\alpha=0$. Then all contributions from genes in the gene set do not take into account the rank metric $s$ in \eqref{eq:2}. In effect, this gives all genes equal weight.

$$\begin{equation} \label{eq:4} F^{G_k}_i = \frac{1}{n_k} \sum\limits_{t=1}^i \mathbb{1}_{\{gene_t \in G_k\}} \end{equation}$$

Under these circumstances, the $S^{GSEA}_k$ calculated using \eqref{eq:4} is the ‘classic’ version of the Kolmogorov-Smirnov (K-S) statistic. It is more common to represent $F$ as an empirical cumulative distribution function (ecdf) of an order statistic ($X_{(\cdot)}$). Without loss of generality, define the order statistic as the rank or index of each gene in $L$ that is, $X_{(1)}=1, X_{(2)}=2, \cdots, X_{(i)}=i, \cdots, X_{(n)}=n$. For ease of notation, we drop the subscript brackets.

Using the order statistic notation, equation \eqref{eq:4} describing the contribution of genes within the gene set can be expressed as

$$\begin{equation} \label{eq:5} \hat{F}_{n}^{G_k}(x) = \frac{1}{n_k} \sum\limits_{t=1}^n \mathbb{1}_{\{X_t \leq x\}} \cdot \mathbb{1}_{\{gene_t \in G_k\}} \end{equation}$$

In the same way, we can express equation \eqref{eq:3} for the genes outside the gene set as

$$\begin{equation} \label{eq:6} \hat{F}_{n}^{\bar{G}_k}(x) = \frac{1}{n-n_k} \sum\limits_{t=1}^n \mathbb{1}_{\{X_t \leq x\}} \cdot \mathbb{1}_{\{gene_t \in \bar{G}_k\}} \end{equation}$$

Rather than plotting a single running sum (Figure 2) we can plot its constituent $\hat{F}_n^{G_k}(x)$ and $\hat{F}_n^{\bar{G}_k}(x)$ on the same axes, such as in Figure 3.

Upon close inspection, it is simple to relate the running sum (bottom Figure 2) with the ecdfs (Figure 3): Increases in the running sum correspond to increments in the ecdf of genes within the gene set (Figure 3, red) and likewise, decreases in the running sum correspond to increments in the ecdf for genes outside (blue). Then, the enrichment score - the largest deviation in the running sum - corresponds to the largest vertical distance between ecdf plots (dotted green).

We are now ready to handle the most important question that gets to the very validity of GSEA: When is the $S^{GSEA}_k$ ‘significant’?.

The Kolmogorov-Smirnov goodness-of-fit test

We have learned that the $S^{GSEA}$ can be represented as the biggest distance between component ecdfs. To ask whether a given $S^{GSEA}$ is significant is equivalent to asking whether the component ecdfs represent different ‘samples’ of the same cdf and whether their differences are attributable to minor sampling errors. For didactic purposes, we’ll first present a simpler case, where we set up a K-S goodness-of-fit test between a single empirical cdf derived from a sample and a theoretical ‘specified’ cdf. We will also define concepts mentioned already in a more rigorous fashion.

One sample K-S goodness-of-fit test

Suppose we have an independent and identically distributed sample $X_{1}, \ldots, X_{n}$ with some unknown cdf $\mathbb{P}$ and we wish to test whether it is equal to a specified cdf $\mathbb{P_0}$. The null hypothesis is

$$\begin{equation*} H_0 : \mathbb{P} = \mathbb{P_0} \end{equation*}$$

Concretely, define the empirical cumulative distribution function (ecdf) that is generated from the data.

$$\begin{equation} \label{eq:7} \hat{F}_n(x) = \sum\limits_{t=1}^{n} \mathbb{1}_{\{X_t \leq x\}} \end{equation}$$

Then the K-S goodness-of-fit test proceeds by assuming that the ecdf in equation \eqref{eq:7} is an estimate of a specified cdf $F(x)$.

Definition The Kolmogorov-Smirnov statistic $D_n$ is the distance between two functions.

$$\begin{equation} \label{eq:8} D_n = \sup_{x} \left|\hat{F}_n(x)-F(x)\right| \end{equation}$$

So given the null hypothesis that the ecdf is a sample from the specified cdf, we want to know how the $D_n$ behaves. In other words, if we calculate a value of $D_n$ then we wish to know if it is a discrepancy that is worthy of further investigation. It turns out that there is an analytic solution to this question, that is, there is an equation for the probability distribution of $D_n$ so that we can derive a p-value under our null hypothesis. To get to this point, we’ll need a few theorems, which we present without proof.

Theorem 1 The Glivenko-Cantelli theorem:

$$\begin{equation} \label{eq:9} D_n \xrightarrow[]{a.s.} 0 \text{ as } n \rightarrow \infty \end{equation}$$

The significance of this theorem is subtle: The sample points we use to construct our ecdf along with all those points in between are sure to be within a specified distance of the target cdf when $n > N$. This is also termed ‘uniform convergence’.

Theorem 2 The distribution-free property states that the distribution of $D_n$ is the same for all continuous underlying distribution functions $F$.

The distribution-free property is a key aspect of the K-S test and pretty powerful result. It says that regardless of whether the $F$ is normal, uniform, gamma or even some completely unknown distribution, they all have the same $D_n$ distribution. This is particularly useful because we won’t have any idea what the distribution of the ecdfs used to construct our GSEA running sum will be.

The Glivenko-Cantelli and the distribution-free properties are nice, but not very useful in practice. Knowing that an ecdf will converge regardless of the form of the distribution is just the start. What we really want to know is how the convergence happens, that is, how $D_n$ is distributed. This leads us right into the next theorem.

Theorem 3 Define $K(x)$ as the Kolmogorov-Smirnov distribution then,

$$\begin{equation} \label{eq:10} P(\sqrt{n} D_n \leq x) \rightarrow K(x) = \sum\limits_{k=-\infty}^{\infty} (-1)^{k}\exp(-2k^2x^2) \end{equation}$$

Theorem 3 is the culmination of an extensive body of knowledge surrounding empirical process theory that, in simple terms, describes the distribution of random walks between two fixed endpoints and bounds (i.e. the interval $[0,1]$ as these are the cdf bounds) otherwise known as a ‘Brownian Bridge’.

The practical outcome of these theorems is that it gives us an equation from which we can calculate the exact probability of our maximum ecdf deviation from a specified cdf. In the K-S hypothesis testing framework, we set an a priori significance level $\alpha$ and calculate the probability of our observed $D_n$ or anything more extreme, denoting this the p-value $P$. If $P \lt \alpha$ then this would suggest a discrepancy between the data and the specified cdf causing us to doubt the validity of $H_0$.

Recall that by definition, the enrichment score for any given candidate gene set depends upon its size. In this case, the scores are not identically distributed and must be normalized prior to calculating the p-values. We reserve this discussion for the next section we when broach the topic of null distributions.

Two sample K-S goodness-of-fit test

We are now ready to describe the setup that occurs in GSEA where we compare ecdfs representing the distribution of genes within and outside the gene set. Suppose we have a sample $X_1,\ldots,X_{a}$ with cdf $F^{G_k}(x)$ and a second sample $Y_1,\ldots,Y_{b}$ with cdf $F^{\bar{G}_k}(x)$. GSEA is effectively a test of the null hypothesis

$$\begin{equation*} H_0 : F^{G_k} = F^{\bar{G}_k} \end{equation*}$$

The corresponding ecdfs are $\hat{F}_a^{G_k}$ and $\hat{F}_b^{\bar{G}_k}$as before and the K-S statistic is

$$\begin{equation} \label{eq:11} D_{ab} = \left(\frac{ab}{a+b}\right)^{1/2} \sup_{x} \left|\hat{F}^{G_k}_a(x)-\hat{F}^{\bar{G}_k}_b(x)\right| \end{equation}$$

and the rest is the same.

Unequal weights: Boosting sensitivity

Consider the case when $\alpha=1$ which is the recommended setting in GSEA. Then the global statistic is weighted by the value of the rank metric $s$ in \eqref{eq:2}. Effectively, this renders the enrichment score more sensitive to genes at the top and bottom of the gene list $L$ compared to the classic K-S case.

The choice to depart from the classic K-S statistic was two-fold. First, Subramanian et al. noticed an unwanted sensitivity of the enrichment score to genes in the middle of the list.

In the original implementation, the running sum statistic used equal weights at every stop, which yielded high scores for the sets clustered near the middle of the ranked list…These sets do not represent biologically relevant correlation with the phenotype.

Second, the authors were unsatisfied with the observation that the unweighted method failed to identify well-known gene expression responses, in particular, the p53 transcriptional program in wild-type cells.

In the examples described in the text, and in many other examples not reported, we found that p = 1 (weighting by the correlation) is a very reasonable choice that allows significant genes with less than perfect coherence, i.e. only a subset of genes in the set are coordinately expressed, to score well.

While weighting the enrichment score calculation with the rank metric can increase the power of the approach to detect valid gene sets, it does have several consequences that must be kept in mind.

Recall that the enrichment score is a function of the size of the gene set $n_k$. This means that enrichment scores must be normalized for gene set size. In the unweighted case, an analytic expression for the normalization factor can be estimate but when terms are weighted, this is no longer the case and once again, we defer to bootstrap measures to derive it.

Another trade-off incurred by departing from the classic K-S statistic is that we no longer have an analytic expression for the null distribution of the enrichment scores as described by equation \eqref{eq:10}. This motivates the empirical bootstrap procedure for deriving the null distribution discussed in the following section concerning significance testing.

The above discussion motivates the next section on how GSEA generates null distributions for candidate gene set enrichment scores.

SAFE Step 3. Significance testing

To recap, GSEA uses the set of rank metrics for a gene list to calculate a set of enrichment scores for candidate gene sets. The primary issue at this point is which scores are indicative of enrichment? In hypothesis testing jargon, we wish to determine the statistical significance of each global statistic. We accomplish this by deriving a p-value $P$ representing the probability of observing a given enrichment score or anything more extreme. To do this, we require some understanding of how statistics are distributed.

Null distributions

From our discussion of the global statistic, using a weighted enrichment score leaves us without an analytic description of their null distribution. That is, weighting the enrichment score $S_k^{GSEA}$ with the local statistic deviates from the classic Kolmogorov-Smirnov statistic which would typically follow a K-S-like distribution.

GSEA employs ‘resampling’ or ‘bootstrap’ methods to derive an empirical sample of the null distribution for the enrichment scores of each gene set. The GSEA software provides a choice of two flavours of permutation methods that underlie the null distribution calculations (Figure 4).

Phenotype permutation

The SAFE and GSEA publications describe ‘phenotype’ permutation approach to sample the null distribution (Figure 4A). For a given gene set $G_k$, this amounts to randomly swapping sample labels and recalculating a permutation enrichment score $S^{GSEA}_{k,{\pi_p}}$. This process is repeated $N_k$ times to generate $S^{GSEA}_{k,{\boldsymbol{\pi(k)}}} = \{S^{GSEA}_{k,{\pi_p}}: p=1,\ldots,N_k\}$ which is a vector of points sufficient to depict the underlying distribution.

From an intuitive standpoint, this generates a sample of the enrichment score distribution under the assumption of no particular association between gene rank and phenotype. In other words, we get a sense of how widely enrichment scores vary and how often when the two groups are effectively the same.

From a statistical standpoint, the authors claim that this provides a more accurate description of the null model.

Importantly, the permutation of class labels preserves gene-gene correlations and, thus, provides a more biologically reasonable assessment of significance than would be obtained by permuting genes.

Indeed, variations on GSEA that purport to simply the methodology rely on an assumption of gene-gene independence (Irizarry 2009). However, empirical comparisons suggest that ignoring these correlations leads to variance inflation - a wider and flatter null distribution - resulting in a much higher risk of false discovery for gene sets (Tamayo 2016).

Gene set permutation

The workflow we recommend uses as input a ‘pre-ranked’ list of genes are ordered by a function of the p-value for differential expression. In GSEA software this is called ‘GSEAPreranked’ and precludes phenotype permutation.

Rather, a gene set permutation approach is used to generate the null distribution (Figure 4B). For each gene set $G_k$ of size $n_k$, the same number genes are randomly selected from the ranked list $L$ and the corresponding enrichment score $S^{GSEA}_{k,{\pi_p}}$ is calculated. This process is repeated $N_k$ times to generate the sample null distribution consisting of the vector $S^{GSEA}_{k,{\boldsymbol{\pi(k)}}}$.

The GSEA team recommends using phenotype permutation whenever possible. This preserves the correlation structure between the genes in the dataset. Gene set permutation creates random gene sets and so disrupts the gene-gene correlations in the data. Thus, gene set permutation provides a relatively weaker (less stringent) assessment of significance.

P-value calculation

Once an empirical null distribution of values is in hand using either of the permutation methods described above, it is straightforward to calculate the p-value $P$ for an enrichment score (Figure 4C). By definition, $P$ is the probability of observing a statistic or anything more extreme under the null hypothesis. In practice, we derive an empirical p-value $P_k$ by calculating the fraction of null values $S^{GSEA}_{k,{\pi_p}}$ greater than or equal to our observed value $S^{GSEA}_{k}$.

$$\begin{equation}\label{eq:12} P_k = \frac{1}{N_k}\sum\limits_{p=1}^{N_k} I\{S^{GSEA}_{k,{\pi_p}} \geq S^{GSEA}_{k}\} \end{equation}$$

On last thing…

Note that throughout this discussion we have chosen to depict the GSEA global statistic as positive deviations of the running sum (Figure 2). Of course, there is no particular reason why the deviations could not be negative. Indeed a more accurate description of the results of permutation would see a bimodal distribution representing global statistics that lie above and below zero (Figure 5). Thus, it is more accurate to say that the p-values are calculated from the positive or negative region of the empirical null distribution.

SAFE Step 4. Multiple testing correction

When we test a family of hypotheses, the chance of observing a statistic with a small p-value increases. When smaller than the significance level, they can be erroneously classified as discoveries or Type I errors. Multiple testing procedures attempt to quantify and control for these.

In GSEA, the collection of gene sets interrogated against the observed data are a family of hypotheses. The recommended procedure for quantifying Type I errors is the false discovery rate (FDR). The FDR is defined as the expected value of the fraction of rejected null hypotheses that are in fact true. In practice, GSEA establishes this proportion empirically.

In general, given a specified threshold $T$ of the global statistic, the FDR is the number of true null hypotheses larger than $T$ divided by the sum of true and false null hypotheses larger than $T$. For GSEA, $T$ would be some value of the enrichment score and a true null hypothesis would represent a gene set that is in fact not enriched. In practice, we won’t directly observe the latter so it is estimated from the values of the empirical null that exceed $T$ (Figure 6).

Suppose there are $n_{obs}$ observed and $n_{null}$ empirical null distribution data points. If the number of observed and null statistics beyond the threshold $T$ are $s_{obs}$ and $s_{null}$, respectively, then the empirical false discovery rate is

$$\begin{equation} \label{eq:13} \hat{FDR} = \frac{s_{null}}{s_{obs}} \cdot \frac{n_{obs}}{n_{null}} \end{equation}$$

For a nice primer on empirical methods for FDR estimation, we refer the reader to a well written piece by William S. Noble entitled ‘How does multiple testing correction work?’ (Noble 2009).

At first glance, this appears to be quite a straightforward exercise: We retrieve our sets of null distributions $S_{k,\boldsymbol{\pi(k)}}^{GSEA}$ and observed enrichment statistics $S_k^{GSEA}$ and calculate the ratio for a given enrichment score threshold. However, there is one glaring problem: The enrichment statistic depends on gene set size which precludes comparison across gene sets.

Normalized enrichment score: Accounting for gene set size

That the enrichment score depends on gene set size can be seen from the definition in equations \eqref{eq:1} - \eqref{eq:3}. Consequently, unless we are in an unlikely scenario where all the gene sets we test are of equal size, the enrichment scores will not be identically distributed and hence cannot be directly compared. This precludes the calculation of an empirical false discovery rate.

GSEA solves this problem by applying a transformation to calculated enrichment scores such that they lie on a comparable scale. In particular, this normalized enrichment score (NES) is the enrichment score divided by the expected value (i.e. average) of the corresponding null distribution. Concretely, for each gene set we derive an enrichment score $S_k^{GSEA}$ and a corresponding null distribution $S_{k,\boldsymbol{\pi(k)}}^{GSEA}$ via gene set permutation.

The normalized enrichment score $\zeta_k^{GSEA}$ is the score divided by the expected value of the corresponding null

$$\begin{equation} \label{eq:14} \zeta_k^{GSEA} = \frac{S_k^{GSEA}}{E[S_{k,\boldsymbol{\pi(k)}}^{GSEA}]} \end{equation}$$

In practice, we will partition the null into positive and negative values and use the average of these to normalize the positive and negative enrichment scores, respectively.

$$\begin{equation*} \zeta_k^{GSEA} = \begin{cases} \frac{S_k^{GSEA}}{E[S_{k,\boldsymbol{\pi(k)}}^{GSEA} \mid S_{k,\boldsymbol{\pi(k)}}^{GSEA} \geq 0]} & \text{for } S_{k}^{GSEA} \geq 0\\ \frac{S_k^{GSEA}}{E[S_{k,\boldsymbol{\pi(k)}}^{GSEA} \mid S_{k,\boldsymbol{\pi(k)}}^{GSEA} \lt 0]} & \text{for } S_{k}^{GSEA} \lt 0\\ \end{cases} \end{equation*}$$

Remember that, in addition to the gene set enrichment scores (Figure 6, blue), we will also need a normalized null distribution $\zeta_{k, \boldsymbol{\pi(k)}}^{GSEA}$ for every gene set (Figure 6, red). Each element of a given null distribution will be determined in a similar fashion

$$\begin{equation*} \zeta_{k,\pi_p}^{GSEA} = \begin{cases} \frac{S_{k,\pi_p}^{GSEA}}{E[S_{k,\boldsymbol{\pi(k)}}^{GSEA} \mid S_{k,\boldsymbol{\pi(k)}}^{GSEA} \geq 0]} & \text{for } S_{k,\pi_p}^{GSEA} \geq 0\\ \frac{S_{k,\pi_p}^{GSEA}}{E[S_{k,\boldsymbol{\pi(k)}}^{GSEA} \mid S_{k,\boldsymbol{\pi(k)}}^{GSEA} \lt 0]} & \text{for } S_{k,\pi_p}^{GSEA} \lt 0\\ \end{cases} \end{equation*}$$

We can now restate the null hypothesis

$$\begin{equation*} \begin{split} H_0^{GSEA} &: \zeta_{k}^{GSEA}, \ldots, \zeta_{k}^{GSEA}, \ldots, \zeta_{K}^{GSEA} \text{ are identically distributed and } \\ \zeta_{k}^{GSEA} &\sim F_0^{permutation}\\ \end{split} \end{equation*}$$

where $F_0^{permutation}$ is derived empirically through $\zeta_{k,\boldsymbol{\pi(k)}}^{GSEA}$ over all gene sets.

---

Given the observed and null normalized enrichment scores, the FDR can be calculated using an approach similar to that described in the previous section (equation \eqref{eq:13}).

---

References

• Ackermann M & Strimmer K.A general modular framework for gene set enrichment analysis. BMC Bioinformatics. 2009 Feb 3;10:47. doi: 10.1186/1471-2105-10-47.
• Barry WT et al. Significance analysis of functional categories in gene expression studies: a structured permutation approach. Bioinformatics. 2005 May 1;21(9):1943-9. doi: 10.1093/bioinformatics/bti260. Epub 2005 Jan 12.
• Damian D & Gorfine M. Statistical concerns about the GSEA procedure. Nat Genet. 2004 Jul;36(7):663; author reply 663. doi: 10.1038/ng0704-663a.
• Irizarry RA et al. Gene set enrichment analysis made simple. Stat Methods Med Res. 2009 Dec;18(6):565-75. doi: 10.1177/0962280209351908.
• Khatri P & Drăghici S. Ontological analysis of gene expression data: current tools, limitations, and open problems. Bioinformatics. 2005 Sep 15;21(18):3587-95. doi: 10.1093/bioinformatics/bti565. Epub 2005 Jun 30.
• Khatri P et al. Ten years of pathway analysis: current approaches and outstanding challenges. PLoS Comput Biol. 2012;8(2):e1002375. doi: 10.1371/journal.pcbi.1002375. Epub 2012 Feb 23.
• Mootha VK et al. PGC-1alpha-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet. 2003 Jul;34(3):267-73. doi: 10.1038/ng1180.
• Noble WS. How does multiple testing correction work?Nat Biotechnol. 2009 Dec;27(12):1135-7. doi: 10.1038/nbt1209-1135.
• Subramanian A et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles.Proc Natl Acad Sci U S A. 2005 Oct 25;102(43):15545-50. doi: 10.1073/pnas.0506580102. Epub 2005 Sep 30.
• Tamayo P et al. The limitations of simple gene set enrichment analysis assuming gene independence. Stat Methods Med Res. 2016 Feb;25(1):472-87. doi: 10.1177/0962280212460441. Epub 2012 Oct 14.
• Creixell P et al. Pathway and network analysis of cancer genomes. Nat Methods. 2015 Jul;12(7):615-621. doi: 10.1038/nmeth.3440.