Bioinformatics Q and A > Beginning Bioinformatics for Biologists

BioPerl이란 이곳에 있는 library를 말합니다. BioPerl이 어떤 역할을 하는 것인지 실제 예를 들어서 한 번 설명을 해봅니다.

여러분이 하고자 하는 생물정보학적인 일은 당연히 어떤 연구에 수반되는 여러 가지 과정의 일부분일 것입니다. 또한, 그 연구에서 필요로 하는 생물정보학적인 일들 그 자체도, 다시 이런저런 다양한 작업들로 구성되어 있을 것입니다. 그 중에서 하나로, 예를 들어 GenBank에서 mammal의 hemoglobin의 subunit들(또는 자신의 연구에 필요한 다른 gene)에 대한 mRNA sequence들을 구해야 하는 일이 있다고 합시다. 이것은 http://www.ncbi.nlm.nih.gov에 가면 나오는 Search 기능에서, database는 Nucleotide로 선택을 하고, for라고 되어 있는 부분 다음에 나오는 입력 박스에 아래를 copy-and-paste 해주면 됩니다. (이처럼 web browser를 마우스 클릭으로 사용하는 방식도 경우에 따라서는 더 편리하고 자연스러운 방식이 되는 경우도 많습니다. 이 방식 그 자체가 문제가 있다는 것이 아니라, 이것만으로는 전혀 충분치가 않다는 점이 문제입니다.)

Mammalia[Organism] AND hemoglobin[Title] AND ((((1900[MDAT] : 3000[MDAT] NOT gbdiv_sts[PROP]) NOT srcdb_tpa_ddbj/embl/genbank[PROP]) NOT gbdiv_htg[PROP]) NOT gbdiv_pat[PROP]) AND biomol_mrna[PROP] AND gene_in_genomic[PROP] AND srcdb_refseq[PROP]

여러분이 하고자 하는 생물정보학적인 일은 당연히 어떤 연구에 수반되는 여러 가지 과정의 일부분일 것입니다. 또한, 그 연구에서 필요로 하는 생물정보학적인 일들 그 자체도, 다시 이런저런 다양한 작업들로 구성되어 있을 것입니다. 그 중에서 하나로, 예를 들어 GenBank에서위의 쿼리문을 어떻게 작성을 했는가 하는 것은 다른 글에서 자세히 적겠습니다. (작성이 된 다음에는 이곳에 링크가 될 것입니다.) 물론 위의 옵션이나 키워드 등을 다 외워야 하는 것이 아니고, 마우스 클릭으로 선택을 해주면 만들어지도록 되어 있습니다. 지금 자세한 세부에 대한 설명은 생략을 하고, 이렇게 하면 mammal에서 hemoglobin의 chain들에 대해 현재 밝혀진 mRNA sequence들 (이 글을 쓰고 있는 현재로서는) 69개가 잡히게 됩니다. 이제 무엇을 해야 할까요? 물론 sequence를 하나씩 클릭을 해서 일일이 파일로 저장을 하고, 그 각각을 일일이 눈으로 들여다보고 하는 방식도 가능할 것입니다. 이 글에서는 당연히 이와는 다른 방식에 대해서 설명을 합니다.

우선, 위의 검색 결과 화면에서 'Show'라고 되어 있는 부분을 500으로 선택을 해줍니다. 그런 다음 각 record의 설명을 훑어봅니다. 분명히 자신이 원하지 않는 것들도 소수 포함이 되어 있을 것입니다. 훑어보면서 원하는 것은 record 번호 앞에 있는 네모를 클릭해서 주욱 선택을 합니다. (사실 "전체 선택" 기능이 있어서 일단 전체를 선택을 하고, 원하지 않는 것만 unselect를 하는 방식도 필요할 터인데, 이 점은 조금 불만이군요.) 그리고는, 'Display"라고 되어 있는 부분을 'GenBank"로 선택을 해줍니다. 이렇게 하면 엄청나게 큰 파일이 웹 브라우저 상으로 전송이 되기 시작하는데, 컴퓨터의 메모리가 적은 경우에는 컴퓨터의 반응 속도가 심하게 느려지는 문제 같은 것까지 생겨나게 될 것입니다. 다 전송되기를 기다리지 말고, 곧장 'Sent to"라고 되어 있는 부분을 File로 선택을 해줍니다. (이렇게 어색하게 하지 않아도 되도록 Interface를 만들어 놓을 수도 있을 터인데, 이점 역시 조금 불만이군요.) 그리고는 파일명을 입력해서 자신의 하드 디스크에 저장을 하면 드디어 mammal들의 hemoglobin의 subunit(또는 연구에 필요한 다른 gene)들에 대한 mRNA sequence들이 구해진 것입니다. (물론 대부분의 경우에는 이대로가 자신이 원하는 것은 아닐 것입니다. 예를 들어 splicing variant는 제거를 하는 등과 같은 여러 가지 고려들이 더 필요할 것입니다.)

