scripts/utilities/bp_remote_blast.pl

   1 #!/usr/bin/perl
   2 #
   3 # BioPerl module for bp_remote_blast.pl
   4 #
   5 # Revived by Evan Weaver for bioperl-1.5.1
   6 # 3/14/2006
   7 #
   8 # Copyright Jason Stajich, Evan Weaver
   9 #
  10 # You may distribute this module under the same terms as perl itself
  11
  12 # POD documentation - main docs after the code
  13
  14 use strict;
  15 use warnings;
  16 use vars qw($USAGE);
  17
  18 use Bio::Tools::Run::RemoteBlast;
  19 use Bio::SeqIO;
  20 use Getopt::Long;
  21
  22 $USAGE = "remote_blast.pl [-h] [-p prog] [-d db] [-e expect] [-mod
  23 Blast] [-f seqformat] -z=\"entrez query\" -v 1 -t output_format -i
  24 seqfile\n";
  25
  26 my ($prog, $db, $expect,$method) = ( 'blastp', 'nr', '10', 'Blast');
  27
  28 my ($sequencefile,$sequenceformat,$help, $entrez, $outputformat,
  29 $verbose) = (undef, 'fasta',undef, undef, undef, 1);
  30
  31 &GetOptions('prog|p=s'               => \$prog,
  32                                 'db|d=s'                 => \$db,
  33                                 'expect|e=s'             => \$expect,
  34                                 'blsmod|module|method=s' => \$method,
  35                                 'input|i=s'              => \$sequencefile,
  36                                 'format|f=s'             => \$sequenceformat,
  37                                 'help|h'                 => \$help,
  38                                 'entrez|z=s'             => \$entrez,
  39                                 'output_format|t=s'      => \$outputformat,
  40                                 'verbose|v=s'            => \$verbose
  41                           );
  42
  43 if( $help ) {
  44     exec('perldoc', $0);
  45     die;
  46 }
  47
  48 if( !defined $prog ) {
  49     die($USAGE . "\n\tMust specify a valid program name ([t]blast[pxn])\n");
  50 }
  51 if( !defined $db ) {
  52     die($USAGE . "\n\tMust specify a db to search\n");
  53 }
  54 if( !defined $sequencefile ) {
  55     die($USAGE . "\n\tMust specify an input file\n");
  56 }
  57
  58 my $blastfactory = new Bio::Tools::Run::RemoteBlast ('-prog' => $prog,
  59                                                                                                                                           '-data'      => $db,
  60                                                                                                                                           '-expect'    => $expect,
  61                                                                                                                                           'readmethod' => $method,
  62                                                                                                                                          );
  63
  64 if ($entrez) {
  65   if ($verbose) {
  66     print "Entrez query (submission side): $entrez\n";
  67   }
  68   #$Bio::Tools::Run::RemoteBlast::RETRIEVALHEADER{ FORMAT_ENTREZ_QUERY} = $entrez;
  69   $Bio::Tools::Run::RemoteBlast::HEADER{ ENTREZ_QUERY } = $entrez;
  70 }
  71 if ($outputformat) {
  72   print "Don't use output format type; it doesn't work.\n";
  73   $Bio::Tools::Run::RemoteBlast::RETRIEVALHEADER{ FORMAT_TYPE } = $outputformat;
  74 }
  75
  76 # submit_blast can only currenly handle fasta format files so I'll
  77 # preprocess outside of the module but I'd rather be sure here
  78
  79 my $input;
  80 if( $sequenceformat !~ /fasta/ ) {
  81         my @seqs;
  82         my $seqio = new Bio::SeqIO('-format' => $sequenceformat,
  83                                                                                 '-file'   => $sequencefile );
  84         while( my $seq = $seqio->next_seq() ) {
  85                 push @seqs, $seq;
  86         }
  87         $input = \@seqs;
  88 } else {
  89         $input = $sequencefile;
  90 }
  91
  92 my $r = $blastfactory->submit_blast($input);
  93 #my $r = $factory->submit_blast(?amino.fa?);
  94
  95 print STDERR "waiting...\n" if( $verbose > 0 );
  96 while ( my @rids = $blastfactory->each_rid ) {
  97   foreach my $rid ( @rids ) {
  98     my $rc = $blastfactory->retrieve_blast($rid);
  99     if( !ref($rc) ) {
 100       if( $rc < 0 ) {
 101         $blastfactory->remove_rid($rid);
 102       }
 103       print STDERR "   checking $rid\n" if ( $verbose > 0 );
 104       sleep 5;
 105     } else {
 106       my $result = $rc->next_result();
 107       #save the output
 108       my $filename = $result->query_name()."\.out";
 109       $blastfactory->save_output($filename);
 110       $blastfactory->remove_rid($rid);
 111       print "\nQuery Name: ", $result->query_name(), "\n";
 112       while ( my $hit = $result->next_hit ) {
 113         next unless ( $verbose > 0);
 114         print "\thit name is ", $hit->name, "\n";
 115         while( my $hsp = $hit->next_hsp ) {
 116           print "\t\tscore is ", $hsp->score, "\n";
 117         }
 118       }
 119     }
 120   }
 121   print STDERR scalar(@rids) . " left\n";
 122 }
 123
 124
 125 __END__
 126
 127 #
 128 # BioPerl module for bp_remote_blast.pl
 129 #
 130 # Cared for by Jason Stajich <jason-at-bioperl-dot-org>
 131 #
 132 # Copyright Jason Stajich
 133 #
 134 # You may distribute this module under the same terms as perl itself
 135
 136 # POD documentation - main docs before the code
 137
 138 =head1 NAME
 139
 140 bp_remote_blast.pl - script for submitting jobs to a remote blast server
 141 (ncbi blast queue at this time)
 142
 143 =head1 SYNOPSIS
 144
 145 % bp_remote_blast.pl -p blastp -d ecoli -e 1e-5 -i myseqs.fa
 146
 147 =head1 DESCRIPTION
 148
 149 This module will run a remote blast on a set of sequences by
 150 submitting them to the NCBI blast queue and printing the output of the
 151 request.
 152
 153 =head1 FEEDBACK
 154
 155 =head2 Mailing Lists
 156
 157 User feedback is an integral part of the evolution of this and other
 158 Bioperl modules. Send your comments and suggestions preferably to
 159 the Bioperl mailing list.  Your participation is much appreciated.
 160
 161   bioperl-l@bioperl.org                  - General discussion
 162   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
 163
 164 =head2 Reporting Bugs
 165
 166 Report bugs to the Bioperl bug tracking system to help us keep track
 167 the bugs and their resolution.  Bug reports can be submitted via the
 168 web:
 169
 170   https://github.com/bioperl/bioperl-live/issues
 171
 172 =head1 AUTHOR - Jason Stajich
 173
 174 Email jason-at-bioperl-dot-org
 175
 176 =cut