利用Bioperl进行序列输入和输出（Bioperl HOWTO翻译13）

利用Bioperl进行序列输入和输出

作者

版权

Ewan Birney拥有此文档版权。基于 Perl Artistic License协议有限共享。

基础

作者假定阅读此文章的读者没有用过BioPerl，也许是一个想要获取一些序列信息的生物学者；也许是一个想学习一下目前流行的所谓“生物新信息学”的IT人士。我们讲解的第一个脚本就是如何从包含一个或多个序列的序列文件中提取相关信息。

一个小建议：坚持使用Bio::SeqIO模块来处理序列。先来看一下完成任务所需要的前几行脚本：

#!/bin/perl

use strict;
use Bio::SeqIO;

my $file         = shift; # get the file name, somehow
my $seqio_object = Bio::SeqIO->new(-file => $file);
my $seq_object   = $seqio_object->next_seq;

为什么要用Bio::SeqIO？部分原因是SeqIO能识别各种不同的序列格式并从中提取然后创建相应的BioPerl对象。有些序列格式很简单，例如FASTA。fasta格式仅包含序列和一个简单的序列标识符，而不包含序列内部的一些详细的特征或信息（Feature-Annotation HOWTO里具体讲了序列的特征）。一般来说，SeqIO能识别不同的序列格式，并创建相应的对象，例如输入fasta，SeqIO会创建一个Bio::Seq对象，而输入Genbank或EMBL这种包含注释信息的序列文件，SeqIO则会创建一个Bio::Seq::RichSeq对象。

至于SeqIO针对不同格式的序列文件，到底创建了哪种对象，我们并不需要关心。SeqIO会把自动这些细节问题处理好。

十秒速览

生物信息学的工作会接触到各种格式的序列。而又有几乎同样多的程序来处理这些序列。SeqIO的优势是可以通过创建和输出序列对象的方式来处理很多不同格式的序列。可以将SeqIO想象成一个智能的文件控制系统。

背景

SeqIO系统可以处理很多复杂的序列格式。提供给Bioperl输入序列（序列可以是普通文件、标准输入输出（STDIN和STDOUT）、变量等等）以及序列格式，Bioperl即可从序列中提取所有可以提取的信息。有时候并不需要说明序列格式。SeqIO可通过文件扩展名或文件内容来猜测出序列格式。如果文件没有扩展名或者文件内容不完整甚至压根输入的就不是一个文件，Bioperl直接默认是fasta格式（译者注：而这样就很可能出错）。除非明确知道所用的序列就是fasta格式，最好还是明确告知Bioperl你所输入序列的格式。

SeqIO可以接受多种形式的输入。唯一的要求是序列必须包含在Perl所规定的“句柄”中（详见IO::Handle）。大多数人使用文件句柄或者标准输入/输出（STDIN/STDOUT）。Perl还提供了如何将一段字符串转为句柄（见下文）。这样一来，几乎所有东西都可以作为SeqIO的输入，只要其中含有序列信息。SeqIO的主要工作流程就是获取一个句柄，然后依次将一条条序列信息转化为SeqI对象，若有需要则创建一个句柄以输出序列文件。SeqIO可以识别多个序列之间的分隔符，如genbank格式中的“// ”，fasta格式中的“>”，也能识别用关键字或符号标识的其他序列信息。

格式

BioPerl的SeqIO模块能识别很多序列格式并且互相转换. 表1是最新版本的SeqIO所支持的格式。（译者注：表格见原文。如前所述，SeqIO将输入的不同格式序列用相应的模块转换成相应的序列对象）

注: Bio::SeqIO需要bioperl-ext包来和Staden包中的io_lib库才能读取scf, abi, alf, pln, exp, ctf,和ztr格式的序列.

也许有些人会对上面提到的不同序列格式对应不同序列对象产生疑问，如“为了转换序列格式，如何才能把一个PrimarySeq对象转换为RichSeq对象？”答案是：这种问题压根就不用担心，SeqIO会自动转换。之所以不同的格式对应不同的对象，因为有些序列有很少信息，而有些格式又包含很多信息。（译者注：可能是避免大炮打蚊子似的浪费）你并不需要关心SeqIO是如何转换的。列出表2，展示一下常见的序列格式和对应的序列对象，也可满足一些人的好奇心。（表格见原文。）

