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

2012/07/27评论4,098

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

发表评论