初品BioPerl(第三篇:从本地文件中获取fasta序列)

十月 16, 2010

上一篇讲的是怎样自己“编造”一条fasta序列,说实话确实没什么实际上的用处,充其量也就是了解fasta序列的构造以及Bio::Seq对象能使用的属性、方法等(注意啦:如果你对这些内容还不是很清楚,请先不要往下看,先回去把上一篇看明白哦!)。在实际工作中,我们并不会真的像上一篇讲得那样,自己一项一项地给$seq_obj“赋值”,正如同我们很少对一个哈希逐项赋值一样(我们一般都是从文件或其他信息来源中读入、提取数据并存放到哈希中),而是通过其他来源(其实主要就是存放着fasta序列的文本文件)来将序列“读”到$seq_obj中。让Perl自己来输入显然比我们自己动手敲键盘快多喽!这样我们就有更多的时间来上网冲浪。 🙂
      嗯,那我们还是拿上一篇提到的fasta序列为例子,现在这条fasta序列是放在一个名字叫做ecorho.fasta的文本文件里的(比上一篇要生动多喽,上一篇这条fasta序列可是我们自己手工输入到程序源代码里的),大家可以自己去NCBI网站下载一条fasta序列来练习,也可以点击这里下载这条ecorho.fasta序列来练习。其实,一个文本文件可以放成百上千条fasta序列,但现在我们先练习一条。
      文件准备好了,我们要干什么呢?还记得fasta序列的三要素吗?我们当然是想知道它的名称、描述和序列内容喽!有了这些信息,我们就可以做别的事情,比如判断它是DNA还是蛋白质啦,看看序列有多长啦,各种碱基或氨基酸比例啦……在学习BioPerl之前,我们一般会这么做:

#!/usr/bin/perl -w
open FH, "<ecorho.fasta";         # 打开ecorho.fasta文件
while (<FH>)
{
    chomp;
    #  如果这一行的开头是>,就说明是描述的一行,可以提取序列的名称和描述
    if(/^>/)
    {    # 从>开始到第一个空白出现为“名称”,之后的内容为“描述”
        ($display_name,$desc)=/>(.*?)\s(.*)$/;
    }
    #  如果这一行的开头不是>,则就是序列行。由于序列可以分为好几行,所以要把每一行的序列都连接起来。
    else
    {
        $seq.=$_;
    }
}
#  现在“三要素”已经提取出来,可以进一步分析了。先来计算序列的长度。
$seq_length =length($seq);
#  我们可以判断序列是DNA还是蛋白质
if ($seq=~/[^atgcATGC]/)
{
    $seq_type="protein";
}
else      {$seq_type="dna";}
#  打印一些信息
print
<<EOF;
sequence name: $display_name
sequence description: $desc
sequence type: $seq_type
sequence length: $seq_length
EOF

      输出的结果将会是这样:

sequence name: gi|147605|gb|J01673.1|ECORHO
sequence description: E.coli rho gene coding for transcription termination factor
sequence type: dna
sequence length: 1880

      看起来完全正确,不是吗?且慢,先别急着下结论呢!其实这个程序仅适用于文件中只有一条fasta序列的情况,如果文件中放了好几条fasta序列(这种情况可是相当普遍呵),该程序就乱套了,因为你很难分析这一行是属于哪条序列的!当然,如果你会使用一些高级技巧(就是小驼书中没有提到的内容,比如,你可以把输入分隔符 $/ 定为“>”而不是默认的"\n",这样Perl就会把两个>之间的内容看成“一行”,不同的fasta序列就分开来了),也可以解决问题,但是我们发现这样一来问题并没有变得简单。如果要分析的是多个庞大的genbank序列,还是会把人搞糊涂。

      好啦,我们来看看用BioPerl该怎么做。

     在传统的思路中,我们一般是通过钻石操作符(<>)从文本文件中一行一行地读取数据(可以是真正意义上的“一行”,也可以是“一块”,只要你设定了输入分隔符 $/就行),把这一行的数据以字符串的形式存放在某个标量变量中(一般就是$_),然后再用正则表达式提取所需的信息。在BioPerl中情况稍有不同,我们是把一整条序列的“信息”存放到某个序列对象里(我喜欢写成$seq_obj),然后通过该对象的一些属性方法来获取我们所需要的信息(而不是正则表达式,因为$seq_obj可不是字符串哦!)。
     上一篇中,我们是自己敲打键盘把序列的信息输入到$seq_obj对象里面的,这样太麻烦了。既然现在fasta序列已经规规矩矩地排在文件里面了,有没有办法能直接把整条序列“输入”到$seq_obj里呢?有的。这个方法很特殊,它要借助另外一个模块生成的对象:Bio::SeqIO对象。
      Bio::SeqIO模块可以打开、读取特定的序列文件,比如fasta文件,得到的结果是一个特殊的对象:SeqIO对象。它的属性和方法与之前提到的Seq对象不大一样。(说实话,关于对象的事情比较纠结:整个Perl世界中只有一种数组,一种哈希,看过小驼书的人都知道怎样操作数组、哈希;然而整个Perl世界中却有无数种对象,它们的操作方法没有统一的规则,只有自己去看模块说明文档)
     首先我们要用use语句来加载需要用的模块:

