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

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

上一篇讲的是怎样自己“编造”一条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来分析这些文本文件。 ^_^

尊重他人的劳动成果,转载请注明出处:Bluesky's blog » 初品BioPerl(第三篇:从本地文件中获取fasta序列)

分享到:更多 ()

评论 抢沙发

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