comparison variant_effect_predictor/Bio/Tools/Primer3.pm @ 0:1f6dce3d34e0

Uploaded

0:1f6dce3d34e0

#
# 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.

  [e-mail]          - 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:

  [e-mail]
  http://bugzilla.bioperl.org/


=head1 AUTHOR - Chad Matsalla

[e-mail]

=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