lib/Bio/Tools/Run/StandAloneNCBIBlast.pm

   1 #
   2 # BioPerl module for Bio::Tools::Run::StandAloneBlast
   3 #
   4 # Copyright Peter Schattner
   5 #
   6 # You may distribute this module under the same terms as perl itself
   7
   8 # POD documentation - main docs before the code
   9
  10 =head1 NAME
  11
  12 Bio::Tools::Run::StandAloneNCBIBlast - Object for the local execution
  13 of the NCBI BLAST program suite (blastall, blastpgp, bl2seq). With
  14 experimental support for NCBI rpsblast.
  15
  16 =head1 SYNOPSIS
  17
  18  # Do not use directly; see Bio::Tools::Run::StandAloneBlast
  19
  20 =head1 DESCRIPTION
  21
  22 See Bio::Tools::Run::StandAloneBlast
  23
  24 =head1 FEEDBACK
  25
  26 =head2 Mailing Lists
  27
  28 User feedback is an integral part of the evolution of this and other
  29 Bioperl modules. Send your comments and suggestions preferably to one
  30 of the Bioperl mailing lists.  Your participation is much appreciated.
  31
  32   bioperl-l@bioperl.org                  - General discussion
  33   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  34
  35 =head2 Support
  36
  37 Please direct usage questions or support issues to the mailing list:
  38
  39 I<bioperl-l@bioperl.org>
  40
  41 rather than to the module maintainer directly. Many experienced and
  42 reponsive experts will be able look at the problem and quickly
  43 address it. Please include a thorough description of the problem
  44 with code and data examples if at all possible.
  45
  46 =head2 Reporting Bugs
  47
  48 Report bugs to the Bioperl bug tracking system to help us keep track
  49 the bugs and their resolution.  Bug reports can be submitted via
  50 the web:
  51
  52   https://github.com/bioperl/bioperl-live/issues
  53
  54 =head1 AUTHOR - Peter Schattner
  55
  56 Email schattner at alum.mit.edu
  57
  58 =head1 MAINTAINER - Torsten Seemann
  59
  60 Email torsten at infotech.monash.edu.au
  61
  62 =head1 CONTRIBUTORS
  63
  64 Sendu Bala  bix@sendu.me.uk (reimplementation)
  65
  66 =head1 APPENDIX
  67
  68 The rest of the documentation details each of the object
  69 methods. Internal methods are usually preceded with a _
  70
  71 =cut
  72
  73 package Bio::Tools::Run::StandAloneNCBIBlast;
  74
  75 use strict;
  76 use warnings;
  77
  78 use base qw(Bio::Tools::Run::StandAloneBlast);
  79
  80 our $AUTOLOAD;
  81 our $DEFAULTREADMETHOD = 'BLAST';
  82
  83 # If local BLAST databases are not stored in the standard
  84 # /data directory, the variable BLASTDATADIR will need to be
  85 # set explicitly
  86 our $DATADIR = $Bio::Tools::Run::StandAloneBlast::DATADIR;
  87
  88 our %GENERAL_PARAMS  = (i => 'input',
  89                         o => 'outfile',
  90                         p => 'program',
  91                         d => 'database');
  92 our @BLASTALL_PARAMS = qw(A B C D E F G K L M O P Q R S W X Y Z a b e f l m q r t v w y z n);
  93 our @BLASTALL_SWITCH = qw(I g J T U n V s);
  94 our @BLASTPGP_PARAMS = qw(A B C E F G H I J K L M N O P Q R S T U W X Y Z a b c e f h j k l m q s t u v y z);
  95 our @RPSBLAST_PARAMS = qw(F I J L N O P T U V X Y Z a b e l m v y z);
  96 our @BL2SEQ_PARAMS   = qw(A D E F G I J M S T U V W X Y a e g j m q r t);
  97
  98 our @OTHER_PARAMS = qw(_READMETHOD);
  99
 100
 101 =head2 new
 102
 103  Title   : new
 104  Usage   : my $obj = Bio::Tools::Run::StandAloneBlast->new();
 105  Function: Builds a newBio::Tools::Run::StandAloneBlast object
 106  Returns : Bio::Tools::Run::StandAloneBlast
 107  Args    : -quiet => boolean # make program execution quiet
 108            -_READMETHOD => 'BLAST' (default, synonym 'SearchIO') || 'blast_pull'
 109                            # the parsing method, case insensitive
 110
 111 Essentially all BLAST parameters can be set via StandAloneBlast.pm.
 112 Some of the most commonly used parameters are listed below. All
 113 parameters have defaults and are optional except for -p in those programs that
 114 have it. For a complete listing of settable parameters, run the relevant
 115 executable BLAST program with the option "-" as in blastall -
 116 Note that the input parameters (-i, -j, -input) should not be set directly by
 117 you: this module sets them when you call one of the executable methods.
 118
 119 Blastall
 120
 121   -p  Program Name [String]
 122         Input should be one of "blastp", "blastn", "blastx",
 123         "tblastn", or "tblastx".
 124   -d  Database [String] default = nr
 125         The database specified must first be formatted with formatdb.
 126         Multiple database names (bracketed by quotations) will be accepted.
 127         An example would be -d "nr est"
 128   -e  Expectation value (E) [Real] default = 10.0
 129   -o  BLAST report Output File [File Out]  Optional,
 130             default = ./blastreport.out ; set by StandAloneBlast.pm
 131   -S  Query strands to search against database (for blast[nx], and tblastx). 3 is both, 1 is top, 2 is bottom [Integer]
 132             default = 3
 133
 134 Blastpgp (including Psiblast)
 135
 136   -j  is the maximum number of rounds (default 1; i.e., regular BLAST)
 137   -h  is the e-value threshold for including sequences in the
 138             score matrix model (default 0.001)
 139   -c  is the "constant" used in the pseudocount formula specified in the paper (default 10)
 140   -B  Multiple alignment file for PSI-BLAST "jump start mode"  Optional
 141   -Q  Output File for PSI-BLAST Matrix in ASCII [File Out]  Optional
 142
 143 rpsblast
 144
 145   -d  Database [String] default = (none - you must specify a database)
 146         The database specified must first be formatted with formatdb.
 147         Multiple database names (bracketed by quotations) will be accepted.
 148         An example would be -d "Cog Smart"
 149   -e  Expectation value (E) [Real] default = 10.0
 150   -o  BLAST report Output File [File Out]  Optional,
 151             default = ./blastreport.out ; set by StandAloneBlast.pm
 152
 153 Bl2seq
 154
 155   -p  Program name: blastp, blastn, blastx. For blastx 1st argument should be nucleotide [String]
 156     default = blastp
 157   -o  alignment output file [File Out] default = stdout
 158   -e  Expectation value (E) [Real]  default = 10.0
 159   -S  Query strands to search against database (blastn only).  3 is both, 1 is top, 2 is bottom [Integer]
 160     default = 3
 161
 162 =cut
 163
 164 sub new {
 165     my ($caller, @args) = @_;
 166     my $self = $caller->SUPER::new(@args);
 167
 168     # StandAloneBlast is special in that "one can modify the name of
 169     # the (ncbi) BLAST parameters as desired as long as the initial letter (and
 170     # case) of the parameter are preserved". We handle this by truncating input
 171     # args to their first char
 172     my %args = @args;
 173     @args = ();
 174     while (my ($attr, $value) = each %args) {
 175         $attr =~ s/^-//;
 176         $attr = substr($attr, 0, 1) unless $attr =~ /^_/;
 177         push(@args, $attr, $value);
 178     }
 179
 180     $self->_set_from_args(\@args, -methods => {(map { $_ => $GENERAL_PARAMS{$_} } keys %GENERAL_PARAMS),
 181                                                (map { $_ => $_ } (@OTHER_PARAMS,
 182                                                                   @BLASTALL_PARAMS,
 183                                                                   @BLASTALL_SWITCH,
 184                                                                   @BLASTPGP_PARAMS,
 185                                                                   @RPSBLAST_PARAMS,
 186                                                                   @BL2SEQ_PARAMS))},
 187                                   -code => { map { $_ => 'my $self = shift;
 188                                                           if (@_) {
 189                                                               my $value = shift;
 190                                                               if ($value && $value ne \'F\') {
 191                                                                   $value = \'T\';
 192                                                               }
 193                                                               else {
 194                                                                   $value = \'F\';
 195                                                               }
 196                                                               $self->{\'_\'.$method} = $value;
 197                                                           }
 198                                                           return $self->{\'_\'.$method} || return;' } @BLASTALL_SWITCH },  # these methods can take boolean or 'T' and 'F'
 199                                   -create => 1,
 200                                   -force => 1,
 201                                   -case_sensitive => 1);
 202
 203     my ($tfh, $tempfile) = $self->io->tempfile();
 204     my $outfile = $self->o || $self->outfile || $tempfile;
 205     $self->o($outfile);
 206     close($tfh);
 207
 208     $self->_READMETHOD($DEFAULTREADMETHOD) unless $self->_READMETHOD;
 209
 210     return $self;
 211 }
 212
 213 # StandAloneBlast is special in that "one can modify the name of
 214 # the (ncbi) BLAST parameters as desired as long as the initial letter (and
 215 # case) of the parameter are preserved". We handle this with AUTOLOAD
 216 # redirecting to the automatically created methods from _set_from_args() !
 217 sub AUTOLOAD {
 218     my $self = shift;
 219     my $attr = $AUTOLOAD;
 220     $attr =~ s/.*:://;
 221
 222     my $orig = $attr;
 223
 224     $attr = substr($attr, 0, 1);
 225
 226     $self->can($attr) || $self->throw("Unallowed parameter: $orig !");
 227
 228     return $self->$attr(@_);
 229 }
 230
 231 =head2  blastall
 232
 233  Title   : blastall
 234  Usage   :  $blast_report = $factory->blastall('t/testquery.fa');
 235         or
 236                $input = Bio::Seq->new(-id=>"test query",
 237                                       -seq=>"ACTACCCTTTAAATCAGTGGGGG");
 238                $blast_report = $factory->blastall($input);
 239         or
 240               $seq_array_ref = \@seq_array;
 241          # where @seq_array is an array of Bio::Seq objects
 242               $blast_report = $factory->blastall($seq_array_ref);
 243  Returns : Reference to a Blast object containing the blast report.
 244  Args    : Name of a file or Bio::Seq object or an array of
 245            Bio::Seq object containing the query sequence(s).
 246            Throws an exception if argument is not either a string
 247            (eg a filename) or a reference to a Bio::Seq object
 248            (or to an array of Seq objects).  If argument is string,
 249            throws exception if file corresponding to string name can
 250            not be found.
 251
 252 =cut
 253
 254 sub blastall {
 255     my ($self, $input1) = @_;
 256     $self->io->_io_cleanup();
 257     my $executable = 'blastall';
 258
 259     # Create input file pointer
 260     my $infilename1 = $self->_setinput($executable, $input1) || $self->throw("$input1 not Bio::Seq object or array of Bio::Seq objects or file name!");
 261     $self->i($infilename1);
 262
 263     my $blast_report = $self->_generic_local_blast($executable);
 264 }
 265
 266 =head2  blastpgp
 267
 268  Title   : blastpgp
 269  Usage   :  $blast_report = $factory-> blastpgp('t/testquery.fa');
 270         or
 271                $input = Bio::Seq->new(-id=>"test query",
 272                                       -seq=>"ACTADDEEQQPPTCADEEQQQVVGG");
 273                $blast_report = $factory->blastpgp ($input);
 274         or
 275               $seq_array_ref = \@seq_array;
 276          # where @seq_array is an array of Bio::Seq objects
 277               $blast_report = $factory-> blastpgp(\@seq_array);
 278  Returns : Reference to a Bio::SearchIO object containing the blast report
 279  Args    : Name of a file or Bio::Seq object. In psiblast jumpstart
 280            mode two additional arguments are required: a SimpleAlign
 281            object one of whose elements is the query and a "mask" to
 282            determine how BLAST should select scoring matrices see
 283            DESCRIPTION above for more details.
 284
 285            Throws an exception if argument is not either a string
 286            (eg a filename) or a reference to a Bio::Seq object
 287            (or to an array of Seq objects).  If argument is string,
 288            throws exception if file corresponding to string name can
 289            not be found.
 290  Returns : Reference to Bio::SearchIO object containing the blast report.
 291
 292 =cut
 293
 294 sub blastpgp {
 295     my $self = shift;
 296     my $executable = 'blastpgp';
 297     my $input1 = shift;
 298     my $input2 = shift;
 299     # used by blastpgp's -B option to specify which
 300     # residues are position aligned
 301     my $mask = shift;
 302
 303     my ($infilename1, $infilename2 ) = $self->_setinput($executable,
 304                                                         $input1, $input2,
 305                                                         $mask);
 306     if (!$infilename1) {$self->throw("$input1 not Bio::Seq object or array of Bio::Seq objects or file name!");}
 307     $self->i($infilename1);     # set file name of sequence to be blasted to inputfilename1 (-i param of blastpgp)
 308     if ($input2) {
 309         unless ($infilename2) {$self->throw("$input2 not SimpleAlign Object in pre-aligned psiblast\n");}
 310         $self->B($infilename2); # set file name of partial alignment to inputfilename2 (-B param of blastpgp)
 311     }
 312
 313     my $blast_report = $self->_generic_local_blast($executable);
 314 }
 315
 316 =head2  rpsblast
 317
 318  Title   : rpsblast
 319  Usage   :  $blast_report = $factory->rpsblast('t/testquery.fa');
 320         or
 321                $input = Bio::Seq->new(-id=>"test query",
 322                                       -seq=>"MVVLCRADDEEQQPPTCADEEQQQVVGG");
 323                $blast_report = $factory->rpsblast($input);
 324         or
 325               $seq_array_ref = \@seq_array;
 326          # where @seq_array is an array of Bio::Seq objects
 327               $blast_report = $factory->rpsblast(\@seq_array);
 328  Args    : Name of a file or Bio::Seq object or an array of
 329            Bio::Seq object containing the query sequence(s).
 330            Throws an exception if argument is not either a string
 331            (eg a filename) or a reference to a Bio::Seq object
 332            (or to an array of Seq objects).  If argument is string,
 333            throws exception if file corresponding to string name can
 334            not be found.
 335  Returns : Reference to a Bio::SearchIO object containing the blast report
 336
 337 =cut
 338
 339 sub rpsblast {
 340     my ($self, $input1) = @_;
 341     $self->io->_io_cleanup();
 342     my $executable = 'rpsblast';
 343
 344     # Create input file pointer
 345     my $infilename1 = $self->_setinput($executable, $input1) || $self->throw("$input1 not Bio::Seq object or array of Bio::Seq objects or file name!");
 346     $self->i($infilename1);
 347
 348     my $blast_report = $self->_generic_local_blast($executable);
 349 }
 350
 351 =head2   bl2seq
 352
 353  Title   : bl2seq
 354  Usage   : $factory-> bl2seq('t/seq1.fa', 't/seq2.fa');
 355         or
 356           $input1 = Bio::Seq->new(-id=>"test query1",
 357                                   -seq=>"ACTADDEEQQPPTCADEEQQQVVGG");
 358           $input2 = Bio::Seq->new(-id=>"test query2",
 359                                   -seq=>"ACTADDEMMMMMMMDEEQQQVVGG");
 360           $blast_report = $factory->bl2seq ($input1,  $input2);
 361  Returns : Reference to a BPbl2seq object containing the blast report.
 362  Args    : Names of 2 files  or 2 Bio::Seq objects containing the
 363            sequences to be aligned by bl2seq.
 364
 365            Throws an exception if argument is not either a pair of
 366            strings (eg filenames) or references to Bio::Seq objects.
 367            If arguments are strings, throws exception if files
 368            corresponding to string names can not be found.
 369
 370 =cut
 371
 372 sub bl2seq {
 373     my $self = shift;
 374     my $executable = 'bl2seq';
 375     my $input1 = shift;
 376     my $input2 = shift;
 377
 378     # Create input file pointer
 379     my ($infilename1, $infilename2 ) = $self->_setinput($executable,
 380                                                           $input1, $input2);
 381     if (!$infilename1){$self->throw(" $input1  not Seq Object or file name!");}
 382     if (!$infilename2){$self->throw("$input2  not Seq Object or file name!");}
 383
 384     $self->i($infilename1);     # set file name of first sequence to
 385                             # be aligned to inputfilename1
 386                             # (-i param of bl2seq)
 387     $self->j($infilename2);     # set file name of first sequence to
 388                             # be aligned to inputfilename2
 389                             # (-j param of bl2seq)
 390
 391     my $blast_report = $self->_generic_local_blast($executable);
 392 }
 393
 394 =head2  _generic_local_blast
 395
 396  Title   : _generic_local_blast
 397  Usage   : internal function not called directly
 398  Returns : Bio::SearchIO
 399  Args    : Reference to calling object and name of BLAST executable
 400
 401 =cut
 402
 403 sub _generic_local_blast {
 404     my $self = shift;
 405     my $executable = shift;
 406
 407     # Create parameter string to pass to Blast program
 408     my $param_string = $self->_setparams($executable);
 409
 410     # run Blast
 411     my $blast_report = $self->_runblast($executable, $param_string);
 412 }
 413
 414 =head2  _runblast
 415
 416  Title   :  _runblast
 417  Usage   :  Internal function, not to be called directly
 418  Function:   makes actual system call to Blast program
 419  Example :
 420  Returns : Report Bio::SearchIO object in the appropriate format
 421  Args    : Reference to calling object, name of BLAST executable,
 422            and parameter string for executable
 423
 424 =cut
 425
 426 sub _runblast {
 427         my ($self, $executable, $param_string) = @_;
 428         my ($blast_obj, $exe);
 429         if (! ($exe = $self->executable($executable)) ) {
 430                 $self->warn("cannot find path to $executable");
 431                 return;
 432         }
 433
 434     # Use double quotes if executable path have empty spaces
 435     if ($exe =~ m/ /) {
 436         $exe = "\"$exe\"";
 437     }
 438         my $commandstring = $exe.$param_string;
 439
 440         $self->debug("$commandstring\n");
 441         system($commandstring) && $self->throw("$executable call crashed: $? | $! | $commandstring\n");
 442
 443     # set significance cutoff to set expectation value or default value
 444         # (may want to make this value vary for different executables)
 445         my $signif = $self->e() || 1e-5;
 446
 447     # get outputfilename
 448         my $outfile = $self->o();
 449
 450     # this should allow any blast SearchIO parser (not just 'blast_pull' or 'blast',
 451     # but 'blastxml' and 'blasttable').  Fall back to 'blast' if not stipulated.
 452     my $method = $self->_READMETHOD;
 453         if ($method =~ /^(?:blast|SearchIO)/i )  {
 454         $method = 'blast' if $method =~ m{SearchIO}i;
 455                 $blast_obj = Bio::SearchIO->new(-file => $outfile,
 456                                         -format => $method);
 457         }
 458     # should these be here?  They have been deprecated...
 459     elsif ($method =~ /BPlite/i ) {
 460                 if ($executable =~ /bl2seq/i)  {
 461                         # Added program info so BPbl2seq can compute strand info
 462                         $self->throw("Use of Bio::Tools::BPbl2seq is deprecated; use Bio::SearchIO modules instead");
 463                 }
 464         elsif ($executable =~ /blastpgp/i && defined $self->j() && $self->j() > 1) {
 465                         $self->throw("Use of Bio::Tools::BPpsilite is deprecated; use Bio::SearchIO modules instead");
 466                 }
 467         elsif ($executable =~ /blastall|rpsblast/i) {
 468                         $self->throw("Use of Bio::Tools::BPlite is deprecated; use Bio::SearchIO modules instead");
 469                 }
 470         else {
 471                         $self->warn("Unrecognized executable $executable");
 472                 }
 473         }
 474     else {
 475                 $self->warn("Unrecognized readmethod $method");
 476         }
 477
 478         return $blast_obj;
 479 }
 480
 481 =head2  _setparams
 482
 483  Title   : _setparams
 484  Usage   : Internal function, not to be called directly
 485  Function: Create parameter inputs for Blast program
 486  Example :
 487  Returns : parameter string to be passed to Blast
 488  Args    : Reference to calling object and name of BLAST executable
 489
 490 =cut
 491
 492 sub _setparams {
 493     my ($self, $executable) = @_;
 494     my ($attr, $value, @execparams);
 495
 496     if    ($executable eq 'blastall') { @execparams = (@BLASTALL_PARAMS,
 497                                                        @BLASTALL_SWITCH); }
 498     elsif ($executable eq 'blastpgp') { @execparams =  @BLASTPGP_PARAMS;  }
 499     elsif ($executable eq 'rpsblast') { @execparams =  @RPSBLAST_PARAMS;  }
 500     elsif ($executable eq 'bl2seq'  ) { @execparams =  @BL2SEQ_PARAMS;    }
 501
 502     # we also have all the general params
 503     push(@execparams, keys %GENERAL_PARAMS);
 504
 505     my $database = $self->d;
 506     if ($database && $executable ne 'bl2seq') {
 507         # Need to prepend datadirectory to database name
 508         my @dbs = split(/ /, $database);
 509         for my $i (0..$#dbs) {
 510             # (works with multiple databases)
 511             if (! (-e $dbs[$i].".nin" || -e $dbs[$i].".pin") &&
 512                 ! (-e $dbs[$i].".nal" || -e $dbs[$i].".pal") ) {
 513                 $dbs[$i] = File::Spec->catdir($DATADIR, $dbs[$i]);
 514             }
 515         }
 516         $self->d('"'.join(" ", @dbs).'"');
 517     }
 518
 519     # workaround for problems with shell metacharacters [bug 2707]
 520     # simply quoting does not always work!
 521     my $tmp = $self->o;
 522     $self->o(quotemeta($tmp)) if ($tmp && $^O !~ /^MSWin/);
 523
 524     my $param_string = $self->SUPER::_setparams(-params => [@execparams],
 525                                                 -dash => 1);
 526
 527     $self->o($tmp) if ($tmp && $^O !~ /^MSWin/);
 528
 529     $self->d($database) if $database;
 530
 531     if ($self->quiet()) {
 532         $param_string .= ' 2> '.File::Spec->devnull;
 533     }
 534
 535     return $param_string;
 536 }
 537
 538 1;