[Bioperl-l] num_hits and next_hit() bug in Bio::Search::Result::BlastResult?

Peng Yu pengyu.ut at gmail.com
Mon Jun 7 06:04:55 UTC 2010


I have the following two sequences which don't have any hits. But the
perl code gives num_hits =1 (see below). Is it a bug in parsing blast
results?

$ make
blastn -query <(echo
CACGACAACCCTGTATGCACACTGGTGTCTTCCAGATCGGAAGGGCGGTTCAGCAGGAATGCCGAGGCCGGTATC)
-subject <(echo
GGAAGACACCAGTGTGCATACAGGGGTGTCGTGAGATCGGAAGAGCGTCGTGGAGGGAAAGAGTGGAGATCTCGG)
> HWI-EAS11X_10097_4_1_1523_15064.txt
$ cat first1.fa
>first1
CACGACAACCCTGTATGCACACTGGTGTCTTCCAGATCGGAAGGGCGGTTCAGCAGGAATGCCGAGGCCGGTATC
$ cat second1.fa
>second1
GGAAGACACCAGTGTGCATACAGGGGTGTCGTGAGATCGGAAGAGCGTCGTGGAGGGAAAGAGTGGAGATCTCGG


$ cat HWI-EAS11X_10097_4_1_1523_15064.txt
BLASTN 2.2.23+


Query=
Length=75

Subject=
Length=75


***** No hits found *****



Lambda     K      H
    1.33    0.621     1.12

Gapped
Lambda     K      H
    1.28    0.460    0.850

Effective search space used: 4761




Matrix: blastn matrix 1 -2
Gap Penalties: Existence: 0, Extension: 0

#####

But the following code and output shows that there is one hit?

$ cat main.pl
#!/usr/bin/env perl

use strict;
use warnings;
use Bio::Tools::Run::StandAloneBlastPlus;
use Bio::Perl;
use Data::Dumper;

my $factory = Bio::Tools::Run::StandAloneBlastPlus->new();

my $seq1 = Bio::Perl::read_sequence('first1.fa');
my $seq2 = Bio::Perl::read_sequence('second1.fa');

my $result=$factory->bl2seq(-method=>'blastn',
  -query=> $seq1,
  -subject=> $seq2
);

$factory->cleanup;

print 'ref($factory): ', ref($factory), "\n";
print 'ref($result): ', ref($result), "\n";

print "num_hits: ",  $result->num_hits, "\n";
while(my $hit = $result->next_hit()) {
  print ref($hit), "\n";
  print Dumper($hit), "\n";
}
$ ./main.pl
ref($factory): Bio::Tools::Run::StandAloneBlastPlus
ref($result): Bio::Search::Result::BlastResult
num_hits: 1
Bio::Search::Hit::BlastHit
$VAR1 = bless( {
                 '_hsps' => undef,
                 '_iterator' => 0,
                 '_description' => '',
                 '_query_length' => '75',
                 '_accession' => 'second1',
                 '_length' => '75',
                 '_name' => 'second1',
                 '_rank' => 1,
                 '_algorithm' => 'BLASTN',
                 '_root_verbose' => 0,
                 '_hsp_factory' => bless( {
                                            'interface' =>
'Bio::Search::HSP::HSPI',
                                            'type' =>
'Bio::Search::HSP::GenericHSP',
                                            '_loaded_types' => {

'Bio::Search::HSP::GenericHSP' => 1
                                                               },
                                            '_root_verbose' => 0
                                          }, 'Bio::Factory::ObjectFactory' )
               }, 'Bio::Search::Hit::BlastHit' );



-- 
Regards,
Peng



More information about the Bioperl-l mailing list