Footsteps on my way!
perl/linux/测序分析
        找回密码

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

2016-12-03 来源：师大生信者 分类：生物信息 / 首页 阅读(1049) 评论(0)

这一篇哪，其实主要是复习性质的。你可以松口气，边欣赏、边回味用世界上最美的程序语言（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的某些方法之后能返回什么。

genbank

当你看了这个表（或这张图）之后可能会非常吃惊：我们调用这些方法所返回的都是些无关紧要的信息。我们是想要这条序列的发表日期、序列号、关键词吗？显然不是。我们更想要物种名以及序列里面包含的基因、mRNA、CDS等元件的位置、名称、序列等信息。这么多内容跑哪儿去了？哈哈，这个要留待以后再说。BioPerl有另外一套机制来获取它们。

处理多个Genbank文件

到目前位置我们使用Bio::SeqIO模块都是只读取一个序列文件。而有时我们经常想读取好几个Genbank文件，当然一个很简单的方法是先用某些Unix命令把要处理的一连串gbk文件都合并为一个，再进行操作就行了。比如可以使用Shell重定向来合并文件：

cat *.gbk >combine.gbk     #把当前目录下面所有的Genbank文件都合并到combine.gbk

当然对于没有相关命令可用的用户来说就太痛苦了。其实BioPerl也是可以读取多个文件的，只是模块要换一下，不是Bio::SeqIO了，改成Bio::SeqIO::MultiFile，

use Bio::SeqIO::MultiFile;
#  读取a.gbk和b.gbk两个文件
$catch_seq = Bio::SeqIO::MultiFile -> new(-files=>['a.gbk','b.gbk'], -format=>'genbank');
while($seq_obj = $catch_seq -> next_seq)  {
     do sth...
}

或者这样：

use Bio::SeqIO::MultiFile;
#  读取当前目录下面所有后缀名为gbk的文件。
$catch_seq = Bio::SeqIO::MultiFile -> new(-files=>[glob "*.gbk"], -format=>'genbank');
while($seq_obj = $catch_seq -> next_seq)  {
     do sth...
}

语法和Bio::SeqIO很像，当然你应该写成-files而不是-file，把要处理的文件的列表传给-files就行了，别忘了加上中括号，因为列表是不能当成哈希的值的哦！引用才可以。

从命令行中传入参数

如果我们想用前面的代码处理另一个Genbank文件（比如seq.gbk），我们必须打开源代码文件，把-file => ‘ecorho.gbk’ 改成 -file => ‘seq.gbk’，这样子太麻烦了。事实上这些要处理的文件的名字应该由用户来决定，而不是由你写在源代码里。有没有办法像“钻石操作符”（<>）那样通过命令行参数读入文件名呢？当然可以，而且还很简单，请看：

use Bio::SeqIO;
$filename = shift @ARGV;       #  从命令行中读取文件名,而不是写在程序里
$catch_seq = Bio::SeqIO -> new(-file => $filename, -format => 'genbank');

呵呵，这样子这几行代码就可以处理世界上所有的Genbank文件了，再也不需要修改源代码了。
如果想读取多个文件，也是可以的，稍微修改一下就行了，换汤不换药的。

use Bio::SeqIO::MultiFile;
@filenames = @ARGV;
$catch_seq = Bio::SeqIO::MultiFile -> new(-files => [@filenames], -format =>'genbank');

尊重他人劳动成果，转载请注明出处:Bluesky's blog » 初品BioPerl（第四篇：从本地文件中获取Genbank序列）

分享到：更多 ()

相关推荐

评论 抢沙发

取消

  • 昵称 (必填)
  • 邮箱 (必填)
  • 网址