举例

先举个最简单的例子，如何从序列文件中提取accession number。

# first, bring in the SeqIO module

use Bio::SeqIO;

# Notice that you do not have to use any Bio:SeqI
# objects, because SeqIO does this for you. In fact, it
# even knows which SeqI object to use for the provided
# format.

# Bring in the file and format, or die with a nice
# usage statement if one or both arguments are missing.
my $usage  = "getaccs.pl file format\n";
my $file   = shift or die $usage;
my $format = shift or die $usage;

# Now create a new SeqIO object to bring in the input
# file. The new method takes arguments in the format
# key => value, key => value. The basic keys that it
# can accept values for are '-file' which expects some
# information on how to access your data, and '-format'
# which expects one of the Bioperl-format-labels mentioned
# above. Although it is optional, it is good
# programming practice to provide > and < in front of any
# filenames provided in the -file parameter. This makes the
# resulting filehandle created by SeqIO explicitly read (<)
# or write(>).  It will definitely help others reading your
# code understand the function of the SeqIO object.

my $inseq = Bio::SeqIO->new(
                            -file   => "<$file",
                            -format => $format,
                            );
# Now that we have a seq stream,
# we need to tell it to give us a $seq.
# We do this using the 'next_seq' method of SeqIO.

while (my $seq = $inseq->next_seq) {
    print $seq->accession_number,"\n";
}

运行此脚本需要两个参数，输入序列的文件名和文件格式。上例显示了最基本的处理genbank文件的方式。其他格式，如fasta,swissprot, ace都是类似的，需要给bioperl提供其所支持的格式名称。（译者注：一般情况下bioperl可以自己猜测出来序列的格式，而不需提供format参数，为保险期间最好说明格式类型。）

需要注意的是，SeqIO默认是处理多个序列的。每一次调用next_seq时，返回的是下一个序列，直到返回一个undef，即文件末尾。脚本是一个一个序列依次读取，相比一次性把所有序列读取后再处理更节省内存。undef的重要性在于，当文件达到末尾的时候，undef可终止while循环。

下例是把EMBL文件中的所有序列信息都放先在一个数组里，然后从数组里提取所需信息。

use strict;
use Bio::SeqIO;

my $input_file = shift;

my $seq_in  = Bio::SeqIO->new(
                              -format => 'embl',
                              -file   => $input_file,
                              );

# loads the whole file into memory - be careful
# if this is a big file, then this script will
# use a lot of memory

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

# now do something with these. First sort by length,
# find the average and median lengths and print them out

@seq_array = sort { $a->length <=> $b->length } @seq_array;

my $total = 0;
my $count = 0;
foreach my $seq ( @seq_array ) {
    $total += $seq->length;
    $count++;
}

print "Mean length ",$total/$count," Median ",$seq_array[$count/2]->length,"\n";

现在来看如何转换序列格式。当你从一个序列文件创建一个Bio::SeqIO对象后，每次调用next_seq，Bioperl都会通过运行很多相关的脚本，来提取下一条序列的信息，并记录到一个SeqI对象中。神奇的是，next_seq一次性提取整条序列的所有信息，因为Bioperl能够理解序列的起始和终止位置，当提取完一条序列信息后自动停止，等待下一次next_seq的调用。同时，Bioperl也能识别出序列的注释信息，例如在genbank序列 LOCUS行中的DIVISION信息。要转换格式，就是把存储序列信息的一个Bio::SeqI对象以另一种格式输出到新文件中。这时就要用Bio::SeqIO中的另一个方法，write_seq。write_seq的工作原理与刚才所说的读取序列的过程正好相反，并通过next_seq将输入和输出的过程连接起来。Bioperl可以将一个SeqI对象按照所需格式输出到新的序列文件中。看下面的示例会更清楚一些：

