Custom kraken DB test

Kraken에서 database란 아주 구체적으로 말하자면 디렉토리 명칭이다. DB 작성을 위해 모은 염기서열 fasta file 모음은 library라 부른다.Kraken DB가 수록한 내용은 library에서 추출한 모든 k-mer에 대해서 taxonomy 정보를 연결해 놓은 것이라고 생각하면 간단하다. 만약 어떤 k-mer가 여러 genome에 공통적으로 존재한다면 어떤 taxonomy 정보를 연결할 것인가? 이 경우에는 lowest common ancestor(LCA)를 연결한다. 보다 상세한 내용은 Kraken 공식 매뉴얼을 참조하라.

Kraken DB를 만들려면 taxonomy 정보가 포함된 서열 자료(fasta files)와 NCBI taxonomy 파일이 필요하다. Standard DB란, NCBI에서 제공하는 bacteria 및 virus의 유전체 서열을 이용한 것이다. 이들 파일에는 서열 ID에 gi 번호가 있어서 taxonomy와 연결되는 키 역할을 한다. 실제로 라이브러리 파일 하나를 열어보면 다음과 같다.

>gi|255767013|ref|NC_000964.3| Bacillus subtilis subsp. subtilis str. 168 chromosome, complete genome

단, 이제는 NCBI에서 bacteria 유전체 서열 전부를 압축한 all.fna.tar.gz를 더 이상 업데이트하지 않고 있으며 보관 위치도 변경되었다. 따라서 Kraken package의 download_genomic_library.sh를 다음과 같이 일부 수정해야 한다. 다행스럽게도 virus에 대해서는 아직 변경되지 않은 위치에서 all.fna.tar.gz를 제공한다. bacteria 라이브러리는 2786건의 유전체(업데이트되지 않음), 그리고 viruses 라이브러리는 2017년 3월 28일 현재 4391 건의 유전체가 존재한다. plamids library는 download_genomic_library.sh에 따르면 $FTP_SERVER/genomes/Plasmids/plasmids.all.fna.tar.gz를 받는 것으로 되어있지만, 이 파일은 현재 보이지 않는다. NCBI ftp 트리를 뒤져보니 ftp://ftp.ncbi.nlm.nih.gov/genomes/archive/old_refseq/Plasmids/plasmids.all.fna.tar.gz가 존재한다.

FTP_SERVER="ftp://ftp.ncbi.bln.nih.gov"
(수정 전) wget $FTP_SERVER/genomes/Bacteria/all.fna.tar.gz
(수정 후) wget $FTP_SERVER/genomes/archive/old_refseq/Bacteria/all.fna.tar.gz

현재의 GitHub 사이트에서 배포하는 소스에는 Bacteria용 download URL이 수정된되었음을 확인하였다. — Haeyoung Jeong 2016/12/13 19:33

Bacteria/all.fna.tar.gz는 5242개의 complete sequence로서 복수의 염색체가 존재하는 경우 별도의 .fna 파일로 분리된 상태이다. kraken-build –standard라고 실행을 하면 fasta file과 taxonomy 파일의 다운로드부터 이루어진다. 실제 standard kraken library를 만드는 명령어를 다음에 소개하였다. 사용한 서버는 256 GB의 메모리가 장착되었으며 16개의 thread를 사용하였다(Intel(R) Xeon(R) CPU E5-2630 v3 @ 2.40GHz x 2ea). 여기에서는 kraken 프로그램이 설치된 디렉토리에서 모든 작업을 하는 것으로 가정하였고, 목표가 되는 DB의 이름은 kraken-standard이다. kraken 설치 디렉토리를 $PATH 환경변수에 선언하고 임의의 위치에서 실행해도 된다. jellyfish 1.x의 실행파일이 $PATH에 반드시 있어야 한다. kraken database의 이름(아래 사례에서는 kraken-standard)은 실제적으로는 디렉토리 명칭으로서, kraken-build를 실행하면 새로 만들어진다. Braken을 할 계획이라면 –clean을 사용하지 말기 바란다. 라이브러리 파일(.fna)과 gi2seqid.map 파일이 필요하기 때문이다.

