[R] ssgsea 의 scoring방법을 실습을 통해 이해해 보기

1. ssgsea 간단 설명

ssgsea는 single sample gene set enrichment analysis 의 약자입니다.

 

GSEA 라는 분석의 경우 두 그룹간의 발현값들을 비교하여 결과를 내는 방식인 반면(GSEA 계산하는 방법도 추후에 알아보고 올리도록 하겠습니다.) ssgsea의 경우 각 샘플 내에서 원하는 Gene set이 얼만큼 발현이 되었는지에 대해서 확인할 수가 있습니다.

 

계산을 어떻게 하는지 궁금해서 찾아봤는데, ssgsea의 경우 다음 연구에서 처음 사용이 되었습니다.

 

https://www.nature.com/articles/nature08460#online-methods

 

2. ssgsea의 계산 방법

 

 위에 그림에나오는 부분이 ssgsea의 식입니다. 저도 수학 전공이 아니다 보니 식을 보면 정말 난감할 때가 많습니다. 그래도 차근 하나씩 풀어보면 다음과 같이 볼 수가 있습니다.

 

위 식에서는 PwG는 gene set에서의 score 방법이고, 다음은 gene set 이외에서 score하는 방법입니다.

 

그리고 Geneset 미포함된 score의 합을 Geneset 포함 score의 합에 빼줍니다. 그러면 score가 나옵니다.

 

말로 표현하기가 힘들기 때문에 예시를 통해서 설명을 드리고자 합니다.

 

3. ssgsea 계산해보기

예시로 우리가 원하는 Geneset을 (Gene1, Gene 3, Gene5)라고 가정합니다.

 

실제로 가지고 있는 데이터는 다음과 같습니다.

 

Gene	Gene 1	Gene 2	Gene 3	Gene 4	Gene 5	Gene 6
Expression value	30	2	4	10	23	13

 

위 식에는 absolute value라고 이야기하는데 이 이유는 원래 microarray기반이기 때문에 상대값으로 -값이 나올 수가 있다고 가정한 것입니다.

 

그리고 위 논문을 읽어보시면 설명을 하는 부분에서 rank로 먼저 normalized를 진행을 하고 Empirical Cumulative Distribution Functions(ECDF)를 이용 하여 Enrichment score를 구한다라고 이야기하고 있습니다.

 

그렇기 때문에 rank normalized를 진행하면 다음과 같고, 위에서 나온 식을 이용하여 진행해보겠습니다.

 

Gene	Gene 1	Gene 2	Gene 3	Gene 4	Gene 5	Gene 6
Expression value	30	2	4	10	23	13
rank	1	6	5	4	2	3

rank로 순서대로 재 배열하였습니다. 그리고 우리가 원하는 Geneset에다 빨간색으로 표시를 하였습니다.

Gene	Gene 1	Gene 5	Gene 6	Gene 4	Gene 3	Gene 2
Expression value	30	23	13	10	4	2
rank	1	2	3	4	5	6

 

3-1. Geneset 포함 부분계산

Gene	Gene 1	Gene 5	Gene 6	Gene 4	Gene 3	Gene 2
Expression value	30	23	13	10	4	2
rank	1	2	3	4	5	6
rank normalized	6	5	4	3	2	1
Step1	6^0.25	5^0.25	0	0	2^0.25	0
Step2	6^0.25/ sum(Step1)	5^0.25/ sum(Step1)	0^0.25/ sum(Step1)	0^0.25/ sum(Step1)	2^0.25/ sum(Step1)	0^0.25/ sum(Step1)
Step2 (계산)	0.3682864	0.3518765	0	0	0.2798371	0
Step3	0.3682864	0.7201629	0.7201629	0.7201629	1	1

먼저 rank normalized를 위해서 rank가 높은 순으로 값을 크게 주도록 하겠습니다.

 

Geneset 포함된 유전자의 rank normalized값에다 0.25를 지수로 곱해줍니다. Geneset에 미포함인 경우 0으로 채워넣습니다. (Step1)

 

그리고 Step1의 합에 대한 각 유전자들의 비율을 확인합니다. (Step2)

 

그리고 누적함수로 값들을 누적시켜줍니다. (Step3)

 

마지막으로 Step3의 값들을 합쳐줍니다. 

 

그렇게 되면 4.528775이 나옵니다.

 

3-2 Geneset 미포함 부분 계산

Gene	Gene 1	Gene 5	Gene 6	Gene 4	Gene 3	Gene 2
Expression value	30	23	13	10	4	2
rank	1	2	3	4	5	6
rank normalized	6	5	4	3	2	1
Step1	0	0	1	1	0	1
Step2	0/ sum(Step1)	0/ sum(Step1)	1/ sum(Step1)	1/ sum(Step1)	0/ sum(Step1)	1/ sum(Step1)
Step2 (계산)	0	0	0.3333333	0.3333333	0	0.3333333
Step3	0	0	0.3333333	0.6666667	0.6666667	1

위와 마찬가지로 계산을 해줍니다. 포함된 부분과의 차이점은 Step1에서 값들을 1로 만든다는 점입니다.

 

결과적으로 미포함부분의 점수는 2.666667이 나옵니다.

 

그리고 포함 - 미포함을 진행해주면 됩니다.

 

값은 1.862108가 나옵니다.

 

3-3 실제로 확인해보기

 

실제로 위와 같은 값을 만들어서 구하면 같은 값이 나옵니다.

mat=matrix(c(30,2,4,10,23,13),ncol=1)
rownames(mat)=c('Gene1','Gene2','Gene3','Gene4','Gene5','Gene6')
geneset=list('test'=c('Gene1','Gene2','Gene5'))

library(gsva)
test_score=gsva(expr=mat,gset.idx.list=geneset,mx.diff=FALSE,
			method='ssgsea',verbose=TRUE,ssgsea.norm=F,min.sz=1,tau=0.25)

 

4. ssgsea 마무리

 

값들을 계산하는 부분을 공부하면서 ssgsea의 경우 geneset에 발현 rank가 높은 유전자들이 많이 포함될수록 높은 score가 만들어지고, 전체 유전자들이 많을수록 score는 작아질수 있다는 생각이 들었습니다(뭐 다들 보시면서 알아차렸을 수도 있지만요 ㅎㅎ).

 

다음 포스팅에서는 GSEA부분에 대해서 공부하고 포스팅 해보겠습니다.

 

감사합니다.