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

评论3,927

从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/

发表评论

匿名网友