MASON 시뮬레이터도구기능분석 Analysis MASON Read Generation Tool 권대건 부산대학교컴퓨터공학과 duskan@pusan.ac.kr Abstract 2000 년대후반부터유전정보를분석하는 NGS 분야에대한연구가활발히진행되고있다. 그중 RNA 에대한연구가활발히진행되고있는데주로 mrna 를이용하여 mrna 에서발현되는특징이전체유전자염기서열에서어느부분에해당하는가를찾는연구가큰비중을차지하고있다. 이러한연구는전체유전자서열에서특정부분이우리몸의어떤부분의발현에관여하는지를알아내는가장중요한연구가된다. 하지만현재유전정보에대한연구가활발히진행됨에도불구하고아직전체유전자서열에대해밝혀지지않은부분이많아, 실제유전데이터를이용하여도구의성능을측정하는것은힘들다. 그러므로이러한도구의성능을측정하기위해서는참조서열로부터가상의리드를생성하고, 매핑도구실행시나오는결과가얼마나정확한지평가하는시뮬레이터도구가필요하다. 본보고서에서는 wgsim 과더불어매핑도구성능평가에주로사용되는시뮬레이터도구인 MASON 에대해서알아보고 MASON 에서제공하는 Sanger, 454, Illumina 플랫폼별리드생성방식을비교한다. 그리고이를실제로수행해보고각플랫폼별생성되는리드를비교해보도록한다. Keywords: RNA, Simulator Tool, MASON, NGS 1 리드생성시뮬레이터의필요성 세계적으로큰관심을불러일으켰던 2004년인간게놈프로젝트 (Human Genome Project) 의종료선언이후개인의유전자정보를얻기위한연구는계속되었으며 Sanger[1] 방법을이용하여 2007년에처음으로개인의유전지도정보를얻을수있었다. 이후 2008년에는 FLX 454를이용하여 2007년에발표한개인유전자지도의 1% 의비용으로염기서열을분석이가능해지며이를이용하여대량의유전자정보를얻을수있게되었으며이후에도이에대한연구가계속되어왔다.[2] 대량의염기서열을얻을수있게되면서염기서열에대한광범위한분석이가능해졌고이를이용해염기서열과발현형태를분석하기위한여러가지매핑도구가연구, 개발되었다.[3] 이러한매핑도구의성능을평가하기위해서는리드가참조서열에어느정도로정확하게매핑되는지확인이되어야하는데, 하지만지속적인연구에도불구하고아직염기서열의많은부분이밝혀지지않고있으며, Indel, SNP과같은자연적인변이와유전자로부터염기서열을읽을때 1
발생하는오류로인해실제유전정보를이용하기에는많은제약이따르기때문에, 실제염기서열데이터를이용하여매핑도구들의정확한성능평가를한다는것은사실상불가능하다. 이러한매핑도구의성능을평가하기위해대부분의논문에서가상의리드서열을생성하는시뮬레이터를사용하고있으며이로인해염기서열을분석하는데있어서매핑도구뿐만아니라시뮬레이터도구또한중요한연구로자리잡고있다. 본보고서에서는 wgsim[4] 과더불어매핑도구의성능평가를하기위해주로사용되는시뮬레이터인 MASON[5] 시뮬레이터도구의기본배경과 MASON에서제공하는기능에대해분석해보고실제염기서열을시뮬레이션해봄으로써현재시뮬레이터의성능을확인하고추후시뮬레이터연구에있어서고려해야될사항을검토하려고한다. 2 MASON 시뮬레이터전체구조 MASON은 wgsim과더불어여러논문에서인용되는대표적인시뮬레이터도구이다. MASON 에서는염기서열의 RNA 발현과정에서자연적으로생기는변이와실제플랫폼에서시퀀스를읽을때오류가발생하는환경과유사한환경을시뮬레이션하여리드를생성한다. 실제플랫폼에서시퀀스를읽을때플랫폼의종류에따라생성되는리드서열이조금식차이가나게된다. 현재 NGS에서주로사용되는플랫폼은 Sanger, 454, Illumina, Solid인데, MASON에서는이중 Solid를제외한나머지 3개의플랫폼에대해오류모델을제공한다. 전체 MASON 기능구조는그림1와같이크게공통 (Global), 에러모델 (Error Model), Quality 모델에관련된옵션으로나뉜다. 공통옵션의경우주로리드생성자체에관련된옵션이나, 리드서열생성시자연적인변이값을설정하는옵션이다. MASON에서리드를생성하는방법은아래순서로진행된다. 1. 참조서열을파일로부터읽어들인다. 2. 참조서열로부터 indel, subtitution 를적용한리드를생성한다. 이때각모델의 Error설정에따라리드를생성한다. 3. 염기서열의 Quality를계산한다. 각모델별로 Quality를계산하는방법에차이가있기때문에, Quality의계산은각모델에맞게수행된다. 4. 리드에대한메타정보를추가한다. RNA 발현시일어나는자연적인변이작용은사용자가입력한 indel, subtitution 비율에따라 랜덤하게적용된다. 이후사용자가입력한플랫폼과에러정보에따라실제플랫폼과유사하게 2
MASON Option ] Sanger 454 Illumina Global ErrorModel Quaility Global ErrorModel Global ErrorModel Quaility 그림 1: MASON 전체구조. Mason 에서는 Sanger, 454, Illumina 플랫폼별기능을나누어제공한다. 시뮬레이션하여리드를생성한다. MASON의각플랫폼모델별특징은아래표1에잘나타나있다. 최근에가장많이사용되는 Illumina는다른플랫폼과달리리드길이가고정적이고매우짧다는단점이있으나비용이저렴하고빠르기때문에개인 DNA 분석과관련된분야에서가장각광받는플랫폼이다. MASON에서는실제플랫폼에서시퀀싱을수행하게될때신호를감지하지못하는경우, 신호를잘못읽는경우, 신호방출이늦게일어나서신호방출전에이미시퀀싱이완료되는경우와같은이유로잘못된데이터가삽입, 삭제, 변환되는것을플랫폼에서리드를읽어들이는단계에따라고려하여시뮬레이션한다. 시뮬레이션을통해리드를생성할때생성된리드가편향되지않도록하지만실제플랫폼에서유전자를읽어들일때화학적방법을사용하여읽기때문에, 리드서열에서에러를발생할시어느정도위치를고려하여리드를생성한다. NGS분야가생겨난후최초로상용화된 454 플랫폼은 Illumina에비해 throughput이낮아 NGS관련분야에많이사용되지는않다. 하지만 Sanger에비해저비용으로염기서열을추출할수있으며 Illumina보다리드길이가상대적으로길기때문에 De novo Sequencing 과같이몇몇특정한용도에서는비중이높게사용된다. 454에서의에러모델은 454 기반으로한시뮬레이터도구인 MetaSIM[6] 의에러모델을이용하였다. 표 1: 플랫폼별특징정리 플랫폼 Quaility 지원 리드길이 리드형태 Error model특징 Illumina O 짧다 고정길이 Position based Model 454 X 길다 가변길이 MetaSIM Sanger O 길다 가변길이 CelSIM Ramp Function 3
Sanger 플랫폼은 NGS이전에사용된고전적인염기서열결정방법으로비용이상대적으로비싸고시퀀싱과정이오래걸리기때문에현재실제산업에는잘사용되지는않다. 하지만 Sanger에서생성되는리드서열은다른플랫폼에비해정확하기때문에연구용으로아직많이사용한다. MASON에서는 1999년발표된 Sanger를기반으로한시뮬레이터도구인 CelSIM[7] 을참고하여에러모델을고려하였다. 3 MASON 시뮬레이터기능 MASON 은 3 개의플랫폼에대해리드생성을지원한다. 플랫폼마다리드의생성방식이나, 에 러발생이차이가있기때문에, MASON 에서는다음과같이세개의플랫폼을각각분리하여 제공한다../MASON sanger [OPTIONS] SEQUENCE./MASON 454 [OPTIONS] SEQUENCE./MASON illumina [OPTIONS] SEQUENCE 위세개의플랫폼에대해 MASON에서는공통적인옵션과각플랫폼의특징을제어하는옵션을제공한다. 공통적인옵션은주로리드서열의생성과일반적인변이에관한옵션들로이루어져있다. 또한앞서언급한것과같이플랫폼별로리드의생성방식이조금씩다르기때문에각플랫폼에맞는세분화된옵션을제공한다. 세분화된옵션은유전체로부터서열을읽을때플랫폼에서생성되는리드의길이나, 리드생성과정에서발생하는오류에대한옵션이다. MASON의옵션에대한자세한내용은다음과같다. 3.1 공통옵션공통적으로들어있는옵션은주로리드의형식을지정하는옵션과, 입출력파일형식을지정하는옵션이다. 대표적인옵션은다음과같다. 1. Main Option - ang : 참조서열중 N(None) 에대해서도리드에포함하도록한다. (default false) - N : 생성할리드서열의개수를결정한다. (Default : 1000000) - s : seed Number를설정한다. (Default : 0) - spa, -spc, -spg : 리드에서 A,C,G가나오는비율을설정한다. T는앞의 A,C,G 비율에의해결정된다.(Default : 0.25) 4
- o : outputfile 을설정한다. 2. Mate-Pair Option - mp : 해당옵션을선택하면 Mate-pair 형태로리드를생성한다. - ll : Mate-pair 의길이범위를의미한다 (Default : 1000) 3. Haplotype Opton - hn : Haplotype 의수를지정한다 (Default : 1) - hs : Haplotype에서 SNP 비율을설정한다. (Default : 0.001) - hi : Haplotype에서 indel 비율을설정한다. (Default : 0.001) - hm : Haplotype에서최소 indel 크기를설정한다. (Default : 1) - hm : Haplotype에서최대 indel 크기를설정한다. (Default : 6) - hnn : 삽입, 치환시 N을허용하지않는다. 위의옵션은리드를시뮬레이션할때주로사용되는옵션으로이외에도리드의 naming 설 정을비롯한 10 여개의옵션이더존재한다. 공통옵션에서는주로리드파일에대한설정이나, 변이와관련된부분에대해설정할수있다. 3.2 플랫폼별옵션플랫폼별옵션에는리드의길이, 리드생성시발생하는에러의모델, 그리고 Quality에대한설정을할수있다. 454는 Quility를지원하지않고, Illumina는고정길이의리드서열을생성하는것과같이각각의모델별로리드생성형식이약간식차가난다. 그렇기때문에각각의 Option 은약간씩의차이가있으며, 그중중요한몇몇 Option들만아래와같이소개한다. 3.2.1 Illumina 옵션 1. Illumina Read Length - n : 리드의길이를설정한다. 모든리드는동일한길이를지니고있다.(Default : 36) 2. Illumina Error Model - pi : insert가발생할확률을설정한다. (Default : 0.001) - pd : delete가발생할확률을설정한다. (Default : 0.001) 5
- pmm : mismatch 의평균비율을설정한다. (Default :0.004) - pmmb, - pmme : 리드의시작, 끝부분에서 mismatch 가생갈확률을설정한다. - nn : 이를설정할경우리드에 N 이나타나지않는다. Illumina 의에러모델에서는 Quality 의전체평균, 분산 mismatch 에서의평균, 분산대해서도 설정할수있다. 또한 Illumina 의경우리드길이가다른플랫폼에비해짧으며평균과분산으로 리드가설정되는것이아닌고정된길이를가지는리드를생성한다. 3.2.2 454 옵션 454는 Quality에대한설정은지원하지않는다.454의리드길이와에러모델에대한설정은다음과같다. 1. 454 Read Lenght Parametars - nu : 해당옵션을설정하면리드의길이를일정하게설정한다. - nm : 리드의길이를설정한다 (Default : 400) - ne : nu 옵션을설정하지않았을경우이값은편차를의미한다. nu가설정되어있다면일정한간격을의미한다. 2. 454 Error Mdoel parameters -nsq 설정할경우오류계산시제곱근연산을사용하지않는다. -k : 에러발생의비례상수의편차를설정한다. (Default : 0.15) 위에서언급한옵션이외에도에러모델에서리드에 noise가생길확률에대한옵션이존재한다. 3.2.3 Sanger 옵션 Sanger의경우리드길이, 에러모델, Quaility에대해설정할수있다. 1. Sanger Read Lenght Parametars Sanger의리드길이설정은 454와동일하여생략. 2. Sanger Error Model Parameters - pmb, -pme : 리드서열의처음, 끝에서 mismatch가발생할확률을설정한다. - pib, pie : 리드서열의처음, 끝에서 insert가발생할확률을설정한다. - pdb, pde : 리드서열의처음, 끝에서 delete가발생할확률을설정한다. 6
3. Sanger Quality Model Parameters 리드의시작, 끝부분의매칭, 오류에대한평균과표준편차를설정한다. Sanger의 Quaility에대한자세한설명은생략하였다. 에러모델의경우리드서열의처음, 끝에서 Error가발생할확률을설정하고전체에러비율을앞서설정한처음, 끝에서발생하는에러확률을이용하여 Ramp function에의해계산하여결정한다. Ramp function은다음과같다. x (x >= 0) y = 0 (x < 0) (1) 4 MASON 구동및 Option 비교 이제까지 MASON의특징을살펴보고 MASON에어떠한옵션이있는지확인하였다. 이제실제로 MASON을실행해보고임의로생성한염기서열로부터몇가지옵션에따른결과를비교해보도록한다. MASON의입력은 Fasta Format형태의시퀀스데이터이다. 출력은콘솔에서는 MASON 전체옵션값이출력되며, File 형태로 Fastq Format의리드시퀀스와리드의실제매핑위치결과인 SAM Format 데이터가출력된다. MASON 실행은 3개의플랫폼중 Illumina 플랫폼에서만테스트하였다. 참조서열은임의로생성한길이 100의시퀀스를사용하였다. MASON 옵션으로는리드생성개수 10개리드길이를 8bp로하고오류율을 10% 로두었으며, 이외에다른 Option들은기본옵션으로사용하였다. 실행문은아래와같다../MASON illumina -N 10 -n 8 -pi 0.1 -rnp read -o sample1.fastq sample1.fasta sample1.result 리드생성결과는그림2와같이출력된다. 리드의실재매핑결과로생성된 SAM File에서 8번리드의 Position과 cigar code를확인해보면 83, 7M1I 즉 83번째위치부터순서대로 7개매핑되고마지막한개는 Insert 발생되었다고표시되어있음을확인할수있다. 실제로생성된참조서열과생성된리드를비교하보면 참조서열 (83-90) : TGCATCAT 생성된 8 번리드 : TGCATCAN 으로 SAM 결과파일과일치함을확인할수있엇다. 7
그림 2: 참조서열길이 =100, N=10, n=8, Insert=0.1 일때 Fastq Format 의리드시퀀스 이후 MASON의리드생성개수와 Insert 비율을조정하며전체 Insert 개수를비교하였다. 표2 에서는 5000의길이를가지는참조서열로부터리드길이가 10 인리드를개수를달리하여 출력하고그결과를비교하였다. 위결과에서보면리드길이를제외한모든옵션을기본설정된 값으로둘때전체리드개수의 3% 정도되는 Insert가생김을확인할수있었다. 또한표3에서와 같이. 생성하는리드의개수를를 10000으로고정하고 Insert의비율을바꾸어가며비교해본 결과각각의 Insert 수에따라차이는있으나, 전체 Insert된리드의개수는어느정도비례하여 올라감을확인할수있었다. 표 2: 리드개수에따른 Insert된리드수 insert 수 -n 100 -n 1000 -n 10000 1 2 25 253 2 1 3 17 3 0 3 20 4 0 1 8 5 0 0 3 표 3: Insert Option에따른 Insert된리드수 insert 수 -pi 0.001 -pi 0.005 -pi 0.010 1 253 578 1071 2 17 17 30 3 20 20 20 4 8 5 1 5 3 1 2 8
5 결론및추후연구과제 매핑도구를검증하기위한시뮬레이터는매핑도구에대한연구와더불어 NGS분야에서앞으로유전자염기서열을연구하는데있어서매우중요한연구중하나이다. 본보고서에서는대표적인시뮬레이터중하나인 MASON의 3가지플랫폼에서리드생성방식의기본배경을알아보고각각의옵션을분석해보았다. MASON에서는각플랫폼에따라리드의길이, Quality, 에러모델에대한옵션을제공하였으며이를통해유사한형태의리드데이터를생성할수있엇다. 또한 MASON에서리드를생성하였을때생성된리드가여러곳에매핑이가능한경우, 실제올바르게매핑이되더라도 MASON에서제공하는 SAM File통해검증할경우잘못된결과출력될수있는가능성을확인할수있엇다. 이후연구에서는 MASON 시뮬레이터도구에서제공하는여러특징과발견된문제점을고려하여새로운시뮬레이터도구에대해연구하려고한다. References [1] Coulson AR Sanger F, Nicklen S, Dna sequencing with chain-terminating inhibitors, Proc Natl Acad Sci, vol. 74, no. 12, pp. 5463 5467, 1977. [2] Olena Morozova and Marco A. Marra, Applications of next-generation sequencing technologies in functional genomics, Genomics, vol. 92, no. 5, pp. 255 264, 2008. [3] Badr A Zhang G Zhang J, Chiodini R, The impact of next-generation sequencing on genomics, Genomics, vol. 38, no. 3, pp. 95 109, 2011. [4] Heng Li, Bob Handsaker, Alec Wysoker, Tim Fennell, Jue Ruan, Nils Homer, Gabor Marth, Goncalo Abecasis, Richard Durbin, and 1000 Genome Project Data Processing Subgroup, The sequence alignment/map format and samtools, Bioinformatics, vol. 25, no. 16, pp. 2078 2079, 2009. [5] Holtgrewe M, A read simulator for second generation sequencing data, Tech Rep, 2010. [6] Daniel C. Richter, Felix Ott, Alexander F. Auch, Ramona Schmid, and Daniel H. Huson, Metasim - a sequencing simulator for genomics and metagenomics, PLoS ONE, vol. 3, no. 10, pp. e3373, 10 2008. [7] Mayers G, A dataset generator for whole genome shotgun sequencing, 1 1999. 9