[Bioperl-l] Errors in a bioperl script, pls help me to have a check, thanks in advance

[Andreas Leimbach]

Hi Lyn,

there's a mistake in your sort subroutine, you're missing one level (the 
hsp level). Thus the module cannot call the method blessed with the hsp 
object. Why don't you use a hash to store each hit and overwrite if 
there's another hit with a higher identity hsp? Something like:

my %best_hit;
while (my $result = $parser->next_result) {
     my $query = $result->query_name;
     while (my $hit = $result->next_hit) {
     my $hit_name = $hit->name;
         while (my $hsp = $hit->next_hsp) { # <- the one you're missing
             my $ident = $hsp->percent_identity;
             if (!$best_hit{$hit_name} || 
$best_hit{$hit_name}->{'perc_identity'} < $ident){
                 $best_hit{$hit_name} = $query;
             }
         }
     }
}
foreach (keys %best_hit) {
     print "$_\t$best_hit{$_}\n";
}
...

HTH,
Andreas


--
Andreas Leimbach
Universität Münster
Institut für Hygiene
Mendelstr. 7
D-48149 Münster
Germany

Tel.: +49 (0)551 39 33843
E-Mail: andreas.leimbach at uni-wuerzburg.de

On 6.6.13 05:22, LynHsiong wrote:
> 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.
> _______________________________________________
> Bioperl-l mailing list
> Bioperl-l at lists.open-bio.org
> http://lists.open-bio.org/mailman/listinfo/bioperl-l
>