BioPerl（二）：使用BioPerl读取fasta文件

Overview

在 BioPerl（一）：安装BioPerl 中我们不只安装了BioPerl，还给出了一个使用BioPerl的手动构造了一个Bio::Seq对象的例子，这个对象中包含了我们手动填入的fasta格式的信息。既然可以手动构造，那么也就可以从fasta文件中读取序列信息，由BioPerl自动填充成Bio::Seq对象。

1. 使用`BioPerl`读取单条fasta序列

我们有一个名字叫做singleEntry.fasta的文件，里面有一条序列内容：

>gi|28898513|ref|NP_798118.1| hypothetical protein VP1739 [Vibrio parahaemolyticus RIMD 2210633]|1
MTSQKNIEEQDEVVVIEERDKRTHIYIAIAAVLGLAFGGLAGSVLTANKWESTYQVLEDKYQALAQDKTALVSQVKTREESLDKEIQEKVDALLAEKDAAHQKALKKLQTQLTEVEKVNLSLESQVKQQNDKLSSSKTENEKLTRQADMQATMFERSREVFQKELKISQDLEALEKEREKLLPKIDKLKKECDVFLEGKSWDVKSDACDKHDEANSRLSQIDQLIEVYKMDLKQIKEITSDMGL

这是一个标准的fasta格式的序列，一共有两行，第一行是序列的各种标识，第二行是序列内容，注意，整个序列是一行，只是因为序列太长显示不下，才会自动换行。

我们写了一个名叫readSingleFastaByBioPerl.pl的脚本，内容如下：

#!/usr/bin/perl -w
use strict;
use warnings;
# Bio::SeqIO类用来读取文件
use Bio::SeqIO;
# Bio::Seq类用来存储fasta序列，一条fasta序列对应一个Bio::Seq对象
use Bio::Seq;

# 读取fasta文件，文件名可以用单引号，也可以用双引号，
# -format指定了文件的内容是fasta格式，注意这里的文件不一定是.fasta后缀名，只要内部序列是fasta格式就可以了，文件名可以是.txt或者.out等后缀。
my $catchseq_seqio_obj = Bio::SeqIO->new(-file=>"singleEntry.fasta", -format=>'fasta');
# Bio::SeqIO读取一个文件，会创建一个Bio::SeqIO对象，这里名字叫$catchseq_seqio_obj，有多少条序列，就会有多少个Bio::Seq对象
#下面这句取Bio::Seq对象的方法看起来像是一个典型的迭代器的遍历方法，只是我们这里只有一条序列，所以取一次就可以了。
my $seq_obj = $catchseq_seqio_obj->next_seq;

# Bio::Seq对象($seq_obj)只能改包含了很多序列相关的属性，使用下面的方式取出来。
#  序列名字
my $display_name = $seq_obj->display_name; 
#  序列的描述
my $desc = $seq_obj->desc; 
#   序列字符串                  
my $seq = $seq_obj->seq;    
#   序列的类型（dna还是蛋白质？）                           
my $seq_type = $seq_obj->alphabet; 
#   序列的长度                   
my $seq_length = $seq_obj->length;                   

#使用双引号时，$display_name会显示为$display_name存储的内容，如果想显示$display_name这几个字符，需要使用\$display_name，去掉$的引用含义。
#使用单引号时，$display_name会显示为这几个字符本身
print "\$display_name:\n$display_name\n";
print "\$desc:\n$desc\n";
print "\$seq:\n$seq\n";
print "\$seq_type:\n$seq_type\n";
print "\$seq_length:\n$seq_length\n";

运行脚本readSingleFastaByBioPerl.pl

perl readSingleFastaByBioPerl.pl

显示如下：

$display_name:
gi|28898513|ref|NP_798118.1|
$desc:
hypothetical protein VP1739 [Vibrio parahaemolyticus RIMD 2210633]|1
$seq:
MTSQKNIEEQDEVVVIEERDKRTHIYIAIAAVLGLAFGGLAGSVLTANKWESTYQVLEDKYQALAQDKTALVSQVKTREESLDKEIQEKVDALLAEKDAAHQKALKKLQTQLTEVEKVNLSLESQVKQQNDKLSSSKTENEKLTRQADMQATMFERSREVFQKELKISQDLEALEKEREKLLPKIDKLKKECDVFLEGKSWDVKSDACDKHDEANSRLSQIDQLIEVYKMDLKQIKEITSDMGL
$seq_type:
protein
$seq_length:
244

可以看出，Bio::Seq对象帮我们解析了序列信息并按照指定的格式存储起来。第一行的序列标志被分成了两部分，以第一个空格为分割，前面是序列名称，存在属性display_name中，后面是属性描述，存在属性desc中。序列被完整地存在了seq属性中，除了这些序列中包含的信息之外，Bio::Seq对象还分析了序列的其他属性，如序列类型，序列长度等。

其实Bio::Seq对象不只有这些属性，参考HOWTO:Beginners，可以找到这个对象所有的属性。

与 BioPerl（一）：安装BioPerl 中构造Bio::Seq对象的例子比较不难发现，每个属性都有一个get和set方法。以display_name属性为例，当我们构造Bio::Seq对象时，我们使用display_name属性的set方法$seq_obj->display_name("gi|147605|gb|J01673.1|ECORHO");为该属性赋值，当我们读取Bio::Seq对象的display_name时，我们使用该属性的get方法$seq_obj->display_name，这是一个通用的方式，可以以JAVA的get/set方式类比这种存取变量的方式。

2. 使用`BioPerl`读取多条fasta序列

在实际过程中，我们读取的文件中通常包含了多条fasta序列。其实前面已经提到了，一个Bio::SeqIO对象对应一个读取的文件，里面有多少fasta序列，Bio::SeqIO对象中就有多少个Bio::Seq对象，因此只要循环取一下就可以了：

#!/usr/bin/perl -w
use strict;
use warnings;
use Bio::SeqIO;
use Bio::Seq;

my $catchseq_seqio_obj = Bio::SeqIO->new(-file=>"multipleEntry.fasta", -format=>'fasta');

#  在这儿处理每个序列的信息

#更像Iterator了，$catchseq_seqio_obj->next_seq返回一个Bio::Seq对象，当读取完所有的Bio::Seq对象之后，$catchseq_seqio_obj->next_seq就会返回undef，循环结束。
while(my $seq_obj = $catchseq_seqio_obj->next_seq)
{
    # 序列名字
    my $display_name = $seq_obj->display_name; 
    # 省略下面的...               
    print "\$display_name:\n$display_name\n";
    # 省略下面的...
}