Bracken: Baysian Reestimation of Abundance with Kraken

Whole-genome shotgun sequencing 방식의 metagenomics에서 read(or contig)의 분류(또는 binning)과 더불어 중요한 목표는 군집을 이루는 각 미생물의 abundance를 추정하는 것이다. 즉 실제 미생물이 각 종에 따라서 어느 정도의 개체수를 유지하고 있는지를 알아내는 것을 의미한다. 이 수치는 직접적으로 측정하는 것은 불가능하고 단지 추정만 가능할 뿐이다.

전통적인 abundance estimation 방법으로는 유전체 안에 한 벌씩만 존재하는 마커에 대하여 서열을 query로 던져서 검색하여 추정하는 MetaPhyler,  PhyloSift 등의 방법이 있다. 이것 말고도 그 기본적인 원리를 나의 얄팍한 지식으로는 완벽하게 이해하기 어려운 방법들이 많이 존재한다.

Abundance estimation의 결과물은 각 taxon에 해당하는 퍼센트 수치이다. 그렇다면 서열 데이터를 각 taxon에 맞게 분류하면 얼마나 쓸모가 많겠는가? 이런 binning software로서 가장 최근에 공부한 것은 다음의 두 가지이다.

MaxBin: an automated binning method to recover individual genomes from metagenomes using an expectation-maximization algorithm.

PhyloPythiaS+: a self-training method for the rapid reconstruction of low-ranking taxonominic bins from metagenomes.

Read classifier로는 그동안 Kraken을 가장 즐겨 사용했었다. 하나의  short read에 하나의 taxonimic label! 매우 빠르고 깔끔한 결과가 나온다. 그러면 이러한 classifier를 사용하여  metagenomic sample의 abundance를 추정할 수는 없을까? Kraken 개발팀(존스홉킨스 대학 전산생물학 센터)에서는 최근에 Bracken이라는 도구를 발표하였다.

Bracken: estimating species abundance in metagenomics data. 프로그램 웹사이트

논문을 읽어보니 성능이 꽤 좋은 것으로 생각된다. Bracken을 사용하려면 Kraken 데이터베이스를 만들 때 사용한 라이브러리 파일(.fna)과 gi2seqid.map 파일이 필요하다. 예전에 만든 standard kraken DB는 --clean 명령을 사용해서 taxonomy 디렉토리와 database.kbd + database.idx 파일만 남기고 싹 지워버려서 Braken에 사용할 수가 없다. 따라서 다시 데이터베이스를 만드는 작업을 시작하였다.

오랜만에 standard kraken DB를 다시 만들다보니 궁금한 점이 생겼다.

1. Kraken 매뉴얼 페이지에서는 low-complexity region을 NCBI BLAST 패키지의 dustfilter로 제거하는 것을 강력히 권장하였다. 그런데 가장 적당한 dustfilter의 옵션값은 무엇인가? Readme 파일과 논문을 봐도 잘 모르겠다. 이 필터 기능이 Kraken에 포함되어 있다면 매우 좋겠지만 그렇지 못하다. Kraken의 저자 중 하나인 Salzberg에 의하면 개발팀 내부적으로는 필터를 쓰고 있다고 하였다. dustfilter의 기본 옵션을 그대로 사용하여 필터링을 하는 사례는 여기에 소개되었다. dustfilter를 거치면 소문자로 전환되므로 이를 N으로 바꾸는 트릭이 필요하다. 실제로 필터가 끝난 파일을 하나 열어보니 최소 일곱개의 염기가 반복되면 마스킹이 된다. 너무 짧은 repeat을 필터링하는 것은 아닌지?
2.  NCBI에서 자동으로 내려받는 파일 중 박테리아의 것은 all.fna.tar.gz(genome) 하나인데 바이러스의 것은 all.fna.tar.gz + all.ffn.tar.gz(gene) 두 가지이다. 왜 ffn 파일까지 받아서 DB 작성을 한단 말인가? 특별히 문제가 될 것은 없겠으나 - 모든 유전자 파일은 유전체 서열 파일에 포함될 것이므로 - 시간은 약간 낭비될 것이다. 다음은 download_genomic_library.sh 스크립트의 해당 라인이다.

 76       wget $FTP_SERVER/genomes/Viruses/all.fna.tar.gz
 77       wget $FTP_SERVER/genomes/Viruses/all.ffn.tar.gz
 78       echo -n "Unpacking..."
 79       tar zxf all.fna.tar.gz

metAMOS를 써 보려고 일년 가까이 노력하였고 이제 조건을 다 잡았다고 생각하였으나. HMP project의 mock sample(NCBI BioProject 48475)을 분석해 본 결과 metAMOS 논문에서 보여준 Figure 3이 재현되지 않았다. Deicoccocus-Thermus phylum만 100%가 나오는 해괴한 결과가 도무지 개선되지를 않았다. metAMOS를 설치하면서 Linuxbrew 및 python  업그레이드 등의 많은 기법을 익힌 것은 분명히 긍정적이었지만, benchmark test 결과를 신뢰할 수 없는 상태에서는 분석 실무에 쓰기가 곤란하다. Assembly에 기반한 abundance estimation이라는 것도 썩 마음에 들지는 않는다.