바이오인포매틱스

Short Read Alignment 그것을 알려주마

hongiiv 2012. 5. 3. 16:03
반응형
지난주 Boston에서 열린 Bio-IT World Conference & Expo를 보고 나서 느낀게 있다면 요즘 이바닥은 점점 상용화쪽으로 흘러가고 있다는 것이다. Alignment Software 하나를 보더라도 예전에는 논문을 쓰고 학술적인것에 의미를 두었다면, 이제는 좀 더 빠르고 정교하게 (정교하다는 표현이 맞을지는 모르겠지만, 다른 Alignment Software가 놓칠 수 있는 Variants까지 잡아낼 수 있는) 만들어 이를 상업적으로 내놓고 있다는 것이다.

그렇다고 무작정 NGS를 하고 난뒤 이러한 상용의 서비스에 맡겨서 분석하기에는 가격이 만만치 않다. 아직 우리에게는 BWA에 있지 않은가! BWA라도 잘 알고 사용한다면 분석에 있어서 좀 더 의미있는 결과를 낼 수 있지 않을까하는 작은 바램에서 본 글을 작성해 본다. 물론 아무것도 모르는 초딩도 이해 할 수 있을 수준의 글이니 난 그런거 몰라도 돼라고 생각되는 분들도 한번쯤 읽어 본다면 유익한 글이 되지 않을까 생각되는 닥치고 한번 지루하더라도 읽어보길 권한다.

Benjamin Franklin Award for Open Access in the Life Sciences

Bio-IT World Conference에서는 Benjamin Franklin Award for Open Access in the Life Sciences의 수상식까지 겸하고 있으며, 매년 생명공학 분야에서 Open Access 즉, Open Source Software나 공공DB를 만드는데에 공헌한 분들에게 주어지는 상이다. 이제까지의 수상자를 보면 James Kent, Licoln D. Stein, Ewan Birney(Ensembl, BioPerl), Michael Ashburner(FlyBase), Rovert Gentleman(BioConductor), Philip E. Bourne(PDB), Alex Bateman(Pfam) 등으로 올해 12년도는 SAMtools, BWA, MAQ 등을 만든 (Sanger에서 포닥하면서...) Heng Li로 현재는 Broad에서 David Reich, David Altshuler와 함께 일하고 있다. - 이분들 이름만 들어도 설레네...

수상식 당일 Heng Li의 모습 

 자 그럼 2009년도 Sanger 시절 발표한 BWA (요즘 제일 잘 나가는 Short Read를 위한 Aligning Software)에 대해서 짧지만 임팩트 있게 살펴보도록 하겠다.

Short Read Alignment using BWA

Sequencing을 나누자면 여러가지로 나눌 수 있겠지만, NGS 전/후 즉, Short Read 이전 시대와 NGS와 Short Read 시대로 나누어 볼 수 있겠다. 그동안의 이분야에서는 BLAST로 대변되는 시대였다면 NGS와 함께 엄청한 Short 다리들이 생산되기 시작하면서 Short 다리를 위한 Aligning Software가 나오기 시작했고 이들을 어느정도 평정한 것이 바로 BWA이다.  Burrows-Wheeler Aligner는 이름에서도 알 수 있듯이 버로스와 휠러가 1994년 캘리포니아 건포도를 씹으면서 팔로알토의 DEC Systems Research Center에서 만든 data compression에 사용되는 알고리즘으로 block-sorting compression이라고 불리는 Burrows-Wheeler Transform(BWT)를 기반으로 하고 있다. - 뭐 이런건 몰라도 되지만 이런거 찾는게 더 재미있다. -.-;;

Heng Li 흉아는 숏다리(200bp 이하의)와 롱다리(200bp~100kbp)를 위해 두개의 Software를 친절히도 만들어 하나는 우리가 흔히 부르는 BWA로 알고 있는 BWA와 (-.-;;), BWA-SW라 칭하셨으며, 두 놈다 앞서 언급한 BWT 알고리즘을 기반으로 하고 있다.