use Bio::SeqIO;
# get command-line arguments, or die with a usage statement
my $usage         = "x2y.pl infile infileformat outfile outfileformat\n";
my $infile        = shift or die $usage;
my $infileformat  = shift or die $usage;
my $outfile       = shift or die $usage;
my $outfileformat = shift or die $usage;

# create one SeqIO object to read in,and another to write out
my $seq_in = Bio::SeqIO->new(
                             -file   => "<$infile",
                             -format => $infileformat,
                             );
my $seq_out = Bio::SeqIO->new(
                              -file   => ">$outfile",
                              -format => $outfileformat,
                              );

# write each entry in the input file to the output file
while (my $inseq = $seq_in->next_seq) {
    $seq_out->write_seq($inseq);
}

可以将$seq_in和$seq_out想象成两个特殊的文件句柄，并且这个文件句柄“知道”序列及其格式。用文件句柄时一般用类似<F>的操作符，而$seq_in和$seq_out则使用next_seq()方法来读取或输出序列对象，如用“$seqio->write_seq($seq_object)”相对于“print F $line”。

注：Bio::SeqIO允许使用一种办法来模仿文件句柄的输入输出（也不知道用这种方法是傻还是聪明），让<F>存储序列，“print F”则是输出序列。但这种方式对很多人来说是一种非主流想法，看起来很别扭，有时候也会引起一些不必要的困惑。

这个可通用与序列格式转换的脚本只比前面例子里提取accession号和计算序列平均长度的脚本多了几行代码。（主要增加的部分是为了从命令行获取参数）这就是Bioperl的优势所在，可以用很少的代码完成很复杂的任务。

接着，我们来看通过简单修改上面例子中的代码来展示一下SeqIO的灵活性。例如，用其他程序输出一段内容，然后用Bioperl来接受其输出的内容，处理后输出到新文件中。这里需要明确亮点，一个和Perl有关，另一个是Bioperl特有的。Perl可以通过在一个句柄的名字前加“*”号将其转换（GLOB）成一个文件句柄，以便于后面使用。另外，Bio::SeqIO可以将前面转换后的文件句柄作为-fh的参数，以代替-file参数。下面展示的例子中，all2y.pl可以从其他程序的输出端获取内容，使用时类似于：

    >cat myseqs.fa | all2y.pl fasta newseqs.gb genbank

其代码如下：

use Bio::SeqIO;
# get command-line arguments, or die with a usage statement
my $usage     = "all2y.pl informat outfile outfileformat\n";
my $informat  = shift or die $usage;
my $outfile   = shift or die $usage;
my $outformat = shift or die $usage;

# create one SeqIO object to read in, and another to write out
# *STDIN is a 'globbed' filehandle with the contents of Standard In
my $seqin = Bio::SeqIO->new(
                            -fh     => \*STDIN,
                            -format => $informat,
                            );
my $seqout = Bio::SeqIO->new(
                             -file   => ">$outfile",
                             -format => $outformat,
                             );

# write each entry in the input file to the output file
while (my $inseq = $seqin->next_seq) {
    $seqout->write_seq($inseq);
}

有时候可能根本就用不着文件。可以直接使用 STDIN和STDOUT，将输出序列做为另一程序的输入，比如：

    cat *.seq | in2out.pl EMBL Genbank | someother program

代码如下：

use Bio::SeqIO;
# get command-line arguments, or die with a usage statement
my $usage     = "in2out.pl informat outformat\n";
my $informat  = shift or die $usage;
my $outformat = shift or die $usage;

# create one SeqIO object to read in, and another to write out
my $seqin = Bio::SeqIO->new(
                            -fh     => \*STDIN,
                            -format => $informat,
                            );
my $outseq = Bio::SeqIO->new(
                             -fh     => \*STDOUT,
                             -format => $outformat,
                             );

# write each entry in the input to the output
while (my $inseq = $seqin->next_seq) {
    $outseq->write_seq($inseq);
}

字符串相关

