Hi monks I'm a first time poster so please bear with me :)

Problem: My bioperl searchIO object is empty, even though the blast output file is not.

I have searched high and low for the past couple of hours for an explanation to the output that I get from my perl script. I hope someone here can shed some light on it.

Details: I am reading a set of nucleotide sequences (solexa reads) from a file in FASTA format. The nucleotide sequences gets translated into protein sequences using all 6 reading frames. If any protein sequences contain stop codons these sequences are discarded. The remaining sequences are then passed on to a Bio::Tools::Run::StandAloneBlast blast object and the method blastall() is called. The method yields a Bio::SearchIO object. From the SearchIO object I should be able to retrieve all the results from my blast using the method next_result() and from that object I should be able to retrieve the blast hits using the method next_hit(). But instead I get this error message telling me that the variable storing the next_result() output is undefined. 
Can't call method "next_hit" on an undefined value at read_blast_stats
+.pl line
    150, <GEN2> line 253 (#1)
    (F) You used the syntax of a method call, but the slot filled by t
+he
    object reference or package name contains an undefined value.  Som
+ething
    like this will reproduce the error:
    
        $BADREF = undef;
        process $BADREF 1,2,3;
        $BADREF->process(1,2,3);
    
Uncaught exception from user code:
    Can't call method "next_hit" on an undefined value at read_blast_s
+tats.pl line 150, <GEN2> line 253.
 at read_blast_stats.pl line 150

[download]
I store the results from the blast search in a file blast.out. After running my script and therby also blasting. This file is not empty, but something is wrong in this file too because only the results from one of the blasted sequences are stored, not all of them as one would have expected.

System:
  • Linux dist. : Red Hat 3.4.6-8
  • Perl v5.8.5 built for x86_64-linux-thread-multi
  • BioPerl v. 1.5.2
  • blastall v. 2.2.18


    • Code: 
    # Creates BioPerl input sequence stream.
my $in_seq = Bio::SeqIO->new(-file   => "<$input_file", -format => $se
+q_format );

# Extracting sequences from stream.
print("- reading and translating sequences...\n");
my $seq_counter = 0;
my $unique_count = 0;
my @translated_seqs;

while (my $seq = $in_seq->next_seq()) 
{
    
    $seq_counter++;
    print("-\tsequences read: $seq_counter\r");
    
    my $temp_seq;
    
    # Translates sequence and stores the translated sequences if they 
+do not
    # contain any stop codons.
    
    # Reading frames 1, 2 and 3.
    $temp_seq = $seq->translate();
    $temp_seq->display_id("test_id-$unique_count");
    $unique_count++;
    push(@translated_seqs, $temp_seq) unless ($temp_seq->seq() =~ /\*+
+/);
    $temp_seq = $seq->translate(-frame => 1);
    $temp_seq->display_id("test_id-$unique_count");
    $unique_count++;
    push(@translated_seqs, $temp_seq) unless ($temp_seq->seq() =~ /\*+
+/);
    $temp_seq = $seq->translate(-frame => 2);
    $temp_seq->display_id("test_id-$unique_count");
    $unique_count++;
    push(@translated_seqs, $temp_seq) unless ($temp_seq->seq() =~ /\*+
+/);

    # Creates the reverse complement sequence.
    my $rc_seq = $seq->revcom();
    
    # Reading frames -1, -2, and -3
    $temp_seq = $rc_seq->translate();
    $temp_seq->display_id("test_id-$unique_count");
    $unique_count++;
    push(@translated_seqs, $temp_seq) unless ($temp_seq->seq() =~ /\*+
+/);
    $temp_seq = $rc_seq->translate(-frame => 1);
    $temp_seq->display_id("test_id-$unique_count");
    $unique_count++;
    push(@translated_seqs, $temp_seq) unless ($temp_seq->seq() =~ /\*+
+/);
    $temp_seq = $rc_seq->translate(-frame => 2);
    $temp_seq->display_id("test_id-$unique_count");
    $unique_count++;
    push(@translated_seqs, $temp_seq) unless ($temp_seq->seq() =~ /\*+
+/);
    
}

# Blasting...
my @params = (-program => 'blastp', -database => '/share/raid1/databas
+e/ftp.ncbi.nih.gov/blast/db/20090602/nr', -expectation => '200000', -
+Word => '2', -Matrix => 'PAM30', -Gapcost => '9.1', -FilterString => 
+'F', _READMETHOD => 'blasttable' , -m => '8', -outfile => 'blast.out'
+);
my $blaster = Bio::Tools::Run::StandAloneBlast->new(@params);
my $blast_report = $blaster->blastall(@translated_seqs);

# Parsing results
my $result = $blast_report->next_result();
my $hit = $result->next_hit();

    [download]


    My thoughts:
  • Because my sequences are created the way that they are, the sequence objects do not contain any other information than the protein sequence. So I thought, maybe blastall is behaving funny because of that. But after giving all sequences unique IDs I still get the same result
  • Blastall will give me the expected output file (with multiple results) if I, after translating the sequences write them all to the to the file blast.in and then calls blastall from the script using that input. But it must be possible to do this without having to write to the disc using BioPerl.
  • I tried taking the appropriate output created without bioperl like described above and process it using Bio::SearchIO. This works fine. So I guess the problem must lie somewhere around StandAloneBlast. I did also try to exchange my array of sequence objects with just one sequence object, but it made no difference.


    • Workaround: The following code is my far-from-optimal workaround so far. I really don't want to make that blast.in file. If you don't have a solution but a better workaround, then I would love to hear it. I have a ton of data and almost unlimited memory so IO operations is bad mmkay ;) 
    # This code is inserted instead of 
