lib/Bio/SeqFeature/Primer.pm

   1 #
   2 # BioPerl module for Bio::SeqFeature::Primer
   3 #
   4 # This is the original copyright statement. I have relied on Chad's module
   5 # extensively for this module.
   6 #
   7 # Copyright (c) 1997-2001 bioperl, Chad Matsalla. All Rights Reserved.
   8 #           This module is free software; you can redistribute it and/or
   9 #           modify it under the same terms as Perl itself.
  10 #
  11 # Copyright Chad Matsalla
  12 #
  13 # You may distribute this module under the same terms as perl itself
  14 # POD documentation - main docs before the code
  15 #
  16 # But I have modified lots of it, so I guess I should add:
  17 #
  18 # Copyright (c) 2003 bioperl, Rob Edwards. All Rights Reserved.
  19 #           This module is free software; you can redistribute it and/or
  20 #           modify it under the same terms as Perl itself.
  21 #
  22 # Copyright Rob Edwards
  23 #
  24 # You may distribute this module under the same terms as perl itself
  25 # POD documentation - main docs before the code
  26
  27 =head1 NAME
  28
  29 Bio::SeqFeature::Primer - Primer Generic SeqFeature
  30
  31 =head1 SYNOPSIS
  32
  33   use Bio::SeqFeature::Primer;
  34
  35   # Primer object with explicitly-defined sequence object or sequence string
  36   my $primer = Bio::SeqFeature::Primer->new( -seq => 'ACGTAGCT' );
  37   $primer->display_name('test_id');
  38   print "These are the details of the primer:\n".
  39         "Name:     ".$primer->display_name."\n".
  40         "Tag:      ".$primer->primary_tag."\n".   # always 'Primer'
  41         "Sequence: ".$primer->seq->seq."\n".
  42         "Tm:       ".$primer->Tm."\n\n";            # melting temperature
  43
  44   # Primer object with implicit sequence object
  45   # It is a lighter approach for when the primer location on a template is known
  46   use Bio::Seq;
  47   my $template = Bio::Seq->new( -seq => 'ACGTAGCTCTTTTCATTCTGACTGCAACG' );
  48   $primer   = Bio::SeqFeature::Primer->new( -start => 1, -end =>5, -strand => 1 );
  49   $template->add_SeqFeature($primer);
  50   print "Primer sequence is: ".$primer->seq->seq."\n";
  51   # Primer sequence is 'ACGTA'
  52
  53 =head1 DESCRIPTION
  54
  55 This module handles PCR primer sequences. The L<Bio::SeqFeature::Primer> object
  56 is a L<Bio::SeqFeature::Subseq> object that can additionally contain a primer
  57 sequence and its coordinates on a template sequence. The primary_tag() for this
  58 object is 'Primer'. A method is provided to calculate the melting temperature Tm
  59 of the primer. L<Bio::SeqFeature::Primer> objects are useful to build
  60 L<Bio::Seq::PrimedSeq> amplicon objects such as the ones returned by
  61 L<Bio::Tools::Primer3>.
  62
  63 =head1 FEEDBACK
  64
  65 =head2 Mailing Lists
  66
  67 User feedback is an integral part of the evolution of this and other
  68 Bioperl modules. Send your comments and suggestions preferably to one
  69 of the Bioperl mailing lists.  Your participation is much appreciated.
  70
  71   bioperl-l@bioperl.org                  - General discussion
  72   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  73
  74 =head2 Support
  75
  76 Please direct usage questions or support issues to the mailing list:
  77
  78 I<bioperl-l@bioperl.org>
  79
  80 rather than to the module maintainer directly. Many experienced and
  81 reponsive experts will be able look at the problem and quickly
  82 address it. Please include a thorough description of the problem
  83 with code and data examples if at all possible.
  84
  85 =head2 Reporting Bugs
  86
  87 Report bugs to the Bioperl bug tracking system to help us keep track
  88 the bugs and their resolution.  Bug reports can be submitted via the
  89 web:
  90
  91   https://github.com/bioperl/bioperl-live/issues
  92
  93 =head1 AUTHOR
  94
  95 Rob Edwards, redwards@utmem.edu
  96
  97 The original concept and much of the code was written by
  98 Chad Matsalla, bioinformatics1@dieselwurks.com
  99
 100 =head1 APPENDIX
 101
 102 The rest of the documentation details each of the object
 103 methods. Internal methods are usually preceded with a _
 104
 105 =cut
 106
 107
 108 package Bio::SeqFeature::Primer;
 109
 110 use strict;
 111 use Bio::PrimarySeq;
 112 use Bio::Tools::SeqStats;
 113
 114 use base qw(Bio::SeqFeature::SubSeq);
 115
 116
 117 =head2 new()
 118
 119  Title   : new()
 120  Usage   : my $primer = Bio::SeqFeature::Primer( -seq => $seq_object );
 121  Function: Instantiate a new Bio::SeqFeature::Primer object
 122  Returns : A Bio::SeqFeature::Primer object
 123  Args    : -seq , a sequence object or a sequence string (optional)
 124            -id  , the ID to give to the primer sequence, not feature (optional)
 125
 126 =cut
 127
 128 sub new {
 129     my ($class, %args) = @_;
 130
 131     # Legacy stuff
 132     my $sequence = delete $args{-sequence};
 133     if ($sequence) {
 134         Bio::Root::Root->deprecated(
 135             -message => 'Creating a Bio::SeqFeature::Primer with -sequence is deprecated. Use -seq instead.',
 136             -warn_version  => '1.6',
 137             -throw_version => '1.8',
 138         );
 139         $args{-seq} = $sequence;
 140     }
 141
 142     # Initialize Primer object
 143     my $self = $class->SUPER::new(%args);
 144     my ($id) = $self->_rearrange([qw(ID)], %args);
 145     $id && $self->seq->id($id);
 146     $self->primary_tag('Primer');
 147     return $self;
 148 }
 149
 150
 151 # Bypass B::SF::Generic's location() when a string is passed (for compatibility)
 152
 153 sub location {
 154     my ($self, $location) = @_;
 155     if ($location) {
 156         if ( not ref $location ) {
 157             # Use location as a string for backward compatibility
 158             Bio::Root::Root->deprecated(
 159                 -message => 'Passing a string to location() is deprecated. Pass a Bio::Location::Simple object or use start() and end() instead.',
 160                 -warn_version  => '1.6',
 161                 -throw_version => '1.8',
 162             );
 163             $self->{'_location'} = $location;
 164         } else {
 165             $self->SUPER::location($location);
 166         }
 167     }
 168     return $self->SUPER::location;
 169 }
 170
 171
 172 =head2 Tm()
 173
 174  Title   : Tm()
 175  Usage   : my $tm = $primer->Tm(-salt => 0.05, -oligo => 0.0000001);
 176  Function: Calculate the Tm (melting temperature) of the primer
 177  Returns : A scalar containing the Tm.
 178  Args    : -salt  : set the Na+ concentration on which to base the calculation
 179                     (default=0.05 molar).
 180          : -oligo : set the oligo concentration on which to base the
 181                     calculation (default=0.00000025 molar).
 182  Notes   : Calculation of Tm as per Allawi et. al Biochemistry 1997
 183            36:10581-10594. Also see documentation at
 184            http://www.idtdna.com/Scitools/Scitools.aspx as they use this
 185            formula and have a couple nice help pages. These Tm values will be
 186            about are about 0.5-3 degrees off from those of the idtdna web tool.
 187            I don't know why.
 188
 189            This was suggested by Barry Moore (thanks!). See the discussion on
 190            the bioperl-l with the subject "Bio::SeqFeature::Primer Calculating
 191            the PrimerTM"
 192
 193 =cut
 194
 195 sub Tm {
 196     my ($self, %args) = @_;
 197     my $salt_conc = 0.05; # salt concentration (molar units)
 198     my $oligo_conc = 0.00000025; # oligo concentration (molar units)
 199     if ($args{'-salt'}) {
 200         # Accept object defined salt concentration
 201         $salt_conc = $args{'-salt'};
 202     }
 203     if ($args{'-oligo'}) {
 204         # Accept object defined oligo concentration
 205         $oligo_conc = $args{'-oligo'};
 206     }
 207     my $seqobj = $self->seq();
 208     my $length = $seqobj->length();
 209     my $sequence = uc $seqobj->seq();
 210     my @dinucleotides;
 211     my $enthalpy;
 212     my $entropy;
 213     # Break sequence string into an array of all possible dinucleotides
 214     while ($sequence =~ /(.)(?=(.))/g) {
 215         push @dinucleotides, $1.$2;
 216     }
 217     # Build a hash with the thermodynamic values
 218     my %thermo_values = ('AA' => {'enthalpy' => -7.9,
 219                                   'entropy'  => -22.2},
 220                          'AC' => {'enthalpy' => -8.4,
 221                                   'entropy'  => -22.4},
 222                          'AG' => {'enthalpy' => -7.8,
 223                                   'entropy'  => -21},
 224                          'AT' => {'enthalpy' => -7.2,
 225                                   'entropy'  => -20.4},
 226                          'CA' => {'enthalpy' => -8.5,
 227                                   'entropy'  => -22.7},
 228                          'CC' => {'enthalpy' => -8,
 229                                   'entropy'  => -19.9},
 230                          'CG' => {'enthalpy' => -10.6,
 231                                   'entropy'  => -27.2},
 232                          'CT' => {'enthalpy' => -7.8,
 233                                   'entropy'  => -21},
 234                          'GA' => {'enthalpy' => -8.2,
 235                                   'entropy'  => -22.2},
 236                          'GC' => {'enthalpy' => -9.8,
 237                                   'entropy'  => -24.4},
 238                          'GG' => {'enthalpy' => -8,
 239                                   'entropy'  => -19.9},
 240                          'GT' => {'enthalpy' => -8.4,
 241                                   'entropy'  => -22.4},
 242                          'TA' => {'enthalpy' => -7.2,
 243                                   'entropy'  => -21.3},
 244                          'TC' => {'enthalpy' => -8.2,
 245                                   'entropy'  => -22.2},
 246                          'TG' => {'enthalpy' => -8.5,
 247                                   'entropy'  => -22.7},
 248                          'TT' => {'enthalpy' => -7.9,
 249                                   'entropy'  => -22.2},
 250                          'A' =>  {'enthalpy' => 2.3,
 251                                   'entropy'  => 4.1},
 252                          'C' =>  {'enthalpy' => 0.1,
 253                                   'entropy'  => -2.8},
 254                          'G' =>  {'enthalpy' => 0.1,
 255                                   'entropy'  => -2.8},
 256                          'T' =>  {'enthalpy' => 2.3,
 257                                   'entropy'  => 4.1}
 258                         );
 259     # Loop through dinucleotides and calculate cumulative enthalpy and entropy values
 260     for (@dinucleotides) {
 261         $enthalpy += $thermo_values{$_}{enthalpy};
 262         $entropy += $thermo_values{$_}{entropy};
 263     }
 264     # Account for initiation parameters
 265     $enthalpy += $thermo_values{substr($sequence,  0, 1)}{enthalpy};
 266     $entropy  += $thermo_values{substr($sequence,  0, 1)}{entropy};
 267     $enthalpy += $thermo_values{substr($sequence, -1, 1)}{enthalpy};
 268     $entropy  += $thermo_values{substr($sequence, -1, 1)}{entropy};
 269     # Symmetry correction
 270     $entropy -= 1.4;
 271     my $r = 1.987; # molar gas constant
 272     my $tm = $enthalpy * 1000 / ($entropy + ($r * log($oligo_conc))) - 273.15 + (12* (log($salt_conc)/log(10)));
 273
 274     return $tm;
 275  }
 276
 277 =head2 Tm_estimate
 278
 279  Title   : Tm_estimate
 280  Usage   : my $tm = $primer->Tm_estimate(-salt => 0.05);
 281  Function: Estimate the Tm (melting temperature) of the primer
 282  Returns : A scalar containing the Tm.
 283  Args    : -salt set the Na+ concentration on which to base the calculation.
 284  Notes   : This is only an estimate of the Tm that is kept in for comparative
 285            reasons. You should probably use Tm instead!
 286
 287            This Tm calculations are taken from the Primer3 docs: They are
 288            based on Bolton and McCarthy, PNAS 84:1390 (1962)
 289            as presented in Sambrook, Fritsch and Maniatis,
 290            Molecular Cloning, p 11.46 (1989, CSHL Press).
 291
 292            Tm = 81.5 + 16.6(log10([Na+])) + .41*(%GC) - 600/length
 293
 294            where [Na+] is the molar sodium concentration, %GC is the
 295            %G+C of the sequence, and length is the length of the sequence.
 296
 297            However.... I can never get this calculation to give me the same result
 298            as primer3 does. Don't ask why, I never figured it out. But I did
 299            want to include a Tm calculation here because I use these modules for
 300            other things besides reading primer3 output.
 301
 302            The primer3 calculation is saved as 'PRIMER_LEFT_TM' or 'PRIMER_RIGHT_TM'
 303            and this calculation is saved as $primer->Tm so you can get both and
 304            average them!
 305
 306 =cut
 307
 308 sub Tm_estimate {
 309
 310     # This should probably be put into seqstats as it is more generic, but what the heck.
 311
 312     my ($self, %args) = @_;
 313     my $salt = 0.2;
 314     if ($args{'-salt'}) {
 315         $salt = $args{'-salt'}
 316     };
 317     my $seqobj = $self->seq();
 318     my $length = $seqobj->length();
 319     my $seqdata = Bio::Tools::SeqStats->count_monomers($seqobj);
 320     my $gc=$$seqdata{'G'} + $$seqdata{'C'};
 321     my $percent_gc = ($gc/$length)*100;
 322
 323     my $tm = 81.5+(16.6*(log($salt)/log(10)))+(0.41*$percent_gc) - (600/$length);
 324
 325     return $tm;
 326 }
 327
 328 =head2 primary_tag, source_tag, location, start, end, strand...
 329
 330 The documentation of L<Bio::SeqFeature::Generic> describes all the methods that
 331 L<Bio::SeqFeature::Primer> object inherit.
 332
 333 =cut
 334
 335 1;