3 # This is the original copyright statement. I have relied on Chad's module
4 # extensively for this module.
6 # Copyright (c) 1997-2001 bioperl, Chad Matsalla. All Rights Reserved.
7 # This module is free software; you can redistribute it and/or
8 # modify it under the same terms as Perl itself.
10 # Copyright Chad Matsalla
12 # You may distribute this module under the same terms as perl itself
13 # POD documentation - main docs before the code
15 # But I have modified lots of it, so I guess I should add:
17 # Copyright (c) 2003 bioperl, Rob Edwards. All Rights Reserved.
18 # This module is free software; you can redistribute it and/or
19 # modify it under the same terms as Perl itself.
21 # Copyright Rob Edwards
23 # You may distribute this module under the same terms as perl itself
24 # POD documentation - main docs before the code
29 Bio::Seq::PrimedSeq - A representation of a sequence and two primers
30 flanking a target region
34 # The easiest way to use this is probably either, (i), get the
35 # output from Bio::Tools::Run::Primer3, Bio::Tools::Primer3, or
36 # Bio::Tools::PCRSimulation:
38 # For example, start with a fasta file
40 use Bio::Tools::Run::Primer3;
42 my $file = shift || die "need a file to read";
43 my $seqin = Bio::SeqIO->new(-file => $file);
44 my $seq = $seqin->next_seq;
46 # use primer3 to design some primers
47 my $primer3run = Bio::Tools::Run::Primer3->new(-seq => $seq);
48 $primer3run -> run; # run it with the default parameters
50 # create a file to write the results to
51 my $seqout = Bio::SeqIO->new(-file => ">primed_sequence.gbk",
52 -format => 'genbank');
54 # now just get all the results and write them out.
55 while (my $results = $primer3run->next_primer) {
56 $seqout->write_seq($results->annotated_seq);
59 # Or, (ii), to create a genbank file for a sequence and its cognate
63 use Bio::Seq::PrimedSeq;
65 # have a sequence file ($file) with the template, and two primers
66 # that match it, in fasta format
68 my $file = shift || die "$0 <file>";
69 my $seqin = Bio::SeqIO->new(-file => $file);
71 # read three sequences
72 my ($template, $leftprimer, $rightprimer) =
73 ($seqin->next_seq, $seqin->next_seq, $seqin->next_seq);
74 # set up the primed sequence object
75 my $primedseq = Bio::Seq::PrimedSeq->new(-seq => $template,
76 -left_primer => $leftprimer,
77 -right_primer => $rightprimer);
78 # open a file for output
79 my $seqout = Bio::SeqIO->new(-file => ">primed_sequence.gbk",
80 -format => 'genbank');
81 # print the sequence out
82 $seqout->write_seq($primedseq->annotated_sequence);
84 # This should output a genbank file with the two primers labeled.
88 This module is a slightly glorified capsule containing a primed sequence.
89 It was created to address the fact that a primer is more than a seqfeature
90 and there need to be ways to represent the primer-sequence complex and
91 the behaviors and attributes that are associated with the complex.
93 The primers are represented as Bio::SeqFeature::Primer objects, and should
94 be instantiated first.
96 A simple way to create a PrimedSeq object is as follows:
98 my $primedseq = Bio::Seq::PrimedSeq->new(
99 -seq => $seq, # Bio::Seq object,
100 -left_primer => $left, # Bio::SeqFeature::Primer object,
101 -right_primer => $right # Bio::SeqFeature::Primer object,
104 From the PrimedSeq object you should be able to retrieve
105 information about melting temperatures and what not on each of the primers
108 This is based on the PrimedSeq.pm module started by Chad Matsalla, with
109 additions/improvements by Rob Edwards.
113 User feedback is an integral part of the evolution of this and other
114 Bioperl modules. Send your comments and suggestions preferably to one
115 of the Bioperl mailing lists. Your participation is much appreciated.
117 bioperl-l@bioperl.org - General discussion
118 http://bioperl.org/wiki/Mailing_lists - About the mailing lists
122 Please direct usage questions or support issues to the mailing list:
124 L<bioperl-l@bioperl.org>
126 rather than to the module maintainer directly. Many experienced and
127 reponsive experts will be able look at the problem and quickly
128 address it. Please include a thorough description of the problem
129 with code and data examples if at all possible.
131 =head2 Reporting Bugs
133 Report bugs to the Bioperl bug tracking system to help us keep track
134 the bugs and their resolution. Bug reports can be submitted via the
137 http://bugzilla.open-bio.org/
141 Rob Edwards, redwards@utmem.edu
143 Based on a module written by Chad Matsalla, bioinformatics1@dieselwurks.com
147 The rest of the documentation details each of the object
148 methods. Internal methods are usually preceded with a _
153 # Let the code begin...
155 package Bio
::Seq
::PrimedSeq
;
158 use Bio
::SeqFeature
::Primer
;
160 use vars qw
($AUTOLOAD @RES %OK_FIELD $ID);
162 use base
qw(Bio::Root::Root Bio::SeqFeature::Generic);
165 @RES = qw(); # nothing here yet, not sure what we want!
166 foreach my $attr (@RES) {$OK_FIELD{$attr}++}
169 $ID = 'Bio::Tools::Analysis::Nucleotide::PrimedSeq';
173 my $attr = $AUTOLOAD;
175 $self->throw("Unallowed parameter: $attr !") unless $OK_FIELD{$attr};
176 $self->{$attr} = shift if @_;
177 return $self->{$attr};
184 Usage : $primed_sequence = Bio::SeqFeature::Primer->new(
186 -left_primer => $left_primer,
187 -right_primer => $right_primer);
188 Function: A constructor for an object representing a primed sequence
189 Returns : A Bio::Seq::PrimedSeq object
190 Args : -seq => a Bio::Seq object (required)
191 -left_primer => a Bio::SeqFeature::Primer object (required)
192 -right_primer => a Bio::SeqFeature::Primer object (required)
194 Many other parameters can be included including all of the output
195 parameters from the primer3 program. At the moment most of these
196 parameters will not do anything.
202 # note, I have cleaned up a lot of the script that Chad had written here,
203 # and I have removed the part where he removed the - before the tags.
206 my($class,%args) = @_;
207 my $self = $class->SUPER::new
(%args);
208 # these are the absolute minimum components required to make
211 foreach my $key (keys %args) {
212 if ($key =~ /^-seq/i) {
213 $self->{target_sequence
} = $args{$key};
217 ($okey = $key) =~ s/^-//;
218 if (($okey eq "left_primer" || $okey eq "right_primer") &&
219 ref($args{$key}) && $args{$key}->isa('Bio::SeqI') ) {
220 # We have been given a Bio::Seq object,
221 # make it a Bio::SeqFeature::Primer object
222 $self->{$okey} = Bio
::SeqFeature
::Primer
->new(-seq
=> $args{$key});
223 push @
{$self->{'arguments'}},$okey;
227 $self->{$okey} = $args{$key};
228 push @
{$self->{'arguments'}},$okey;
231 # and now the insurance - make sure that things are ok
232 if (!$self->{target_sequence
} || !$self->{left_primer
} ||
233 !$self->{right_primer
} ) {
234 $self->throw("You must provide a -seq, -left_primer, and -right_primer to create this object.");
237 if (! ref($self->{target_sequence
}) ||
238 ! $self->{target_sequence
}->isa('Bio::SeqI') ) {
239 $self->throw("The target_sequence must be a Bio::Seq to create this object.");
241 if (! ref($self->{left_primer
}) ||
242 ! $self->{left_primer
}->isa("Bio::SeqFeature::Primer") ||
243 ! ref($self->{right_primer
}) ||
244 ! $self->{right_primer
}->isa("Bio::SeqFeature::Primer")) {
245 $self->throw("You must provide a left_primer and right_primer, both as Bio::SeqFeature::Primer to create this object.");
248 # now we have the sequences, lets find out where they are
249 $self->_place_seqs();
256 Title : get_primer();
257 Usage : $primer = $primedseq->get_primer(l, left, left_primer,
258 -left_primer) to return the left primer or
259 $primer = $primedseq->get_primer(r, right, right_primer,
260 -right_primer) to return the right primer or
261 $primer = $primedseq->get_primer(b, both, both_primers,
263 to return the left primer, right primer array
264 Function: A getter for the left primer in thie PrimedSeq object.
265 Returns : A Bio::SeqFeature::Primer object
266 Args : Either of (l, left, left_primer, -left_primer) to get left
268 Either of (r, right, right_primer, -right_primer) to get
270 Either of (b, both, both_primers, -both_primers) to get
272 Note that this is plural. [default]
277 my ($self, $arg) = @_;
278 if (! defined $arg ) {
279 return ($self->{'left_primer'}, $self->{'right_primer'});
280 } elsif( $arg =~ /^l/ || $arg =~ /^-l/) {
281 # what a cheat, I couldn't be bothered to write all those or statements!
282 # Hah, now you can write leprechaun to get the left primer.
283 return $self->{'left_primer'}
285 elsif ($arg =~ /^r/ || $arg =~ /^-r/) {return $self->{'right_primer'}}
286 elsif ($arg =~ /^b/ || $arg =~ /^-b/) {return ($self->{'left_primer'}, $self->{'right_primer'})}
291 =head2 annotated_sequence
293 Title : annotated_sequence
294 Usage : $annotated_sequence_object = $primedseq->annotated_sequence()
295 Function: Get an annotated sequence object containg the left and right
297 Returns : An annotated sequence object or 0 if not defined.
299 Note : Use this method to return a sequence object that you can write
300 out (e.g. in GenBank format). See the example above.
304 sub annotated_sequence
{
306 if (exists $self->{annotated_sequence
}) {return $self->{annotated_sequence
}}
313 Usage : my $amplicon = $primedseq->amplicon()
314 Function: Retrieve the amplicon as a sequence object
315 Returns : A seq object. To get the sequence use $amplicon->seq
322 my ($self,@args) = @_;
323 my $id = $self->{'-seq'}->{'id'};
324 unless ($id) {$id=""}
325 # this just prevents a warning when $self->{'-seq'}->{'id'} is not defined
326 $id = "Amplicon from ".$id;
328 my $seqobj=Bio
::Seq
->new(-id
=>$id, seq
=>$self->{'amplicon_sequence'});
336 Usage : my $seqobj = $primedseq->seq()
337 Function: Retrieve the target sequence as a sequence object
338 Returns : A seq object. To get the sequence use $seqobj->seq
346 return $self->{target_sequence
};
352 Usage : $self->_place_seqs()
353 Function: An internal method to place the primers on the sequence and
354 set up the ranges of the sequences
357 Note : Internal use only
364 # we are going to pull out the target sequence, and then the primer sequences
365 my $target_sequence = $self->{'target_sequence'}->seq();
368 my $left_seq = $self->{'left_primer'}->seq()->seq();
370 my $rprc = $self->{'right_primer'}->seq()->revcom();
372 my $right_seq=$rprc->seq();
374 # now just change the case, because we keep getting screwed on this
375 $target_sequence=uc($target_sequence);
376 $left_seq=uc($left_seq);
377 $right_seq=uc($right_seq);
379 unless ($target_sequence =~ /(.*)$left_seq(.*)$right_seq(.*)/) {
380 unless ($target_sequence =~ /$left_seq/) {$self->throw("Can't place left sequence on target!")}
381 unless ($target_sequence =~ /$right_seq/) {$self->throw("Can't place right sequence on target!")}
384 my ($before, $middle, $after) = ($1, $2, $3); # note didn't use $`, $', and $& because they are bad. Just use length instead.
386 # cool now we can figure out lengths and what not.
387 # we'll figure out the position and compare it to known positions (e.g. from primer3)
389 my $left_location = length($before). ",". length($left_seq);
390 my $right_location = (length($target_sequence)-length($after)-1).",".length($right_seq);
391 my $amplicon_size = length($left_seq)+length($middle)+length($right_seq);
393 if (exists $self->{'left_primer'}->{'PRIMER_LEFT'}) {
394 # this is the left primer from primer3 input
395 # just check to make sure it is right
396 unless ($self->{'left_primer'}->{'PRIMER_LEFT'} eq $left_location) {
397 $self->warn("Note got |".$self->{'left_primer'}->{'PRIMER_LEFT'}."| from primer3 and |$left_location| for the left primer. You should email redwards\@utmem.edu about this.");
401 $self->{'left_primer'}->{'PRIMER_LEFT'}=$left_location;
404 if (exists $self->{'right_primer'}->{'PRIMER_RIGHT'}) {
405 # this is the right primer from primer3 input
406 # just check to make sure it is right
407 unless ($self->{'right_primer'}->{'PRIMER_RIGHT'} eq $right_location) {
408 $self->warn("Note got |".$self->{'right_primer'}->{'PRIMER_RIGHT'}."| from primer3 and |$right_location| for the right primer. You should email redwards\@utmem.edu about this.");
412 $self->{'right_primer'}->{'PRIMER_RIGHT'}=$right_location;
415 if (exists $self->{'PRIMER_PRODUCT_SIZE'}) {
416 # this is the product size from primer3 input
417 # just check to make sure it is right
418 unless ($self->{'PRIMER_PRODUCT_SIZE'} eq $amplicon_size) {
419 $self->warn("Note got |".$self->{'PRIMER_PRODUCT_SIZE'}."| from primer3 and |$amplicon_size| for the size. You should email redwards\@utmem.edu about this.");
423 $self->{'PRIMER_PRODUCT_SIZE'} = $amplicon_size;
426 $self->{'amplicon_sequence'}= lc($left_seq).uc($middle).lc($right_seq); # I put this in a different case, but I think the seqobj may revert this
428 $self->_set_seqfeature;
431 =head2 _set_seqfeature
433 Title : _set_seqfeature
434 Usage : $self->_set_seqfeature()
435 Function: An internal method to create Bio::SeqFeature::Generic objects
439 Note : Internal use only. Should only call this once left and right
440 primers have been placed on the sequence. This will then set
441 them as sequence features so hopefully we can get a nice output
447 sub _set_seqfeature
{
449 unless ($self->{'left_primer'}->{'PRIMER_LEFT'} &&
450 $self->{'right_primer'}->{'PRIMER_RIGHT'}) {
451 $self->warn("hmmm. Haven't placed primers, but trying to make annotated sequence");
454 my ($start, $length) = split /,/, $self->{'left_primer'}->{'PRIMER_LEFT'};
455 my $tm=$self->{'left_primer'}->{'PRIMER_LEFT_TM'} || $self->{'left_primer'}->Tm || 0;
457 my $seqfeatureL=Bio
::SeqFeature
::Generic
->new(
458 -start
=> $start+1, -end
=> $start+$length, -strand
=> 1,
459 -primary_tag
=> 'left_primer', -source
=> 'primer3',
460 -tag
=> {new
=> 1, author
=> 'Bio::Seq::PrimedSeq', Tm
=> $tm}
463 ($start, $length) = split /,/, $self->{'right_primer'}->{'PRIMER_RIGHT'};
464 $tm=$self->{'right_primer'}->{'PRIMER_RIGHT_TM'} || $self->{'right_primer'}->Tm || 0;
466 my $seqfeatureR=Bio
::SeqFeature
::Generic
->new(
467 -start
=> $start-$length+2, -end
=> $start+1, -strand
=> -1,
468 -primary_tag
=> 'right_primer', -source
=> 'primer3',
469 -tag
=> {new
=> 1, author
=> 'Bio::Seq::PrimedSeq', Tm
=> $tm}
472 # now add the sequences to a annotated sequence
473 $self->{annotated_sequence
} = $self->{target_sequence
};
474 $self->{annotated_sequence
}->add_SeqFeature($seqfeatureL);
475 $self->{annotated_sequence
}->add_SeqFeature($seqfeatureR);