Bioperl:从本地文件中获取fasta序列

从NCBI上下载一个fasta格式的文件,20条WRKY家族基因的蛋白序列,wrky.fasta

文件准备好了,我们想知道它的名称、描述和序列内容!有了这些信息,我们就可以做别的事情,比如判断它是DNA还是蛋白质啦,看看序列有多长,各种碱基或氨基酸比例:

其实这个程序仅适用于文件中只有一条fasta序列的情况

#!/usr/bin/perl -w

open FH, "<wrky.fasta";
while (<FH>)
{ 
    chomp;

    if(/^>/)  #  如果这一行的开头是>,就说明是描述的一行,可以提取序列的名称和描述
    {    # 从>开始到第一个空白出现为“名称”,之后的内容为“描述”
        ($display_name,$desc)=/>(.*?)\s(.*)$/;
    }
    #  如果这一行的开头不是>,则就是序列行。由于序列可以分为好几行,所以要把每一行的序列都连接起来。
    else
    {
        $seq.=$_;
    }
} 
$seq_length =length($seq);      #  计算序列的长度。 
if ($seq=~/[^atgcATGC]/)        #  判断序列是DNA还是蛋白质 
{ 
    $seq_type="dna";
} 
else      {$seq_type="protein";}
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

BioPerl做法:

在BioPerl中,把一整条序列的“信息”存放到某个序列对象里($seq_obj),通过该对象的一些属性和方法来获取所需要信息。

直接把整条序列“输入”到$seq_obj里,它要借助另外一个模块生成的对象:Bio::SeqIO对象

Bio::SeqIO模块可以打开、读取特定的序列文件,比如fasta文件,得到的结果是一个特殊的对象:SeqIO对象。

use Bio::SeqIO;     # 注意大小写
use  Bio::Seq;        #  其实这句话可以不写,因为SeqIO模块“包括”了Seq模块

然后就可以创建一个SeqIO对象来读取文件(注意,SeqIO对象是用来读写文件使用的,而Seq对象是用来存放序列使用的),创建模块对象(Bio::SeqIO)+瘦箭头(->)+方法名(new),所有面向对象的模块(Bio::Seq,Bio::SeqIO)都有一个叫做new的方法,它的功能简单而且重要:创建一个新的对象。我们直接把SeqIO对象的创建与“赋值”合而为一。因为SeqIO对象是用来读取文件的(或写入文件,我们以后再说),所以你需要告诉Perl两件事:(1)你要读取的序列文件叫什么名字?然后把它“赋值”给-file!(2)你要读取的序列文件是什么格式?然后把它“赋值”给-format!
对于本例来说,我们要读取的序列文件叫做wrky.fasta,文件格式是fasta,所以我们要创建SeqIO对象可以这么写:

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

我们希望把完整的一条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的属性和方法还很多呢。比如,可以从里面截取一段序列,还可以求它的互补序列。

处理多条序列:

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

 $seq = $seq_obj->seq;                                  

 $seq_type = $seq_obj->alphabet;                      

 $seq_length = $seq_obj->length;                    

 }

本文来自:http://blog.163.com/guojinyan2008@126/blog/static/3074706720111018112829613/

  • 文章来源: 未知。文章来源待更新,请等待。
  • 版权说明: 除非特殊说明,本站文章版权归于文章来源网站或投稿作者。未标记来源文章,请原作者联系管理员更新版权信息

发表评论

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen: