Bracken: estimating species abundance in metagenomic data

Kraken은 매우 빠르고 능률적인 read classifier이다. 즉 read 하나에 대해서 taxonomic label을 부여하는 것이 고유 기능이다. 그러나 metagenomics 연구에서 가장 근본적인 질문 중 하나는 샘플 안에 존재하는 미생물이 각 종별로 얼마나 존재하는가, 즉 abundance를 추정하는 것이다. 이를 위해서는 각 종의 genome size 등 고려할 것들이 많다. 미국 존스홉킨스 대학에서는 Kraken의 성능을 그대로 이용하되 베이즈 추론을 사용하여 species abundance를 추정하는 새로운 프로그램인 Bracken(Baysian Reestimation of Abundance after Classification with KrakEN)을 2017년 1월에 발표하였기에 이를 소개하고자 한다.

논문: https://peerj.com/articles/cs-104/ (PeerJ Comput Sci)
프로그램 홈페이지: http://ccb.jhu.edu/software/bracken/index.shtml

사용 방법

Bracken 웹사이트의 매뉴얼이 약간 혼동을 초래하는 면이 없지 않아서, 프로그램 패키지에 포함된 README 파일을 기준으로 사용 방법을 정리한다. Bracken analysis는 다음의 네 단계로 이루어진다.

DNA read를 Kraken으로 분류하기
Kraken database(library, 즉 .fna file)를 분류하기
Kmer 분포를 계산하기
Abundance estimation

Bracken 분석에 필요한 Kraken 결과물은 kraken-report가 생성하는 파일이다. kraken 명령어가 만들어내는 각 read마다 한줄씩 classification 정보가 들어있는 파일이 아님에 유의하라. 또한 두번째 과정을 수행하려면 Kraken DB 작성에 투입한 FASTA genome sequence file이 남아있어야 한다. 따라서 직접 DB를 만든 뒤 해당 디렉토리를 청소하면 안된다.

Step 1. Kraken classification of DNA reads

Kraken에 대한 상세한 설명 자료는 소개 및 매뉴얼을 참조한다.

$ KRAKEN_DB=/usr/local/apps/kaken/standard_kraken # 각자 상황에 맞게
$ kraken --db $KRAKEN_DB --threads=16 READS.fastq > READS.fastq.kraken
$ kraken-report --db $KRAKEN_DB READS.fastq.kraken > READS.fastq.kraken.report

Step 2. Classifying a full kraken database

라이브러리는 ${KRAKEN_DB}/library 디렉토리 아래에 .fna라는 확장자로 존재해야 한다. 다음의 사례는 .fa와 .fasta 파일도 작업 대상에 포함시킨다.

$ kraken --db=${KRAKEN_DB} --fasta-input --threads=10 <( find -L library -name "*.fna" -o -name "*.fa" -o -name "*.fasta" -exec cat {} + )  > database.kraken
$ perl count-kmer-abundances.pl --db=${KRAKEN_DB} --read-length=75 database.kraken  > database75mers.kraken_cnts

Step 3. Generation of Kmer distribution file

$ python generate_kmer_distribution.py -i database75mers.kraken_cnts -o KMER_DISTR.TXT

Step 4. Abundance estimation

이제 최종 단계만 남았다. 필요한 파일은 query reads에 대한 Kraken report, 그리고 바로 직전에 만든 KMER_DISTR.TXT 파일이다.

$ python estimate_abundance.py -i READS.fastq.kraken.report -k KMER_DISTR.TXT -l CLASSIFICATON_LEVEL -t THRESHOLD -o OUTPUT_FILE.TXT

CLASSIFICATION_LEVEL: Default = 'S' for Species. K (kingdom level), P (phylum), C (class), O (order), F (family), and G (genus).
THRESHOULD: Default = 10. Species classification을 할 때 이 수치 이하의 read를 가진 종은 abundance estimation 과정에서 상위 taxonomy level의 read를 추가적으로 배정받지 못한다.

이 명령을 실행하는 동안 표준 에러로 다음과 같은 메시지가 출력된다.

PROGRAM START TIME: 03-29-2017 23:50:00
BRACKEN SUMMARY (Kraken report: report)
    >>> Threshold: 10 
    >>> Number of genuses in sample: 332 
        >> Number of genuses with reads > threshold: 92 
        >> Number of genuses with reads < threshold: 240 
    >>> Total reads in sample: 6562065
        >> Total reads kept at genuses level (reads > threshold): 5904658
        >> Total reads discarded (genuses reads < threshold): 615
        >> Reads distributed: 101669
        >> Reads not distributed (eg. no genuses above threshold): 487
        >> Unclassified reads: 554636
BRACKEN OUTPUT PRODUCED: OUTPUT_FILE.txt
PROGRAM END TIME: 03-29-2017 23:50:00

Genome Informatics Laboratory at KRIBB

Table of Contents