# the Blasting and parsing part of the original code.

# Creating blast.in
open(OUTPUT_FILE, "> blast.in") or die "Could not create the output fi
+le \"blast.in\"\n";

# Writing blast.in sequence file in FASTA format
foreach(@translated_seqs){
    print OUTPUT_FILE (">".$_->id()."\n".$_->seq()."\n");}

# Performing blast
`blastall -i "blast_beta.in" -d /share/raid1/database/ftp.ncbi.nih.gov
+/blast/db/20090602/nr -p blastp -e 200000 -W 2 -M PAM30 -G 9.1 -F F -
+m 8 -o "blast.out"`;

# Creating SearchIO
my $blast_report = Bio::SearchIO->new(-format => "blasttable", -file  
+=> "blast.out");

# Testing SearchIO for data.
my $result = $blast_report->next_result();
my $hit = $result->next_hit();
print("Number of hits: ".$hit->length()."\n");

    [download]


    Even if you couldn't help me, I thank you for taking the time to read this. Best wishes :)

    Edit: Changed the filehandle from $temp_file to OUTPUT_FILE.

    In reply to BioPerl StandAloneBlast is returning unexpected undefined SearchIO object by WonkoTheSane

    Title:
    Use:  <p> text here (a paragraph) </p>
    and:  <code> code here </code>
    to format your post, it's "PerlMonks-approved HTML":


  • Posts are HTML formatted. Put <p> </p> tags around your paragraphs. Put <code> </code> tags around your code and data!
  • Titles consisting of a single word are discouraged, and in most cases are disallowed outright.
  • Read Where should I post X? if you're not absolutely sure you're posting in the right place.
  • Please read these before you post! —
  • Posts may use any of the Perl Monks Approved HTML tags:
    a, abbr, b, big, blockquote, br, caption, center, col, colgroup, dd, del, details, div, dl, dt, em, font, h1, h2, h3, h4, h5, h6, hr, i, ins, li, ol, p, pre, readmore, small, span, spoiler, strike, strong, sub, summary, sup, table, tbody, td, tfoot, th, thead, tr, tt, u, ul, wbr
  • You may need to use entities for some characters, as follows. (Exception: Within code tags, you can put the characters literally.)
            For:     Use:
    & &amp;
    < &lt;
    > &gt;
    [ &#91;
    ] &#93;
  • Link using PerlMonks shortcuts! What shortcuts can I use for linking?
  • See Writeup Formatting Tips and other pages linked from there for more info.