Bioperl:使用和解析 BLAST


转自:Perlchina.org

使用和解析 BLAST
在生物信息学中最流行的数据库搜索程序就是 BLAST 了,和前面获取序列数据的原则一样,也有两种方法可以让你使用BLAST:在 NCBI 等站点上运行它并得到返回结果,或者下载到本地运行。

BLAST 程序的用法超出本文范围,这里就不介绍了,请自行参照 BLAST 相关文档。现在我们假设你已经有了本地的 BLAST 程序,并且已经使用 formatdb 程序为本地序列数据库 db.fa 建立好了 index。

现在我们来看 Bioperl 如何做 BLAST,它使用Bio::Tools::Run::StandAloneBlast 模块:

use Bio::Seq;
use Bio::Tools::Run::StandAloneBlast;

@params = (program => ‘blastn’,
database => ‘db.fa’ );
$blast_obj = Bio::Tools::Run::StandAloneBlast->new(@params);
$seq_obj = Bio::Seq->new(-id =>”test_query”,
-seq => “TTTAAATATATTTTGAAGTATAGATTATATGTT”);

$report_obj = $blast_obj->blastall($seq_obj);
$result_obj = $report_obj->next_result;
print $result_obj->num_hits;

简单解释一下上面的例子。@params 储存了运行blast程序需要的参数,然后利用它建立了 $blast_obj 对象,以调用 blastall 方法。 $blast_obj 对象储存了用来做 BLAST 搜索的序列, $result_obj 储存了 BLAST 的每一个返回结果,最后用 num_hits 方法 print 出了 BLAST 得到的结果数。

你可能会对 $report_obj 和$result_obj 这两个对象的使用有些迷惑,实际上它们是来自 SearchIO 模块,可以参照 searchIO 的文档进一步了解它们的详细用法。

前面讲了如何调用 BLAST 程序,现在来看如何解析 BLAST 的结果。

use Bio::SearchIO;
my $cutoff = '0.001'; #BlAST 搜索的域值
my $file = 'BOSS_Ce.BLASTP', #结果文件
my $in = new Bio::SearchIO(-format => 'blast',
-file => $file);
while( my $r = $in->next_result ) { #读取每一个 result
print "Query is: ", $r->query_name, " ",
$r->query_description," ",$r->query_length," aa\n";
print " Matrix was ", $r->get_parameter('matrix'), "\n";
while( my $h = $r->next_hit ) { #每个 result 中有多个 hit
last if $h->significance > $cutoff;
print "Hit is ", $h->name, "\n";
while( my $hsp = $h->next_hsp ) { #每个 hit 有它自己的属性
print " HSP Len is ", $hsp->length('total'), " ",
" E-value is ", $hsp->evalue, " Bit score ",
$hsp->score, " \n",
" Query loc: ",$hsp->query->start, " ",
$hsp->query->end," ",
" Sbject loc: ",$hsp->hit->start, " ",
$hsp->hit->end,"\n";
}
}
}

它打印的结果是像下面这个样子的(如果你熟悉 BLAST 就很容易理解每一项的意思):

Query is: BOSS_DROME Bride of sevenless protein
precursor. 896 aa
Matrix was BLOSUM62
Hit is F35H10.10
HSP Len is 315 E-value is 4.9e-11 Bit score 182
Query loc: 511 813 Sbject loc: 1006 1298
HSP Len is 28 E-value is 1.4e-09 Bit score 39
Query loc: 508 535 Sbject loc: 427 454

可以看到 Bio::SearchIO 有很丰富的方法和对象可以处理 BLAST 的结果文件,同样它还可以解析 FASTA,HMMER,UCSC等等许多格式的文件,使用如下命令可以查看它的文档:

perldoc Bio::SearchIO


6条回应:“Bioperl:使用和解析 BLAST”

  1. 对生物信息学不太了解,想问下,我看你写了好几篇解析BLAST的,不知BLAST解析后是用来作什么的?谢谢

  2. 您好,我想使用bioperl是批量提取一个给定物种名(或蛋白质gi号)的rank或者Lineage(如Lineage (full): root; cellular organisms; Eukaryota; Fungi/Metazoa group; Metazoa; Eumetazoa; Bilateria; Coelomata; Deuterostomia; Chordata; Craniata; Vertebrata; Gnathostomata; Teleostomi; Euteleostomi; Sarcopterygii; Tetrapoda; Amniota; Mammalia; Theria; Eutheria; Euarchontoglires; Primates; Haplorrhini; Simiiformes; Catarrhini; Hominoidea; Hominidae; Homininae; Homo)呢?