Gene Set Enrichment Analysis

Original post

Gene Set Enrichment Analysis

Gene Set Enrichment Analysis (GSEA) is a statistical analysis method used to determine the degree of association between a gene set and a phenotype of interest. This method is particularly useful when we have prior knowledge of well-established gene sets, such as pathway information.

The analysis involves calculating the Enrichment Score (ES) using the modified Kolmogorov-Smirnov (KS) statistic for each gene set. The modification is a straightforward approach where weights, such as fold-change or correlation coefficient between gene and trait of interest, are assigned to individual genes within a specific gene set before calculating the KS statistic.

The figure below illustrates the overall workflow of GSEA.

From the paper, “ Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles”, 2005

Figure A represents an expression matrix, while B represents a pre-existing known gene set. B is sorted based on its correlation with the current phenotype of interest through A. This sorted gene set, with rankings or correlations as weights, is used to calculate the KS statistic. The subset of the gene set before reaching the peak is referred to as the leading edge subset.

The specific calculation of the score can be seen on the far right of the figure. Based on the ordered gene set, values are added in order from left to right. A running-sum statistic is calculated by adding + if a gene is in the gene set (with correlation as a weight) or - if it's not. This distribution is shown in the bottom right diagram.

Accordingly, if the specific gene set is consistently placed towards the lower ranks, the graph will be drawn in a pattern where the values keep decreasing and then increasing. In this case, a peak will be formed predominantly below the halfway mark, causing more negative peak to fall below zero compared to positive peaks. Consequently, the ES will be negative.

The specific equation for calculating ES is as follows.

 

$$ P_{hit}(S,i)=\sum_{g_j\in S, j\leq i}\frac {|r_j|^p} {N_R}, where N_R=\sum_{g_j\in S}|r_j|^p $$

$$ P_{miss}(S,i)= \sum_{g_j\in S, j\leq i}\frac {1} {N-N_H} $$

 

The ES is the maximum deviation from zero of $P_{hit} – P_{miss}$. In this equation, $S$ represents the Signature set, $r_j$ is the weight of the $j^{th}$ gene, $p$ is the power term that adjusts the weights, $N_R$ is the absolute sum of the adjusted weights, and $N_H$ is the total number of signature genes.

First, the genes are sorted in descending order, and if they are present in the signature, they are considered a hit; otherwise, they are considered a miss. The sorted index of each gene corresponds to the $i$ in the equation, and for each $i$, the value of $P_{hit} - P_{miss}$ is calculated to determine the ES value when it is farthest from 0. Weight and penalty are normalized separately and added or subtracted from each other. Ultimately, the maximum value is 1 (when there are no misses and all hits are continuously at the front), and the minimum value is -1 (when there are no hits and only misses are continuously at the front). This is simply a formulation of a running sum.

The example below illustrates this concept easily.

Interpreting this concept is not particularly difficult, but if you don't grasp it well now, it may become challenging to understand further analyses involving the ES, such as the Connectivity score, in the future.

The significance is investigated using a permutation test. With this method, the labels of the traits are shuffled to create a null distribution by calculating the ES, enabling the calculation of a p-value for each gene set. The important point here is that the permutation of class labels preserves the gene-gene correlation. Therefore, it can provide a more biologically valid assessment than permuting genes.

 

Issue: R package, Original, Permutation test version

https://support.bioconductor.org/p/9148303/#9148312

Implementation of the "original" GSEA algorithm in R

 

Original paper

https://www.pnas.org/doi/10.1073/pnas.0506580102

 

정리

Gene Set Enrichment Analysis

  Gene Set Enrichment Analysis (GSEA)는 pathway 분석 기법으로, gene set이 관심있는 표현형과 얼마나 연관성이 있는지 확인하는 통계적 기법이다. 전통적인 gene set과 표현형간의 통계적 기법은 DAVID와 같이 초기하분포를 가정하는 검정부터 PAEA와 같은 새로운 방법의 검정법이 있다. 이러한 여러 분석 기법 중 GSEA는 독보적인 위치에 있다. 그 원리를 한번 확인해보자.

 

  GSEA는 굉장히 단순하며, 이해하기 쉬움으로부터 오는 해석력까지 갖추고 있다. 우선 유전자를 표현형에 관하여 줄세운다. logFC & correlation과 같은 기준을 사용함! 그리고 관심있는 유전자 집합이 그 중 어디에 위치해 있는지 살펴본다. 아래와 같이 앞쪽에 몰려있으면, gene set이 phenotype과 양의 상관성이 있을것으로 예상 가능하다. 그렇다면 이를 어찌 수치적으로 표현 가능할까?

  가장 먼저 줄을 세워보자. 관심있는 phenotype과 유전자 마다의 연관성을 기준으로 줄세우는 경우를 생각해보자.

