다중목적함수진화알고리즘을이용한마이크로어레이프로브디자인 501 다중목적함수진화알고리즘을이용한마이크로어레이프로브디자인 (Microarray Probe Design with Multiobjective Evolutionary Algorithm) 이인희 신수용 조영민 양경애 장병탁 (In-Hee Lee) (Soo-Yong Shin) (Youngmin Cho) (Kyung-Ae Yang) (Byoung-Tak Zhang) 요약프로브 (probe) 디자인은성공적인 DNA 마이크로어레이 (DNA microarray) 실험을위해서필수적인작업이다. 프로브가만족시켜야하는조건은마이크로어레이실험의목적이나방법에따라다양하게정의될수있는데, 대부분의기존연구에서는각각의조건에대하여각각독립적으로정해진한계치 (threshold) 값을넘지않는프로브를탐색하는방법을취하고있다. 그러나, 본연구에서는프로브디자인을두가지목적함수를지닌다중목적함수최적화문제 (multiobjective optimization problem) 로정의하고, ε - 다중목적함수진화알고리즘 (ε -multiobjective evolutionary algorithm) 을이용하여해결하는방법을제시한다. 제시된방법은 19 종류의고위험군인유두종바이러스 (Human Papillomavirus) 유전자들에대한프로브디자인과 52 종류의애기장대칼모듈린유전자군 (Arabidopsis Calmodulin multigene family) 에대한프로브디자인에각각적용되었다. 제안한방법론을사용하여기존의공개프로브디자인프로그램인 OligoArray 및 OligoWiz 에비해목표유전자에더적합한프로브를찾을수있었다. 키워드 : 진화연산, 마이크로어레이프로브디자인, 다중목적함수최적화, ε - 다중목적함수진화알고리즘 Abstract Probe design is one of the essential tasks in successful DNA microarray experiments. The requirements for probes vary as the purpose or type of microarray experiments. In general, most previous works use the simple filtering approach with the fixed threshold value for each requirement. Here, we formulate the probe design as a multiobjective optimization problem with the two objectives and solve it using ε -multiobjective evolutionary algorithm. The suggested approach was applied in designing probes for 19 types of Human Papillomavirus and 52 genes in Arabidopsis Calmodulin multigene family and successfully produced more target specific probes compared to well-known probe design tools such as OligoArray and OligoWiz. Key words : Evolutionary computation, Microarray probe design, Multiobjective optimization, ε -Multiobjective evolutionary algorithm 본연구는교육인적자원부 BK21-IT, 산업자원부차세대신기술개발사업의분자진화컴퓨팅 (MEC) 과제및과학기술부국가지정연구실 (NRL) 사업, 그리고 2006년정부 ( 교육인적자원부 ) 의재원으로한국학술진흥재단의지원을받아수행된연구임 (KRF-2006-214-D00140). 또한이연구를위해장비를지원하고공간을제공한서울대학교컴퓨터연구소에도감사드린다. 학생회원 : 서울대학교컴퓨터공학부 ihlee@bi.snu.ac.kr 비회원 : 서울대학교병원의료정보센터 syshin@snuh.org 비회원 : Department of Computer Science and Engineering, University of California, San Diego yoc002@cs.ucsd.edu 비회원 : 서울대학교컴퓨터연구소연구원 kayang@bi.snu.ac.kr 종신회원 : 서울대학교컴퓨터공학부교수 btzhang@bi.snu.ac.kr 논문접수 : 2008년 2월 22일심사완료 : 2008년 6월 25일 Copyright@2008 한국정보과학회ː개인목적이나교육목적인경우, 이저작물의전체또는일부에대한복사본혹은디지털사본의제작을허가합니다. 이때, 사본은상업적수단으로사용할수없으며첫페이지에본문구와출처를반드시명시해야합니다. 이외의목적으로복제, 배포, 출판, 전송등모든유형의사용행위를하는경우에대하여는사전에허가를얻고비용을지불해야합니다. 정보과학회논문지 : 소프트웨어및응용제35권제8호 (2008.8)
502 정보과학회논문지 : 소프트웨어및응용제 35 권제 8 호 (2008.8) 1. 서론 DNA 마이크로어레이 (microarray) 는많은수의유전자의발현정도를동시에측정할수있는실험도구로서다양한생물학실험에널리사용되고있다. 마이크로어레이의표면은미세하게구획이나누어져있으며, 각각의구획내부에는발현정도를관찰하고자하는개별유전자에특이적으로결합할수있는 DNA 분자 ( 프로브, probe) 가부착되어있어서, 세포샘플과반응시킨후, 각구획별로프로브와유전자가결합한정도를측정하여해당유전자들의발현정도를알수있게된다. 이때, 표면에부착되는프로브의생성방식에따라 cdna 마이크로어레이와올리고뉴클레오타이드 (oligonucleotide) 마이크로어레이의두가지로나눠진다. 이중 cdna 마이크로어레이는서열에대한제어를하기힘들고, 서열생성과정상오류를완전히배제하기어렵기때문에정확한결과를요구하는실험에는잘쓰이지않고있다. 반면에올리고뉴클레오타이드마이크로어레이는 DNA 합성기술의발달로인해프로브서열을제어하기쉽고, 사용자의용도에맞게프로브의서열을디자인하여사용할수있다는장점이있어서널리사용되고있다. 그러므로본연구에서는올리고뉴클레오타이드마이크로어레이에사용될프로브의최적화방법에한정하여논의하고자한다. 마이크로어레이실험에서목표유전자의발현정도는프로브와목표유전자와의결합정도로부터추정되기때문에, 개별프로브가각각의목표유전자와만반응하도록디자인하는것이중요하다. 만약, 한프로브가둘이상의유전자와결합가능할경우, 해당프로브가원래목표했던유전자의발현정도를정확히알수없게되어, 부정확한분석결과를유도하게될수있다. 따라서목표유전자에특이적인프로브를디자인하는문제는마이크로어레이를사용한연구에서가장중요하고근본적인문제중의하나라고할수있으며, 이점은관련문제를다룬연구논문들의방대함에서도확인할수있다 [1-8]. 대부분의관련연구에서공통적으로언급되는좋은프로브의조건은다음의세가지이다 : 1) 목표유전자와특이적으로강하게상보결합을이루되, 2) 프로브자체적인 2차구조의형성가능성이낮아야하며, 3) 녹는점 (Tm, melting temperature) 등의상호결합반응에영향을줄수있는화학적인성질이균질해야한다. 이러한조건을만족시키는프로브를디자인하기위해다양한연구그룹으로부터여러가지프로브평가방법이제시되어왔는데, 대부분 DNA의결합에필요한자유에너지 (free energy), 목표유전자가아닌서열과프로 브사이의 BLAST 검색점수, 프로브의 2차구조형성에필요한자유에너지와프로브자체의녹는점등을각각검사하여사용자가정한기준치를모두만족시키는프로브를찾아내는방식을취하고있다 [1,2]. 이외에도프로브서열의빈도수에따른방식 [3], Shannon entropy에기반한정보이론적인접근방식 [4], 유전체정보에기반한방식 [5] 등이있다. 그런데프로브의평가방법의다양함과는대조적으로적당한프로브를찾기위한탐색방식에있어서는대개의연구에서공통적인방식 ( 그림 1의 (a)) 을취하고있다 : 먼저프로브를디자인하고자하는목표유전자의서열로부터슬라이딩윈도우 (sliding window) 방식등으로후보프로브서열을생성시킨다. 그후, 앞서언급한여러평가방법들을사용하여후보프로브서열들에점수를부여하되, 미리정의된한계치 (threshold) 값을만족시키지못하는후보서열은제거하고마지막으로남은서열들중에서각각의평가항목에대한점수를통합한결과가가장우수한후보서열을해당목표유전자의프로브로선정한다. 이과정을각각의목표유전자마다수행하여전체목표유전자집합에대한프로브집합을구성하게된다. 그러나이와같은탐색방식에서는몇가지문제점을찾을수있다. 우선프로브평가점수와실제마이크로어레이실험결과에서의형광강도값의연관관계가아직불명확하기때문에평가점수에대한명확한한계치를설정하기가어렵고, 개별실험환경에따라여러조절이필요하다는점을들수있다. 그리고각각의목표유전자에대한프로브탐색을독립적으로진행하기때문에개별목표유전자에대한프로브들사이에서발생할수있는비특이적결합이고려되지않을수있다. 또한최근의연구결과에의하면 [6-8], 우수한프로브는한두가지기준으로정의될수없기때문에근본적으로다중목적함수최적화 (multiobjective optimization) 를필요로한다. 이와같은관찰을바탕으로다중목적함수진화알고리즘 (MOEA, multiobjective evolutionary algorithm) 을올리고뉴클레오타이드마이크로어레이를위한프로브디자인에적용하는방법이제안된바있으며 [9], 본연구에서는보다효율적인탐색과사용상의편의를위하여추가의과정을도입하고, 기존의프로그램과비교분석을통하여다중목적함수진화알고리즘을이용한프로브디자인의유용성을검증하였다. 그리고온라인인터페이스 (EvoOligo, http://cbit.snu.ac.kr/evooligo) 를제공하여관련연구자의편의를도모하였다. 프로브최적화를다중목적함수최적화문제 (MOP, multiobjective optimization problem) 로보고이를다중목
다중목적함수진화알고리즘을이용한마이크로어레이프로브디자인 503 그림 1 프로브탐색방법의비교. (a) 기존의유전자별프로브탐색방법. (b) 다중목적함수진화알고리즘을이용한프로브조합탐색방법. 적함수진화알고리즘으로해결할때의장점으로는우선기존의방법에서처럼다양한평가항목들을하나의수치로표현하는데서발생할수있는왜곡현상을줄일수있고, 적절한한계점값이나각각의평가항목들을하나의수치로표현하기위한가중치 (weight) 조절에따른반복적인작업을피할수있다는것을들수있다. 그리고진화연산의확률적 (stochastic) 이고병렬적인탐색특성으로인하여여러평가항목들사이의다양한타협점 (trade-off point) 들을찾아서한번에여러가지의프로브디자인이가능하다 ( 그림 1의 (b)). 또한, 사용가능한사전정보 (prior knowledge) 또는특정목적에따른사용자의기호등을목적함수의변형또는추가를통하여손쉽게반영할수있다. 본논문은다음과같이구성되었다. 2장에서프로브디자인을다중목적함수최적화문제의관점에서살펴보고, 3장에서다중목적함수진화알고리즘을이용한프로브최적화에대해구체적으로설명한다. 구체적인적용사례로서고위험군인유두종바이러스 (Human Papillomavirus) 유전자군과애기장대칼모듈린유전자군 (Arabidopsis Calmodulin multigene family) 에대한실험결과가 4장에서제시되며, 마지막으로 5장에서요약및결론을내리고자한다. 2. 다중목적함수최적화를통한프로브디자인 2.1 다중목적함수최적화문제다중목적함수최적화문제 (MOP) 란여러개의목적함수 (objective) 를한꺼번에최적화할수있는해를찾는문제를말하는데, 주로목적함수들사이에상충하는관계가존재하여하나의목적함수에대한최적화방향이다른목적함수에대한최적화방향과일치하지않아서동시에최적화시키기어렵고, 또한목적함수들사이의우선순위 (priority) 가주어지지않아서우선순위에따른사전적순서 (lexicographic order) 대로최적화하는것도불가능한경우의문제를의미한다. MOP의일반적인형태는다음과같다 :,,. 여기에서 는개별목적함수 (objective), 는개별제한조건 (constraint), 은목적함수의개수, 은제한조건의개수를각각의미하며, 와 는변수 가가질수있는최소값과최대값을의미한다. MOP에서변수공간 (variable space) 에서하나의해 에는이에상응하는목적함수공간
504 정보과학회논문지 : 소프트웨어및응용제 35 권제 8 호 (2008.8) (objective space) 상의벡터 가존재하며, 변수공간에서서로다른해사이의우열은이에상응하는목적함수공간의벡터를비교하여가리게된다. 일반적으로두해의비교에는 dominance 관계를사용하는데, dominace 관계의정의는다음과같다. 를각각최소화해야하고 ( 최대화의경우부등호의방향이반대이다 ), 제한조건이모두동등할경우, 변수공간의서로다른해 와 및각각에상응하는목적함수공간의벡터 와 에대하여, 의두조건을모두만족시킬때, 는 보다더좋은해로평가받고, 이를 가 를 dominate한다 또는 로표현한다. 이것은여러개의목적함수사이의우선순위가존재하지않기때문에하나의해가다른해보다낫다고말하기위해서는모든목적함수값의면에서못하지않고, 또한하나이상의목적함수에서더나은값을가져야함을의미한다. 따라서 MOP에서최적해는변수공간의어떤해에도 dominate되지않는것을의미하며이를 Pareto-최적 (Pareto-optimal) 해라고한다. 또한두해중어느쪽도다른하나를 dominate하지못하는경우가있을수있는데, 이때두해는 서로 non-dominate한다 고표현하고서로동등한수준의해로간주된다. 일반적인 MOP에서 Pareto-최적해는목적함수간의 trade-off 관계때문에서로 non-dominate 인여러개의타협점들로구성된다. 2.2 다중목적함수진화알고리즘실제 MOP의응용에서는최적화알고리즘을사용하여다양한타협점을찾아내고, 마지막에이들타협점중에서해당응용문제및사용자의요구에따라적당한것을선택하게된다. 이때가능한한 Pareto-최적해에가까운타협점들을찾아내기위하여목적함수공간의경사면 (gradient) 을이용하거나목적함수들에임의의가중치를부여하여합산한값을최적화하는등의여러가지방법들이제시되어왔다 [10]. 그러나이러한방법들은대개특정 MOP의형태또는최초시작위치에의존적이거나, 지역최적해 (local optimal solution) 에서벗어나지못한다는단점을가지고있다. 또한, 한번에하나의타협점만을찾도록되어있어서여러개의다양한타협점을찾기위해서는여러번의최적화를반복해야한다. 고전적인다중목적함수최적화알고리즘의이러한단점을해결할수있는방안중의하나로주목받고있는것이진화알고리즘이다 [11]. 일반적인진화알고리즘은초기에일정한수의해를무작위로생성하여주어진문 제에따라적합도를평가한후이에따라좋은해를선택하여변형시키는과정을반복하는방식으로단계적으로성능이좋은해를찾아나가게된다. 이때적합도에따라해들을비교하여좋은해를선택하는대신 dominance 관계에따라선택하는방식으로 MOP를해결하도록할수있다. 이러한진화알고리즘의한갈래를다중목적함수진화알고리즘 (MOEA) 이라한다 [10]. MOEA 의장점으로는, 1) 문제의형태에종속되지않는일반적인알고리즘이므로어떠한문제에도적용할수있으며, 2) 해집합을통하여여러개의타협점들을한번에찾을수있고, 3) 진화알고리즘의확률적인 (stochastic) 탐색특성으로인하여지역최적해에서보다쉽게빠져나올수있다는점을들수있다. 2.3 다중목적함수최적화문제로서의프로브디자인일반적으로좋은프로브가갖추어야할조건은앞서 1장에서언급한것과같이크게다음의세가지를들수있다. 1) 목표유전자이외의유전자와의상보결합최소화. 2) 프로브자체의 2차구조형성최소화. 3) 녹는점등의화학적인성질의균질성최대화. 이상의조건이외에도기본적으로프로브의정의상프로브와상보적인서열이목표유전자이외에존재해서는안된다는제한조건이추가된다. 이조건은첫번째조건에포함될수있으나, 프로브디자인문제의정의를보다명확히하기위하여제한조건으로추가하였다. 그런데이중에서세번째조건은대개다른두조건과동등한목적함수라기보다는별도의제한조건으로취급되는경우가많다. 따라서본연구에서세번째조건은목적함수가아니라 MOEA를통해서생성한여러타협점중에서선택할때의척도로사용하였다. 따라서프로브디자인은하나의제한조건과두목적함수를가진 MOP로정의될수있고, 이를보다명확하게수식화하면다음과같다. 목표유전자의집합 와이에대한프로브의집합 가있을때, 각각의목표유전자 에대한프로브 는 의부분서열에상보적이고 모두일정한길이 인문자열 로가정한다. 따라서프로브의정의에따른제한조건과두목적함수는다음과같다.,,. 여기서 는두서열 와 에서 에상보
다중목적함수진화알고리즘을이용한마이크로어레이프로브디자인 505 적인서열 가 의부분서열로포함되어있을때만 1 의값을가지고그외는 0의값을갖는지시함수 (indicator function) 이며, 는두서열 와 가가장안정적으로결합하기위해필요한자유에너지를나타낸다. 마찬가지로 는한서열 가이루는 2차구조중에서가장안정적인구조에필요한자유에너지를나타낸다. 위의목적함수값은 OligoArray[2] 에사용된변형된 Mfold[12] 프로그램을이용하여계산하였고 NN 모델 (nearest neighbor model)[13] 에기반하였다. 이들자유에너지의값이작을수록안정적인결합이라할수있는데, 와 는목표유전자가아닌유전자와의결합이나스스로의 2차구조등프로브의기능을저해할수있는결합의안정성을표현하므로이들함수의값이클수록우수한프로브이다. 이상에서와같이프로브디자인은 이라는제한조건을만족시키면서두목적함수 와 를최대화시켜야하는 MOP로간주할수있다. 또한이두목적함수는실제유전자서열을대상으로테스트한결과어느정도상충되는관계에있음이관찰되어 [9], 프로브디자인이 MOP의정의에잘맞음을보여준다. 3. 다중목적함수진화알고리즘을이용한프로브디자인 2.3절에서의정의를바탕으로 MOEA를프로브디자인문제에적용하였다. MOEA를이용한프로브디자인과정은그림 2에서와같이크게 3단계를거치는데, 첫번째단계는적절한프로브탐색영역을찾는전처리과정이고, 두번째단계는찾아진영역에서 ε -다중목적함수진화알고리즘 (ε -MOEA, ε -multiobjective evolu- 그림 2 프로브디자인과정 tionary algorithm) 를이용해서프로브를탐색하는최적화과정이며, 마지막단계는 MOEA의결과로나온다양한프로브집합들중에서우선적으로사용될수있는프로브를추천해주는후처리과정이다. 또한이러한전체과정을 web상에서처리할수있는온라인인터페이스 (EvoOligo, http://cbit.snu.ac.kr/evooligo) 를제공한다. 3.1 프로브영역선택을위한전처리과정목표유전자서열의어느부분이나프로브가될수있지만, 효율적인탐색을위하여보다적합한프로브서열을찾을수있을것으로예상되는영역내로탐색을제한하였다. 이를위하여목표유전자들의서열들을다중정렬 (multiple alignment) 하여공통적으로보존되어있는영역은제외하고, 각각의목표유전자만의특색이있는영역을찾는전처리과정을도입하였다. 이것은목표유전자들사이에연관성이있어서서열상에공통점이있다는것을전제로하는데, 현재유전체전체의발현분석이가능한상용마이크로어레이서비스가존재하는것과는별개로, 특정유전자군에대한추가분석을위한프로브디자인이필요한경우도많다는것을감안한가정이다. 각각의목표유전자별로특색이있을것으로예상되는영역을찾기위하여전처리과정에서는중합효소연쇄반응 (PCR, polymerase chain reaction) 에필요한프라이머 (primer) 디자인도구인 Primer3[14] 를사용하여후보 PCR 영역들을찾은다음이를 ClustalW[15] 에서다중정렬을통해생성한대표서열과비교하여적절한영역을선택하여이를프로브탐색영역으로사용하였다. 여기서 PCR 영역을사용하는것은 PCR 프라이머선택과정에서프라이머사이의영역이본래의서열이나프라이머와결합하는일을최소화하도록선택되므로프로브탐색에도적절하기때문이다. 그러나프로브탐색영역은다중정렬된서열과비교하여가능한한여러유전자에공통적인부분을피하도록선정하였다 ( 그림 3의 3)). 구체적인전처리과정은그림 3에자세히설명되어있다. 전처리단계에서이후의최적화단계에서탐색할영역이결정되므로, 필요한인자들 (,, ) 은목표유전자서열들사이의유사성을고려하여결정되어야한다. 3.2 ε - 다중목적함수진화알고리즘을이용한프로브최적화프로브최적화단계에서는다른 MOEA와의비교에서탁월한성능을보인 ε -MOEA[16] 를적용하였다. ε -MOEA는 steady-state 구조의진화알고리즘에기반한 MOEA로서 ε -dominance를이용한일종의군집화 (clustering) 효과와현재세대까지 non-dominate인해들을별도의집합 (archive) 으로보관하는방법을통하여
506 정보과학회논문지 : 소프트웨어및응용제 35 권제 8 호 (2008.8) 그림 3 프로브탐색영역선택을위한전처리과정 우수한성능을보였다 [17]. 여기서 ε -dominance[17] 란 ε -MOEA의가장큰특징으로서하나의해가다른해를 dominate하기위해서최소한 ε 만큼더나은목적함수값을가져야한다는것을의미한다. 가 를 ε -dominate한다는것은 로표현하고, 수식으로는다음과같다.,. 따라서이 ε -dominance 개념을이용하여 non-dominate하는해들을선택하면각각의해사이의거리가최소 ε 이상이되어목적함수공간에서지나치게가까운해들을제거해주는효과를볼수있다. 또한이렇게해서찾아진해들을별도의집합에보관하여여러세대에걸쳐진화에참여하게함으로써좋은형질을계속해서유전시킬수있도록하였다. ε -MOEA를사용한프로브최적화과정은그림 4와같다. 각세대마다개체군과 archive에서각각선택된부모해로부터생성 ( 그림 4의 3)) 된새로운해를기존의해들과비교하여더좋을경우에만 archive나개체군에반영 ( 그림 4의 5) 와 6)) 되도록되어있다. 또한프로브최적화에는제한조건이포함되어있으므로 5) 의단계에서 archive를갱신할때제한조건을적게위반하거나만족시키는해를우선적으로선택하도록하였다. 또한제한조건을만족시키는해들사이에서는 dominate되지않는해를우선시하였다. 따라서 archive의해들은제한조건을만족시키면서차츰 Pareto-최적해에가까워지게되며 ε -dominance의정의상목적함수공간의특정구간에밀집되지도않게된다. 3.3 최종프로브추천을위한후처리과정한번에다양한해를찾아준다는것이 MOEA의장점이지만, 비용이나시간등의제약으로인하여해집합의일부만필요할경우를위하여실제실험에우선적으로사용해볼만한프로브를추천해주는단계도도입하였다. EvoOligo에서는세가지기준을두고이에따라서로다른세프로브집합을선택하여사용자에게추천하도록하였다. 첫번째기준은 BLAT[18] 으로측정한비- 목표유전자와의서열유사도 ( ) 의최소화이고, 두번째는 NACST/Sim[19] 으로측정한비-목표유전자와의결합안정성 ( ) 을최소화하는것이다. 결합안정성의기준으로는프로브와비-목표유전자가결합했을때의녹는점을측정하여 보다높은지의여부를사용하였다. 마지막으로프로브집합내의프로브들의녹는점의균일도 ( ) 의최소화를기준으로이용하였다. 보다명확하게정의하면다음과같다.,,. 여기서 는 -MOEA의수행결과최종 archive의 번째프로브집합에서 번째프로브를의미하고, 는 BLAT을이용해서계산한두서열 와 의서열유사도로서큰값일수록높은유사도를의미한다. 와 는각각지시함수와표준편차를나타내며, 는두서열 와 가가장안
다중목적함수진화알고리즘을이용한마이크로어레이프로브디자인 507 그림 4 ε -MOEA 를사용한프로브최적화과정 정적으로결합했을때의가장높은녹는점을나타낸다. 4. 프로브디자인적용결과제시된방법의검증을위하여두가지프로브디자인문제에적용해보았다. 하나는고위험성인유두종바이러스 (HPV, Human Papillomavirus) 에대한프로브들을디자인하는것인데, ( 주 ) 바이오메드랩의협조를통하여관련연구 [20] 에서자궁암을일으킬확률이높은것으로분류된 19종의고위험성인유두종바이러스를선별하여사용하였다. 다른하나는 52종류의애기장대칼모듈린유전자군 (AtCaMs/AtCMLs, Arabidopsis Calmodulin multigene family) 에대한프로브디자인이다. 애기장대칼모듈린유전자군은 6종류의애기장대칼모듈린 (Arabidopsis thaliana Calmodulins) 과 46종류의애기장대유사칼모듈린 (Arabidopsis thaliana Calmodulin like proteins) 으로구성되어있는데 [21], 기능이나발현상의특징에차이가있지만유전자서열의유사도는높은편이어서이들의분석을위해서는프로브디자인에주의를기울여야한다. 프로브의길이는 HPV의경우 30 염기, AtCaMs/ AtCMLs의경우에는 24 염기로지정하였다. ε -MOEA 에사용된인자로는교차와변이연산의확률을각각통 상적인값인 0.9와 0.01로놓았으며, ε 값으로는 1을사용하였다. 그리고개체군의크기와최대세대수는각각 50과 300으로정하였다. 전처리과정에사용된인자,, 값은여러실험을거쳐가장좋은성능을보이는값을선택하였는데, HPV의경우각각 0.05, 0.5, 0.9, AtCaMs/AtCMLs의경우 0.25, 0.0, 0.7로정하였다. 그외에 Primer3, ClustalW, Mfold, NACST/Sim, BLAT, BLAST에필요한인자는기본값을사용하였다. 성능평가를위하여널리알려진프로브디자인프로그램인 OligoArray[2] 와 OligoWiz[5] 를사용하여디자인한프로브와 EvoOligo에서최적화한프로브를비교분석하고, 전처리과정의효과에대하여살펴보았다. 4.1 프로브의비교평가를위한척도각기다른두방법으로생성된두개의프로브집합을정확하게비교하여평가하는것은어려운문제이다. 실제마이크로어레이실험을통하여평가하는것이가장신뢰성있는방법이겠으나, 실험시의환경조건등에영향을받기쉽다는단점도있다. 그렇지만프로브의길이가 50 또는 70 염기인경우 [6] 와 20 염기인경우 [7] 에대한경험적인지침을다룬문헌에따르면목표유전자와의결합, 비-목표유전자와의결합및 2차구조에필요한자유에너지와서열유사도가마이크로어레이상에
508 정보과학회논문지 : 소프트웨어및응용제 35 권제 8 호 (2008.8) 서프로브의신호세기와어느정도연관을갖는것으로나타났다. 본연구에서는이러한결과를참조하여표 1과같은척도를사용하였다. BLAST match와 (non-target) 은모두프로브의특이성 (specificity) 에관한척도이다. BLAST match는비-목표유전자와프로브의서열유사도 ( ) 를계산하여특이성을판단하는데반하여, (non-target) 은프로브를디자인하고자하는유전자중프로브의목표유전자가아닌것과의결합에필요한자유에너지를계산하여특이성을측정하게된다. BLAST match의계산에는 BLAST 프로그램 [22] 을이용하였는데, HPV의경우전역데이타베이스로 NCBI(http://www.ncbi.nlm.nih.gov/) 의 "nr" 데이타를사용했으며, AtCaMs/AtCMLs의경우에는 TAIR[23] 을사용했다. (self structure) 는프로브자체의안정적인 2차구조형성정도를측정함으로서프로브의목표유전자에대한반응감수성 (sensitivity) 을알아본다. 마지막의 Tm( C) 은프로브의반응조건의유사도에대한척도이다. BLAST match는작은값일수록, (non-target) 과 (self structure) 는모두의도하지않은결합에필요한자유에너지값의평균을의미하므로값이클수록좋은프로브집합을의미한다. Tm( C) 에서녹는점의평균은높고표준편차는작을수록목표유전자와안정적으로결합하면서프로브들의반응조건이많이유사한좋은프로브집합을의미한다. 4.2 기존프로그램과의비교평가비교평가대상으로는널리잘알려진프로그램인 OligoArray[2] 와 OligoWiz[5] 를사용하였다. Oligo- Array에서는서열유사도와자유에너지계산후단순한필터링방법을통하여적당한프로브를찾아내는반면, OligoWiz에서는 dynamic programming 방법을사용하였다. 디자인할프로브의길이와프로브의녹는점범위등은 EvoOligo와같도록하고, 그외의설정은두프로그램모두기본값을사용하도록하였다. EvoOligo에서생성된최종프로브집합중, NACST/Sim으로측정한비-목표유전자와의결합안정성이가장낮은 ( 즉, 를 최소화하는 ) 프로브집합과비교한결과만제시되었으나, 다른기준으로선택된프로브집합과의비교역시같은양상을보였다. 또한 HPV의경우 ( 주 ) 바이오메드랩에서전문가가디자인하여상용바이오칩에사용된프로브집합 [9] 도비교대상에포함하였다. 표 2의비교평가결과와같이 HPV나 AtCAMs/ AtCMLs의두경우모두에서 EvoOligo를이용해서디자인한프로브가프로브의반응특이성이나감수성측면에서다른두프로그램이나전문가가디자인한프로브 (Biomedlab) 보다더좋은것을알수있다. 유전자서열의차이가보다큰편인 HPV의경우, 비-목표유전자와의결합또는 2차구조에필요한자유에너지뿐만아니라 BLAST를통한서열비교에서도 EvoOligo와다른프로그램이나전문가가디자인한프로브와의차이가확연히남을알수있다. 반면에유전자서열의유사도가높은편인 AtCaMs/AtCMLs의경우, EvoOligo에서생성한프로브가여전히가장좋은것으로평가되기는했지만, 다른프로그램과의차이는 HPV의경우보다적었는데, 이는 AtCaMs/AtCMLs에대한프로브디자인이쉽지않은문제임을보여준다. 그러나녹는점의유사도면에서는 EvoOligo 쪽의표준편차가더크게나타나는데, 다른두프로그램에서녹는점을균일하게맞추기위하여우선필터링을한다음프로브최적화를한다는것과, EvoOligo의프로브최적화과정에서녹는점의균일도가목적함수에포함되지않았다는점을감안하면우수한성능이라고할수있다. 또한자유에너지값이나녹는점은어느정도서열의길이에비례하므로, 프로브의길이가 30인 HPV의경우에서길이가 24인 AtCaMs/AtCMLs의경우보다전반적으로자유에너지값이더낮고녹는점이더높게나오는경향이관찰되었다. 4.3 전처리단계에대한분석마지막으로프로브탐색영역을결정하는전처리단계의영향을알아보기위하여전처리과정을거쳐서디자인한프로브와그렇지않은프로브를비교하였다 ( 표 3). 또한, HPV의경우, 전문가들이경험적으로선택한프로브탐색영역 (L1 영역 ) 이존재하는데, ( 주 ) 바이오메 표 1 프로브의비교평가에사용된척도 이름정의설명 BLAST match BLAST를이용하여계산한전역데이타베이스에서의 비-목표유전자와프로브의서열유사도. (non-target) Mfold 를이용하여계산한비 - 목표유전자와프로브의안정적인결합정도. (self structure) Mfold 를이용하여계산한프로브의안정적인 2 차구조형성정도. Tm ( C) 과 프로브의녹는점들의평균과표준편차.
다중목적함수진화알고리즘을이용한마이크로어레이프로브디자인 509 표 2 디자인된프로브의비교결과. Biomedlab으로표시된항목은 [9] 에서비교에사용된전문가가디자인한프로브 집합을의미한다. 프로브의녹는점은나트륨이온농도 0.5M를가정하고 NN 모델을사용하여계산한값이다. EvoOligo OligoArray OligoWiz Biomedlab 2 3 4 4 BLAST match HPV -21.2801-25.9314-27.6639-21.8649 (non-target) (kcal/mol) -0.8545-2.4670-2.2126-2.0208 (self structure) (kcal/mol) 67.49±3.22 73.59±2.72 72.93±1.27 73.35±5.00 Tm( C) 4 4 5 BLAST match AtCaMs/AtCMLs -19.5756-20.5116-21.2913 (non-target) (kcal/mol) -0.7454-1.3851-0.8898 (self structure) (kcal/mol) 68.33±4.00 69.73±2.57 70.63±2.10 Tm ( C) 표 3 프로브탐색영역의변화에따른프로브성능비교. Evo는 EvoOligo의전처리과정에서선정한영역에서디자인한프로브, WS는전처리과정을거치지않고전체유전자서열에서프로브탐색을수행한경우, L1은 HPV 서열의 L1 영역에서디자인한프로브를각각의미한다. HPV AtCaMs/AtCMLs EvoOligo Evo EvoOligo L1 EvoOligo WS Biomedlab L1 2 7 0 4 BLAST match -21.2801-22.6492-22.9594-21.8649 (non-target) (kcal/mol) -0.8545-2.4497-1.0747-2.0208 (self structure) (kcal/mol) 67.49±3.22 71.30±2.44 67.78±3.85 73.35±5.00 Tm ( C) 4 3 BLAST match -19.5756-18.5346 (non-target) (kcal/mol) -0.7454-0.9746 (self structure) (kcal/mol) 68.33±4.00 66.81±4.49 Tm ( C) 드랩의상용바이오칩에서이영역을대상으로하여디자인한프로브 ( 표 3의 Biomedlab L1) 와도비교하였다. 표 3에서프로브탐색영역이어떻게정해지는지에따라최종적으로디자인된프로브의성능에많은차이가있음을알수있다. HPV에서전문가가선정한영역 (L1 영역 ) 에서는, EvoOligo를이용하여디자인한프로브가녹는점의편차는더작지만, 반응특이성이나감수성측면에서약간뒤떨어지는것으로나타났다. 그러나 EvoOligo 자체의전처리과정에서선택한영역에서디자인한프로브가다른프로브보다모든면에서더좋은값을나타내는것을알수있다. 그리고 AtCaMs/ AtCMLs의경우, 전체유전자서열을사용한경우가반응특이성이더좋지만, 감수성이나녹는점의편차가더큰것으로나타났다. AtCaMs/AtCMLs의경우유전자서열의전체적인다양성이 HPV의경우보다커서전처리단계에서좋은프로브탐색영역을찾기가어려웠는데, 이런특성때문에좋은프로브서열을찾을수있는영역을놓쳐서전처리단계를거친프로브의반응특이성이오히려낮아진것으로생각된다. 또한, 전처리단계에서프로브탐색영역을제한함으로인해서 ε -MOEA의탐색공간이줄어들어수렴속도를빠르게하는효과를관찰할수있었다 ( 그림 5). 그 그림 5 전처리단계에따른수렴속도의변화비교림 5는 AtCaMs/AtCMLs의경우각각의세대에서 의최소값의변화를관찰하여전처리단계를사용했을때와그렇지않을때를비교한것인데, 전처리를사용한경우가그렇지않은경우보다더우수한프로브를더빨리발견할수있음을알수있다. 이나 의경우에도비슷한경향성이나타났으며, HPV에서도마찬가지의결과를보였다. 5. 결론본연구에서는생물정보학분야에서중요한문제의
510 정보과학회논문지 : 소프트웨어및응용제 35 권제 8 호 (2008.8) 하나인마이크로어레이실험에사용될프로브의최적화문제를다중목적함수최적화문제로정의하고, 이를다중목적함수진화알고리즘인 ε -MOEA를사용하여해결하는방안을제시하였으며, 19종류의고위험군인유두종바이러스유전자와 52종류의애기장대칼모듈린유전자군에대한프로브디자인문제에서널리알려진다른프로브최적화도구와비교하여제시된방법의성능을입증하였다. 또한제시된방법은온라인도구인 EvoOligo(http://cbit.snu.ac.kr/EvoOligo) 로구현되어관련연구자가편리하게접근하여사용할수있도록하였다. 진화알고리즘의특성상프로브의최적화에약간의시간이걸리기때문에디자인된프로브는사용자에게 e-메일로전달된다. 또한, EvoOligo는마이크로어레이의통합분석 / 디자인플랫폼인 DNAChipBench의한부분인 ProbeBench로도사용되었다 (http://cbit.snu.ac. kr/~dnachipbench/). 현재 EvoOligo에서는전처리단계에서여러인자를사용자가지정해줄필요가있는데, 만약이러한인자를입력된유전자서열의분석을통해서자동적으로지정할수있다면더욱편리할것이다. 따라서이후추가연구의방향은프로브디자인에필요한여러인자의자동화방안을연구하는것이될것이다. 참고문헌 [1] Gordon, P. M. K. and Sensen, C. W., "Osprey: a comprehensive tool employing novel methods for the design of oligonucleotides for DNA sequencing and microarrays," Nucleic Acids Research, Vol.30, No.17, pp.e133, 2004. [ 2 ] Rouillard, J.-M., Zuker, M. and Gulari, E., "OligoArray 2.0: design of oligonucleotide probes for DNA microarrays using a thermodynamic approach," Nucleic Acids Research, Vol.31, No.12, pp.3057-3062, 2003. [3] Drmanac, S., Stravropoulos, N. A., Labat, I., Vonau, J., Hauser, B., Soares, M. B. and Drmanac, R., "Gene-representing cdna clusters defined by hybridization of 57,419 clones from infant brain libraries with short oligonucleotide probes," Genomics, Vol.37, No.1, pp.29-40, 1996. [4] Herwig, R., Schmitt, A. O., Steinfath, M., O'Brien, J., Seidel, H., Meier-Ewert, S., Lehrach, H. and Radelof, U., "Information theoretical probe selection for hybridisation experiments," Bioinformatics, Vol.16, No.10, pp.890-898, 2000. [ 5 ] Wernersson, R. and Nielsen, H., "OligoWiz 2.0- integrating sequence feature annotation into the design of microarray probes," Nucleic Acids Research, Vol.33, Web Server issue, pp.w611-w615, 2005. [6] He, Z., Wu, L., Li, X., Fields, M. W. and Zhou, J., "Empirical establishment of oligonucleotide probe design criteria," Applied and Environmental Microbiology, Vol.71, No.7, pp.3753-3760, 2005. [ 7 ] Matveeva, O. V., Shabalina, S. A., Nemtsov, V. A., Tsodikov, A. D., Gesteland, R. F. and Atkins, J. F., "Thermodynamic calculations and statistical correlations for oligo-probes design," Nucleic Acids Research, Vol.31, No.14, pp.4211-4217, 2003. [8] Wu, C., Carta, R. and Zhang, L., "Sequence dependence of cross-hybridization on short oligo microarrays," Nucleic Acids Research, Vol.33, No.9, pp.e84, 2005. [9] Shin, S.-Y., Lee, I.-H. and Zhang, B.-T., "Microarray probe design using ε -multi-objective evolutionary algorithms with thermodynamic criteria," Lecture Notes in Computer Science (EvoBio 2006), Vol.3907, pp.184-195, 2006. [10] Deb, K., Multi-Objective Optimization using Evolutionary Algorithms, John Wiley & Sons, Ltd., 2001. [11] B ck, T., Evolutionary Algorithms in Theory and Practice, Oxford University Press, 1996. [12] Zuker, M., "Mfold web server for nucleic acid folding and hybridization prediction," Nucleic Acids Research, Vol.31, No.13, pp.3406-3415, 2003. [13] SantaLucia, J. Jr., "A unified view of polymer, dumbbell, and oligonucleotide DNA nearest-neighbor thermodynamics," Proceedings of the National Academy of Sciences of the United States of America, Vol.95, No.4, pp.1460-1465, 1998. [14] Rozen, S. and Skaletsky, H., "Primer3 on the WWW for general users and for biologist programmers," Methods in Molecular Biology, Vol.132, pp.365-386, 2000. [15] Chenna, R., Sugawara, H., Koike, T., Lopez, R., Gibson, T. J., Higgins, D. G. and Thompson, J. D., "Multiple sequence alignment with the Clustal series of programs," Nucleic Acids Research, Vol.31, No.13, pp.3497-3500, 2003. [16] Deb, K., Mohan, M. and Mishra, S., "A fast multi-objective evolutionary algorithm for finding well-spread Pareto-optimal solutions," Kanpur Genetic Algorithm Laboratory, Indian Institute of Technology Kanpur, KanGAL Report 2003002, 2003. [17] Laumanns, M., Thiele, L., Deb, K. and Zitzler, E., "Combining convergence and diversity in evolutionary multiobjective optimization," Evolutionary Computation, Vol.10, No.3, pp.263-282, 2002. [18] Kent, W. J., "BLAT-the BLAST-like alignment tool," Genome Research, Vol.12, No.4, pp.656-664, 2002. [19] Shin, S.-Y., Jang, H.-Y., Tak, M.-H. and Zhang, B.-T., "Simulation of DNA hybridization chain reaction based on thermodynamics and artificial chemistry," Preliminary Proceedings of 9th Inter-
다중목적함수진화알고리즘을이용한마이크로어레이프로브디자인 511 national Meeting on DNA Based Computer, pp.451, 2004. [20] Walboomers, J. M. M., Jacobs, M. V., Manos, M. M., Bosch, F. X., Kummer, J. A., Shah, K. V., Snijders, P. J. F., Peto, J., Meijer, C. J. L. M. and Munoz, N., "Human papillomavirus is a necessary cause of invasive cervical cancer worldwide," The Journal of Pathology, Vol.189, No.1, pp.12-19, 1999. [21] McCormack, E., Tasi, Y.-C. and Braam, J., "Handling calcium signaling: Arabidopsis CaMs and CMLs," Trends in Plant Science, Vol.10, No.8, pp.383-389, 2005. [22] Altschul, S. F., Madden, T. L., Schaffer, A. A., Zhang, J., Zhang, Z., Miller, W. and Lipman, D. J., "Gapped BLAST and PSI-BLAST: a new generation of protein database search programs," Nucleic Acids Research, Vol.25, No.17, pp.3389-3402, 1997. [23] Rhee, S.Y., Beavis, W., Berardini, T. Z., Chen, G., Dixon, D., Doyle, A., Garcia-Hernandez, M., Huala, E., Lander, G., Montoya, M., Miller, N., Mueller, L. A., Mundodi, S., Reiser, L., Tacklind, J., Weems, D. C., Wu, Y., Xu, I., Yoo, D., Yoon, J. and Zhang, P., "The Arabidopsis Information Resource (TAIR): a model organism database providing a centralized, curated gateway to Arabidopsis biology, research materials and community," Nucleic Acids Research, Vol.31, No.1, pp.224-228, 2003. 조영민 2001 년서울대학교컴퓨터공학부학사. 2006 년 ~ 현재 University of California, San Diego Computer Science 박사과정 양경애 1998 년제주대학교 ( 학사 ). 2000 년제주대학교 ( 석사 ). 2005 년경상대학교 ( 박사 ). 2006 년 ~ 현재서울대학교바이오정보기술연구센터 ( 박사후연수생 ) 장병탁정보과학회논문지 : 소프트웨어및응용제 35 권제 6 호참조 이인희 2001년 2월서울대학교컴퓨터공학부학 사. 2001년 3월~현재서울대학교컴퓨 터공학부석박사통합과정. 관심분야는 진화연산, 생물정보학, 기계학습, DNA 컴퓨팅 신수용 1998년 2월서울대학교컴퓨터공학부학사. 2000년 2월서울대학교컴퓨터공학부석사. 2005년 8월서울대학교전기, 컴퓨터공학부박사. 2006년 4월~2008년 3 월 NIST 연구원. 2008년 4월~현재서울대학교병원의료정보센터연구교수. 관심분야는진화연산, data mining, 생물정보학, 의료정보학, 기계학습, DNA 컴퓨팅