use Bio::SeqIO;
use  Bio::Seq;        #  透露个小秘密:其实这句话可以不写,因为SeqIO模块“包括”了Seq模块

     注意大小写,不要写成Seqio了,BioPerl的难点就在这里:你必须熟记很多模块和方法的名称。虽然记不住的时候可以翻阅文档,但还是记几个比较常用的为好。
     然后就可以创建一个SeqIO对象来读取文件(注意,SeqIO对象是用来读写文件使用的,而Seq对象是用来存放序列使用的,别混淆啦!),还记得上一篇提过的怎样创建一个新的对象吗?其实熟练了以后是很简单的哦:模块名(Bio::SeqIO)+瘦箭头(->)+方法名(new),上一篇也说过,所有面向对象的模块(Bio::Seq,Bio::SeqIO)都有一个叫做new的方法,它的功能简单而且重要:创建一个新的对象。在这里,我们直接把SeqIO对象的创建与“赋值”合而为一。那么怎样为SeqIO对象“赋值”呢?因为SeqIO对象是用来读取文件的(或写入文件,我们以后再说),所以你当然需要告诉Perl两件事:(1)你要读取的序列文件叫什么名字?然后把它“赋值”给-file!(2)你要读取的序列文件是什么格式?然后把它“赋值”给-format!
      对于本例来说,我们要读取的序列文件叫做ecorho.fasta,文件格式是fasta,所以我们要创建SeqIO对象可以这么写:

$catchseq_seqio_obj = Bio::SeqIO->new(-file=>'ecorho.fasta', -format=>'fasta');   #  注意,现在“打开“文件不是写成 open 了!

      我在这里创建的SeqIO对象的名字取得很复杂:$catchseq_seqio_obj,其实是有原因滴!后缀名seqio_obj表示这是一个SeqIO对象,以便于在后续调用中让自己易于识别(Perl本身并不需要为创建的对象取任何后缀名,但是我自己需要呵!),前面的名称catchseq意思是我创建的这个SeqIO对象是用来读取文件的。(因为SeqIO对象还有写入文件的作用,但语法有所不同。以后我们会遇到这样的问题:从一个大的序列文件中读取几条序列,再以另一种格式写入一个新的文件中,如果都叫做$seqio_obj,Perl不一定会混淆,但你自己肯定会混淆)     
      好了,文件ecorho.fasta已被“打开”,怎样读取里面的内容呢?相信你现在应该已经熟悉BioPerl的思维(其实是面向对象的思维)了吧?我们不是希望像以前那样一行一行地读取文本字符串,而是希望把完整的一条fasta序列读进Seq对象里(就是上一篇常写的$seq_obj啦!)。方法是:调用刚刚建立的SeqIO对象的next_seq方法

$seq_obj = $catchseq_seqio_obj->next_seq;     #   next_seq方法将返回一个Seq类型的对象,这里写作$seq_obj

      现在,这条fasta序列的各种信息都已经存在$seq_obj里了,于是我们就可以像上一篇那样通过调用Seq对象的方法来获得需要的信息。

