Footsteps on my way!
perl/linux/测序分析

初品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的某些方法之后能返回什么。

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序列)

分享到:更多 ()

评论 抢沙发

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