[Statistics] FDR과 Benjamini-Hochberg procedure을 통한 adjusted p-value

- 전사체 분석을 하다보면서 adjusted p-value의 개념을 공부할 일이 있어 관련 개념을 간략하게 정리

 

1. multiple comparison problem과 FDR

- 대량의 데이터 (ex) 전사체 프로파일)에서 비교대상 간 통계적 유의성 판단할 때 필요한 개념

- 비교 집단이 2개를 초과하는 경우 (ex) 약물 처리에 따른 유전자 10000개의 transcriptome data, 약물처리 후 시간에 따른 transcriptome data), 그만큼 type 1 error*가 발생할 확률이 높아짐

 * Type 1 error - 귀무가설이 참인데 그걸 기각하는 오류 = 음성인 걸 양성으로 판정하는 경우 = False positives

 

- Type 1 error 발생확률이 높아진다는 것 풀어서 설명

  ex) 전사체 분석의 예시

   - RNA-seq의 불완전성 + 동일 실험군 내 실험 샘플 간의 필연적 차이 등으로 특정 유전자의 발현량은 매우 조금씩이라도 다를 수 밖에 없음

   - technical & biological replication이 매우 많아지면, 특정 gene의 발현량은 정규분포의 형태를 띠게 될 것

   - 그러나, 보통 제한된 수로 실험을 하기때문에 매우 낮은 확률로 같은 특성을 지닌 개체들의 유전자 발현량이 다르다고 나올 수 있음 -> 보통 그 기준을 p-value=0.05로 잡음 -> 즉, 하나의 gene에 대해서 5% 확률로 false positive가 발생할 수 있다는 것

   - 이게 한 번의 통계검정이면 5%가 큰 문제가 아닐텐데..

    gene 갯수는 eukaryote 등에서 10000개 이상 넘어가는 큰 숫자이기 때문에,

    결국 그만큼 모든 gene들에 대한 통계 결과가 true positive만 나올 확률은

    (0.95) * (0.95) * ... * (0.95)로 전체 gene 갯수 n이 커짐에 따라 점점 더 작아짐 (반대로 false positive 발생 확률은 커지게 됨)

    - 전사체 분석을 하나의 큰 통계검정으로 본다면 10000개 유전자 중 500개는 false positive이라는 뜻이기도 한데, 이는 매우 위험한 결과 해석으로 이어질 수 있음 

  -> 즉, Type 1 error 발생 증가의 결과인 False positive를 쳐내고 True positive만을 남기기 위해선 기존 p-value가 아닌 특정한 보정을 거친 adjusted p-value가 필요함 

 

- Type 1 error 보정 방법 (adjusted p-value)

 : Type 1 error 보정의 conceptualization 방법에는 Familywise error rate, False discovery rate가 있음

 : 전자는 Bonferroni correction으로 이를 보정하고, 후자는 Benjamini-Hochberg procedure(B-H method)을 사용

 : False positives를 쳐내려는 과정에서 필연적으로 True positive까지 손실될 위험(Type 2 error 증가, false negative의 증가)이 있는데, 위 2방법 중 B-H method가 손실의 정도가 더 적기 때문에 보통 이를 많이 활용함

 -> Bonferroni correction = (각 p-value 결과값) / (전체 p-value 갯수)
 -> 유전자 10000개면은 새로운 adjusted p-value cut-off = 0.05 / 10000 = 0.000005 -> 남는 true positive가 너무 없어짐 (=too conservative)

2. FDR 감소를 위한 B-H method 원리와 adjusted p-value 구하는 방법

 - FDR: False positive / Total positive (Total = False + True)

 - FDR=0.05 -> 10000개 중 500개는 거짓이라는 뜻 -> 앞서 말했듯 유전자 500개가 잘못 해석되는 불상사 발생 가능

 - 이를 보정해주기 위한 B-H method

 

 2.1. B-H method 원리

 : 그림 1처럼 유전자 중 발현량 차이없는 것들은 위 그래프처럼 동일한 p-value의 분포를 보일 것이고,

   유전자 발현량이 차이가 나는 것들은 0에 치우쳐서 분포할 것이다

그림1. 유의성 여부에 따라 다른 p-value의 분포 (ref: StatQuest 유투브영상 중 발췌)

 

 : 그러면 전체 유전자에 대한 p-value 분포는 각 bin별 p-value의 합을 y값으로 가지는 그래프가 될 것이다 (그림 2)

그림 2. 전체 유전자에 대한 p-value 분포 그림은 각 구간별로 p-value의 갯수 값을 더한 형태이다 (ref: StatQuest 유투브영상 중 발췌)

 

 - 그리고, uniform distribution의 y값 평균선을 확인할 수 있고, 이를 통해 "True positive"가 어느 정도 분포하는지 눈대중으로 알 수 있다 (eyeball method, 그림 3)

그림 3. uniform distribution의 y값 평균선을 통한 true positive의 대략적인 도출 (eyeball method) (ref: StatQuest 유투브영상 중 발췌)

 - 이런 원리를 기반으로 B-H method는 "유의성이 있다고 나온 값들 중 false positive를 빼내는 방식"으로써, 이 때 그 기준을 눈대중으로 하는 대신 FDR의 threshold를 두고 거기에 맞게 샘플들을 선별하고, 이 정보를 담고 있는 새로운 p-value를 제시해줌

  ex) FDR<0.05로 기준을 잡을 경우,

     기존 0.05 이하의 p-value들 중에서 false positive의 비율이 5% 이내가 되도록 일부만 빼냄 (그림 4)

     비록 5%의 false positive가 있지만 95%의 true positive를 남기면서 꽤 많은 수의 false positive를 제거할 수 있기 때문에 의미가 있음

   

그림 4. B-H method를 통해 true positive를 많이 살려두면서 false positive를 많이 줄일 수 있다 (ref: StatQuest 유투브영상 중 발췌)

 

 

2.2. B-H method를 통한 adjusted p-value 계산 방법

 - 생각보다 매우 간단

 - 모든 p-value를 오름차순으로 정렬 -> 다음의 식 계산 (아래 식, 그림 5) (보면 알겠지만, 엑셀로 p-value 리스트 만들어서 정렬해놓고 쉽게 계산해서 얻을 수 있을 것이다 ㅎㅎ)

 

adjusted p-value = p-value값 * (전체 p-value 갯수 / 현재 보정 진행하는 p-value의 오름차순 시 순서)

 

그림 5.  B-H method에 의한 adjusted p-value 계산 방법 (ref: StatQuest 유투브영상 중 발췌)

3. 요약

- 수많은 유전자들의 발현량을 검증하는 경우 type 1 error가 늘어나고, 이를 회피해주기 위해서 FDR threshold을 이용한 B-H method로 p-value를 보정해준다

- B-H method에서, 각 p-value를 오름차순으로 정렬 후 p-value * (전체 갯수/해당 p-value의 순서값)으로 adjusted p-value를 얻을 수 있다

- 여러분 StateQuest 채널 구독하세요 두번하세요 진짜짱짱임.. (오늘의 진짜 결론)

References
[1] Youtube, StateQuest with Josh Starmer, False Discovery Rates,FDR,clearly explained

 

Copyright 2021. komok’s sight All Rights Reserved.