初品BioPerl（第二篇：构造一条fasta序列）

接下来我们就来讨论BioPerl的用法。根据一句经典名言，Perl的用途“90%与文本处理有关，10%与其它事务有关”（改编自小驼书），Perl语言的强项就在于文本处理（当然主要是纯文本，和许多Unix工具如grep, awk, sed等工具类似哦！），而恰好大多数生物信息学数据都是以纯文本的形式存放的（包括蛋白质、核酸等序列文件，序列比对文件，进化树文件等），所以BioPerl当初设计的初衷就是为了分析、处理这些文本数据文件。当然随着模块越来越多，BioPerl的功能也扩展了，现在也有人喜欢用BioPerl作为下载工具来下载序列，或者用来运行blast程序等。虽然不是不可以，但我并不建议这么做，因为通过BioPerl来运行的blast程序肯定没有直接运行的blast程序来得快，来得灵活。正如同想要删除一个文件，我们一般都会执行rm file.txt，而不会另外编一个Perl程序说：perl -e 'unlink "file.txt"' 。
进入正题啦！首先，生物信息学处理最多的问题是什么？当然就是蛋白质和核酸序列喽！生物大分子序列的书写有好多种格式，如Fasta, Genbank, EMBL, SwissProt等。其中Fasta是最简单的序列格式，所以我们先来学习使用BioPerl来构造fasta序列。当然这在实际工作上意义不大，因为大部分序列应该从数据库中下载得到（或者通过程序运算出来），而不是自己构造出来的。我们在下一篇会学习如何从文件中提取fasta序列。现在还是先打基础吧。 ^_^

      fasta格式的序列是这样的：

也就是说，上面这条序列的“名称”是：    gi|147605|gb|J01673.1|ECORHO
“描述”是：    E.coli rho gene coding for transcription termination factor
序列则是下面的一大段
当然，这条序列的名称比较复杂，如果是自己“构造”序列，自然没必要用这么复杂的名称，可以只叫ECORHO或者叫Ecolirho都行。但千万注意。由于BioPerl一般是把从大于号开始一直到第一个空格出现的字段认为是“名称”，所以序列的名字中不要出现空格。如果你把序列的名字叫做E.coli rho，把头字段写成这样：

BioPerl会认为“名称”是E.coli（从大于号开始到第一个空格出现），而“描述”是 rho    E.coli rho gene coding for transcription termination factor，结果会把你搞得稀里糊涂的，要想不引发岐义，可以这么写：

呵呵，只加上一条下划线，就对了。如果我们的序列是从NCBI这样的数据库里下载下来的，无须担心名称和描述的问题，但有时序列是通过其他程序运算出来，需转入下一道工序加工，这时取个好听的名字就很重要。
知道了fasta序列的格式之后，我们可以用一般的Perl程序“编”出这样一条序列：

#!/usr/bin/perl -w 
%a_fasta_seq=(
"display_name" => "gi|147605|gb|J01673.1|ECORHO", #这是序列的名称
"desc" => "E.coli rho gene coding for transcription termination factor", #这是序列的描述
"seq" =>"AACCCTAGCACTGCGCCGAAATATGGCATCCGTGGTATCCCGACTCTGCTGCTGTTCAAAAACGGTGAAG",)  

print <<EOF;    #这是根据fasta格式依次打印名称、描述和序列，使用了Here文档的写法，也可以拆分成多条print语句
> $a_fasta_seq{dispaly_name}\t$a_fasta_seq{desc}
$a_fasta_seq {seq}
EOF

这种方法的优点是步骤十分清晰，只是太麻烦了一点。如果要构造一个GenBank序列，则得写N长的print。还有，如果要构造许多条fasta序列，也不方便。总不能写%fastaseq1, %fastaseq2, %fastaseq3吧？

      下面我们就来看看用BioPerl该怎么写：

#!/usr/bin/perl -w
use Bio::Seq;              #加载Bio::Seq模块。
$seq_obj = Bio::Seq->new;  #调用Bio::Seq模块的new方法，可以建立一个序列对象，命名为$seq_obj。
#给这个新建的$seq_obj对象赋予三个属性的值：dispaly_name（序列的名称），desc（序列的描述）以及seq（序列的内容） 
$seq_obj->display_name("gi|147605|gb|J01673.1|ECORHO");
$seq_obj->desc("E.coli rho gene coding for transcription termination factor");
$seq_obj->seq("AACCCTAGCACTGCGCCGAAATATGGCATCCGTGGTATCCCGACTCTGCTGCTGTTCAAAAACGGTGAAG");

第一句是模块的加载，会把相关的新函数（其实是“方法”）加载进来。
第二句的语法其实在小驼书中已提到过：调用了Bio::Seq模块中的new方法。对于调用面向对象的模块的方法（以前说“函数”，在OO术语中称“方法”）要使用“全名”的方式来调用，语法是：模块名（Bio::Seq）+瘦箭头（->）+方法名（new）。顺便透露一个小秘密：所有面向对象的模块都有一个叫做“new”的方法，它的功能简单且重要：创建一个新的对象。
会一点Perl的人都应该会创建一些所谓的“数据结构”，比如说：
创建一个标量：  $scalar;
创建一个数组：  @array;
创建一个哈希：  %hash;
因为标量、数组和哈希都太“简单”啦，而且是Perl内置的数据结构，在标量、数组或哈希里以什么样的方式来组织数据都是Perl自己已经定好了的（注意是组织数据的方式而不是数据的内容哦！）。而对象则可以看成是自己定义的一种数据结构，我们可以给对象赋予一些属性来描述它（有点类似于在哈希中放入键值对来描述你创建的哈希）。但是创建的对象不像标量、数组或哈希那样有前缀（$, @和%），在Perl中对象是用存放它的地址来表示的，而恰恰凑巧的是，地址（一般叫做“引用”）的前缀就是$。
所以，创建一个Bio::Seq对象就是这样：

