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";
    # 省略下面的...
}

3.参考资料