BWA는 다시 두개의 part로 나뉘어졌는데 1) bwa aln 과 2) bwa samse/sampe 가 바로 그것이다. 덴장 Align 하는데 무슨 두개의 part로 나누고 지랄이삼이라고 칭얼거릴 수 도 있겠지만은 헹리흉아 나름 깊은 뜻이 있는 것이니 그 깊은 뜻 조금이나마 이해하려면 앞서 계속 걍 넘어가려고 했던 알고리즘 부분을....하지만 난 초간단히 "Align하려고 하는 Reference에 대해서 BWT 알고리즘을 이용해서 색인을 만들고 Suffix Array를 통해서 read들이 Reference에 잘 맞았는지에 대한 suffix array(SA)를 찾게 된다"

즉, bwa aln 과정에서는 SA에 대한 정보만 존재하기 때문에 single read냐 paired read냐에 따라서 bwa samse/sampe 명령을 통해서 실제 genomic상의 위치(coordinates)로 바꾸어 주게 된다. 따라서, sampe 명령을 통해 sam  파일을 생성할 때 "bwa sampe <in.db.fasta> <in1.sai> <in1.fa>"와 같이 reference, fastq, aln의 생성물인 sai를 모두 넣어주는 것이 SA를 genomic coordinates로 바꾸어 주기 위해서 이렇듯 모든 파일을 필요로 하는 것이다.

Paired-End Read에 대한 처리

BWA는 Single-End와 Paired-End에 대해서 각각 처리가 가능하며, Paired-End의 경우 read간의 거리(insert size)를 Library Construction 과정에서 지정하게 된다. nsert size에 대해서 지정을 해주면 좀더 빠르고 accurary하게 alignment를 수행하게 된다. 

Maximum distance between paired read (default 500)
paired read에 대해서 최대 insert size를 지정하게 되며, 일반적으로 Illumina의 paried-end의 경우 500정도이기 때문에 default값은 500이 지정된다.

BWA에서의 Gap에 대한 처리하기

기본적으로 BWA는 gap을 지원한다. gap이란 reference에 대해 reads를 align할때 빈 공간(gap)을 인식해서 이를 지원할것이냐는 것인데, 다음과 같은 ref, read가 있다고 할때

(ref)      AAGGCCTT
(read)   AACCTG


이를 gap을 지원하도록 align을 수행하면,

(ref)   AAGGCCTT
(read)   AA  -  -  CCTG


요렇게 '-'로 표시된 즉, gap을 고려해서 align하게 된다. sequencing에 문제가 있어서 또는 ref와는 달리 sequencing된 sample의 경우 'G'가 deletion되었을 수 도 있기 때문에 이러한 gap을 허용하는 것은 추후 variation calling에 있어서 중요한 역활을 하게 되지만, 이러한 gap을 허용함에 있어서 무조건적인 허용이 아닌 일정한 penalty를 부가하거나 read상에서 무한정 gap을 허용하는 것이 아니라 일정한 갯수의 gap만을 허용한다던지 gap의 크기를 위의 예에서는 1(하나의 염기)로 설정했지만 몇개나 허용할지에 대해서 설정이 가능하다.

(1) Maximum number of gap opens (default 2)
하나의 read에서  gap이 몇번 발생할지에 대한 옵션으로 지정된 값 이상의 gap이 발생한다면 해당 read는 align되지 않는다. 위의 예에서 하나의 read에서 총 gap은 1번 발생.

(2) Maximum gap size
하나의 read에서 gap의 크기까지 허용할 것인가에 대한 옵션 만약 10으로 설정하고, (1)번을 10으로 설정했다고 할 경우 read의 총 length가 90bp일 경우 10x10=100bp의 구멍까지 허용되기 때문에 왠만한 read는 ref에 붙게 되는 동시에 align 시간도 오래 걸리게 될 것이다. 위의 예에서 gap의 size는 2가 되겠다.

(3) Gap penalty (default 11)
각 read가 ref에 align시 mapping quality score를 통해 해당 read가 ref에 얼마나 잘 align된 건지를 표현하게 되며, gap이 많이 발생할 경우 당근 mapping quality score는 낮아져야 한다. 이때 gap이 발생할 경우 이 gap에 대해서 얼만큼의 패널티를 부여할 것인가를 지정한다. 

(4) Extension penalty (default 4)
gap의 크기가 커질수록 패널티를 부여하게 되는데, 이에 대한 패널티 값을 설정한다.

BWA Alignment시 mismatches처리하기

