다만 아쉬운 점이 있다면 나는 실무에서 상용 tool인 CLC Genomics Workbench를 써서 대부분의 분석 작업을 실시하지만, 교육에서는 그럴 수가 없다는 것이다. 따라서 사용이 다소 불편하더라도 공개된 소프트웨어를 이용하는 방식으로 교육을 하지 않을 수 없다. 상용 툴이 모든 것을 다 해주는 것은 아니기에 공개 툴을 조금씩 쓰기는 하지만 이를 남에게 전수하려면 더욱 많은 사전 경험이 필요하므로 실수가 없도록 철저하게 매뉴얼을 숙지하고 점검에 점검을 거듭하고 있다.
어제는 미생물 NGS 데이터를 이용한 variant calling 작업을 하다가 오류를 만났다.
$ samtools mpileup -uf ref.fa BL21.sorted.bam | bcftools view -bvcg - > var.raw.bcf
[mpileup] 1 samples in 1 input files
Set max per-file depth to 8000
[bcf_sync] incorrect number of fields (0 != 5) at 0:0
[afs] 0:0.000
빨강색으로 표시한 것이 오류 메시지이다. 구글을 열심히 뒤져 보았으나 이것과 딱 맞는 사례는 없었다. 도무지 이유를 알 수가 없었다. 필드의 수가 5개가 아니고 6개라는 오류 메시지를 보고한 경우는 매우 많았으나 나처럼 아예 아무것도 나오지 않은 사례는 보이질 않았다. 문제 해결을 위해 이것저것 한참을 알아보다가 내 시스템에 몇 가지 프로그램이 버전이 맞지 않은 상태로 뒤섞여 존재함을 알 수 있었다.
문제는 바로 de novo assembly pipeline인 A5-miseq 때문이었다. 내가 사용하는 A5-miseq에는 0.1.99-4428cd 버전의 samtools, bcftools 및 vcfutils.pl이 포함되어 있다. 그런데 별도로 설치한 samtools 1.2와 뒤섞여 실행이 되면서 오류가 난 것이다. samtools 0.1.18에는 bcftools/vcftools.pl이 같이 포함되어 있었지만 1.* 대의 버전으로 올라오면서 각자 따로 배포되고 있다. 이에 대해서는 최신 배포판 사이트를 참조하면 된다. Scaffolding tool인 SSPACE에도 bwa와 bowtie가 같이 들어있다. 이를 설치할 때에는 시스템에 먼저 깔려있던 것들과 서로 혼동되지 않도록 주의가 필요하다. 아무 생각 없이 bwa를 실행하였을 때 도대체 어느 경로에 설치된 바이너리가 돌아가고 있는지는 알고 있어야 한다는 뜻이다.
그렇다면 /usr/local/bin에 위치해서 나를 혼동스럽게 했던 samtools 1.2를 지우고 A5-miseq에 포함된 samtools와 bcftools만을 이용하여 variant calling 작업을 해 보자. 이제는 아무런 오류 없이 성공적으로 수백개의 variant가 출력되었다. 이런 성가신(?) 작업이 CLC Genomics Workbench에서는 일관적인 GUI 환경 안에서 편리하게 수행되지만 편리한 만큼 적지 않은 돈이 드는 선택이라 강요할 수는 없는 노릇이다. Command line interface 위주의 공개형 도구는 비록 사용은 불편할지라도 분석 과정의 기본을 익히기에는 좋으니 이것과 담을 쌓고 살 수는 없는 노릇이다.
Variant Calling 과정의 상세
글을 작성한 김에 변이 추출(variant calling) 과정에 대하여 좀 더 자세히 기술해 보고자 한다. 변이 추출은 유전체 resequencing의 가장 기본적인 목적이다. copy number variation이나 structural variation처럼 까다로운 것은 다음으로 미루고, SNV의 추출에 대한 것만을 정리하겠다.
매핑 도구는 bowtie2를 쓰기로 한다. 따라서 bowtie2에 포함된 공식 매뉴얼을 참조하는 것이 바람직할 것이다.
- Using SAMtools/BCFTools downstream
- Fast gapped-read alignment with Bowtie 2 (Nature Methods 2012년)
bowtie2로 만들어진 SAM 파일을 정렬 및 인덱싱을 하여 samtools pileup 명령으로 처리하는 것이 기본이다. ref.fa 파일은 미리 samtools faidx 명령으로 인덱스를 생성해 두어야 한다.
samtools mpileup -u -f ref.fa sample.sorted.bam > sample.bcf
bcftools view -v -c -g sample.bcf > sample.raw.vcf
파이프를 사용하면 한 줄의 명령어로 줄여서 실행할 수도 있다.
samtools mpileup -uf ref.fa sample.sorted.ban | bcftools view -bvcg - > sample.raw.bcf
이쯤에서 SAMtools의 공식 문서를 참조해 보자.
위 단계에서 추출한 변이를 100% 전부 믿고 받아들일 수는 없는 노릇이다. 이제 vcfutils.pl로 적절히 필터링하여 정확한 결과만을 남기도록 한다. read depth가 100 이상인 곳에 대해서만 필터링하려면 다음과 같이 한다.
bcftools view sample.raw.bcf | vcfutils.pl varFilter -D100 > sample.flt.vcfShort read aligner의 최강자는 누구인가?
여기에 대해서는 꽤 많은 논문이나 포스팅이 있으니 몇 개를 간추려서 소개하는 것으로 갈음하도록 하자.
- (쉬운 글) BWA, SOAP, MAQ, BLAT...OMG Alignment Methods! (2014년)
- Benchmarking short sequence mapping tools (BMC Bioinformatics 2013년 7월 7일)
- DNA aligner accuracy: BWA, Bowtie, Soap and SubRead tested with simulated reads (2014년 11월 4일)
댓글 없음:
댓글 쓰기