clusterProfiler를 사용한 KEGG enrichment analysis

전통적인 microarray 실험이나 RNA-seq data 분석을 통해서 얻어지는 가장 중요한 중간 자료는 차등발현유전자(DEG, differentially expressed gene)이다. 가령 대장균 배양체에 열을 갑자기 가한 뒤 발현된 유전자를 찾아서 그 기능을 알아보면 heatshock protein gene이 많은 비율을 차지할 것이다. 그런데 여기에서 중요한 것은 대장균이 원래 갖고 있는 유전자 중 heatshock 관련 유전자의 수 자체가 얼마인지를 아는 것이다. 고등학교 시절에 배운 통계와 확률 수준에서 이야기한다면, 흰 공과 검은 공이 각각 50개씩 섞인 주머니에서 5개의 공을 꺼냈는데 흰 공이 3개, 검은 공이 2개 나왔다면 이것이 정말 무작위로 꺼낸 것인지 아니면 흰 공이 편향되게 많이 나온 것인지를 결정해야 하는 문제와 비슷한 것이다. 즉 초기하분포(hypergeometric distribution)의 영역이 된다.

물론 내가 아주 정확한 예를 든 것은 아니다. 모집단과 샘플의 크기, 그리고 복원 추출이냐 아니냐 등의 문제를 언급하지 않았기 때문이다. 일단 대충 넘어가도록 하자^^

다시 (gen)omics의 영역으로 넘어오자. 연구 대상이 되는 생명체의 전체 유전자에 대한 functional classification이 아주 잘 되어있는 모델 생명체는 이런 분석을 하기에 매우 적합하다. 단지 gene ID만 있으면 DEG set에 어떤 기능이 특출나게 더 존재하는지를 알려주는 도구가 많이 있기 때문이다. 여기서는 대조군에 대해서 발현량이 더 많거나 적은지 자체가 중요하지 않다. 오직 '차이'가 있는 유전자 전체를 뭉뜽그려서 하는 이야기이다.

그런데 세균 유전체를 다루는 나 같은 사람에게는 이 분석 작업이 오히려 불편하다. 전체 gene set에 대한 functional classification을 스스로 해결해야 하기 때문이다. Blast2GO가 이런 용도로는 매우 편리하다. 그리고 일단 각 유전자에 대한 GO term이 확정되면, Blast2GO(How to perform a gene set enrichment analysis with Blast2GO 링크) 혹은 CLC Genomics Workbench에서 DEG set이 더 많이 갖고 있는 기능들에 대한 통계적 분석을 할 수 있다.

오늘 발견한 R package는 clusterProfiler라는 것이다. 잘 알려진 모델 생명체가 아니더라도 KEGG에 등록된 유전체라면 functional enrichment analysis가 가능하게 만드는 패키지이다. 참조한 사례는 R-Bloggers.com에 패키지 저자들이 올린 포스팅이다.

KEGG enrichment analysis with latest online data using clusterProfiler

업무용 컴퓨터에 설치된 기존 R 버전(3.1.2)에서는 지원이 되지 않아서 최신 버전(3.3.2)를 깔고 DOSE와 clusterProfiler 패키지를 설치하여 사용하였다. 너무 오랫만에 R을 쓰려니 텍스트 파일에 수록된 gene list를 scan() 함수를 입력하는 것도 힘들다. 입력하려는 데이터가 숫자가 아니라고 불평을 하는 것 아닌가? 해결책은 다음과 같았다.

> gene = scan("gene_list.txt", what="string")

나오는 결과는 대충 다음과 같다.

패키지의 저자인 Guangchuang Yu와 Li-Gen Wang에게 감사를! 남은 숙제는 갖가지 기능으로 장식된 Blast2GO Pro 최신 버전을 더욱 알뜰하게 활용하는 것이다.