一个很常见的问题是：“如果我有一个字符串，其中包含了一系列的某格式的序列，我该如何才能把这些序列转换成Bio::Seq对象？”比如有些时候，需要用户在网页表格内输入一段序列，然后对这段序列做一些处理。可以通过把用户输入的一段字符串转化成符合Perl标准的句柄(Perl 5.8.0以后可以直接用open来完成），然后-fh参数就可设置成这个句柄，而不用设置成文件名。

下面的两个例子都不是完整的程序，只是给个基本的示例。假设通过CGI脚本获取到用户输入的序列和序列格式，并分别放到两个字符串变量中。刚才提到Bio::seqIO可通过文件后缀名识别序列格式，这里输入的序列不是文件，所以必须告知Bioperl文件格式。

use IO::String;   # only needed for Perl versions previous to 5.8.0
use Bio::SeqIO;

## get a string into $string somehow, with its format in $format, say from a web form.
my $string   = ">SEQ1\nacgt\n>revseq1\ntgca\n";
my $format   = "fasta";

my $stringfh = IO::String->new($string);                                         # Use this for Perl BEFORE 5.8.0
open($stringfh, "<", \$string) or die "Could not open string for reading: $!";   # Use this for Perl AFTER 5.8.0 (inclusive)

my $seqio = Bio::SeqIO-> new(
                             -fh     => $stringfh,
                             -format => $format,
                             );

while( my $seq = $seqio->next_seq ) {
    # process each seq
    print $seq->id . ' = '.$seq->seq()."\n";
}

同样道理，也可以将一个序列对象以一字符串的形式输出。示例（注意open函数内“<”和“>”的区别，“<”用于输入，“>”用于输出）：

use IO::String;   # only needed for Perl versions BEFORE 5.8.0
use Bio::SeqIO;

my $string;
my $stringfh = IO::String->new(\$string);                                        # Use this for Perl BEFORE 5.8.0
open($stringfh, ">", \$string) or die "Could not open string for writing: $!";   # Use this for Perl AFTER 5.8.0 (inclusive)

my $seqOut = Bio::SeqIO->new(
                             -format => 'swiss',
                             -fh     => $io,
                             );
$seqOut->write_seq($seq_obj);
print $string;

更多示例

Bio::SeqIO中的-file参数可以是多个文件。也可以是一段字符串，用于表明从管道输入的内容。格式是

'-file' => '命令 |'。注意后面的“|”符号。当需要处理很大的压缩文件而又没有足够空间解压缩的时候，这种办法很管用。这里给出一个从gzip压缩的genbank文件中提取序列并转换成新的fasta格式的文件，使用方法如下：

     gzip2fasta.pl gbpri1.seq.gz Genbank gbpri1.fa

下面是gzip2fasta.pl的代码：

use Bio::SeqIO;
# get command-line arguments, or die with a usage statement
my $usage    = "gzip2fasta.pl infile informat outfile\n";
my $infile   = shift or die $usage;
my $informat = shift or die $usage;
my $outfile  = shift or die $usage;

# create one SeqIO object to read in, and another to write out
my $seqin = Bio::SeqIO->new(
                            -file   => "/usr/local/bin/gunzip -c $infile |",
                            -format => $informat,
                            );

my $seqout = Bio::SeqIO->new(
                             -file   => ">$outfile",
                             -format => 'Fasta',
                             );

# write each entry in the input to the output file
while (my $inseq = $seqin->next_seq) {
    $seqout->write_seq($inseq);
}

当然Bioperl也可一“管道输出”。使用格式是：'-file' => "| 命令"。这个时候“|”变到了命令的前面。下面举个例子，用Bioperl将一个genbank序列文件转为fasta格式的序列，但不输出文件，直接传递给WashU Blast，让其建立序列数据库。用法：

    any2wublastable.pl myfile.gb Genbank mywublastable p

any2wublastable.pl的代码：

use Bio::SeqIO;

# get command-line arguments, or die with a usage statement
my $usage     = "any2wublastable.pl infile informat outdbname outdbtype\n";
my $infile    = shift or die $usage;
my $informat  = shift or die $usage;
my $outdbname = shift or die $usage;
my $outdbtype = shift or die $usage;

# create one SeqIO object to read in, and another to write out
my $seqin = Bio::SeqIO->new(
                            -file   => "<$infile",
                            -format => $informat,
                            );
my $seqout = Bio::SeqIO->new(
                             -file => "| /usr/local/bin/xdformat -o $outdbname -${outdbtype} -- -",
                             -format => 'Fasta',
                             );

# write each entry in the input to the output
while (my $inseq = $seqin->next_seq) {
    $seqout->write_seq($inseq);
}

可能一些有经验的读者已经意识到new方法返回的是一个引用变量（reference），引用的可以是任意格式的数据。现在来看一个用引用变量存储序列信息的例子。任务是要从一个genbank序列文件中提取出人类的序列并输出到human.gb文件，剩下的序列输出到other.gb文件。这里需要用到两个哈希数组（hash），并分别以“human”和“other”作为下标（key）。

     splitgb.pl inseq.gb

splitgb.pl的代码:

use Bio::SeqIO;

# get command-line argument, or die with a usage statement
my $usage  = "splitgb.pl infile\n";
my $infile = shift or die $usage;

my $inseq = Bio::SeqIO->new(
                            -file   => "<$infile",
                            -format => 'Genbank',
                            );

my %outfiles = (
                'human' => Bio::SeqIO->new(
                                           -file   => '>human.gb',
                                           -format => 'Genbank',
                                           ),
                'other' => Bio::SeqIO->new(
                                           -file   => '>other.gb',
                                           -format => 'Genbank',
                                           ),
                );

while (my $seqin = $inseq->next_seq) {
    # here we make use of the species attribute, which returns a
    # species object, which has a binomial attribute that
    # holds the binomial species name of the source of the sequence
    if ($seqin->species->binomial =~ m/Homo sapiens/) {
        $outfiles{'human'}->write_seq($seqin);
    } else {
        $outfiles{'other'}->write_seq($seqin);
    }
}

现在来看用更复杂的多维hash来存储Genbank和fasta的序列信息。

use Bio::SeqIO;
# get command-line argument, or die with a usage statement
my $usage  = "splitgb.pl infile\n";
my $infile = shift or die $usage;

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);
    }
}