A: 유전자를 관심있는 phenotype에 연관성이 있는 순서로 정렬한다. B: 세로의 바가 A의 정렬된 유전자와 대응되는 전체의 gene 이고, 가로의 줄 여러개가 전체 중에서 관심있는 gene들에 해당한다.

 

앞의 그림에서 바를 가로로 눞혀놓았다. 가로의 바가 전체 gene 이고 correlation with phenotype으로 정렬되어있는 상태이다. 모든 줄의 set이 결국 관심있는 gene set이다.

 

  줄 세워진 정보를 사용하여 아래 그래프에서 빨간 선으로 표시되는 random walk를 구한다. Random walk는 전체의 gene set에서 왼쪽으로부터 오른쪽으로 진행되며 구해지는데, 관심있는 gene을 만날때(hit)마다 어떤 값을 양의 방향으로 더한다.

hit 할때마다 random walk가 위로 솟고, miss 할때마다 꺼진다.

 

  그 값은 보통 표현형에 대한 유전자의 가중치를 정규화 한 값이며, 이 경우에는 전부 더하는경우 1이 되는 정규화 된 correlation with phenotype이 된다. 반면, 등장하지 않는 경우(miss), 음의 방향으로 값을 더하며, 이 값 역시 다 더하면 -1이 나온다. Miss에는 가중치를 사용하지 않고, 빈도수를 기준으로 유전자마다 균일한 값을 준다. 이렇게 계산된 random walk는 그림과같이 그려지게 된다. 

 

  GSEA에서 구하고자하는 Enrichment Score (ES) 라는 수치는, 0으로 부터 random walk 그래프의 peak까지 거리의 최댓값이 된다. 아래의 경우는 위로의 peak가 가장 큰 값을 갖으므로, 0부터 양의 방향에 있는 peak 까지의 거리가 ES 값이 된다. Hit를 앞단이나, 뒷단에서 몰아서 만나면 큰 ES가 나오게 된다. 중간에 몰려있다면, 아래 위 비슷한 0으로부터의 거리를 지닌 peak가 두개 생기게 되어 mild한 ES가 나올 것이다. 몰려있지 않고 섞여 있다면 말할것도 없이 낮은 ES가 나온다. (한번 상상해보세요!)

ES = 위쪽의 peak까지의 거리

 

  ES는 정규화 된 값을 계산하여 1과 -1 사이의 값을 갖게되며, 아래와 같이 식으로 표현 가능하다.

 

$$ P_{hit}(S,i)=\sum_{g_j\in S, j\leq i}\frac {|r_j|^p} {N_R}, where N_R=\sum_{g_j\in S}|r_j|^p $$

$$ P_{miss}(S,i)= \sum_{g_j\in S, j\leq i}\frac {1} {N-N_H} $$

 

The ES is the maximum deviation from zero of $P_{hit} – P_{miss}$.

 

  이러한 ES의 통계량은 cumulative density를 사용한 일종의 Kolmogorov-Smirnov (KS) statistic이다. 만약 $P_{hit} $에 weight를 사용하지 않고, 단순빈도를 사용한다면 original KS statistic과 값이 같아진다.

 

  그렇다면, 유의성은 어떻게 파악하는가? GSEA에서는 permutation test를 기반으로 통계적 유의성을 판단한다. 여기서 두가지 방법이 있는데, 표본의 수가 충분하지 않다면 유전자의 위치를 섞음으로 영분포를 생성한다.

 

 

반면 표본의 수가 충분하다면, class label을 섞음으로 영분포를 생성하는데, 이는 유전자간의 연관성을 깨뜨리지 않아 앞의 방식보다 더 biologically sense한 결과가 나오게 된다. 

---