Metaphyler 결과물로 krona chart 그리기

kraken read classifier의 결과물을 krona chart로 그리는 방법은 이미 잘 알고 있다. 그 과정을 복습해 보자. kraken의 출력 포맷은 매뉴얼 웹 페이지를 보면 나온다. 결과 파일의 두번째와 세번째 컬럼, 즉 sequence ID와 taxonomy ID를 뽑아서 새로운 파일을 하나 만든다.

cut -f2,3 kraken-results > kraken-results.krona

값을 지니는 옵션의 경우 이렇게 붙여 쓸 수도 있구나! 간혹 긴 옵션을 --longOption=value처럼 표준 형태로 쓰지 않고 --longOption value 또는 -longOption value처럼 비표준적인 방법으로 표기하는 명령어들이 있어서 가끔 짜증이 난다.

그 다음에는 이 명령어 한 줄이면 끝난다.

ktImportTaxonomy kraken-results.krona -o kraken-results.krona.html

그러면 Metaphyler의 결과 파일(*.classify.tab)을 보자. 어떻게 구성되어 있을까?

OUTPUT:
       "prefix".classify.tab
           column 1: query sequence id
           column 2: phylogenetic marker gene name
           column 3: best reference gene hit
           column 4: % similarity with best hit
           column 5: classification rule
           column 6-10: taxonomic label at genus,family,order,class,phylum level

아이쿠, 눈으로 보기에는 좋은데 NCBI tax ID는 직접적으로 보여주고 있지를 못하다. 두번째 컬럼이 tax ID를 탐색할 수 있는 유일한 키가 되는데, 유전체의 accession이 아니라 Metaphyler를 구성한 phylogenetic marker 정보이다. 실제의 형식은 NZ_AAXV01000023_61651_61292이다. 즉 genomeAccession_start_stop인 것이다. 가장 근본적인 방법은 이 컬럼을 조작하여 genome accession을 뽑은 뒤 NCBI의 efetch를 이용하거나 taxonomy 파일을 전부 다운로드한 뒤 BioPerl 등을 이용하여 accession에서 tax id를 추출하는 스크립트를 짜는 것이다.

그러나 Metaphyler가 이미 계층적으로 분류한 결과물을 제공하고 있는데, 뭐하러 이런 미련한 방법을 쓰겠는가? 계층적 구조를 갖는 텍스트가 준비된 상태라면 당연히 krona tool은 이를 파싱하여 챠트를 그려낼 것이다. 그런 기능을 하는 명령어는 ktImportText이다. 사용법은 여기를 참조하라. 이제 해결 방법을 찾았다. Metaphyler의 *.classify tab에서 6번 - 10번 컬럼을 역순으로 추출하여 파일을 작성한 뒤, ktImportText에 공급하면 그만이다.

awk 'OFS="\t"{print 1, $10, $9, $8, $7, $6}' classify.tab > test.result
ktImportText test.result -o test.html

나머지 컬럼의 값이 같은 것을 전부 세어서 첫번째 컬럼의 숫자로 삼으면 test.result 파일이 훨씬 간단해질 것이다. 자, 그럼 firefox를 띄워 보자. 휼륭하다!

내가 쓸 줄 아는 스크립트는 Perl이 전부이니 언젠가는 Bio::DB::Taxonomy를 이용한 간단한 스크립트를 작성해 봐야 한다. 맨날 쓰는 것이라고는 Bio::Seq(IO)가 전부이니 원...