[Bioperl-l] Errors in a bioperl script, pls help me to have a check, thanks in advance
LynHsiong
LynHsiong at hotmail.com
Thu Jun 6 03:22:04 UTC 2013
Dear friends:
I have a blast report, and want to extract following information for each
result(query): the Query_name, hit_number, name and description of the
hit(HSP) with the highest identity.
my bioperl script is:
#!/usr/bin/perl -w
use strict;
use warnings;
use Bio::SearchIO;
if (@ARGV != 2){die "Usage: $0 <BLAST-report-file> <output-file>\n"};
my ($infile,$outfile) = @ARGV;
print "Parsing the BLAST result ...";
my $blast = new Bio::SearchIO(
-format => 'blast',
-file => $infile);
open (OUT,">$outfile") or die "Cannot open $outfile: $!";
print OUT "query\tNumber of hits\tGreatest identity %\tAccession (identity
%)\tDescription (identity %)\n";
while (my $result = $blast->next_result){
print OUT $result->query_name . "\t";
print OUT $result->num_hits. "\t";
if (my $hit = &sort_hit($result)){
if (my $hsp=$hit->hsp){
print OUT $hsp->percent_identity.
"\t";
print OUT $hit->accession. "\t";
print OUT $hit->description. "\n";
}
}
}
close OUT;
print " DONE!!!\n";
sub sort_hit{
my $result = shift;
my @hits;
while (my $hit = $result->next_hit) {
push @hits, $hit;
}
my @hit_sort = sort { $b-> hsp ->percent_identity * 100
<=> $a-> hsp ->percent_identity * 100 } @hits;
$result->rewind;
my $flag=0;
my $temp_hit;
while (my $hit=$result->next_hit){
if($flag==0){
$temp_hit=$hit;
$flag++;
}
return $temp_hit
}
}
error information:
-------------EXCEPTION: Bio::Root::Exception -------------
MSG: Can't get HSPs: data not collected.
STACK: Error::throw
STACK: Bio::Root::Root::throw
/home/xiongtl/Bioperl/dag/lib/perl5/Bio/Root/Root.pm:368
STACK: Bio::Search::Hit::GenericHit::hsp
/home/xiongtl/Bioperl/dag/lib/perl5/Bio/Search/Hit/GenericHit.pm:688
STACK: main::sort_hit xtl_new.pl:37
STACK: xtl_new.pl:20
-----------------------------------------------------------
Parsing the BLAST result ...
please give me some clues, all your words are welcome!
--
View this message in context: http://bioperl.996286.n3.nabble.com/Errors-in-a-bioperl-script-pls-help-me-to-have-a-check-thanks-in-advance-tp16986.html
Sent from the Bioperl-L mailing list archive at Nabble.com.
More information about the Bioperl-l
mailing list