Bio/Tools/Primer3.pm

   1 #
   2 # BioPerl module for Bio::Tools::Primer3
   3 #
   4 # Copyright (c) 2003 bioperl, Rob Edwards. All Rights Reserved.
   5 #           This module is free software; you can redistribute it and/or
   6 #           modify it under the same terms as Perl itself.
   7 #
   8 # Copyright Rob Edwards
   9 #
  10 # You may distribute this module under the same terms as perl itself
  11 # POD documentation - main docs before the code
  12
  13 =head1 NAME
  14
  15 Bio::Tools::Primer3 - Create input for and work with the output from
  16 the program primer3
  17
  18 =head1 SYNOPSIS
  19
  20  # parse primer3 output to get some data
  21  # this is also called from Bio::Tools::Run::Primer3
  22  use Bio::Tools::Primer3;
  23
  24  # read a primer3 output file
  25  my $p3 = Bio::Tools::Primer3->new(-file=>"data/primer3_output.txt");
  26
  27  # how many results were there?
  28  my $num = $p3->number_of_results;
  29  print "There were $num results\n";
  30
  31  # get all the results
  32  my $all_results = $p3->all_results;
  33  print "ALL the results\n";
  34  foreach my $key (keys %{$all_results}) {
  35     print "$key\t${$all_results}{$key}\n";
  36  }
  37
  38  # get specific results
  39  my $result1 = $p3->primer_results(1);
  40  print "The first primer is\n";
  41  foreach my $key (keys %{$result1}) {
  42     print "$key\t${$result1}{$key}\n";
  43  }
  44
  45  # get the results as a Bio::Seq::PrimedSeq stream
  46  my $primer = $p3->next_primer;
  47  print "The left primer in the stream is ",
  48    $primer->get_primer('-left_primer')->seq->seq, "\n";
  49
  50 =head1 DESCRIPTION
  51
  52 Bio::Tools::Primer3 creates the input files needed to design primers using
  53 primer3 and provides mechanisms to access data in the primer3 output files.
  54
  55 This module provides a bioperl interface to the program primer3. See
  56 http://www-genome.wi.mit.edu/genome_software/other/primer3.html
  57 for details and to download the software.
  58
  59 This module is based on one written by Chad Matsalla
  60 (bioinformatics1@dieselwurks.com)
  61
  62 I have ripped some of his code, and added a lot of my own. I hope he
  63 is not mad at me!
  64
  65 This is probably best run in one of the two following ways:
  66
  67   i. To parse the output from Bio::Tools::Run::Primer3.
  68      You will most likely just use next_primer to get the results from
  69      Bio::Tools::Run::Primer3.
  70   ii. To parse the output of primer3 handed to it as a file name.
  71
  72 =head1 FEEDBACK
  73
  74 =head2 Mailing Lists
  75
  76 User feedback is an integral part of the evolution of this and other
  77 Bioperl modules. Send your comments and suggestions preferably to one
  78 of the Bioperl mailing lists.  Your participation is much appreciated.
  79
  80   bioperl-l@bioperl.org                  - General discussion
  81   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  82
  83 =head2 Support
  84
  85 Please direct usage questions or support issues to the mailing list:
  86
  87 I<bioperl-l@bioperl.org>
  88
  89 rather than to the module maintainer directly. Many experienced and
  90 reponsive experts will be able look at the problem and quickly
  91 address it. Please include a thorough description of the problem
  92 with code and data examples if at all possible.
  93
  94 =head2 Reporting Bugs
  95
  96 Report bugs to the Bioperl bug tracking system to help us keep track
  97 the bugs and their resolution.  Bug reports can be submitted via the web:
  98
  99   https://github.com/bioperl/bioperl-live/issues
 100
 101 =head1 AUTHOR -
 102
 103   Rob Edwards
 104
 105   redwards@utmem.edu
 106
 107   Based heavily on work of
 108
 109   Chad Matsalla
 110
 111   bioinformatics1@dieselwurks.com
 112
 113 =head1 CONTRIBUTORS
 114
 115   Brian Osborne bosborne at alum.mit.edu
 116
 117 =head1 APPENDIX
 118
 119 The rest of the documentation details each of the object methods.
 120 Internal methods are usually preceded with a _
 121
 122 =cut
 123
 124 # Let the code begin...
 125
 126 package Bio::Tools::Primer3;
 127
 128 use strict;
 129 use Bio::Seq;
 130 use Bio::Seq::PrimedSeq;
 131 use Bio::SeqFeature::Primer;
 132
 133 use vars qw($AUTOLOAD @PRIMER3_PARAMS %OK_FIELD $ID);
 134
 135 BEGIN {
 136     @PRIMER3_PARAMS = qw(results seqobject);
 137     foreach my $attr (@PRIMER3_PARAMS) {$OK_FIELD{$attr}++}
 138 }
 139
 140
 141 use base qw(Bio::Root::Root Bio::Root::IO);
 142
 143
 144 sub AUTOLOAD {
 145     my $self = shift;
 146     my $attr = $AUTOLOAD;
 147     $attr =~ s/.*:://;
 148     $self->throw("Unallowed parameter: $attr !") unless $OK_FIELD{$attr};
 149     $self->{$attr} = shift if @_;
 150     return $self->{$attr};
 151 }
 152
 153
 154 =head2 new
 155
 156   Title   : new()
 157   Usage   : my $primer3 = Bio::Tools::Primer3->new(-file=>$file);
 158   Function: Parse primer3 output
 159   Returns : Does not return anything. If called with a filename will
 160             allow you to retrieve the results
 161   Args    : -file (optional) file of primer3 results to parse -verbose
 162             (optional) set verbose output
 163   Notes   :
 164
 165 =cut
 166
 167 sub new {
 168     my($class,%args) = @_;
 169     my $self = $class->SUPER::new(%args);
 170     if ($args{'-file'}) {
 171         $self->_readfile($args{'-file'});
 172     }
 173     if ($args{'-verbose'}) {
 174         $self->{'verbose'} = 1;
 175     }
 176     return $self;
 177 }
 178
 179
 180 =head2 number_of_results
 181
 182   Title   : number_of_results()
 183   Usage   : my $count = $primer3->number_of_results();
 184   Function: Retrieve the number of primers returned from Primer3.
 185   Returns : A scalar
 186   Args    : None
 187   Notes   : This returns the count of the primers returned by Primer3
 188              (aka how many of them there are).
 189              This is one more than the maximum offset into the zero
 190              based list of primers that is accessed by primer_results().
 191
 192 =cut
 193
 194 sub number_of_results {
 195     my $self = shift;
 196     return $self->{'maximum_primers_returned'} + 1;
 197 }
 198
 199
 200 =head2 all_results
 201
 202   Title   : all_results()
 203   Usage   : my $results = $primer3->all_results();
 204                or
 205             my $results = $primer3->all_results('primer3 result name', 'other results');
 206   Function: Retrieve the results returned from Primer3.
 207   Returns : A reference to a hash
 208   Args    : Optional array of specific results to retrieve
 209
 210 =cut
 211
 212 sub all_results {
 213     my ($self, @results) = @_;
 214     my %hash;
 215     if (@results) {
 216         # we only want a few things
 217         foreach my $result (@results) {
 218             $hash{$result} = $self->{'results'}->$result;
 219         }
 220     } else {
 221         foreach my $result (keys %{$self->{'results'}}) {
 222             $hash{$result}=$self->{'results'}->{$result};
 223         }
 224     }
 225
 226     return \%hash;
 227 }
 228
 229
 230 =head2 primer_results
 231
 232   Title   : primer_results()
 233   Usage   : my $results = $primer3->primer_results(2); # results for third primer
 234   Function: Retrieve the results returned from Primer3 for specific primer pairs.
 235   Returns : A reference to a hash
 236   Args    : A number between 0 and the maximum number of primers to retrieve
 237
 238 =cut
 239
 240 sub primer_results {
 241     my ($self, $toget) = @_;
 242     if ($toget > $self->{'maximum_primers_returned'}) {
 243         $self->warn("Didn't get any results for $toget");
 244         return 0;
 245     } else {
 246         return \%{$self->{'results_by_number'}->{$toget}};
 247     }
 248 }
 249
 250
 251 =head2 _readfile
 252
 253   Title   : _readfile()
 254   Usage   : $self->_readfile();
 255   Function: An internal function that reads a file and sets up the results
 256   Returns : Nothing.
 257   Args    : None
 258   Notes   :
 259
 260 =cut
 261
 262 sub _readfile {
 263     my ($self, $file) = @_;
 264     $self->_initialize_io(-file=>$file);
 265     my $line;
 266     my $id='primer 3 parsed results'; # hopefully we'll get this, but we can set a temp id in case not.
 267     while (defined($line = $self->_readline()) ) {
 268         chomp $line;
 269         next unless ($line);
 270         my ($return, $value) = split /=/, $line;
 271         if (uc($return) eq "SEQUENCE") {
 272             $self->{seqobject} = Bio::Seq->new(-seq=>$value, $id=>$id);
 273             next;
 274         }
 275         if (uc($return) eq "PRIMER_SEQUENCE_ID") {
 276             if ($self->{seqobject}) {$self->{seqobject}->id($value)} else {$id=$value}
 277         }
 278
 279         $self->{'results'}->{$return} = $value;
 280     }
 281
 282     # convert the results to individual results
 283     $self->_separate();
 284 }
 285
 286
 287 =head2 next_primer
 288
 289   Title   : next_primer()
 290   Usage   : while (my $primed_seq  = $primer3->next_primer()) {
 291   Function: Retrieve the primed sequence and a primer pair, one at a time
 292   Returns : Returns a Bio::Seq::PrimedSeq object, one at a time
 293   Args    : None
 294   Notes   : Use $primed_seq->annotated_seq to get an annotated sequence
 295             object you can write out.
 296
 297 =cut
 298
 299 sub next_primer {
 300     my $self = shift;
 301     # here we are going to convert the primers to Bio::SeqFeature::Primer objects
 302     # and the primer/sequence to Bio::Seq::PrimedSeq objects
 303     # the problem at the moment is that PrimedSeq can only take one sequence/primer pair, and
 304     # yet for each sequence we can have lots of primer pairs. We need a way to overcome this.
 305     # at the moment we can do this as a stream, I guess.
 306
 307     $self->warn("No primers were found for: ".$self->{'seqobject'}->{'primary_id'})
 308       if (! $self->number_of_results);
 309
 310     $self->{'next_to_return'} = 0 unless ($self->{'next_to_return'});
 311     return if ($self->{'next_to_return'} >= $self->number_of_results);
 312     my $results = $self->primer_results($self->{'next_to_return'});
 313
 314     $self->throw("No left primer sequence") unless (${$results}{'PRIMER_LEFT_SEQUENCE'});
 315     $self->throw("No right primer sequence") unless (${$results}{'PRIMER_RIGHT_SEQUENCE'});
 316     $self->throw("No target sequence") unless ($self->{'seqobject'});
 317
 318     my $left_seq = Bio::SeqFeature::Primer->new(
 319          -id         => 'left_primer',
 320          -seq        => ${$results}{'PRIMER_LEFT_SEQUENCE'},
 321          -display_id => ($self->{'next_to_return'} + 1),
 322     );
 323
 324     my $right_seq = Bio::SeqFeature::Primer->new(
 325              -id         => "right_primer",
 326              -seq        => ${$results}{'PRIMER_RIGHT_SEQUENCE'},
 327              -display_id => ($self->{'next_to_return'} + 1) );
 328
 329     # add data to the Primer objects
 330     for my $key (%$results) {
 331         # skip the primer sequence data, already added above
 332         next if ($key =~ /PRIMER_(LEFT|RIGHT)_SEQUENCE/i );
 333         if ($key =~ /PRIMER_LEFT/i) {
 334             $left_seq->add_tag_value($key, $$results{$key});
 335         } elsif ($key =~ /PRIMER_RIGHT/i) {
 336             $right_seq->add_tag_value($key, $$results{$key});
 337         }
 338     }
 339
 340     my $primed_seq = Bio::Seq::PrimedSeq->new(
 341          -target_sequence => $self->{'seqobject'}->clone,
 342          -left_primer     => $left_seq,
 343          -right_primer    => $right_seq,
 344     );
 345
 346     # add data to the the PrimedSeq object that's not specific to the Primers
 347     for my $key (%$results) {
 348         next if ($key =~ /PRIMER_(LEFT|RIGHT)/i );
 349             $primed_seq->add_tag_value($key, $$results{$key});
 350     }
 351
 352     $self->{'next_to_return'}++;
 353     return $primed_seq;
 354 }
 355
 356
 357 =head2 primer_stream
 358
 359   Title   : primer_stream()
 360   Usage   : while (my $primed_seq  = $primer3->primer_stream()) {
 361   Function: Retrieve the primer/sequences one at a time
 362   Returns : Returns a Bio::Seq::PrimedSeq object, one at a time
 363   Args    : None
 364   Notes   : Deprecated, just a link to next_primer
 365
 366 =cut
 367
 368 sub primer_stream {
 369     my $self = shift;
 370     my $primedseq = $self->next_primer;
 371     return $primedseq;
 372 }
 373
 374
 375 =head2 _separate
 376
 377   Title   : _separate()
 378   Usage   : $self->_separate();
 379   Function: An internal function that groups the results by number
 380             (e.g. primer pair 1, etc)
 381   Returns : Nothing.
 382   Args    : None
 383   Notes   :
 384
 385 =cut
 386
 387 sub _separate {
 388     my $self = shift;
 389     my %results; # the results that we find
 390     my $maxlocation = -1; # the maximum number of primers returned
 391     foreach my $key (keys %{$self->{'results'}}) {
 392         next if (${$self->{'input_options'}}{$key}); # don't process it if it is an input key
 393
 394         my $location; # the number of the primer pair
 395         # names will have values like
 396         # PRIMER_RIGHT_SEQUENCE, PRIMER_RIGHT_2_SEQUENCE, PRIMER_PRODUCT_SIZE, and
 397         # PRIMER_PRODUCT_SIZE_3 hence we need to find and remove the number
 398         my $tempkey = $key;
 399         if ($tempkey =~ s/_(\d+)//) {
 400             $location = $1;
 401             if ($location > $maxlocation) {$maxlocation = $location}
 402         } elsif ( $tempkey =~ /PRIMER_(RIGHT|LEFT)_SEQUENCE/ ) {
 403             # first primers reported without a number, therefore set $location to 0
 404             $location = 0;
 405             if ($location > $maxlocation) {$maxlocation = $location}
 406         } else {
 407             $location = 0;
 408         }
 409         # we will hash the results by number, and then by name
 410         ${$results{$location}}{$tempkey}=${$self->{'results'}}{$key};
 411     }
 412     $self->{'results_by_number'} = \%results;
 413     $self->{'maximum_primers_returned'} = $maxlocation;
 414 }
 415
 416
 417 =head2 _set_variable
 418
 419   Title   : _set_variable()
 420   Usage   : $self->_set_variable('variable name', 'value');
 421   Function: An internal function that sets a variable
 422   Returns : Nothing.
 423   Args    : None
 424   Notes   : Used to set $self->{results} and $self->seqobject
 425
 426 =cut
 427
 428 sub _set_variable {
 429     my ($self, $name, $value) = @_;
 430     next unless ($name);
 431     $self->{$name} = $value;
 432 }
 433
 434 1;
 435
 436 __END__
 437