最后要说的是单行程序。Perl允许在命令行里用单行程序的方式编写脚本，而不需要放在脚本文件中。如下所示，-e后面用单引号括住脚本代码，-M后面跟所调用的模块。（译者注：-M后面没有空格，-e后面可有可无，多个模块则写多个-M。）当脚本内部也需要使用单引号时必须用q()来代替，防止混淆。下面举的例子是要找出在一个压缩的序列文件中有多少GSS序列。注：为了便于阅读，下面的代码有分行，但在命令行中一般没有写完代码不分行。（译者注：有单引号存在的时候，也可以分行，直到单引号结束再回车才会运行命令。）

perl -MBio::SeqIO -e 'my $gss=0;
my $in=Bio::SeqIO->new(q(-file)=>q(/usr/local/bin/gunzip -c gbpri1.seq.gz |), q(-format)=>q(Genbank));
while (my $seq = $in->next_seq) { $gss++ if ($seq->keywords =~ m/GSS/);}
 print "There are $gss GSS sequences in gbpri1.seq.gz\n";'

警告

因为BioPerl用一个单一而通用的数据接口来存取所有支持格式的文件，有时候Bioperl就会或多或少改变原本文件的数据结构。例如，从GenBank数据库中下载到一个genbank格式的序列，然后用bioperl转换成一个新的genbank格式的序列。用”diff origfile newfile“后会惊奇地发现两者有所不同，这时你就应该明白我第一句说的是什么了。需要记住Bioperl并不提供“往返行程”，它只是通过一个通用的数据接口来获取多种多样序列的信息，并用于其他脚本或程序从中提取所需信息。

纠错

