下面是一个bioperl小脚本, 用于解析blast的结果.
#!/usr/bin/perl -w use strict; use Bio::SearchIO; my $filename = 'MET30-vs-swissprot.3.WUBLASTP'; my $input = Bio::SearchIO->new(-format => 'blast', -file => $filename); #GOAL: # Print out the query, hit names and percent identity, percent query aligned while( my $r = $input->next_result ) { my $qname = $r->query_name; my $qlen = $r->query_length; while( my $hit = $r->next_hit ) { my $hname = $hit->name; while( my $hsp = $hit->next_hsp ) { my $percent_id = $hsp->frac_identical * 100; my $query_aln_len = $hsp->length('query'); # another way to do this # my $query_aln_len = $hsp->query->length; print join("\t", $qname, $hname, sprintf("%.1f",$percent_id), sprintf("%.2f",100 * ( $query_aln_len / $qlen))), "\n"; } } }
《“Bioperl:解析blast结果”》 有 1 条评论
你好,我现在想通过blast结果判断基因的拷贝数多少。
请问在score和e-value上应该怎样去cutoff值作为单拷贝和多拷贝的界限。谢谢~