[Bioperl-l] remoteblast

Roopa Raghuveer rtbio.2009 at gmail.com
Sun Mar 7 13:11:54 UTC 2010


Hello Mark and everybody,

I have been trying to connect to remote blast to retrieve similar sequences
to a given sequence. But my program is unable to retrieve the sequences from
BLAST, i.e., it is getting executed till the remote blast ids, but it is not
entering the else loop after collecting the rid. Please check this problem
and help me in this regard. I think the problem is in getting the sequence
and going to the 'else' part. i.e.,

else {

              open(OUTFILE,'>',$blastdebugfile);   # I think the problem is
in else part, i.e., it is not taking the next result.#
              print OUTFILE "else entered";
              close(OUTFILE);

          my $result = $rc->next_result();

            #save the output

Please give me your reply.

Thanks and regards,
Roopa.

My code is as follows.

#!/usr/bin/perl

#path for extra camel module
use lib "/srv/www/htdocs/rain/RNAi/";
use rnai_blast;


use Bio::SearchIO;
use Bio::Search::Result::BlastResult;
use Bio::Perl;
use Bio::Tools::Run::RemoteBlast;
use Bio::Seq;
use Bio::SeqIO;
use Bio::DB::GenBank;

$serverpath = "/srv/www/htdocs/rain/RNAi";
$serverurl = "http://141.84.66.66/rain/RNAi";
$outfile = $serverpath."/rnairesult_".time().".html";
$nuc = $serverpath."/nuc".time().".txt";
$debugfile = $serverpath."/debug_".time().".txt";
$blastdebugfile = $serverpath."/blastdebug_".time().".txt";

my $outstring ="";

&parse_form;

print "Content-type: text/html\n\n";
print "<HTML>\n";
print "<head><title>RNAi Result</title>";
print "<META HTTP-EQUIV=\"Refresh\" CONTENT=\"30;
URL=$serverurl/rnairesult_".time().".html\"> \n";
print "</head>\n";
print "<body>\n";
print " Your results will appear <a
href=$serverurl/rnairesult_".time().".html>here</a><br>";
print " Please be patient, runtime can be up to 5 minutes<br>";
print " This page will automatically reload in 30 seconds.";
print "</BODY>\n";
print "</HTML>\n";

defined(my $pid = fork) or die "Can't fork: $!";
exit if $pid;
open STDIN, '/dev/null' or die "Can't read /dev/null: $!";
open STDOUT, '>/dev/null' or die "Can't write to /dev/null: $!";
open STDERR, '>&STDOUT' or die "Can't dup stdout: $!";



open(OUTFILE, '>',$outfile);

print OUTFILE "<HTML>\n
 <head><title>RNAi Result</title>
 <META HTTP-EQUIV=\"Refresh\" CONTENT=\"30;
URL=$serverurl//rnairesult_".time().".html\"> \n
 <meta http-equiv=\"expires\" content=\"0\">
 </head>\n
 <body>\n
  Your results will appear <a
href=$serverurl/rnairesult_".time().".html>here</a><br>
  Please be patient, runtime can be up to 5 minutes <br>
 This page will automatically reload in 30 seconds  <br>
 </BODY>\n
 </HTML>\n";

close(OUTFILE);

@compseqs = blastcode($in{'Inputseq'},$in{'Organism'});

$in{'Inputseq'} =~ s/>.*$//m;
$in{'Inputseq'} =~ s/[^TAGC]//gim;
$in{'Inputseq'} =~ tr/actg/ACTG/;

@out = similar($in{'Inputseq'}, \@compseqs, $in{'Windowsize'},
$in{'Threshold'});


sub blastcode
{

$inpu1= $_[0];

$organ= $_[1];

open(NUC,'>',$nuc);
print NUC $inpu1,"\n";
close(NUC);

 my $prog = 'blastn';
 my $db   = 'refseq_rna';
 my $e_val= '1e-10';
 my $organism= $organ;

$gb = new Bio::DB::GenBank;

 my @params = ( '-prog' => $prog,
         '-data' => $db,
         '-expect' => $e_val,
         '-readmethod' => 'SearchIO',
         '-Organism'   => $organism );

             open(OUTFILE,'>',$blastdebugfile);
               print OUTFILE @params;
              close(OUTFILE);


my $factory = Bio::Tools::Run::RemoteBlast->new(@params, -ENTREZ_QUERY =>
"$organ\[ORGN]");

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

  #change a paramter

 #$Bio::Tools::Run::RemoteBlast::HEADER{'ENTREZ_QUERY'} = 'Trypanosoma
Brucei[ORGN]';

#change a paramter
# $Bio::Tools::Run::RemoteBlast::HEADER{'ENTREZ_QUERY'} = '$input2[ORGN]';

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

 my $str = Bio::SeqIO->new('-file' => $nuc , '-format' => 'fasta' ,
'-organism' => "$organ\[ORGN]");

while (my $input = $str->next_seq())
{
   #Blast a sequence against a database:
    #Alternatively, you could  pass in a file with many
    #sequences rather than loop through sequence one at a time
    #Remove the loop starting 'while (my $input = $str->next_seq())'
    #and swap the two lines below for an example of that.
             open(OUTFILE,'>',$debugfile);
             print OUTFILE $input;
              close(OUTFILE);

#submits the input data to BLAST#

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

                open(OUTFILE,'>',$debugfile);
                print OUTFILE $r;
                close(OUTFILE);


   print STDERR "waiting...." if($v>0);

  while ( my @rids = $factory->each_rid ) {
      open(OUTFILE,'>',$debugfile);
          #    print OUTFILE "while entered";
              close(OUTFILE);
     foreach my $rid ( @rids ) {

               open(OUTFILE,'>',$debugfile);
           #    print OUTFILE "foreach entered";
              close(OUTFILE);
#Retrieving the result ids#

        my $rc = $factory->retrieve_blast($rid);

        if( !ref($rc) )
        {
        if( $rc < 0 )
        {
        $factory->remove_rid($rid);
        }
         open(OUTFILE,'>',$debugfile);
             #  print OUTFILE "if entered";
              close(OUTFILE);
         print STDERR "." if ( $v > 0 );
         sleep 5;
        }

       else {

              open(OUTFILE,'>',$blastdebugfile);   # I think the problem is
in else part, i.e., it is not taking the next result.#
              print OUTFILE "else entered";
              close(OUTFILE);

          my $result = $rc->next_result();

            #save the output
          $blastdebugfile = $serverpath."/blastdebug_".time().".txt";

          open(BLASTDEBUGFILE,'>',$blastdebugfile);
          print BLASTDEBUGFILE $result->next_hit();
          close(BLASTDEBUGFILE);
#saving the output in blastdata.time.out file#

      # $random=rand();

        my $filename = $serverpath."/blastdata_".time()."\.out";
   # open(DEBUGFILE,'>',$debugfile);
         # open(new,'>',$filename);
         # @arra=<new>;
         # print DEBUGFILE @arra;
         # close(DEBUGFILE);
         # close(new);

         $factory->save_output($filename);

       # open(BLASTDEBUGFILE,'>',$debugfile);
       # print BLASTDEBUGFILE  "Hello $rid";
       # close(BLASTDEBUGFILE);

       $factory->remove_rid($rid);

       open(BLASTDEBUGFILE,'>',$blastdebugfile);
     #  print BLASTDEBUGFILE  $organism;
        close(BLASTDEBUGFILE);

    # open(OUTFILE,'>',$outfile);
    # print OUTFILE "Test2 $result->database_name()";
    # close(OUTFILE);

#$hit = $result->next_hit;
#open(new,'>',$debugfile);
#print $hit;
#close(new);
$dummy=0;
   while ( my $hit = $result->next_hit ) {

            next unless ( $v >= 0);

          #     open(OUTFILE,'>',$debugfile);
           #    print OUTFILE "$hit in while hits";
            #  close(OUTFILE);

       my $sequ = $gb->get_Seq_by_version($hit->name);
           my $dna = $sequ->seq(); # get the sequence as a string
        $dummy++;
             open(OUTFILE,'>',$debugfile);
          #     print OUTFILE $dna;
              close(OUTFILE);
          push(@seqs,$dna);
         }
        }
      }
    }
  }

$warum=@seqs;
              open(OUTFILE,'>',$debugfile);
             #  print OUTFILE $warum;
               print OUTFILE @seqs;
              close(OUTFILE);


return(@seqs);   #returning the sequences obtained on BLAST#
}



More information about the Bioperl-l mailing list