[Bioperl-l] Regarding blast in Bioperl

Smithies, Russell Russell.Smithies at agresearch.co.nz
Mon Jan 25 16:14:15 EST 2010


That's a fair mix of incomplete code you've supplied!!
Did you read the documentation for RemoteBlast? The example there will do 99% of what you want.
http://search.cpan.org/~cjfields/BioPerl-1.6.1/Bio/Tools/Run/RemoteBlast.pm

I'm not entirely sure what you're trying to do (as you've left out a bit of your code) but I assume you're trying to retrieve and print the sequence for each hit.

Here's something that works, not sure exactly what/why you want to print but it should get you a bit further.

--Russell


================================
#!perl -w

use Bio::Tools::Run::RemoteBlast;
use Bio::DB::GenBank;

use CGI ':standard';

use strict;

my $q = new CGI;

my @params = (
               -prog         => 'blastn',
               -data         => 'nr',
               -expect       => '1e-30',
               -entrez_query => 'Homo sapiens [ORGN]',
               -readmethod   => 'SearchIO'
);

my $gb = Bio::DB::GenBank->new;

my $factory = Bio::Tools::Run::RemoteBlast->new(@params);

#$v is just to turn on and off the messages
my $v = 1;

my $str = Bio::SeqIO->new( -file => 'test.faa', -format => "fasta" );

while ( my $input = $str->next_seq() ) {

  my $r = $factory->submit_blast($input);

  print STDERR "waiting..." if ( $v > 0 );
  while ( my @rids = $factory->each_rid ) {
    foreach my $rid (@rids) {
      my @seqs = ();
      my $rc   = $factory->retrieve_blast($rid);
      if ( !ref($rc) ) {
        if ( $rc < 0 ) {
          $factory->remove_rid($rid);
        }
        print STDERR "." if ( $v > 0 );
        sleep 5;
      }
      else {
        my $result = $rc->next_result();

        #save the blast output
        my $filename = $result->query_accession . '.out';
        $factory->save_output($filename);
        $factory->remove_rid($rid);
        print "\nQuery Name: ", $result->query_name(), "\n";
        while ( my $hit = $result->next_hit ) {

          # store the hit sequences
          push @seqs, $gb->get_Seq_by_version( $hit->name );

          next unless ( $v > 0 );
          print "\thit name is ", $hit->name, "\n";
          while ( my $hsp = $hit->next_hsp ) {
            print "\t\tscore is ", $hsp->score, "\n";
          }
        }

        ## print the seqs you've retrieved??
        open( OUTFILE, '>', $result->query_accession . '.htm' );
        print OUTFILE $q->start_html('RNAi Result'),
          $q->h1('RNAi Result'),
          $q->h2('Input'),
          $q->pre( toString($input) ),
          $q->h2('Output');

        foreach (@seqs) {

          #there's probably a better way of printing the seq
          print OUTFILE $q->pre( toString($_) );
        }
        print OUTFILE $q->end_html;
        close OUTFILE;
      }
    }
  }
}

sub toString {
  my $s = shift;
  return '>' . $s->display_id . " " . $s->desc . "\n" . $s->seq;
}


=======================================================================
Attention: The information contained in this message and/or attachments
from AgResearch Limited is intended only for the persons or entities
to which it is addressed and may contain confidential and/or privileged
material. Any review, retransmission, dissemination or other use of, or
taking of any action in reliance upon, this information by persons or
entities other than the intended recipients is prohibited by AgResearch
Limited. If you have received this message in error, please notify the
sender immediately.
=======================================================================



More information about the Bioperl-l mailing list