Raw content of Bio::Tools::Primer3 # # Copyright (c) 1997-2001 bioperl, Chad Matsalla. All Rights Reserved. # This module is free software; you can redistribute it and/or # modify it under the same terms as Perl itself. # # Copyright Chad Matsalla # # You may distribute this module under the same terms as perl itself # POD documentation - main docs before the code =head1 NAME Bio::Tools::Primer3 - Create input for and work with the output from the program primer3 =head1 SYNOPSIS Chad will put synopses here by the end of the second week of october, 2002. =head1 DESCRIPTION Bio::Tools::Primer3 creates the input files needed to design primers using primer3 and provides mechanisms to access data in the primer3 output files. =head1 FEEDBACK =head2 Mailing Lists User feedback is an integral part of the evolution of this and other Bioperl modules. Send your comments and suggestions preferably to one of the Bioperl mailing lists. Your participation is much appreciated. bioperl-l@bioperl.org - General discussion http://www.bioperl.org/MailList.html - About the mailing lists =head2 Reporting Bugs Report bugs to the Bioperl bug tracking system to help us keep track the bugs and their resolution. Bug reports can be submitted via email or the web: bioperl-bugs@bio.perl.org http://bugzilla.bioperl.org/ =head1 AUTHOR - Chad Matsalla bioinformatics1@dieselwurks.com =head1 APPENDIX The rest of the documentation details each of the object methods. Internal methods are usually preceded with a _ =cut # Let the code begin... package Bio::Tools::Primer3; use vars qw(@ISA); use strict; use Bio::Seq; use Bio::SeqFeature::Primer; use Bio::Seq::PrimedSeq; use Bio::Seq::SeqFactory; use Bio::Root::Root; use Bio::Root::IO; use Dumpvalue; @ISA = qw(Bio::Root::Root Bio::Root::IO); # Chad likes to use this to debug large hashes. my $dumper = new Dumpvalue; # this was a bunch of the seqio things, now deprecated. delete it soon. # sub _initialize { # my($self,@args) = @_; # $self->SUPER::_initialize(@args); # if( ! defined $self->sequence_factory ) { # $self->sequence_factory(new Bio::Seq::SeqFactory # (-verbose => $self->verbose(), # -type => 'Bio::Seq')); # } # } =head2 new() Title : new() Usage : Function: Returns : Args : Notes : =cut sub new { my($class,@args) = @_; my $self = $class->SUPER::new(@args); my($filename) = $self->_rearrange([qw(FILE)],@args); if (!$filename) { print("Ahh grasshopper, you are planning to create a primer3 infile\n"); return $self; } $self->{filename} = $filename; # check to see that the file exists # I think that catfile should be used here. if (!-f $filename) { print("That file doesn't exist. Bah.\n"); } $self->_initialize_io( -file => $filename ); return $self; } =head2 null Title : Usage : Function: Returns : Args : Notes : =cut =head2 next_primer() Title : next_primer() Usage : $primer3 = $stream->next_primer() Function: returns the next primer in the stream Returns : Bio::Seq::PrimedSeq containing: - 2 Bio::SeqFeature::Primer representing the primers - 1 Bio::Seq representing the target sequence - 1 Bio::Seq representing the amplified region Args : NONE Notes : =cut sub next_primer { my $self = shift; my $fh = $self->_fh(); my ($line,%primer); # first, read in the next set of primers while ($line = $self->_readline()) { chomp ($line); last if ($line =~ /^=/); $line =~ m/(^.*)\=(.*$)/; $primer{$1} = $2; } # then, get the primers as SeqFeature::Primer objects my ($left,$right) = &_create_primer_features(\%primer); # then, create the sequence to place them on my $sequence = Bio::Seq->new(-seq => $primer{SEQUENCE}, -id => $primer{PRIMER_SEQUENCE_ID}); # print("Sequence is ".$primer{SEQUENCE}." and id is ".$primer{PRIMER_SEQUENCE_ID}."\n"); my $primedseq = new Bio::Seq::PrimedSeq( -target_sequence => $sequence, -left_primer => $left, -right_primer => $right, -primer_sequence_id => $primer{PRIMER_SEQUENCE_ID}, -primer_comment => $primer{PRIMER_COMMENT}, -target => $primer{TARGET}, -primer_product_size_range => $primer{PRIMER_PRODUCT_SIZE_RANGE}, -primer_file_flag => $primer{PRIMER_FILE_FLAG}, -primer_liberal_base => $primer{PRIMER_LIBERAL_BASE}, -primer_num_return => $primer{PRIMER_NUM_RETURN}, -primer_first_base_index => $primer{PRIMER_FIRST_BASE_INDEX}, -primer_explain_flag => $primer{PRIMER_EXPLAIN_FLAG}, -primer_pair_compl_any => $primer{PRIMER_PAIR_COMPL_ANY}, -primer_pair_compl_end => $primer{PRIMER_PAIR_COMPL_END}, -primer_product_size => $primer{PRIMER_PRODUCT_SIZE} ); return $primedseq; } =head2 _create_primer_features() Title : _create_primer_features() Usage : &_create_primer_features() Function: This is an internal method used by next_seq() to create the Bio::SeqFeature::Primer objects necessary to represent the primers themselves. Returns : An array of 2 Bio::SeqFeature::Primer objects. Args : None. Notes : This is an internal method. Do not call this method. =cut sub _create_primer_features { my $rdat = shift; my (%left,%right,$updir,$downdir,$var,$trunc); my @variables = qw( PRIMER_DIRECTION PRIMER_DIRECTION_END_STABILITY PRIMER_DIRECTION_EXPLAIN PRIMER_DIRECTION_GC_PERCENT PRIMER_DIRECTION_PENALTY PRIMER_DIRECTION_SELF_ANY PRIMER_DIRECTION_SELF_END PRIMER_DIRECTION_SEQUENCE PRIMER_DIRECTION_TM PRIMER_FIRST_BASE_INDEX ); # create the hash to pass into the creation routine # I do it this way because the primer3 outfile variables are exactly the same for each of # left and right. I create two hashes- one for the left and one for the right primer. foreach $updir (qw(LEFT RIGHT)) { my %dat; foreach (@variables) { ($var = $_) =~ s/DIRECTION/$updir/e; # should you truncate the name of each variable? # for example, should the value be: PRIMER_RIGHT_PENALTY or PENALTY? # i think it should be the second one if (/^PRIMER_DIRECTION$/) { $trunc = "PRIMER"; } elsif (/^PRIMER_FIRST_BASE_INDEX/) { $trunc = "FIRST_BASE_INDEX"; } else { ($trunc = $_) =~ s/PRIMER_DIRECTION_//; } $dat{"-$trunc"} = $rdat->{$var}; } if ($updir eq "LEFT") { %left = %dat; $left{-id} = $rdat->{PRIMER_SEQUENCE_ID}."-left"; } else { %right = %dat; $right{-id} = $rdat->{PRIMER_SEQUENCE_ID}."-right"; } } my $primer_left = new Bio::SeqFeature::Primer(%left); my $primer_right = new Bio::SeqFeature::Primer(%right); return($primer_left,$primer_right); } =head2 get_amplified_region() Title : get_amplified_region() Usage : $primer->get_amplified_region() Function: Returns a Bio::Seq object representing the sequence amplified Returns : (I think) A Bio::Seq object Args : None. Notes : This is not implemented at this time. Note to chad: implement this simple getter. Developer notes: There obviously isn't a way for a single primer to know about its amplified region unless it is paired with another primer. At this time these object will generally be created with another so I will put in this method. If there is no sequence null is returned. THIS DOES NOT BELONG HERE. Put this into something else. =cut sub get_amplified_region { my ($self) = @_; } # end get_amplified_region =head2 get_amplification_error() Title : get_amplification_error() Usage : Function: Returns : Args : Notes : Developer Notes: THIS DOES NOT BELONG HERE. Put this into something else. =cut sub get_amplification_error { my $primer = $_[1]; my $error = $Primer3::primers{$primer}{PRIMER_ERROR}; if ($error) { return $error; } else { return "Some error that primer3 didn't define.\n"; } } =head2 _set_target() Title : _set_target() Usage : &_set_target($self); Function: Returns : Args : Notes : Developer Notes: Really I have no idea why I put this in here. It can is referenced by new_deprecated and by run_primer3 =cut sub _set_target { my $self = shift; my ($sequence,$primer,$primer_left,$primer_right,$position_left,$position_right,$boggle); $boggle = 1; foreach $primer (sort keys %{$self->{primers}}) { $sequence = $self->{primers}{$primer}{SEQUENCE}; $primer_left = $self->{primers}{$primer}{PRIMER_LEFT}; $primer_right = $self->{primers}{$primer}{PRIMER_RIGHT}; if (!$primer_left) { $self->{primers}{$primer}{design_failed} = "1"; } else { $primer_left =~ m/(.*)\,(.*)/; $position_left = $1+$2-1; $primer_right =~ m/(.*)\,(.*)/; $position_right = $1-$2; $self->{primers}{$primer}{left} = $position_left; $self->{primers}{$primer}{right} = $position_right; $self->{primers}{$primer}{amplified} = substr($sequence,$position_left,$position_right-$position_left); } } } =head2 _read_file($self,$filename) Title : _read_file($self,$filename) Usage : Function: Returns : A scalar containing the contents of $filename Args : $self and the name of a file to parse. Notes : Developer notes: Honestly, I have no idea what this is for. =cut sub _read_file { # my ($self,$filename) = @_; # set this to keep track of things.... # $self->{outfilename} = $filename; # to make this better for bioperl, chad should really be using catfile and things. # # my $fh = new FileHandle; # open($fh,$filename) or die "I can't open the primer report ($filename) : $!\n"; # # _parse_report(); # # my %Primer3::primers; # my ($output,$line); # while ($line=<$fh>) { # # print("Adding $line\n"); # $output .= $line; # } # end while # # print("\$output is $output\n"); # return $output; } =head2 _parse_report() Title : _parse_report() Usage : &_parse_report($self,$filename); Function: Parse a primer3 outfile and place everything into an object under {primers} with PRIMER_SEQUENCE_ID being the name of the keys for the {primers} hash. Returns : Nothing. Args : $self and the name of a file to parse. Notes : =cut sub _parse_report { # old # my ($self,$filename) = @_; my ($self,$outputs) = @_; # print("\$self is $self, \$outputs are $outputs\n"); # print("in _parse_report, \$self is $self\n"); # set this to keep track of things.... my ($sequence_name,$line,$counter,$variable_name,$variable_value); my @output = split/\n/,$outputs; foreach $line (@output) { # print("Reading line $line\n"); next if ($line =~ /^\=/); if ($line =~ m/^PRIMER_SEQUENCE_ID/) { $line =~ m/(\S+)=(.*$)/; $variable_name = $1; $sequence_name = $2; $variable_value = $2; } else { $line =~ m/(\S+)=(.*$)/; $variable_name = $1; $variable_value = $2; } # print("$sequence_name\t$variable_name\t$variable_value\n"); $self->{primers}{$sequence_name}{$variable_name} = $variable_value; } # end while <> } # end parse_report =head2 _construct_empty() Title : _construct_empty() Usage : &_construct_empty($self); Function: Construct an empty object that will be used to construct a primer3 input "file" so that it can be run. Returns : Args : Notes : =cut sub _construct_empty { my $self = shift; $self->{inputs} = {}; return; } =head2 add_target(%stuff) Title : add_target(%stuff) Usage : $o_primer->add_target(%stuff); Function: Add an target to the infile constructor. Returns : Args : A hash. Looks something like this: $o_primer2->add_target( -PRIMER_SEQUENCE_ID => "sN11902", -PRIMER_COMMENT => "3831", -SEQUENCE => "some_sequence", -TARGET => "513,26", -PRIMER_PRODUCT_SIZE_RANGE => "100-500", -PRIMER_FILE_FLAG => "0", -PRIMER_LIBERAL_BASE => "1", -PRIMER_NUM_RETURN => "1", -PRIMER_FIRST_BASE_INDEX => "1", -PRIMER_EXPLAIN_FLAG => "1"); The add_target() method does not validate the things you put into this parameter hash. Read the docs for Primer3 to see which fields do what and how they should be used. Notes : To design primers, first create a new CSM::Primer3 object with the -construct_infile parameter. Then, add targets using this method (add_target()) with the target hash as above in the Args: section. Be careful. No validation will be done here. All of those parameters will be fed straight into primer3. Once you are done adding targets, invoke the function run_primer3(). Then retrieve the results using something like a loop around the array from get_primer_sequence_IDs(); =cut sub add_target { my ($self,%args) = @_; my ($currkey,$renamed,$sequence_id,$value); if (!$args{-PRIMER_SEQUENCE_ID}) { print("You cannot add an element to the primer3 infile without specifying the PRIMER_SEQUENCE_ID. Sorry.\n"); } else { $sequence_id = $args{-PRIMER_SEQUENCE_ID}; foreach $currkey (keys %args) { # print("\$currkey is $currkey\n"); next if ($currkey eq "-PRIMER_SEQUENCE_ID"); ($renamed = $currkey) =~ s/-//; # print("Adding $renamed to the hash under $sequence_id\n"); $value = $args{$currkey}; # print("\$value is $value\n"); if ($renamed eq "SEQUENCE") { $value =~ s/\n//g; } $self->{infile}{$sequence_id}{$renamed} = $value; } } } =head2 get_primer_sequence_IDs() Title : get_primer_sequence_IDs() Usage : $o_phred->get_primer_sequence_IDs(); Function: Return the primer sequence ID's. These normally correspond to the name of a sequence in a database but can be whatever was used when the primer3 infile was constructed. Returns : An array containing the names of the primer sequence ID's Args : None. Notes : This would be used as the basis for an iterator to loop around each primer that was designed. =cut sub get_primer_sequence_IDs { my $self = shift; return sort keys %{$self->{primers}}; } # end get keys =head2 dump_hash() Title : dump_hash() Usage : $o_primer->dump_hash(); Function: Dump out the CSM::Primer3 object. Returns : Nothing. Args : None. Notes : Used extensively in debugging. =cut sub dump_hash { my $self = shift; my $dumper = new Dumpvalue; $dumper->dumpValue($self); } # end dump_hash =head2 dump_infile_hash() Title : dump_infile_hash() Usage : $o_primer->dump_infile_hash(); Function: Dump out the contents of the infile hash. Returns : Nothing. Args : None. Notes : Used for debugging the construction of the infile. =cut sub dump_infile_hash { my $self = shift; my $dumper = new Dumpvalue; $dumper->dumpValue($self->{infile}); } 1; __END__ =head2 placeholder Title : This is a place holder so chad can cut and paste Usage : Function: Returns : Args : Notes : =cut =head1 SEE ALSO perl(1). =cut