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:
error information:
NOTE:
BTW, my bioperl version is 1.6.1, I use blastx 2.2.27+.
My blast report(-outfmt 0) have many queries and each query have many hits(against nr database), so each hit maybe also have several HSPs.
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:
Code:
#!/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 } }
-------------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 ...
I don't know what does it mean. Please give me some clues, all your words are welcome!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 ...
NOTE:
BTW, my bioperl version is 1.6.1, I use blastx 2.2.27+.
My blast report(-outfmt 0) have many queries and each query have many hits(against nr database), so each hit maybe also have several HSPs.
Comment