[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