从NCBI等数据库获取序列数据

有两种方法： 1, 使用 Bio::DB::GenBank （使用Web Interface获取序列数据，当需要获取大量数据的时候不建议使用，否则 ip 有可能被 ban） 2, 下载整个 NCBI 数据库到本地，使用 Bio::DB::Flat 对它建立 index。由于这样以后每次的操作都可以在本地进行，所以对大规模的操作来说，这是比较好的方法。

简单的获取序列数据的示例：

use Bio::Perl; 

my $seq = get_sequence('genbank',$acc); #$acc 是 该序列的 accesion number 

print “I got a sequence $seq for $acc\n”;

以下这个例子可以利用 accesion number 获取一个蛋白质的fasta格式序列（或核苷酸序列）并将它输出到STDOUT：

#!/usr/bin/perl -w 

use strict; 

use Bio::DB::GenPept; 

use Bio::DB::GenBank; 

use Bio::SeqIO; 

my $db = new Bio::DB::GenPept(); # my $db = new Bio::DB::GenBank(); 

# if you want NT seqs # use STDOUT to write sequences 

my $out = new Bio::SeqIO(-format => 'fasta'); 

my $acc = 'AB077698'; 

my $seq = $db->get_Seq_by_acc($acc); 

if( $seq ) 

{ $out->write_seq($seq); } 

else { print STDERR "cannot find seq for acc $acc\n"; } 

$out->close();

如果你使用的是第二种方法，并将序列数据库下载到了本地的话，获取一个序列的方法如下：

注意：序列获取使用的是 LWP::Simple 模块，它会自动获取 HTTP_PROXY 环境变量，因此如果需要使用代理，那么设定你的 HTTP_PROXY 即可。

use Bio::DB::Flat; 

my $db = new Bio::DB::Flat(-directory => '/tmp/idx', -dbname => 'swissprot', -write_flag => 1, \

　　　　　　　　　　　　　　　　　　-format => 'fasta', -index => 'binarysearch'); 

$db->build_index('/data/protein/swissprot'); 

my $seq = $db->get_Seq_by_acc('BOSS_DROME');

许多时候你并不知道你要的序列的accesion number是多少，而只是想搜索某些序列。比如，你想从 NCBI 中找到所有 Arabidopsis 的拓扑异构酶的序列，可以用一下方法实现：

use Bio::DB::GenBank; 

use Bio::DB::Query::GenBank;  

$query = "Arabidopsis[ORGN] AND topoisomerase[TITL]";

#建立一个 $query_obj 用来进行查询操作

$query_obj = Bio::DB::Query::GenBank->new(-db => 'nucleotide', -query => $query ); 

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

$stream_obj = $gb_obj->get_Stream_by_query($query_obj); 

while ($seq_obj = $stream_obj->next_seq) 

{ # 循环对每个 sequence object 进行一些操作 

　　print $seq_obj->display_id, "\t", $seq_obj->length, "\n"; 

}

注意：next_***方法是十分常用的一个方法，在本例中，查询会得到多个结果，使用 next_seq 可以取得下一个 seq，用在 which 循环中用于对每一条查询结果进行操作。

===========================================

#!/usr/local/bin/perl -w
use strict;

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

my $db = new Bio::DB::GenPept;
my $out = new Bio::SeqIO(-format => 'fasta');
my $accfile = shift;

open(ACCFILE, $accfile) || die("could not open $accfile");
while(<ACCFILE>) {
chomp;
my $acc = $_;
my $seq = $db->get_Seq_by_acc($acc);
if( $seq ) {
$out->write_seq($seq);
} else {
print STDERR "could not find sequence $acc\n";
}
}