SeqIO

生物学でよく利用するファイル形式のほとんどは Bio::SeqIO モジュールで読み込める。読み込みだけでなく、他のファイル形式に変換して書き出しを行うフォーマット変換も可能である。

ファイルの読み込み

SeqIO->new でファイルの読み込みを行い、ファイルポインタを $inseq に保存し、以降 while 文を利用して $inseq からデータをひとつずつ取り出すことができる。

use strict;
use Bio::SeqIO;

#ファイルを開く
my $inseq = Bio::SeqIO->new(
                 -file => '<complete.1.genomic.gbff',
                 -format =>  'genbank'
            );


#マルチエントリーのファイルを1エントリーずつ読み込む
while (my $seq = $inseq->next_seq) {
	print $seq->accession_number, "\n";
}
#NT_112066
#NT_108597
#NT_086517
#NT_113343
# ...

ファイル中のデータをすべて配列に保存する場合は以下のようにする。

my $inseq = Bio::SeqIO->new(
                 -file => '<complete.1.genomic.gbff',
                 -format => 'genbank'
            );

my @seq_array = ();
while (my $seq = $inseq->next_seq) {
 push(@seq_array, $seq);
}

ファイルの書き出し

$seqio->write_seq のように利用してファイルの書き出しを行う。書き出しを行う前にファイル形式を指定する必要がある。

use strict;
use Bio::SeqIO;

my $inseq = Bio::SeqIO->new(
                 -file => '<complete.1.genomic.gbff',
                 -format => 'genbank'
            );

#書き出しするファイルのフォーマットとファイルの名前を設定する
my $outseq = Bio::SeqIO->new(
                 -file => '<complete.1.genomic.fa',
                 -format => 'fasta'
            );

while (my $seq = $inseq->next_seq) {
	#print $seq->accession_number;
	$outseq->write_seq($seq);
}

標準入力と標準出力

SeqIO では標準入力からもデータ情報を取り込むことができる。gzip などで圧縮されているファイルの場合、解凍しながら読み込めるので非常に役に立ちます。

use strict;
use Bio::SeqIO;

#標準入力からデータを取り込む
my $inseq = Bio::SeqIO->new(
                 -fh => \*STDIN,
                 -format =>  'genbank'
            );

#標準出力にデータを出力
my $outseq = Bio::SeqIO->new(
                 -fh => \*STDOUT,
                 -format =>  'fasta'
            );
while (my $seq = $inseq->next_seq) {
	$outseq->write_seq($seq);
}

コマンドラインで次のように実行すると、gzip で圧縮された Genbank 形式のファイルを解凍しながら FASTA 形式に変更できる。

zcat complete.1.genomic.gbff.gz | perl seqio.pl | grep ">"
>NT_112066 Callithrix jacchus genomic sequence, ENCODE region ENr231.
>NT_108597 Papio anubis genomic sequence, ENCODE region ENm002.
>NT_086517 Callithrix jacchus genomic sequence, ENCODE region ENm014.

SeqIOによるファイル作成

次のサンプルのように、取り込んだデータを FASTA 形式や GenBank 形式に変更し、生物種ごとにまとめることもできる。

use strict;
use Bio::SeqIO;

my $inseq = Bio::SeqIO->new(
                 -file   => "<$infile",
                 -format => 'Genbank',
            );
 
my %outfiles = (
	human => {
        	Genbank => Bio::SeqIO->new(
                                -file   => '>human.gb',
                                -format => 'Genbank',
                           ),
                Fasta   => Bio::SeqIO->new(
                                -file   => '>human.fa',
                                -format => 'Fasta',
                           ),
        },
        other => {
                Genbank => Bio::SeqIO->new(
                                -file   => '>other.gb',
                                -format => 'Genbank',
                           ),
                Fasta   => Bio::SeqIO->new(
                                -file => '>other.fa',
                                -format => 'Fasta',
                           ),
        }
);

 
while (my $seqin = $inseq->next_seq) {
    if ($seqin->species->binomial =~ m/Homo sapiens/) {
        $outfiles{'human'}->{'Genbank'}->write_seq($seqin);
        $outfiles{'human'}->{'Fasta'}->write_seq($seqin);
    } else {
        $outfiles{'other'}->{'Genbank'}->write_seq($seqin);
        $outfiles{'other'}->{'Fasta'}->write_seq($seqin);
    }
}

CDS のコドン出現頻度を計算うする方法

NC_020043.gb ファイルを読み込んで、feature が CDS となっている部分の塩基配列を取り出し、コドンの出現頻度を計算するサンプル。

SeqFeatures メソッドを利用することで、feature を一つずつ取得できる。また、SeqFeature オブジェクトに対して、splice_seq メソッドを呼び出すと、その feature の location 情報に基づいて、その部分だけを切り出しすることができる。

use strict;
use Bio::SeqIO;

my $counts = [];
my $i = 0;

my $seqio = Bio::SeqIO->new(
              -file => 'NC_020043.1.gb', 
              -format => 'genbank'
          );

my $seq = $seqio->next_seq;

for my $ft ($seq->get_SeqFeatures) {
  if ($ft->primary_tag eq 'CDS') {
    $counts->[$i++] = &count($ft->spliced_seq->seq);
  }
}



# 実際のカウントを行う関数(3文字ずつ)
sub count {
  my $seq = shift;
  my $cnt = {};
  for (my $i = 0; $i < length($seq); $i += 3) {
    my $pattern = substr($seq, $i, 3);
    if (exists $cnt->{$pattern}) {
      $cnt->{$pattern}++;
    } else {
      $cnt->{$pattern} = 1;
    }
  }
  return $cnt;
}