'Display' 선택 박스를 보면, GenBank의 sequence 데이터를 상당히 다양한 format으로 가져올 수 있다는 것을 알 수 있을 것입니다. 이 중에서 XML이라고 되어 있는 것이 사실은 가장 현대적인 것이고 전산적으로도 가장 타당한 format인데, 여전히 GenBank format도 많이 쓰입니다. GenBank format은 소위 flat file format이란 것으로, 마치 노트 정리를 항목 구분을 깔끔하게 해서 적어놓는 것과 유사한 방식으로 기록을 해놓은 것입니다. 여기에 대한 자세한 설명은 역시 다른 곳에 적어놓도록 하겠습니다. 이렇게 하면 여러 개의 GenBank 파일들이 연결이 되어서 하나의 상당히 큰 파일이 얻어지게 됩니다. 검색 결과에서 상위 6개의 record만 선택을 해서 이곳에 저장을 해놓았습니다.

자, 이제 문제는 이 파일을 가지고 무엇을 어떻게 할 것인가 하는 점일 터인데요. 일단 "무엇을"에 해당하는 것은, 잘 작동하는 것으로 이미 증명이 되어 있는 PCR primer를 고른다거나 (STS라고 되어 있는 부분이 바로 이 목적으로 사용할 수 있는 정보입니다.) 3' UTR sequence들만 좍 뽑아내고 싶다거나 하는 일이 가능할 것입니다. 이 글에서는 전자를 어떻게 하는지에 대해서 설명을 해봅니다. 이 일을 위해서, 눈으로 일일이 들여다보면서 손으로 원하는 정보를 선택을 해서 또 다른 파일에 copy-and-paste 등으로 옮기고 하는 방식도 물론 가능할 것입니다. 지금 우리는 당연히 이렇게 하자는 것이 아닌 것이고요.

무언가 제대로 컴퓨터의 능력을 활용하는 것이 되어야 할 터인데, 이를 위해서 가장 우선적으로 필요한 것은 무엇일까요? 컴퓨터가 이 파일의 내부를 제대로 항목을 구분해 가면서 (마치 사람이 읽는 것처럼) 읽을 수 있게 해야겠죠? 이렇게 하는 것을 전산 용어로는 "parsing을 한다"라고 합니다. 저도 아주 오래 전에 (심지어) C로 이걸 짜본 적이 있었는데, 정말 엄청난 노가다에다가 한 마디로 보통 일이 아닙니다. 온갖 까다로운 점들이 다 생기게 되죠. 설마 이런 것을 지금 이 글을 읽고 있는 사람들이 각자 다 만들어내야 한다는 소리는 결코 아니겠죠? 그렇다면 해결책은 뭘까요? 누가 대표로 나서서 노가다를 해놓고, 그걸 다들 이용을 하면 되겠죠? 바로 이곳에 나오는 사람들이 그 대표주자들인 사람들입니다. 우린 이 분들이 만들어 놓은 것을 늘 감사하는 마음으로 이용만 하면 되는 것이죠.

일단 BioPerl이라는 강력한 무기를 손에 든 상태에서는, 그 다음부터의 일은 상당히 수월합니다. 아래 Perl script는 위에서처럼 저장한 GenBank 파일에서, 그 안에 들어 있는 모든 sequence들에 대해서, 어떤 feature를 가진 subsequence들을 좍 뽑아주는 것입니다. 여기에 있는 output 파일을 살펴보면, 뭘 하는 것인지 감이 올 것입니다.

#!/usr/bin/perl -w
use strict;
use Bio::SeqIO;
use Bio::Location::Simple;

my $gb_file = shift;
my $feature = shift;

my $in = new Bio::SeqIO(-format => 'genbank', -file => $gb_file);

while( my $seq = $in->next_seq ) {
    my $sequence = $seq->seq();
    foreach my $f ( $seq->top_SeqFeatures ) {
	if ($f->primary_tag eq $feature) {
	    my $desc = $f->primary_tag . " " .
	      ($f->strand < 0 ? 'complement' : '')
		. "(" . $f->start . ".." . $f->end . ")";
	    my $location = Bio::Location::Simple->new(-start  => $f->start,
						      -end   => $f->end,
						      -strand => $f->strand);
	    my $out = Bio::Seq->new(-seq => $seq->subseq($location),
				    -id => $seq->id(),
				    -desc => $desc,
				    -alphabet => "dna");
	    my $seqio_obj = Bio::SeqIO->new(-format => 'fasta' );
	    $seqio_obj->write_seq($out);
	}
    }
}