$ PATH=/home/test/metAMOS-1.5rc3/Utilities/cpp/Linux-x86_64/jellyfish/bin:$PATH
$ ./kraken-build --standard --threads 16 --db kraken-standard
(다운로드 및 압축 해제 과정 진행)
Kraken build set to minimize disk writes.
Creating k-mer set (step 1 of 6)...
Found jellyfish v1.1.6
Hash size not specified, using '11635839817'
K-mer set created. [15m55.827s]
Skipping step 2, no database reduction requested.
Sorting k-mer set (step 3 of 6)...
K-mer set sorted. [1h3m36.646s]
Creating GI number to seqID map (step 4 of 6)...
GI number to seqID map created. [22m51.293s]
Creating seqID to taxID map (step 5 of 6)...
215968 sequences mapped to taxa. [1m10.253s]
Setting LCAs in database (step 6 of 6)...
Finished processing 216310 sequences                   
Database LCAs set. [27m52.924s]
Database construction complete. [Total: 2h11m26.988s]
$ ls -lt
합계 150000908
-rw-rw-r-- 1 hyjeong hyjeong           0 2016-12-13 10:11 lca.complete
-rw-rw-r-- 1 hyjeong hyjeong 71974784200 2016-12-13 10:11 database.kdb
-rw-rw-r-- 1 hyjeong hyjeong    10601108 2016-12-13 09:43 seqid2taxid.map
-rw-rw-r-- 1 hyjeong hyjeong    11149110 2016-12-13 09:42 gi2seqid.map
-rw-rw-r-- 1 hyjeong hyjeong  8589934608 2016-12-13 08:18 database.idx
-rw-rw-r-- 1 hyjeong hyjeong 71974784200 2016-12-13 08:16 database.jdb
drwxrwxr-x 4 hyjeong hyjeong          35 2016-12-13 07:59 library
drwxrwxr-x 2 hyjeong hyjeong        4096 2016-12-13 07:54 taxonomy
$ kraken-build --clean --db kraken-standard
$ ls -lt
합계 79691780
drwxrwxr-x 2 hyjeong hyjeong          38 2016-12-13 10:15 taxonomy
-rw-rw-r-- 1 hyjeong hyjeong 71974784200 2016-12-13 10:11 database.kdb
-rw-rw-r-- 1 hyjeong hyjeong  8589934608 2016-12-13 08:18 database.idx

LCA setting 단계에서 어떻게 215968개의 서열이 처리되었는지는 잘 모르겠다. library에 존재하는 개별 fasta file의 수와는 분명히 다르다. 작업이 끝난 DB는 다른 곳으로 복사하여 사용해도 된다. 아니면 처음부터 원하는 디렉토리를 생성하여 적절한 환경변수(예: $DBNAME)으로 선언하여 작업을 해도 좋다. 2015년 2월 16 thread, 122 GB RAM이 설치된 컴퓨터로 standard library를 구축하는데 걸리는 시간은 다음과 같다(출처: Kraken manual site). 추출한 k-mer set를 정렬하고 LCA value를 설정하는데 가장 많은 시간이 걸린다.

2016년 9월을 기하여 NCBI가 제공하는 GenBank, GenPept 및 FASTA 포맷에서는 더 이상 GI 번호를 쓰지 않게 되었다(공지문). 이를 해결하기 위한 방안이 opinomics 사이트에 게시되었다.

$ kraken --db $DBNAME sequences.fa > sequences.kraken
$ kraken-translate --db $DBNAME sequences.kraken > sequences.labels

3520.notrim.pe.fq - 3119470042 bp / 30885842 seqs; 101.0 average length

30885842 sequences (3119.47 Mbp) processed in 601.197s (3082.4 Kseq/m, 311.33 Mbp/m).
1063743 sequences classified (3.44%)
29822099 sequences unclassified (96.56%)

30885842 sequences (3119.47 Mbp) processed in 559.800s (3310.4 Kseq/m, 334.35 Mbp/m).
501601 sequences classified (1.62%)
30384241 sequences unclassified (98.38%)

custom DB를 사용한 것이 더 나은 read coverage를 보였다. 최소한 kraken DB 자체에는 문제가 없다는 뜻이다. 그런데 왜 metAMOS test script에서는 제대로 작동을 하지 않는 것일까?