转自: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
《 “Bioperl:使用和解析 BLAST” 》 有 6 条评论
对生物信息学不太了解,想问下,我看你写了好几篇解析BLAST的,不知BLAST解析后是用来作什么的?谢谢
对大批量序列的blast比较有用嘛~~这是blast后,提取有用的信息啊~~不用你一个个去整理.
哦,明白,谢谢
請問要如何抓自己建的資料庫裡的序列做比對
没看明白你想问什么?
请到yunbio.com详细提问~ [握手]
您好,我想使用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)呢?