Bioperl:解析blast结果
下面是一个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.000)
- 用BioPerl解析blast的代码 (1.000)
- Bioperl:对序列数据进行操作 (0.860)
- Bioperl:Genbank Seq parsing script (0.860)
- Bioperl:把Genbank格式的序列转换为基因结构图 (0.860)
- AWK详细参考(转载整理) (RANDOM - 0.500)
你好,我现在想通过blast结果判断基因的拷贝数多少。
请问在score和e-value上应该怎样去cutoff值作为单拷贝和多拷贝的界限。谢谢~