$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的属性和方法还很多呢。比如,你可以从里面截取一段序列,还可以求它的互补序列。所有的这些细节都不需要你自己计算,放心地写->就行了!
#   至于这些方法叫什么名字,你可以查看BioPerl网站的教程:初学者起步

      你或许会认为:这样子操作好像与普通的方法相比简单不了多少。别急!我们现在在文件中只放了一条序列,如果有两条以上的序列,我们想分别取出他们的名称、长度等信息,这时你用常规方法就很麻烦了。但是使用BioPerl却相对要简单的多,其秘诀就在SeqIO对象的next_seq方法里。

      文件中有多条fasta序列,怎么办?

      我们假设你在sequence.fasta里放了两条fasta序列。首先,还是老方法,“打开”这个文件(千万不要写open)。

$catchseq_seqio_obj = Bio::SeqIO->new(-file=>'sequence.fasta', -format=>'fasta');

      然后,调用SeqIO对象的next_seq方法,你将得到第一条序列,可以把它放在$seq1_obj里:

$seq1_obj = $catchseq_seqio_obj->next_seq;   #  接下来可以提取它的各种信息了,序列名称、序列长度等

      那么第二条序列怎么办呢?你可以再调用一次SeqIO对象的next_seq方法,这时它返回的将是第二条序列,可以把它放在$seq2_obj里:

$seq2_obj = $catchseq_seqio_obj->next_seq;   #  接下来可以提取它的各种信息了,序列名称、序列长度等

      再往下,其实下面已经没有序列了。如果你的好奇心无可救药,可以再调用一次SeqIO对象的next_seq方法试试看:

$seq_obj = $catchseq_seqio_obj->next_seq;   #  因为下面没有序列了,所以返回来的是undef

      所以,与其要不厌其烦的写$seq1_obj,$seq2_obj……不如来个while循环来得快捷:

while($seq_obj = $catchseq_seqio_obj->next_seq)
{
    #  在这儿处理每个序列的信息
    $display_name = $seq_obj->display_name;
    #  以下省略
}

      每次循环处理一条序列,多么简洁清晰!(当然你要把每次处理的结果保存起来,比如打印到表格里,不然下次循环不就覆盖了上次么?)不管两三条序列,还是成千上万条序列,都不怕了!
      哦哦!这一讲就到这里为止啦。顺便再提一句,在进行这一讲的练习时,需要你自己先从诸如NCBI这样的网站上手动下载一些fasta序列的文本文件(相信这应该难不倒大家吧?),再使用BioPerl来分析这些文本文件。事实上,BioPerl还可以作为一个下载工具从远程数据库下载序列文件。欲知后事如何,请听以后分解。 ^_^

posted in Perl by billzt

Follow comments via the RSS Feed | Leave a comment | Trackback URL

15
说点什么

11 Comment threads
4 Thread replies
1 Followers
 
Most reacted comment
Hottest comment thread
10 Comment authors
Recent comment authors

This site uses Akismet to reduce spam. Learn how your comment data is processed.

  Subscribe  
最新 最旧 得票最多
提醒
游客
xiongjianLuo

Very useful knowledge...hope see more in future~~`~

游客
sglf

你好,我想请教个小问题了。。
就是怎么写个程序从b文件(序列)中获取一些序列,这些序列的gi号在a文件中。。?

游客
billzt

@sglf, 你可以利用b文件建立一个gi序列号-序列的对应关系(比如构建一个哈希),然后依次根据a文件的gi号在哈希中寻找对应的编号(键)就可以了。如还有问题可以发我邮件或msn

游客
L

写的真好

游客
billzt

@L, 多谢鼓励,敬请关注后续内容

游客

太受用了,非常感谢啊。

游客
billzt

@wal, 过奖,敬请关注后续内容 😛

游客
fanglu

有报错信息

游客
fanglu

第一个例子,说找不到EOF

游客
fanglu

我粘贴不上那条错误信息

游客
hzla

请问一下,有fasta文件,但是每行的序列长度不同(不是最后一行),怎样将它们弄成相同的?

游客
cherish

写的太好了!老师

游客
zy

华东师大生信的路过,博文不错!

游客
billzt

@zy, 谢谢!

游客
soybean

相见恨晚~

 

Copyright © 2010-2018 | Powered by Wordpress and MySQL. Theme by Shlomi Noach, openark.org