初品BioPerl（第四篇：从本地文件中获取Genbank序列）

这一篇哪，其实主要是复习性质的。你可以松口气，边欣赏、边回味用世界上最美的程序语言（Perl）是如何解析复杂的、神秘的生物大分子序列的。

BioPerl比普通的Perl更简单吗？

如果你只是解析简单的fasta序列的话，BioPerl未必很简单，甚至刚开始还比普通的Perl更加难懂呢~但是如果你准备解析一些很复杂的序列格式（或程序输出结果），比如Genbank，SwissProt，EMBL，blast……这时BioPerl威力就显现了。

什么时候使用BioPerl？

看一下两个例子，你就明白啦~

例子一：小明、小王、小红的各科成绩。放在文件score.txt里

name        math     English     computer
Xiaoming      90          84           91
Xiaowang      83          88           82
Xiaohong      89          94           85

例子二：三条fasta序列。放在文件sequence.fasta里

>xiaoming
CGCTATTCAACAAGGCATTCCCCCAAGTTTGAATTCTTTGAAATAGATTGCTATTAGCTA
AGATAAGAACGAAAAGGAAGGATATGGCTAAAGAAACATCTGCAATGCTTTTATTAAAAA
>xiaowang
ACCTAGCGACTCTCTCCACCGTTTGACGAGGCCATTTACAAAAACATAACGAACGACAAT
>xiaohong
GCAGACTCGGTGGAAGAGGTGATATATGGTGCAGTGTAGCTAGCAGACTCGGTGGAAGAC
AGATAAGAACGAAAAGGAAGGATATGGCTAAAGAAACATCTGCAATGCTTTTATTAAAAA

看出区别来了吗？第一个文件score.txt，每行的信息都是独立的（第1行：标题；第2行：小明成绩；第3行：小王成绩；……），相信不用我说大家也该知道了：这样的文件很适合用Perl来处理（如果你还想懒一些，甚至连Perl程序都不用写，直接使用awk之类的Unix命令行工具都行）。第二个文件sequence.fasta，不同行的信息却是有关联的（小明的在1~3行，小王的在4~5行，……）由于Perl只能一行一行地读入信息，就比较麻烦。

Genbank序列与fasta序列有什么不同？

从前面两篇中可以看出，fasta序列“非常简单”，只有序列的名称、描述以及序列字符串。当然你也可以在名称或描述中写入其它对你来说有用的信息（比如，基因的边界在哪儿？内含子的边界在哪儿？），但是BioPerl只能摘取整个描述，还是不认得你在里面写了什么。
相比而言，Genbank序列能描述的内容就非常丰富（类似的还有SwissProt，EMBL等格式的序列），除了名称、描述和序列字符串以外还有序列号、形状（线状还是环状？）、发布日期、所属物种以及序列内包含的基因、RNA、CDS、内含子（如果有的话）……What’s more，这些信息并不是你想怎么写就怎么写的，而都是有一定的规范格式，故BioPerl可以准确地识别并提取这些信息。所以，我强烈建议你不要自行修改下载得到的Genbank序列，否则BioPerl会迷路。

如何从本地文件中提取Genbank序列？

差不多该问这个问题了。如果你已经熟记了前面两篇的内容，应该会依葫芦画瓢了（先到NCBI网站上下载一条Genbnak序列来联系吧！但要短一点的那种）。假设我们下载了一条Genbank序列保存在一个叫做ecorho.gbk 的文本文件里（扩展名不是必需的，我只是为了便于辩认），就能很快地写出如下代码：

use Bio::SeqIO;
$catch_seq = Bio::SeqIO -> new(-file => 'ecorho.gbk', -format => 'genbank');
$seq_obj = $catch_seq -> next_seq;

我在这里偷懒了一点：第二句代码中创建的SeqIO对象我直接写成$catch_seq，把后缀名seqio_obj省掉了。如果你一直不放心$catch_seq是什么的话，还是把后缀加上去比较好。也可以使用ref函数来检查一下$catch_seq究竟是什么。

ref($catch_seq);

在这个例子中该函数会返回Bio::SeqIO::genbank，所以说明$catch_seq是一个Bio::SeqIO类型的对象！如果$catch_seq只是一个普通的标量的话，不会有返回值。
可以接着往下说了，因为GenBank文件的信息很丰富，所以在绝大多数情况下从数据库中下载得到一个Genbank的文件一般都只有一条序列，所以我们还可以偷懒一些，把文件的读取和序列的读取合并起来：

$seq_obj = Bio::SeqIO -> new(-file => 'ecorho.gbk', -format => 'genbank') ->next_seq;

一行语句中出现两个瘦箭头，不习惯吧？但这确实是Perl高手们很习惯的写法。注意：如果一个文件中有多条序列，那是不能这么写的。请老老实实按照上一篇所给的格式写。
最后的问题是，-format参数该怎么写？读取fasta文件时我们写-format => ‘fasta’，读取Genbank文件时我们写-format => ‘genbank’，那么其它的格式呢？比如SwissProt格式，是不是写-format => ‘swissprot’呢？不是的，其实应该写成-format => ‘swiss’，关于Bio::SeqIO模块究竟支持哪些序列格式，以及格式的写法，请到http://www.bioperl.org/wiki/HOWTO:SeqIO看一看。

如何获取序列的信息？

怎样从$seq_obj对象里把genbank序列的信息提取出来呢？我们还是依葫芦画瓢，照着上一篇提到的几项内容，代进去试试看：

$display_name = $seq_obj->display_name;    #  序列的名称 
$desc = $seq_obj->desc;                    #   序列的描述
$seq = $seq_obj->seq;                      #   序列字符串
$seq_type = $seq_obj->alphabet;            #   序列的类型（dna还是蛋白质？）
$seq_length = $seq_obj->length;            #   序列的长度 

你就会愉快地发现它们都能返回一些有趣的信息，$seq_obj能调用的方法还有一些，在http://www.bioperl.org/wiki/HOWTO:Beginners的Table3中有详细的描述。下面这张图列举了调用$seq_obj的某些方法之后能返回什么。