初品BioPerl(第六篇:从远程数据库下载序列之一)

十一月 30, 2010

相信大家都等不及了!以前的序列都是自己手动从网上下载下来的,万一序列很长而国际网速又很慢的话,你的浏览器就会被折腾得死去活来。所以你一定很想知道是否使用BioPerl下载起来能方便一些。
      但在这之前,我还是要先唠叨几句:
     (1)Perl的主要任务是处理文本文件。远程下载只是作为一个扩展功能,所以能不用的话尽量别用,不要让你的Perl太吃力。
     (2)你仍然需要有一个稳定的国际网络连接,否则下载到一半即终止是很头疼的事情。
     (3)想清楚你要下载什么物种的序列。如果是人、酵母这样的标准模式生物,还是先去问一下系统管理员吧!说不定他早已经帮大家下载好并放到某个公共的目录里了,你再下载一遍的话又浪费空间又浪费带宽(存放人类全基因组的gbk文件估计有几个G),浪费就是极大的犯罪。如果你发现好心的系统管理员已经帮你下载好了序列,你只需将它们链接到自己的目录下面即可(注意是链接而不是复制,否则主机又白白浪费了一倍的空间)。
     (4)到ftp://ftp.ncbi.nih.gov/genomes上(或者你知道的其他ftp服务器上)找一下你需要的物种基因组是不是已经在上面了。请找仔细一点。如果找到了,就不要麻烦Perl了,还是使用lftp命令连接到ftp上,使用get(或mget)命令直接下载文件吧,多方便!
     (5)如果很不幸没有找到,那么可以让Perl来帮你做这件事。

-----   开始进入正题   -----

      使用BioPerl下载序列需要用到特殊的模块,比如,下载Genbank序列,就会用到Bio::DB::GenBank模块(注意大小写,别写错了);下载EMBL序列,则用Bio::DB::EMBL模块;下载SwissProt序列就是Bio::DB::SwissProt。先检查一下有没有安装这些模块(一般都应该安装了)。
      因为我们通常说的下载“序列”其实在BioPerl中是指下载“序列对象”,也就是说下载一个包含序列字符串、序列名称以及N多注释的集合(这是典型的面向对象思维),所以,Bio::Seq或Bio::SeqIO模块也会用到。我们先把use语句写上。假设我们下载上一篇提到的那条Genbank序列(ECORHO),就是:

use Bio::DB::GenBank;
use Bio::SeqIO;

      下面要开始“下载”了,第一句话是这样的:

$db_obj = Bio::DB::GenBank -> new;

      这句话看起来有点莫名其妙。用通俗的话解释就是:我要下载Genbank序列啦,请给我一个下载工具吧!没错,这句话就是创建了一个“下载工具”(用行话来说叫做“BioPerl数据库对象”),告诉Perl你要下载Genbank序列。$dbobj就是我们的下载工具(不妨把它想象成迅雷^^),接下来就可以通过调用这个对象的某些方法来下载我们所需要的序列了。
      那么,我们怎样告诉BioPerl下载什么序列呢?最好的方法是告诉BioPerl你想要下载的序列的编号。打个比方,想在学校注册的学生信息数据库中搜寻我,当然是输入我的学号061140119来查找最方便。因为这是唯一的!只要你没有输错学号,就一定能得到唯一且正确的信息。如果你希望输入我的名字Zhu Tao来查找,虽然也可能找到我(很可能排在第一个),但却会找到很多其他不相干的信息,比如任何名字中含有Zhu或Tao的人都会被找出来。
      对于Genbank序列来说,常使用它的accession编号来查找(检索)序列,比如ECORHO这个基因的accession号是J01673,那么可以通过调用Bio::DB::GenBank对象(其实就是$dbobj)的getSeqbyacc方法来检索序列:

$seq_obj = $db_obj -> get_Seq_by_acc('J01673');   #   检索accession号为J01673的序列

      下载得到的序列放到哪儿了?从上面这行代码可以看出,它事实上返回了Bio::Seq对象($seq_obj),它并没有变成一个文本文件存放在我们磁盘上!(事实上是存放在内存里了,一旦Perl程序退出,它就会消失得无影无踪!)这显然不是我们所希望的。我们希望的是序列能放在一个看得见、摸得着的文件里。怎么实现呢?嘿嘿!上一篇不是刚刚讲过把序列“写入”文件的方法吗,就看你的记性好不好喽!

$filename = $seq_obj -> accession;
$output_seq = Bio::SeqIO -> new(-file => ">$filename.gbk", -format => 'genbank');
$output_seq -> write_seq($seq_obj);   # 文件总得有名字吧?我采用了最时髦的做法:以它的accession编号作为文件名,然后以序列格式(gbk)为扩展名

      大家差不多可以把这些代码组合起来。运行之后(需要花一点时间,就看你的网速有多快),你会在当前工作目录下看到一个名为J01673.gbk的文件,大功告成!赶快打开它看一看吧!
      你是不是觉得这样子好像与直接用浏览器下载序列相比差别不大呢?别着急,这只是下载一条序列;下一篇我们会学习批量下载序列,这时BioPerl的优势才会真正体现出来。

由于目前网速问题,不推荐此做法,建议采用新的方法:dbfetch

posted in Perl by billzt

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

说点什么

您将是第一位评论人!

提醒
wpDiscuz
 

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