如果所给定的文件不存在，脚本会运行出错并给出相应的错误信息。这用编程术语叫做“抛出异常”。例如：

      user@localhost ~/src/bioperl-live> perl t.pl bollocks silly
      ------------- EXCEPTION  -------------
      MSG: Could not open bollocks for reading: No such file or directory
      STACK Bio::Root::IO::_initialize_io Bio/Root/IO.pm:259
      STACK Bio::SeqIO::_initialize Bio/SeqIO.pm:441
      STACK Bio::SeqIO::genbank::_initialize Bio/SeqIO/genbank.pm:122
      STACK Bio::SeqIO::new Bio/SeqIO.pm:359
      STACK Bio::SeqIO::new Bio/SeqIO.pm:372
      STACK toplevel t.pl:9
      --------------------------------------

这种错误信息非常有用。因为可以从中可以找出出错的地方，也叫“堆栈追踪”。如上例中给出了出错的文件和与错误相关的行号。（译者注：上例中t.pl是用户自己的脚本，出错是在t.pl文件的第9行。其他都是bioperl包含的一些相关模块文件，所列行号只是和错误相关，但一般不会是bioperl本身出错。）

Bioperl会自动识别出类似的错误并终止程序运行。但最好是自己就在脚本中设置一些纠错机制。可以通过类似下例来“捕获异常”。

use strict;
use Bio::SeqIO;

my $input_file  = shift;
my $output_file = shift;

# we have to declare $seq_in and $seq_out before
# the eval block as we want to use them afterwards

my $seq_in;
my $seq_out;

eval {
    $seq_in   = Bio::SeqIO->new(
                                -format => 'genbank',
                                -file   => $input_file,
                                );

    $seq_out  = Bio::SeqIO->new(
                                -format => 'fasta',
                                -file   => ">$output_file",
                                );
};
if( $@ ) { # an error occurred
    print "Was not able to open files, sorry!\n";
    print "Full error is\n\n$@\n";
    exit(-1);
}
my $seq;
while( $seq = $seq_in->next_seq() ) {
    $seq_out->write_seq($seq);
}

eval是perl的一种纠错机制，不止是BioPerl才有。如果eval {BLOCK}中的一段程序出错后，回将错误信息返回到$@变量。需要注意的是：必须提前声明变量 $seq_in和$seq_out。如果在eval {}内用my声明变量，此变量只被设定在eval{}中，当在eval{}之后调用此变量是就会出错（用use strict时会提示错误并终止运行）。

提速, Bio::Seq::SeqBuilder

如果需要处理大批量的序列，且只需要从中提取很少部分信息（如只提取genbank文件中的序列），可以使用Bio::Seq::SeqBuilder来设定Bio::SeqIO只提取需要的部分信息。这样能够节省很多时间从而提升脚本的运行速度。

例如，只需要genbank文件中三类注释信息，可提升大约6倍的速度。

#!/usr/bin/perl

use strict;
use Bio::SeqIO;
use Benchmark qw(:all);

my $file = "gbbct10.seq";

timethis(1, sub {
   my $in = Bio::SeqIO->new(-file => $file, -format => "genbank");
   for (1..1000) {
      my $seq = $in->next_seq;
   }
});

timethis(1, sub {
   my $in = Bio::SeqIO->new(-file => $file, -format => "genbank");
   my $builder = $in->sequence_builder();
   $builder->want_none();
   $builder->add_wanted_slot('display_id','desc','seq');
   for (1..1000) {
      my $seq = $in->next_seq;
   }
});

结果：

timethis 1: 10 wallclock secs ( 9.64 usr +  0.02 sys =  9.66 CPU) @  0.10/s (n=1)
            (warning: too few iterations for a reliable count)
timethis 1:  1 wallclock secs ( 1.63 usr +  0.00 sys =  1.63 CPU) @  0.61/s (n=1)
            (warning: too few iterations for a reliable count)

更多相关内容见HOWTO:Feature-Annotation

Bioops

Bioinformatics=(ACGAAG->AK)+(#!/bin/sh)+(P(A|B)=P(B|A)*P(A)/P(B))

利用Bioperl进行序列输入和输出（Bioperl HOWTO翻译13）