시퀀스에서 Di-nucleotide frequency 계산 문제.

지난해 pung96님의 블로그에 올라왔던 Di-nucleotide composition frequency 세기 문제( http://perlog.pung96.net/20 ) 를 다시 한번 언급해 본다.

문제를 다시 정리해 보면,

AGATAGCGATAGCG
AGATGACGATAGAG
...

위 처럼 DNA 서열이 담긴 input 파일이 있을 때, DNA base 두개씩을 끊은 dinucleotide ( e.g., AG, CT ) 와 같은 16개 조합( DNA base가 4개 이므로 2개의 조합 가지 수는 16개 )의 빈도를 계산하는 것이 원글의 문제다. 이때, 앞에서 순서대로 2개씩 읽되, window size는 1로 하여 1칸씩 옮겨가며 2개씩 읽어 조합을 계산해야 한다.

구체적으로 AGAA 네개 DNA 서열인 경우, AG, GA ,AA 세개의 dinucleotide 조합으로 계산해야 한다는 것이다.

TIMTOWTDI, perl의 특성답게 이를 해결하는 방법은 상당히 다양한데, pung96님의 글에 언급되지 않은 내용을 추가로 정리해 보고자 한다.

substr을 이용하는 경우 ( DNA 서열은 $seq에 저장되어 있는 경우 )

for( my $i=0; $i< length( $seq ) - 1 ;$i++ ){
   $count{ substr($seq,$i,2) }++;
}

정규식의 lookahead 이용하는 경우,

my @di=$seq=~/(?=(\w{2}))/g;
$count{$_}++ for @di;

이걸 한줄로 줄이면,

$count{$_}++ for $seq=~/(?=(\w{2}))/g;

반복을 while로 돌리면,

$count{$1}++ while  $seq=~/(?=(\w{2}))/g;

map 에 넣으면,

map{ $count{$_}++ } $seq=~/(?=(\w{2}))/g; 

와 같은 방식으로도 문제를 해결할 수 있다. 포스팅의 요지는 lookahead 를 이용한 정규식 사용!!

사족을 좀 덧붙여 보면, 이 문제는 연세대 생명공학과 대학원 이인석 교수의 bioinformatics 수업 과제 문제다. 지난해와 올해 연구소에 연대 대학원 학생분들이 있어, 내게 이 문제에 대한 조언을 구해와 본의 아니게 2년 연속으로 이 문제에 대해 생각해볼 기회가 생겨, 포스팅해 본다.