# set up a single primer that can be used in a PCR reaction

 use Bio::SeqFeature::Primer;

 # initiate a primer with raw sequence
 my $primer=Bio::SeqFeature::Primer->new(-seq=>'CTTTTCATTCTGACTGCAACG');

 # get the primery tag for the primer # should return Primer
 my $tag=$primer->primary_tag;

 # get or set the location that the primer binds to the target at
 $primer->location(500);
 my $location=$primer->location(500);

 # get or set the 5' end of the primer homology, as the primer doesn't 
 # have to be the same as the target sequence
 $primer->start(2);
 my $start=$primer->start;

 # get or set the 3' end of the primer homology
 $primer->end(19);
 my $end = $primer->end;

 # get or set the strand of the primer. Strand should be 1, 0, or -1
 $primer->strand(-1);
 my $strand=$primer->strand;

 # get or set the id of the primer
 $primer->display_id('test_id');
 my $id=$primer->display_id;

 # get the tm of the primer. This is calculated for you by the software.
 # however, see the docs.
 my $tm = $primer->Tm;

 print "These are the details of the primer:\n\tTag:\t\t$tag\n\tLocation\t$location\n\tStart:\t\t$start\n";
 print "\tEnd:\t\t$end\n\tStrand:\t\t$strand\n\tID:\t\t$id\n\tTm:\t\t$tm\n";

 @RES=qw(); # nothing here yet, not sure what we want!

 foreach my $attr (@RES) {$OK_FIELD{$attr}++

 my $self = shift;
 my $attr = $AUTOLOAD;
 $attr =~ s/.*:://;
 $self->throw("Unallowed parameter: $attr !") unless $OK_FIELD{$attr};
 $self->{$attr} = shift if @_;
 return $self->{$attr};

  # I have changed some of Chad's code. I hope he doesn't mind. Mine is more stupid than his, but my simple mind gets it.
  # I also removed from of the generic.pm things, but we can put them back....

  my ($class, %args) = @_;  
  my $self = $class->SUPER::new(%args);
  # i am going to keep an array of the things that have been passed
  # into the object on construction. this will aid retrieval of these
  # things later
  foreach my $argument (keys %args) {
   if ($argument eq "-SEQUENCE" || $argument eq "-sequence" || $argument eq "-seq") {
    if (ref($args{$argument}) eq "Bio::Seq") {$self->{seq} = $args{$argument}}
    else {
     unless ($args{-id}) {$args{-id}="SeqFeature Primer object"}
     $self->{seq} = new Bio::Seq( -seq => $args{$argument}, -id => $args{-id});
    }
    $self->{$argument} = $self->{seq};
    push (@{$self->{arguments}}, "seq");
   }
   else {
    $self->{$argument} = $args{$argument};
    push (@{$self->{arguments}}, $argument); # note need to check the BioPerl way of doing this.
   }
  }

  # now error check and make sure that we at least got a sequence
  if (!$self->{seq}) {$self->throw("You must pass in a sequence to construct this object.")}

   # a bunch of things now need to be set for this SeqFeature
   # things like:
   # TARGET=513,26
   # PRIMER_FIRST_BASE_INDEX=1
   # PRIMER_LEFT=484,20

   # these can be added in, and we won't demand them, but provide a mechanism to check that they exist

   $self->Tm();
   return $self;

     my $self = shift;
     return $self->{seq};

     return "Primer";

     my ($self,$insource) = @_;
     if ($insource) { $self->{source} = $insource; }
     return $self->{source};

     my ($self, $location) = @_;
     if ($location) {$self->{location}=$location}
     if ($self->{location}) {return $self->{location}}
     else {return 0}

     my ($self,$start) = @_;
     if ($start) {$self->{start_position} = $start}
     if ($self->{start_position}) {return $self->{start_position}}
     else {return 0}

     my ($self,$end) = @_;
     if ($end) {$self->{end_position} = $end}
     if ($self->{end_position}) {return $self->{end_position}}
     else {return 0}

     my ($self, $strand) = @_;
     if ($strand) {
      unless ($strand == -1 || $strand == 0 ||$strand == 1) {$self->throw("Strand must be either 1, 0, or -1 not $strand")}
      $self->{strand}=$strand;
     }
     if ($self->{strand}) {return $self->{strand}}
     else {return 0}

     my ($self,$newid) = @_;
     if ($newid) {$self->seq()->display_id($newid)}
     return $self->seq()->display_id();

     my ($self, %args) = @_;
     my $salt_conc = 0.05; #salt concentration (molar units)
     my $oligo_conc = 0.00000025; #oligo concentration (molar units)
     if ($args{'-salt'}) {$salt_conc = $args{'-salt'}} #accept object defined salt concentration
     if ($args{'-oligo'}) {$oligo_conc = $args{'-oligo'}} #accept object defined oligo concentration
     my $seqobj = $self->seq();
     my $length = $seqobj->length();
     my $sequence = uc $seqobj->seq();
     my @dinucleotides;
     my $enthalpy;
     my $entropy;
     #Break sequence string into an array of all possible dinucleotides
     while ($sequence =~ /(.)(?=(.))/g) {
         push @dinucleotides, $1.$2;
     }
     #Build a hash with the thermodynamic values
     my %thermo_values = ('AA' => {'enthalpy' => -7.9,
                                   'entropy'  => -22.2},
                          'AC' => {'enthalpy' => -8.4,
                                   'entropy'  => -22.4},
                          'AG' => {'enthalpy' => -7.8,
                                   'entropy'  => -21},
                          'AT' => {'enthalpy' => -7.2,
                                   'entropy'  => -20.4},
                          'CA' => {'enthalpy' => -8.5,
                                   'entropy'  => -22.7},
                          'CC' => {'enthalpy' => -8,
                                   'entropy'  => -19.9},
                          'CG' => {'enthalpy' => -10.6,
                                   'entropy'  => -27.2},
                          'CT' => {'enthalpy' => -7.8,
                                   'entropy'  => -21},
                          'GA' => {'enthalpy' => -8.2,
                                   'entropy'  => -22.2},
                          'GC' => {'enthalpy' => -9.8,
                                   'entropy'  => -24.4},
                          'GG' => {'enthalpy' => -8,
                                   'entropy'  => -19.9},
                          'GT' => {'enthalpy' => -8.4,
                                   'entropy'  => -22.4},
                          'TA' => {'enthalpy' => -7.2,
                                   'entropy'  => -21.3},
                          'TC' => {'enthalpy' => -8.2,
                                   'entropy'  => -22.2},
                          'TG' => {'enthalpy' => -8.5,
                                   'entropy'  => -22.7},
                          'TT' => {'enthalpy' => -7.9,
                                   'entropy'  => -22.2},
                          'A' =>  {'enthalpy' => 2.3,
                                   'entropy'  => 4.1},
                          'C' =>  {'enthalpy' => 0.1,
                                   'entropy'  => -2.8},
                          'G' =>  {'enthalpy' => 0.1,
                                   'entropy'  => -2.8},
                          'T' =>  {'enthalpy' => 2.3,
                                   'entropy'  => 4.1}
                         );
     #Loop through dinucleotides and calculate cumulative enthalpy and entropy values
     for (@dinucleotides) {
        $enthalpy += $thermo_values{$_}{enthalpy};
        $entropy += $thermo_values{$_}{entropy};
     }
     #Account for initiation parameters
     $enthalpy += $thermo_values{substr($sequence, 0, 1)}{enthalpy};
     $entropy += $thermo_values{substr($sequence, 0, 1)}{entropy};
     $enthalpy += $thermo_values{substr($sequence, -1, 1)}{enthalpy};
     $entropy += $thermo_values{substr($sequence, -1, 1)}{entropy};
     #Symmetry correction
     $entropy -= 1.4;
     my $r = 1.987; #molar gas constant
     my $tm = ($enthalpy * 1000 / ($entropy + ($r * log($oligo_conc))) - 273.15 + (12* (log($salt_conc)/log(10))));
     $self->{'Tm'}=$tm;
     return $tm;

 # note I really think that this should be put into seqstats as it is more generic, but what the heck.

 my ($self, %args) = @_;
 my $salt=0.2;
 if ($args{'-salt'}) {$salt=$args{'-salt'}}
 my $seqobj=$self->seq();
 my $length=$seqobj->length();
 my $seqdata = Bio::Tools::SeqStats->count_monomers($seqobj);
 my $gc=$$seqdata{'G'} + $$seqdata{'C'};
 my $percent_gc=($gc/$length)*100;


 my $tm= 81.5+(16.6*(log($salt)/log(10)))+(0.41*$percent_gc) - (600/$length);

 # and now error check compared to primer3
 # note that this NEVER gives me the same values, so I am ignoring it
 # you can get these out separately anyway

# if ($self->{'PRIMER_LEFT_TM'}) {
#  unless ($self->{'PRIMER_LEFT_TM'} == $tm) {
#   $self->warn("Calculated $tm for Left primer but received ".$self->{'PRIMER_LEFT_TM'}." from primer3\n");
#  }
# }
# elsif ($self->{'PRIMER_RIGHT_TM'}) {
#  unless ($self->{'PRIMER_RIGHT_TM'} == $tm) {
#    $self->warn("Calculated $tm for Right primer but received ".$self->{'PRIMER_RIGHT_TM'}." from primer3\n");
#  }
# }

 $self->{'Tm'}=$tm;
 return $tm;

  bioperl-l@bioperl.org                  - General discussion
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists

  http://bugzilla.open-bio.org/

	The rest of the documentation details each of the object
	methods. Internal methods are usually preceded with a _