$seq_obj = Bio::Seq->new;

你可能还会为两个问题而困惑：（1）问：为什么要使用new函数来返回一个对象，而不是 像创建标量那样直接写$seqobj就行了呢？ 答：因为对象是需要我们自定义的，所以我们需要写一个子程序来定义对象的结构（其实是它的属性和方法，只是空的白纸，还没有填入内容），而碰巧这个子程序的名称就叫new，并且别人已经帮我们写好且内嵌在Bio::Seq模块（以及以后我们要用的其他模块）里了。 （2）问：对象的名称一定要以obj结尾吗？ 答：不需要。你完全可以写 $seq =Bio::Seq->new; ，这样也是创建了一个名字叫做seq的对象，但是到后来你会搞不清$seq究竟是一个存放序列字符串的标量标量呢还是一个存放（指向）整个序列信息的对象呢？所以，既然对象没有”前缀“，我们就自己给它加一个“后缀”上去。 😛
创建了一个seq_obj对象之后，我们可以为它“赋值”。先来回顾一下其他数据结构是怎么赋值的：
给标量赋值：  $scalar=”Hello!\n”;
给列表赋值：  @array=(1, 2, 3, 4, 5);
给哈希赋值：  %hash=(“key1″=>”vlaue1”, “key2″=>”value2”);
对象“赋值”的方式比较特殊。对于这个Bio::Seq（以及许多BioPerl对象）对象来说，是采取调用方法来“赋值”。

$seq_obj->display_name("gi|147605|gb|J01673.1|ECORHO");
$seq_obj->desc("E.coli rho gene coding for transcription termination factor");
$seq_obj->seq("AACCCTAGCACTGCGCCGAAATATGGCATCCGTGGTATCCCGACTCTGCTGCTGTTCAAAAACGGTGAAG");

而且，和标量、数据、哈希一样，对象的创建和“赋值”可以合并为一步：

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

$seq_obj = Bio::Seq->new( -display_name => "gi|147605|gb|J01673.1|ECORHO",
-desc => "E.coli rho gene coding for transcription termination factor",
-seq => "AACCCTAGCACTGCGCCGAAATATGGCATCCGTGGTATCCCGACTCTGCTGCTGTTCAAAAACGGTGAAG"
);

第一眼看去，这个程序有点像哈希赋值，但“胖箭头”左边的“键”有点不同寻常，带了短横线。
我们来看看前面一大堆操作创建出来的$seqobj对象到底是怎样的。现在我们已经赋予了它三个属性：displayname, desc, seq，即所谓fasta序列的“三要素”：名称、描述和序列。赋予了这些“值”之后，我们可以查看它们，当然直接打印$seq_obj是行不通的（正如同我们不能自己打印%hash一样），而是要通过调用“方法”来返回可打印的字符串：

my ($fasta_name,$fasta_desc,$fasta_seq); #创建了三个普通的标量变量来存放三个属性值
$fasta_name = $seq_obj->display_name; #调用对象的display_name方法来得到名称
$fasta_desc = $seq_obj->desc;         #调用对象的desc方法来得到描述
$fasta_seq = $seq_obj->seq;           #调用对象的seq方法来得到序列
#现在可以打印三个属性值了 
print "NAME:\t$fasta_name\nDESCRIBE:\t$fasta_desc\nSEQUENCE:\t$fasta_seq\n";

结果将是这样：

NAME:    gi|147605|gb|J01673.1|ECORHO
DESCRIBE:    E.coli rho gene coding for transcription termination factor
SEQUENCE:    AACCCTAGCACTGCGCCGAAATATGGCATCCGTGGTATCCCGACTCTGCTGCTGTTCAAAAACGGTGAAG

如果仅此而已，那BioPerl是没有任何用途的（和直接创建一个哈希%afastaseq一模一样嘛！）。这里面当然是大有文章啦。比如，当我们给一个哈希（或数组）赋值之后，我们除了能查到刚赋值的内容（若干项或所有项）之外，还能查到好多其他内容呢！比如，可以查看数组元素的个数啊（使用scalar @array），哈希的键列表啊（使用keys %hash）等等。对于这个$seq_obj对象，功能就更加强大了。
我们可以直接数出序列的长度：

$fasta_seq_length = $seq_obj->length;  # 得到序列的长度，70

我们可以查阅这段序列究竟是DNA还是蛋白质。

$fasta_property = $seq_obj->alphabet;   # 返回dna，看来Perl猜对啦

由此可以看出来，电脑还是很聪明的，但提供电脑程序的人绝对比电脑更聪明。 😛

尊重他人劳动成果，转载请注明出处:Bluesky's blog » 初品BioPerl（第二篇：构造一条fasta序列）