일루미나 플랫폼으로 라이브러리를 만들어서 시퀀싱을 할 때에는 PhiX174 박테리오파지 DNA가 control로 들어가게 되는데, 이것이 최종 결과물까지 남아서 NCBI에 유전체 정보를 등록할 때에도 그대로 따라가는 경우가 종종 있다. 논문으로도 나올만큼 잘 알려진 문제이다.
Large-scale contamination of microbial isolate genomes by Illumina PhiX control. Standards in Genomic Sciences volume 10, Article number: 18 (2015)
ZGA 파이프라인에는 유전체 조립 결과물에서 PhiX174 박테리오파지 염기서열을 점검하는 기능이 있다. 평소에는 활성화되지 않으므로 다음과 같이 옵션을 주어야 한다.
--check-phix Check genome for presence of PhiX control sequence.
이 기능을 시험하다가 실제 명령어가 어떻게 이루어졌는지 궁금하여 zga.log 파일을 열어 보았다. 단순하에 blastn을 이용하고 있었다.
2022-03-24 10:15:29,851 - DEBUG - Running: blastn -query /opt/miniconda3/envs/zga/lib/python3.7/site-packages/zga/data/phiX174.fasta -subject /data/project/AA37_KRIBB_Edu_2022-Mar/02_data/_check_phiX/assembly/assembly.fasta -outfmt 6 sseqid pident slen length -evalue 1e-6
그런데 옵션 중 하나가 눈길을 끌었다. '-db'가 아니고 '-subject'라고? makeblastdb로 blast DB를 만들지 않고 단지 FASTA file을 제공해도 된단 말인가? Legacy blast를 즐겨 사용하던 나는 충격을 받았다. 왜 나는 지금까지 legacy blast를 쓰고 있었나? 단지 익숙하기 때문이다. 최신 blast, 즉 blast+는 아직도 옵션을 다 외우질 못해서 매번 도움말을 참조해야 한다. 명령행에서 참조를 하거나 아니면 NCBI 웹사이트를 통해서 말이다.그런데 데이터베이스를 만들지 않고도 blast를 실행할 수 있었단 말인가? 물론 legacy blast에서 bl2seq를 통해서 서열 두 개를 비교할 때에는 DB를 쓰지 않지만, 파일 두 개가 아니라 서열 두 개라는 점이 중요하다.
Blast DB를 만들어서 사용하는 이점이 당연히 있을 것이다. 검색 속도가 빠를 것임은 자명하다. 미생물 유전체 32개를 모아서 하나의 파일(all.fna)을 만든 뒤, 두 가지 방법으로 전부 blastn 검색을 해 보았다.
$ time blastn -query query -subject all.fna -out 1.out real 0m2.959s user 0m2.734s sys 0m0.225s $ makeblastdb -in all.fna -dbtype nucl -out DB -title custom_DB Building a new DB, current time: 03/24/2022 11:13:37 New DB name: /data/project/AA37_KRIBB_Edu_2022-Mar/02_data/09_ANI/fasta/DB New DB title: custom_DB Sequence type: Nucleotide Deleted existing Nucleotide BLAST database named /data/project/AA37_KRIBB_Edu_2022-Mar/02_data/09_ANI/fasta/DB Keep MBits: T Maximum file size: 1000000000B Adding sequences from FASTA; added 812 sequences in 0.993721 seconds. $ time blastn -query query -db DB -out 2.out real 0m0.192s user 0m0.149s sys 0m0.043s
DB를 만드는 쪽이 훨씬 빠르다. 너무나 당연한 이야기겠지만... 결과 파일에서는 database를 표시하는 곳에서 다음과 같은 차이가 있다. 먼저 파일 자체를 이용한 경우(-subject all.fna)는 이러하다.
Database: User specified sequence set (Input: all.fna). 812 sequences; 131,110,653 total letters
DB를 만들어서 검색을 한 경우는 이렇다.
Database: custom_DB
812 sequences; 131,110,653 total letters
makeblastdb를 실행할 때 '-title' 옵션의 값으로 제공한 DB의 제목이 여기에 표시된다. 위의 사례에서는 custom_DB였다. 검색 단계에서 DB 타이틀 명칭을 'blastn -db custom_DB'처럼 입력해서는 안 된다. 여기에 들어갈 값은 'formatdb -out' 명령에서 지정한 DB의 이름이다. DB title은 blast 검색 결과 파일에나 표시될 것이다. DB의 타이틀과 이름을 별도로 지정하는 것이 헷갈렸던 것도 최신 blast를 적극적으로 쓰지 않은 구차한 이유 중 하나다. 오늘부터는 그럴 일이 없을 것이다.
댓글 없음:
댓글 쓰기