생물정보학및 RNA-Sequence 매핑도구소개 Introduction of Bioinformatics & RNA-Sequence Mapping Tools 권대건 부산대학교컴퓨터공학과 duskan@pusan.ac.kr Abstract Frederick Sanger 에의해서시퀀싱기술이개발된이후오래동안시퀀싱과관련된연구가계속되었고 2003 년에는 13 년간의연구끝에한성인남성한명의 genome 을분석하기도하였다. 이후 2007 년유전자서열을분석하는 NGS(Next Genaration Sequencing) 분야가생기고꾸준히발달하게되면서유전자서열을분석하기위한수많은기술이발달하였다. NGS 분야의발달로 DNA, RNA 염기서열데이터가증가하게되었고, 늘어난데이터를분석하기위해유전자서열분석에관한수많은연구가진행되었다. 현재컴퓨터를이용한수많은염기서열분석방법에대한연구가활발히진행되었으나아직염기서열의많은부분이미지의영역으로남아있으며염기서열에대한새로발견되는부분과, 미지의부분으로남아있는염기서열을분석하기위한도구의개발은여전히필요하다. 본보고서에서는앞서언급한 NGS 분야에대해설명하고현재발표된분석도구에서가장기본이되는부분인 BWT(Burrows Wheeler Transform) 알고리즘을응용한방법과해시테이블을이용한매핑방법을 BWT 기반의분석도구인 Bowtie 와해시태이블기반의분석도구인 mrfast 의예를통해설명하고이후연구방향에대해논의하고자한다. Keywords: Burrows-Wheeler Transform, Hash Table, Next Generation Sequencing 1 Next Generatrion Sequencing 소개 NGS(Next Generation Sequencing) 는 2007년유전자서열을분석하는대표적인회사인 Sanger 사와 Illumina사가합병되면서사용되기시작한용어이다. 세계적으로큰관심을불러일으켰던 2004년인간게놈프로젝트 (Human Genome Project) 의종료선언이후개인의유전자정보를얻기위한연구는계속되었으며 Sanger방법을이용하여 2007년에처음으로개인의유전지도정보를얻을수있었다. 이후 2008년에는 FLX 454를이용하여 2007년에발표한개인유전자지도의 1% 의비용으로염기서열을분석이가능해졌다. 이후에도저비용으로염기서열지도를생성하는방법에대한연구는계속되었고현재에는인터넷을통해대량의유전자정보를쉽게접할수있게되었다. 염기서열정보량이급증하고대량의데이터로부터기존에알수없던여러가지 1
염기서열의특징을발견하게되면서 NGS에서는쏟아지는데이터를분석하기위한도구가필요하게되었다. SNP(Single Nucleotide Polymorphism), MNP(Multi Nucleotide Polymorphism), Indel(Insertion and Deletion) 등과같은수많은유전체변이를발견하였고이러한유전적변이를분석하기위해많은연구가진행되고있다. 하지만염기서열을장치로부터읽어들이는과정에서오류가발생할수있으며, 앞에서언급한변이현상에대한연구도부족한부분이존재하므로염기서열을읽고분석하는데에는많은연구가필요하다. 본보고서에서는이러한분석도구의연구에앞서분석도구에사용되는대표적인알고리즘인 BWT와해시테이블기반알고리즘에대해현재발표된도구를예로들어기존에발표된알고리즘에대해소개한다. 이후 BWT와해시테이블기반알고리즘으로동작하는여러매핑도구를비교분석하도록한다. 2 매핑알고리즘 2.1 Burrows-Wheeler Transform 기반알고리즘 BWT(Burrows-Wheeler Transform)[1] 알고리즘은 1994년 Burrows와 Wheeler가처음제안한방법으로 BWT를이용하여 BW Matrix을생성성하고 BW Matrix를 SuffixArray, FM-indexing과같이응용하여리드와매칭이되는지확인한다. Bowtie, BWA, Bowtie2와같은가장대표적인분석도구역시이러한 BWT 알고리즘을이용한도구이다. 본논문에서는 BWT 기반의매핑도구인 Bowtie 매핑도구를예로들어 BWT알고리즘과 Bowtie[2] 에서 BWT를활용한매핑과정을설명한다. 2.1.1 BW Matrix 0 BWT를수행하기위해서는원본서열을이용하여정렬된 rotation Matrix을만들어야한다. 본논문에서는이를 BW Matrix이라고정의한다. BW Matrix을만드는과정은아래와같이진행된다. 1) 그림 1의 a에서처럼먼저 Suffixes Matrix을만든다. 시작과끝을구분하기위해맨앞에원본서열맨앞에 $ 기호를표시하고모든접미사서열을만든다. 이때접미사서열은접미사를쓰고뒤에원본서열에서남은부분을쓰도록한다. 이렇게만들어진접미사서열은 Matrix 형태로저장한다. 2) 접미사 Matrix에서첫번째행을기준으로오름차순으로 Sorting 한다. 만약첫글자가같은경우짧은접미사를가지는경우에우선순위를둔다. 이렇게만들어진 Matrix를 BW Matrix 이라고한다. 3) BW Matrix을저장할때에는 BW Matrix의처음행과끝행만저장한다. 처음행과끝행만가지고있으면전체원본서열을복구할수있다. 2
그림 1: BWT(Burrows Wheeler Transform) 에서 BW Matrix 생성과정 이렇게만들어진 BW Matrix은특이한성질을가지고있는데 BW Matrix 마지막행에서 i번째로등장하는문자는 BW Matrix의첫번째행에서 i번째로등장하는문자와동일한문자이다. 즉그림 1의 b에서 4열마지막행의 G는마지막행에서 2번째로등장하는문자이므로첫행에서두번째로나타나는 G인 6열의 G와동일하다는점이다. 이러한 BW Matrix의특징때문에맨처음행과맨마지막행정보만가지고있어도전체 Sequence를추출할수있으며, 이론적으로 O(n) 의시간안에검색이가능하다. 2.1.2 Burrows-Wheeler Transform 매칭 - bowtie 그림 2은 BWT의가장대표적인도구인 Bowtie에서사용되는매칭방법을나타낸것으로저장된 BW Matrix를이용하여원본서열인 AGCTCAT 에서리드조각인 TCA 를찾는과정을보여준다. 매칭방법은다음과같이진행된다. 1) 리드서열의맨마지막글자 A 로시작하는열을찾는다. 2) 찾은열의맨마지막행의글자가 n-1번째글자즉 C 인지확인한다. 3) 확인이완료되면, 첫번째행에서 n-1번째글자와동일한글자를찾는다. 4) 모든글자에대해매칭이완료될때까지위 1 3과정을반복한다. 모든글자에매칭될때까지위과정을반복함으로써정확하게매칭이됨을확인할수있다. 3
그림 2: Bowtie 에서 BW Matrix 를이용한매칭과정 실제 Bowtie 도구에서는 BW Matrix 에각행에참조서열에대한위치정보를추가하여어느곳에 정확히매핑되는지도확인가능하다. 2.1.3 에러율을고려한매칭리드는리드서열을읽는과정에서오류가발생하거나, mrna 생성과정에서자연적으로변이가일어날수있으므로리드에는약간의오류가존재한다. 이럴경우 mrna로부터읽어온리드임에도참조서열에매칭되지않는경우가발생하기도한다. 이렇기때문에 Bowtie를포함안여러매핑도구는사용자가필요한경우입력한값에따라리드를매핑할때약간의오류율의고려하여수행할수있어야한다. Bowtie에서는여러가지에러들중 Subtitution( 치환 ) 에대해서만고려하여수행된다. 매핑수행도중정확히매칭되는서열이없는경우가발생하면 Bowtie는해당리드의각각의글자에대해스코어를계산하여가장낮은스코어를가지는글자를남은세개의글자로변환하면서매칭되는문자열이있는지확인한다. 만약하나를치환하였을때매칭이된다면이를 1-error match 라고한다. Bowtie이외에 BWT기반도구중에는에러의가능성이있는리드를매핑할때 Subtitution이외에도 Indel,Gap을고려하는도구들도있다. 2.2 해시기반알고리즘해시기반알고리즘은참조서열을 k-mer과같이분할한후해시테이블에분할한서열을키값으로하여해당서열의위치를저장한다. 이후매핑과정을수행할때에리드역시 k-mer로나누어해시테이블을통해비교한다. 해시기반도구의경우이론적으로는한번해시테이블을생성하고나면이후에는 O(1) 의수행시간만에리드를매핑하는것이가능하지만실제도구에서는 4
Collision과리드의에러율에대해고려해야하므로 O(1) 의시간내에수행하기는것은약간어렵다. 그러므로해시기반도구는앞의두가지문제점을해결하는방법에따라알고리즘과도구의성능에서차이가발생한다. 본논문에서는여러해시기반도구중 mrfast[3] 의알고리즘에대해분석한다. 2.2.1 해시기반알고리즘 - mrfast mrfast에서는참조서열과리드를 k-mer 로분할하여매칭을수행한다. mrfast에서는매핑과정을수행하기에앞서참조서열을읽어그림3에서와같이해시테이블을생성한다. 해시테이블의키값은 k-mer에서나올수있는모든경우의수를키로설정하며, 값에는키에해당되는염기서열조각의참조서열에서의위치를저장한다. 리드의크기는적게는 100pb부터많게는 1kbp 까지매우다양한길이의리드가나타나는데,Collison을방지하기위해매우긴길이의키값을가지도록해시테이블을구성하기에는현재하드웨어의메모리의한계로인해해시테이블이제대로생성할수없다. 이러한점때문에 mrfast에서는작은길이의키를설정되어있으며어쩔수없이많은 Collision이생기게된다. mrfast에서는그림3에서처럼 list형태로값을저장하여 Collision 이발생하는문제를해결하였다. 2.2.2 매핑과정 리드서열을매핑하기위해다음과같은과정을거친다. 1) 리드서열을해시테이블의키값의길이만큼 k-mer로나눈다. 2) 나누어진서열을해시테이블에대입하여입력한리드조각에대해매핑가능성이있는후보위치리스트들을받아온다. 3) 후보리스트에있는위치가올바른위치인지 seed-and-extend 방식을이용하여검증하고검증에통과한위치정보를매핑결과로출력한다. 다시말해서리드들을 k-mer로자른후해시테이블에대입하여대략적인위치정보리스트를받는다. 이후 seed-and-extend를통해해당위치의앞뒤의염기서열들이리드와일치하는지확인하고일치한다면매핑결과로출력하게된다. 이때후보위치들이많아검증시간이매우길어지기때문에 mrfast에서는 AF를이용하여검증전에전처리작업을통해가능성이없는후보위치들을제거함으로써속도를향상시켰다. 5
그림 3: 해시테이블을이용한 mrfast Indexing 2.2.3 Adjacency Filtering AF는 리드에서자른리드조각들은참조서열에도가까이있을것이다 라는가정하에가능성이낮은리드조각을사전에제거하는방법이다. 그림 4의 a에서처럼리드를조각으로나누게되면이조각들이참조서열에서도비슷한곳에위치하게된다. k=12 일때 324, 459, 535 와같이떨어진경우같은리드라고보기는어렵다. mrfast에서는 AF필터링을통해위와같이검증하기전에미리제거가능한부분을제거하여도구의성능을증가시켰다. 2.2.4 에러율을고려한매핑해시태이블기반알고리즘의경우시퀀스중한개의에러만발생하더라도해시값이전혀다른값으로출력되기때문에 error를찾기힘들다. mrfast에서는허용가능한 edit distance e를두고 e+1개의리드조각에대해 seed-and-extend방식을이용하여검증함으로써매핑되는위치를찾는다. 만약하나의리드에 2개의 e를허용한다고할때 3개의리드조각에대해서검증하면 2개의리드조각이에러를가지고있더라도나머지하나는원래매핑되어야할위치에매핑된다. 이때 CKS(Cheep K-mer Selection) 알고리즘 [4] 을이용하여에러가발생한경우에대한검증과정을최소화하였다. 그림 4의 b에서보면 e = 1일때최소 2개의리드조각에대해검증해야하는데, 이때 1,3번째리드를검증하게되면 1004번의검증을거쳐야하지만가장저비용이드는 1,2번째리드조각을검색하면 6번만검증과정을거치면된다. 즉상대적으로적게매핑된리드조각을선정하여수행함으로써보다빠르게매핑이가능하다. mrfast에서는앞에서설명한것처럼 CKS 알고리즘을이용하여도구의성능을증가시켰다. 6
그림 4: (a)af(adjacency Filtering 에서인접서열선택과정과 (b)cks(cheep K-mer Selection) 에서 Cheep K-mer 를찾기위한과정 3 Burrows-Wheeler Transform 과해시테이블을이용한매핑도구 위에서언급한 Bowtie, mrfast 이외에도표 1 과같이 BWT 와해시테이블을이용한많은매핑 도구가존재한다. 앞에서언급하지못한다른도구들을간단히소개하고이도구들의장단점을 분석하고자한다. 3.1 Burrows-Wheeler Transform을이용한매핑도구 3.1.1 BWA BWA[5] 는 BWT와접미사배열을이용하여정렬을수행하는알고리즘이다. BWT변환에의해생성된 BW-Matrix의접미사가 SA(Surffix Array) 의어느구간에해당하는지알면, BW-Matrix 의매핑결과로 SA구간을찾음으로써원본서열에서의위치를알수있다. BWA는 subtitutioin 뿐만아니라 indel에대해서도 mismatch를수행한다. 백트래킹을통해 mismatch를수행하는데, 이론적으로는 BWA를이용하여모든 k-mismatch를찾을수있으나비용이매우커지므로제한된범위내에서 k-mismatch를수행한다. 3.1.2 SOAP2 SOAP2[6] 는매핑속도에중점을둔 BWT기반도구이다. 전체적인진행과정은 BWA와유사하지만, 속도를올리기위해최대 2개의 mismath만허용한다. 즉탐색도중더이상탐색이불가능한구간에도달하면그위치의문자를다른문자로 subtitution하여수행한다. SOPA2는매우빠른속도를보이지만, mismatch수가제한적이며 error가많을경우속도가저하된다. 7
표 1: 알고리즘별매핑도구특징정리 기반알고리즘 매핑도구 특징 Bowtie BW행렬과 LP mapping 법칙을이용하는알고리즘. Subtitution에대해수행한다. BWT BWA 접미사배열과 BWT를이용하여 SA구간을계산하여정렬을수행. Subtitution, Indel 전부수행 SOAP2 BWA와동일한구조, 속도증가를위해 Subtitution에대해서만수행. mrfast k-mer로나올수있는모든시퀀스에대해해싱한후참조서열의주소를해시테이블에저장. 리드를 kbp 단위로 k-mer 로나누어수행. 해시테이블여러개의탬플릿쌍을이용하여리드를해싱한뒤참조서 MAQ 열을정렬. Stampy 참조서열을 k-mer로해싱하고리드에서 mismatch를허용하는범위안의모든 k-mer를생성한후비교. 3.2 해시테이블을이용한매핑도구 3.2.1 Stampy mrfast에서적은길이의 k-mer를이용하여 k-mer로나오는시퀀스전체를해시테이블로사용한것과달리 Stampy[7] 는참조서열을 k-mer로나누어해시태이블로생성한다. 이렇게해시테이블을만들게될경우메모리비용이매우커지게될수있으므로몇개의 bp씩만겹치게하여해시테이블의크기를줄이게한다. 또한과도하게반복되는 k-mer의경우일정수가넘어갈경우에는반복시퀀스로간주하고더이상해시테이블에저장하지않는다. 해시테이블이생성되고나면리드에서가능한모든 k-mer을생성하여해시테이블과비교한다. 이때최대한개의불일치까지허용된다. Stampy는리드로부터최대한많은 k-mer를생성하여비교하기때문에높은정확도를보이나, 계산량이많아속도가느려지는단점이있다. 3.2.2 MAQ MAQ[8] 는템플릿을이용하여리드서열로부터여러개의해시값 ( 템플릿쌍 ) 을만들어해시테이블을구축한후, 참조서열에비교하여리드의위치를찾는다. 이후하나의탬플릿쌍을선택하여전체참조서열을확인하며템플릿쌍에매칭되는지비교하고다음템플릿쌍으로이를반복한다. MAQ는 k개이하의모든 mismatch쌍을비교할수있으나 mismatch수가증가할수록템플릿수가증가하며템플릿수만큼전체서열을비교해야한다. 또한매번리드를수행할때마다해시테이블을생성해야하는단점이있다. 8
4 결론및향후연구과제 현재 NGS 기술의발달로대량의유전자정보를얻을수있게되면서, 이전에는알지못했던새로운 DNA, RNA와관련된현상들이발견되고있다. 기존에는알지못했던많은양의정보가들어오게되면서유전자를분석을하기위한도구에대한연구도매우중요해졌다. 본논문에는 Bowtie와 mrfast 도구를예로들어 BWT와해시기반매핑방법에대해자세히알아보고 BWT 와 mrfast를이용하는다른도구들을확인하였다. BWT기반의알고리즘은해시테이블기반도구보다유전적변이와오류에대해강인하면서도비교적빠른도구를개발할수있다. 반면에해시테이블에서는 matching의경우에는빠른속도로리드들을참조서열에매핑이가능하며, 정확도가매우높으나오류나변이가많이일어난리드일수록속도가떨어지는단점이존재한다. 하지만두도구의특성이서로다르기때문에두도구간의비교가쉽지않았다. 이후연구에서는두도구를포함한여러가지도구를비교하기위해 Simulator에대해조사하고조사한 Simulator 를이용하여두도구를포함한최신에발표된여러도구들에대해비교하고자한다. References [1] M. Burrows, M. Burrows D. J. Wheeler, and D. J. Wheeler, A block-sorting lossless data compression algorithm, 1994. [2] B Langmead, C Trapnell, M Pop, and SL Salzberg, Ultrafast and memory-efficient alignment of short dna sequences to the human genome, Genome Biol, vol. 10, no. 3, pp. R25, 2009. [3] Marques-Bonet T Aksay G Antonacci F Hormozdiari F Kitzman JO Baker C Malig M Mutlu O Sahinalp SC Gibbs RA Eichler EE Alkan C, Kidd JM, Personalized copy number and segmental duplication maps using next-generation sequencing, Nat Genet, vol. 41, pp. 1061 1067, 2009. [4] Hongyi Xin, Donghyuk Lee, Farhad Hormozdiari, Samihan Yedkar, Onur Mutlu, and Can Alkan, Accelerating read mapping with fasthash, BMC Genomics, vol. 14, no. 1, pp. 1 13, 2013. [5] Durbin R. Li H, Fast and accurate short read alignment with burrows-wheeler transform, Bioinformatics, vol. 25, no. 14, pp. 1754 1760, 2009. [6] R. Li, Soap2: an improved ultrafast tool for short read alignment transform, Bioinformatics, vol. 25, no. 15, pp. 1966 1967, 2009. 9
[7] M. Goodson G. Lunter, Stampy: A statistical algorithm for sensitive and fast mapping of illumina sequence reads, Bioinformatics, vol. 21, pp. 936 939, 2011. [8] R. Durbin H. Li, J. Ruan, Mapping short dna sequencing reads and callingvariants using mapping quality scores, Bioinformatics, vol. 18, pp. 1851 1858, 2008. 10