R로 q-value 구하기

   R로 q-value를 구하기 위해서는 이미 구현되어 있는 package 를 사용할 수 있다. false discovery rate 을 구하기 위한 R package들은 이 곳(strimmer 랩의 FDR에 관한 R 패키지 페이지)에 잘 나와있다. 이 글은 내가 주로 사용하는 package를 중심으로 설명한다.

qvalue : John Storey가 제시한 알고리즘을 구현해 놓은 패키지. feature간의 dependence를 고려하지 않아도 될 때. GUI 가 제공된다.

fdrtool : p-value 뿐만이 아니라 t-score, z-score나 correlation 으로부터 fdr 을 계산할 수 있다.

multtest : bioconductor 하부에 있는 것으로, 꽤 여러 방법이 구현되어 있다. 나는 dependence 의 경우에 사용하기 위해 이것을 사용한다.


   q-value는 간단히 말하면 multiple test를 할 때 false discovery (false positive, Type I error)를 고려한 p-value로 생각할 수 있다. 일반적인 통계에서는 한 두 개의 가설을 검정하기 때문에 p-value 만으로 충분하지만 microarray 와 같이 1만개 이상의 가설을 동시에 검정할 경우, type I error 가 0.05 라 하더라도 잘못된 rejection^[각주:1]을 하는 것은 너무 많은 오류를 범하는 것이기 때문에 많이 사용하는 p-value인 0.05 나 0.01 도 너무 크게 되며, 따라서 이 값을 보정할 필요가 있게 된다. 이런 상황에 맞게 p-value를 보정한 값이 q-value. 보다 자세한 내용은 나중에 별도의 글로 작성하고, 이 글은 q-value를 구하는 방법을 살펴 본다.

QVALUE 사용

   R을 설치한 후, 적당한 패키지를 설치한다. 이 때 패키지 미러 싸이트로 미국 쪽을 택한다^[각주:2]. qvalue는 GUI 까지 제공되기 때문에 사용하기 매우 편리하다. 간단히 그 사용 예를 살펴 보면 다음과 같다.

일단 R을 실행시키고, qvalue library를 부른 후 gui를 실행시킨다.

그러면 다음과 같은 창이 나타난다.

위의 순서대로 실행을 하면 된다. Browse 로 파일을 찾고(1), Load로 파일 내용을 불러 들인다(2). 그 후, true null의 p-value의 분포에 대한 추정을 하는 방법을 선택한다(3). 주로 interpolation 을 하는 방법을 사용하는데, 원론적으로는 bootstrap 으로 하는 것이 더 좋다고 개인적으로 생각하기 때문에 난 주로 bootstrap 으로 한다. 그 후 robust method 를 선택할 수도 있고(4), 아니라면 이제 Execute를 눌러서 q-value를 계산한다(5). 계산된 q-value는 Save Output 으로 저장할 수 있다(6). p-value나 q-value의 분포 등은 Plots에 있는 메뉴를 통해 적당히 볼 수 있다. Q-plots만을 예로 보자면 다음과 같다.

pi_0^[각주:3]의 추정치가 0.27인 것으로 보아 충분히 다른 두 sample 이라 할 수 있어 보인다.




fdrtool 사용

   fdrtool 은 gui가 없긴 한데, 어쨌든 사용은 다음과 같이 할 수 있다.

> library(fdrtool)
> p <- scan("D:\\adnoctum_desktop\\pvalue\\GSE5364.pvalue.txt", what=numeric(0))
Read 6685 items
> pval.estimate.eta0(p, method="bootstrap")
[1] 0.2722513

이 때, GSE5364.pvalue.txt 파일은 p-value 값이 한 줄에 하나씩 나열되어 있는 파일이다. 위까지 하면 출력된 0.27 값이 역시 pi_0의 추정치로 qvalue와 거의 유사한 것을 알 수 있다. 또한 fdrtool 은 이렇게 했을 때 p-value의 distribution 을 그려 준다. 그 후,

> q = fdrtool(p, statistic="pvalue");
Step 1... determine cutoff point
Step 2... estimate parameters of null distribution and eta0
Step 3... compute p-values and estimate empirical PDF/CDF
Step 4... compute q-values and local fdr
Step 5... prepare for plotting

위처럼 fdrtool 로 적당한 option 을 입력하여 실행시키면 값이 계산되며, 객체가 하나 반환되는데, 그 객체의 qval 값이 q-value 값이다. 따라서 이것을 다시 저장해 보자.

> write(q$qval, "GSE5364.qvalue.txt", sep="\n");

multtest 사용
   multtest 는 꽤 여러 알고리즘이 구현되어 있으며, 특히 dependence 의 경우에도 사용할 수 있기 때문에 사용하게 되었다. 실제 사용 예는 다음과 같다.

> library(multtest);
... // 생략.
> q <- mt.rawp2adjp(p, "TSBH");
> length(q$index);
[1] 6685
> write(q$adjp[,2], "D:\\adnoctum_desktop\\test.q.txt", sep="\t");

manual 에 보면 여러 기능이 제공되는데 우리에게 필요한 것은 p-value를 보정하는 것이므로 rawp2adjp 를 사용하면 된다. 또한, 저장된 q-value는 p-value의 오름차순으로 정렬된 값에 대한 q-value가 나열되어 있는 것이라고 한다. 즉, 반환된 값 q$adjp 는 N * 2 의 matrix로, 첫 번째 열이 원래의 p-value, 두 번째 열의 값이 그 q-value라고 manual에 나와 있다. 현재 p-value가 6,685개인데, two-step 단계 중 첫 번째 단계에서 찾아진 true null의 개수는

> q$h0.TSBH
h0.TSBH_0.05
        3825

3,825개 라고 한다. 따라서 이 경우 pi_0의 추정치가 대략 0.57 정도 된다고 할 수 있다.


각 패키지의 보다 자세할 설명은 각 패키지의 메뉴얼을 참고한다.

1. true null을 reject 하는 것을 우리는 discovery 라 하게 되는데, 따라서 이런 것은 false positive. [본문으로]
2. 많은 경우 미러 싸이트는 한국보다는 일본이나 미국을 선택하는 편이 좋다. 심지어 일본을 선택하는 것이 더 빠를 때도 있다. [본문으로]
3. PNAS 논문에 있는데, 이 값은 전체 중에 true null의 비율이다. 즉, 유전자 경우에 대해 얘기하자면, 정상과 암에서 다르지 않다고 할 수 있는 유전자의 비율은 대략 27% 라는 얘기. [본문으로]