(1) Maximun number of read mismatches
Alignment를 수행할때 reference와 다른 서열(아까 gap이 아예 없는 거라면, 이번에는 다른것) mismatches에 대해서 어떻게 처리할 것이냐에 대한 부분이다. 위의 예에서 ref는 'T'이지만 read는 'G'로 서로 사맏디 아니한 경우라 할 수 있겠다.

BWA에서는 이러한 mismatche가 하나의 read내에서 몇개 혹은 read내 몇 %까지 허용할 것인지에 대해서 설정이 가능하며, 양의 정수일때는 갯수를 의미하며 float인 경우는 %를 의미한다. read 하나의 길이가 총 90bp이고 옵션을 0.04로 설정하게 되면 총 3.6개의 mismatches를 허용하게 되는 것이다.

(2) Mismatch penalty
Mismatch도 gap과 마찬가지로 alignment score를 지정할때 감점(Penalty)를 부여하게 된다.

BWA Alignment에서의 Seed 처리

(1) Seed Mismatches
ref에 대해서 read를 align하는데에 전체 read를 사용해서 align하는 경우와 seed를 사용해서 하는 두가지 방법이 있다. seed란 read의 길이보다 작은 read의 일부분을 가지고 ref에 align함으로서 속도면에서 빠른 align을 수행하게 된다. 이렇게 read의 일부분인 seed를 가지고 우선 빠르게 align한 후 seed를 기준으로 확장해 가면서 align을 수행하는데 이 seed(read의 일부분)를 ref에 align할때 몇개까지의 mismatches를 허용해서 align할 것인가를 지정하게 된다. 이부분은 0~3개(기본은 2개정도)가 적당하며, read의 mismatch보다 크게 잡아버리면 이 read의 align은 헛수고가 되는 꼴이 되겠다.

(2) Number of bases in seed region
이건 seed의 길이를 얼마로 지정할지에 대한 부분이며, seed의 길이가 read의 길이와 같거나 혹은 더 큰 값을 설정하게 된다면, seed 옵션을 사용하지 않는다는 것을 의미하게 된다. 일반적으로 short read에서의 seed의 길이는 25~35(기본은 30)으로 주어 align을 수행하게 된다. 

BWA 문서에 따르면 70bp의 read에서 32 bp의 seed와 maximum two differences를 2로 설정할 경우 seeding 없이 수행한 BWA에 비해서 2.5x 빠르게  align을 수행한다고 한다. 

기타 BWA의 Options

(1) Thread 갯수 (-t)
요즘에는 Multi-Core, Multi-CPU가 대부분으로  BWA는 multi-threading mode를 지원한다. 따라서 하나의 BWA Align 작업시 다수의 Core를 할당해서 좀 더 빠른 align이 가능하도록 할 수 있다.

(2) Illumina 1.3 quality score (-I)
FASTQ  포맷에서 사용하는 Quality값을 표현하는 방법은 다양한데, 이는 이전 블로그 (http://hongiiv.tistory.com/693) 글에서 언급한적이 있으므로 참고하기 바란다. 현재 대부분의 데이터는 Sanger (Phred+33, fastq-sanger, sanger) 포맷을 사용하는데 예전 Illumnina 데이터의 경우 Illumina 1.3+ (Phred+64, fastq-illumina) 포맷을 사용하는데 이때 사용하는 옵션이다.

(3) 특정 read group 지정 (-r)
samse/sampe 과정에서 특정 read group을 지정할 수 있도록 한다. '@RG\tID:foo\tSM:bar'의 형태로 read group을 지정해준다. GATK와 samtools의 경우 이 read group이 지정되지 않은 경우 실행되지 않는 경우가 종종 발생하기도 하기 때문에 read group을 지정해주는 습관이 필요하다. 

마지막으로 간단히 BWA 수행에 대해서 정리한다면, ref에 read를 align하는데에 있어서 주어진 길이의 seed를 이용해서 (seed의 mismatches를 허용하는 범위내에서) gap을 고려(gap이 나오는 횟수, 최대길이, 패널티)과 mismatches를 고려하여 read들을 aling하며 이때 align된 read들은 index된 ref의 SA의 값을(SAI포맷) 가지며, 이것은 samse/sampe를 가지고 실제 ref의 위치(coordinate)로 변환 (SAM포맷)되게 된다.

본 내용은 잘못된 부분이 다소 포함되어 있을 수 있으니 이점 고려하기 바라며, 댓글로 언급해준다면 고맙겠습니다.
반응형