biosequence의 전환: gb_to_gff.pl과 Readseq

주로 세균 유전체만을 다루어 오다가 효모 유전체로 범위를 넓힌지도 꽤 되었다. 내가 즐겨사용하는 효모 유전체 자동 주석화 서비스로는 Yeast Genome Annotation Pipeline(YGAP)이 있다. Contig 서열만 넣어주면 매우 신속하게 분석 결과를 제공해 준다. 잘 정리된 효모 유전체만을 비교 대상으로 하기에 속도가 빠를 수밖에 없고, 현재와 같은 효모 속/종 분화기 이루어지기 전의 원시 유전자에 대한 대응 정보도 산출하여 준다.

여기서는 자체 format의 annotation 정보 파일 및 GenBank/EMBL 파일도 결과물로서 제공한다. 문제는 tophat 등의 후속 분석에 쓰려면 GTF 혹은 GFF 파일이 필요한데, GenBank에서 이를 만들어 내려면 상당히 번거로운 스크립트를 짜야만 한다.

BioPerl을 활용하여 내가 만든 스크립트는 GenBank에서 CDS 혹은 translated sequence를 FASTA 형태로 출력하는 것, 혹은 FASTA 형식의 서열 파일을 여러가지로 조작하는 것이 전부였다. GenBank 파일을 GFF로 변환하는 도구는 무엇이 있을까? artemis와 같은 GUI tool에서 Save An Entry As 기능을 쓰면 안되는 것은 아니지만 한번에 여러개의 파일을 처리하기가 곤란하다. command line interface에서 일괄작업을 돌릴 수 있는 유틸리티가 필요하다.

내가 가장 즐겨쓰는(솔직하게 말하자면 쓸 줄 아는 유일한) 도구인 BioPerl을 사용한 방법이 있을 것으로 생각하고 검색을 해 보았다. bioperl-live/exmple/tools에 gb_to_gff.pl이라는 예제가 있다.

#!/usr/bin/perl
use strict;
use Bio::Tools::GFF;
use Bio::SeqIO;
my ($seqfile) = @ARGV;
die("must define a valid seqfile to read") unless ( defined $seqfile && -r $seqfile);
my $seqio = new Bio::SeqIO(-format => 'genbank',
                           -file   => $seqfile);
my $count = 0;
while( my $seq = $seqio->next_seq ) {
    $count++;
    # defined a default name
    my $fname = sprintf("%s.gff", $seq->display_id || "seq-$count");
    my $gffout = new Bio::Tools::GFF(-file => ">$fname" ,
                                     -gff_version => 1);
    foreach my $feature ( $seq->top_SeqFeatures() ) {
        $gffout->write_feature($feature);
    }
}

대단히 단순한 스크립트다. 실행을 해 본 결과 YGAP가 생성한 GenBank 파일에 대해서 잘 작동한다. $seq->top_SeqFeatures() 메서드가 핵심 기능을 하는 것으로 생각된다. 그동안 내가 사용하던 $seq->SeqFeatures()와는 무엇이 다른지 공부를 해야 되겠다. 

이보다는 좀 더 세련된 형태의 전환기는 없을까? 조금 더 검색을 해 보니 Readseq라는 유틸리티가 나온다. Java로 만들어졌으며 GUI 환경도 갖추고 있으나 개발된지 너무 오래되어 GFF3는 지원하지 않는 것으로 보인다. 완벽한 공용 소프트웨어로서 사용에 제한이 없다.

Biological sequence에 수반되는 feature를 기술하고 교환하는 포맷이 정해져 있다 하여도 실제로는 많은 변형(방언? 사투리?)이 존재하여 완벽한 호환을 어렵게 한다. 이를 감안하여 되도록 표준안을 준수하는 버릇을 들여야 될 것이다.