[Bioperl-l] Runnig bl2seq throws error + Failure to parse BlastResult

[Jason Stajich]

You should  continue to ask your question to the mailing list too.

I presume when query == subject there is in fact an alignment which is why you don't get an error.

However, I don't understand why you are calling rewind_results -- this is what is causing the # of hits == 1 even when there is no hit. 

=head2 rewind

 Title   : rewind
 Usage   : $searchio->rewind;
 Function: Allow one to reset the Result iterator to the beginning, so that
           next_result() will subsequently return the first result and so on.

           NB: result objects are not cached, so you will get new result objects
           each time you rewind. Also, note that result_count() counts the
           number of times you have called next_result(), so will not be able
           tell you how many results there were in the file if you use rewind().


I don't think the parser handles the new Blast+ bl2seq output properly is my guess, it is still creating hit objects even when there are no results, I'm entirely sure but I think has to do with workaround for that format not providing the Query= at the beginning of the output in some cases.  blast+ may be better but staring at the output format suggests it is different.

 
You can try the older bl2seq with StandAloneBlast I suppose.

You may want to reconsider another alignment approach anyways if you are trying to:
"Find number of pairwise similar sequences from a set of sequences."

Usually one searches a sequence against a database if they are looking for all hits, using BLAST or FASTA on a database for each query, not enumerating all pairwise as you seem to be doing in your code.

You can always stop having messages displaying on the console by not having perl run with -w  as for this script that will turn off the error message.


On Jun 28, 2011, at 3:10 PM, Das, Abhiram wrote:

> I wanted to stop further checking of hits and then hsps to decide if there was a match. 
> 
> The message is annoying when I am running 100s of sequences. Is there any way I can stop the message printing to the console?
> 
> The message do not show up when the 'query' is same as 'subject'.
> 
> Thanks again.
> -Abhi
> 
> 
> ----- Original Message -----
> From: "Jason Stajich" <jason.stajich at gmail.com>
> To: "Abhiram Das" <abhiram.das at gatech.edu>
> Cc: bioperl-l at bioperl.org
> Sent: Tuesday, June 28, 2011 5:39:00 PM
> Subject: Re: [Bioperl-l] Runnig bl2seq throws error + Failure to parse BlastResult
> 
> To answer your question "there might be a bug"
> 
> But there won't be any hsps or hits so why do you care if the counter is wrong? 
> 
> 
> On Jun 28, 2011, at 2:32 PM, Das, Abhiram wrote:
> 
>> Jason, thanks for your reply.
>> 
>> I have two issues:
>> 
>> 1. While running the program it shows up the following message:
>> 
>> "Use of uninitialized value in numeric le (<=) at /Library/Perl/5.8.8/Bio/SearchIO/IteratedSearchResultEventBuilder.pm line 367, <GEN4> line 26."
>> 
>> 2. Even though there is no hit; $results->num_hits get's me 1 hit.
>> 
>> 
>> What I am trying to do:
>> 
>> "Find number of pairwise similar sequences from a set of sequences."
>> 
>> Please find attached the sequence file.
>> 
>> Thanks again
>> -Abhi
>> 
>> ----- Original Message -----
>> From: "Jason Stajich" <jason.stajich at gmail.com>
>> To: "Abhiram Das" <abhiram.das at gatech.edu>
>> Cc: bioperl-l at bioperl.org
>> Sent: Tuesday, June 28, 2011 5:02:39 PM
>> Subject: Re: [Bioperl-l] Runnig bl2seq throws error + Failure to parse BlastResult
>> 
>> I'm sorry but I don't understand your question, are you worried that the num_hits is 1 simply? I think this is a function of a forced result from bl2seq -- can you better explain what you are trying to do or provide your sequence file so the issue you are concerned with can be replicated.
>> 
>> Are you trying to construct an all-vs-all pairwise distances based on bl2seq?
>> 
>> On Jun 28, 2011, at 1:26 PM, Das, Abhiram wrote:
>> 
>>> Hi, 
>>> 
>>> Running the following code throws the message: 
>>> 
>>> "Use of uninitialized value in numeric le (<=) at /Library/Perl/5.8.8/Bio/SearchIO/IteratedSearchResultEventBuilder.pm line 315, <GEN6> line 26." 
>>> 
>>> Even though the output file shows no hit's, the BlastResult->num_hits returns 1. 
>>> 
>>> 
>>> 
>>> CODE: 
>>> 
>>> use strict; 
>>> use Bio::Seq; 
>>> use Bio::SeqIO; 
>>> use Bio::DB::GenBank; 
>>> use Bio::Tools::Run::StandAloneBlastPlus; 
>>> use Bio::Search::Result::BlastResult; 
>>> use Bio::Search::Hit::HitI; 
>>> 
>>> my $fac = Bio::Tools::Run::StandAloneBlastPlus->new(); 
>>> 
>>> my $seq_obj = Bio::SeqIO->new(-file => "all_reads.fasta", -format => "fasta", -alphabet => "dna"); 
>>> 
>>> #loop through each seq in the seq-obj and blast it with the next sequence 
>>> my @seq_list; 
>>> while(my $seq = $seq_obj->next_seq()){ 
>>> push(@seq_list, $seq); 
>>> } 
>>> 
>>> my $eval = 0.000001; 
>>> my $word = 16; 
>>> my $match = 0; 
>>> no strict; 
>>> for(my $i = 0; $i < @seq_list; $i++){ 
>>> for(my $j = $i+1; $j < @seq_list; $j++){ 
>>> $fac->bl2seq(-method=>'blastn', -query => $seq_list[$j], -subject => $seq_list[$i], -outfile=>'test.out'); 
>>> $fac->rewind_results; 
>>> if($result = $fac->next_result){ 
>>> print "# hits found:: ", $result->num_hits,"\n"; 
>>> if(my $hit = $result->next_hit){ 
>>> $match++; 
>>> } 
>>> } 
>>> } 
>>> } 
>>> $fac->cleanup(); 
>>> 
>>> 
>>> Appreciate any help. 
>>> 
>>> 
>>> 
>>> 
>>> 
>>> 
>>> 
>>> Thanks 
>>> -Abhi 
>>> 
>>> 
>>> _______________________________________________
>>> Bioperl-l mailing list
>>> Bioperl-l at lists.open-bio.org
>>> http://lists.open-bio.org/mailman/listinfo/bioperl-l
>> 
>> <all_reads.fasta>
>