~ubuntu-branches/ubuntu/raring/bioperl/raring

« back to all changes in this revision

Viewing changes to Bio/Seq/PrimedSeq.pm

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2008-03-18 14:44:57 UTC
  • mfrom: (4 hardy)
  • mto: This revision was merged to the branch mainline in revision 6.
  • Revision ID: james.westby@ubuntu.com-20080318144457-1jjoztrvqwf0gruk
* debian/control:
  - Removed MIA Matt Hope (dopey) from the Uploaders field.
    Thank you for your work, Matt. I hope you are doing well.
  - Downgraded some recommended package to the 'Suggests' priority,
    according to the following discussion on Upstream's mail list.
    http://bioperl.org/pipermail/bioperl-l/2008-March/027379.html
    (Closes: #448890)
* debian/copyright converted to machine-readable format.

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
# $Id: PrimedSeq.pm,v 1.23.4.1 2006/10/02 23:10:27 sendu Exp $
1
2
#
2
3
# This is the original copyright statement. I have relied on Chad's module
3
4
# extensively for this module.
25
26
 
26
27
=head1 NAME
27
28
 
28
 
Bio::Seq::PrimedSeq - A representation of a sequence and two primers flanking a
29
 
  target region for amplification
 
29
Bio::Seq::PrimedSeq - A representation of a sequence and two primers 
 
30
flanking a target region
30
31
 
31
32
=head1 SYNOPSIS
32
 
  # The easiest way to use this is probably as one of the following:
33
 
  # (i) to get the output from Bio::Tools::Run::Primer3, Bio::Tools::Primer3,
34
 
  # or Bio::Tools::PCRSimulation
35
 
 
36
 
  #    For example, start with a fasta file
37
 
 
 
33
 
 
34
The easiest way to use this is probably either, (i), get the output 
 
35
from Bio::Tools::Run::Primer3, Bio::Tools::Primer3, or 
 
36
Bio::Tools::PCRSimulation:
 
37
 
 
38
      # For example, start with a fasta file
38
39
      use Bio::SeqIO;
39
40
      use Bio::Tools::Run::Primer3;
40
41
 
41
42
      my $file = shift || die "need a file to read";
42
 
      my $seqin = Bio::SeqIO->new(-file=>$file);
 
43
      my $seqin = Bio::SeqIO->new(-file => $file);
43
44
      my $seq = $seqin->next_seq;
 
45
 
44
46
      # use primer3 to design some primers
45
 
      my $primer3run = Bio::Tools::Run::Primer3->new(-seq=>$seq);
46
 
      $primer3run -> run; # we'll just run it with the default parameters
 
47
      my $primer3run = Bio::Tools::Run::Primer3->new(-seq => $seq);
 
48
      $primer3run -> run; # run it with the default parameters
47
49
 
48
50
      # create a file to write the results to
49
 
      my $seqout=Bio::SeqIO->new(-file=>">primed_sequence.gbk", -format=>'genbank');
 
51
      my $seqout = Bio::SeqIO->new(-file => ">primed_sequence.gbk", 
 
52
                                   -format => 'genbank');
50
53
 
51
54
      # now just get all the results and write them out.
52
 
      while (my $results=$primer3run->next_primer) {
53
 
       $seqout->write_seq($results->annotated_seq);
 
55
      while (my $results = $primer3run->next_primer) {
 
56
         $seqout->write_seq($results->annotated_seq);
54
57
      }
55
58
 
56
 
   #(ii) to create a genbank file for a sequence and its cognate primers
57
 
 
58
 
     #For example:
 
59
Or, (ii), to create a genbank file for a sequence and its cognate primers:
59
60
 
60
61
     use Bio::SeqIO;
61
62
     use Bio::Seq::PrimedSeq;
63
64
     # have a sequence file ($file) with the template, and two primers
64
65
     # that match it, in fasta format
65
66
 
66
 
     my $file=shift || die "$0 <file>";
67
 
     my $seqin=new Bio::SeqIO(-file=>$file);
 
67
     my $file = shift || die "$0 <file>";
 
68
     my $seqin = new Bio::SeqIO(-file => $file);
68
69
 
69
70
     # read three sequences
70
71
     my ($template, $leftprimer, $rightprimer) =
71
 
           ($seqin->next_seq, $seqin->next_seq, , $seqin->next_seq);
 
72
           ($seqin->next_seq, $seqin->next_seq, $seqin->next_seq);
72
73
     # set up the primed sequence object
73
 
     my $primedseq = Bio::Seq::PrimedSeq->new(-seq=>$template, 
74
 
                                              -left_primer=>$leftprimer,
75
 
                                              -right_primer=>$rightprimer);
 
74
     my $primedseq = Bio::Seq::PrimedSeq->new(-seq => $template, 
 
75
                                              -left_primer => $leftprimer,
 
76
                                              -right_primer => $rightprimer);
76
77
     # open a file for output
77
 
     my $seqout=Bio::SeqIO->new(-file=>">primed_sequence.gbk",
78
 
                                -format=>'genbank');
 
78
     my $seqout = Bio::SeqIO->new(-file => ">primed_sequence.gbk",
 
79
                                  -format => 'genbank');
79
80
     # print the sequence out
80
81
     $seqout->write_seq($primedseq->annotated_sequence);
81
82
 
82
 
  # This should output a genbank file with the two primers labeled.
 
83
This should output a genbank file with the two primers labeled.
83
84
 
84
85
=head1 DESCRIPTION
85
86
 
86
 
This module is a slightly glorified capsule containg a primed sequence. It was
87
 
created to address the fact that a primer is more the a seqfeature and there
88
 
need to be ways to represent the primer-sequence complex and the behaviors and
89
 
attributes that are associated with the complex.
 
87
This module is a slightly glorified capsule containing a primed sequence. 
 
88
It was created to address the fact that a primer is more than a seqfeature 
 
89
and there need to be ways to represent the primer-sequence complex and 
 
90
the behaviors and attributes that are associated with the complex.
90
91
 
91
92
The primers are represented as Bio::SeqFeature::Primer objects, and should
92
 
be instatiated first.
93
 
 
94
 
The simplest way to initiate a PrimedSeq object is as follows:
95
 
 
96
 
  my $primedseq=Bio::Seq::PrimedSeq->new(
97
 
  -seq=>Bio::Seq object,
98
 
  -left_primer=>Bio::SeqFeature::Primer object,
99
 
  -right_primer=>Bio::SeqFeature::Primer object,
 
93
be instantiated first.
 
94
 
 
95
A simple way to create a PrimedSeq object is as follows:
 
96
 
 
97
  my $primedseq = Bio::Seq::PrimedSeq->new(
 
98
          -seq          => $seq,  # Bio::Seq object,
 
99
          -left_primer  => $left, # Bio::SeqFeature::Primer object,
 
100
          -right_primer => $right # Bio::SeqFeature::Primer object,
100
101
  );
101
102
 
102
103
From the PrimedSeq object you should be able to retrieve
103
 
information about Tm's and what not of each of the primers 
 
104
information about melting temperatures and what not on each of the primers 
104
105
and the amplicon.
105
106
 
106
107
This is based on the PrimedSeq.pm module started by Chad Matsalla, with 
112
113
Bioperl modules. Send your comments and suggestions preferably to one
113
114
of the Bioperl mailing lists.  Your participation is much appreciated.
114
115
 
115
 
    bioperl-l@bioperl.org          - General discussion
116
 
    http://bio.perl.org/MailList.html             - About the mailing lists
 
116
  bioperl-l@bioperl.org                  - General discussion
 
117
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
117
118
 
118
119
=head2 Reporting Bugs
119
120
 
120
121
Report bugs to the Bioperl bug tracking system to help us keep track
121
 
the bugs and their resolution.  Bug reports can be submitted via email
122
 
or the web:
 
122
the bugs and their resolution.  Bug reports can be submitted via the
 
123
web:
123
124
 
124
 
    bioperl-bugs@bio.perl.org
125
 
    http://bugzilla.bioperl.org/
 
125
  http://bugzilla.open-bio.org/
126
126
 
127
127
=head1 AUTHOR
128
128
 
142
142
 
143
143
package Bio::Seq::PrimedSeq;
144
144
 
145
 
 
146
145
use strict;
147
 
use Bio::RangeI;
148
 
use Bio::SeqFeature::Generic;
149
146
use Bio::SeqFeature::Primer;
150
147
 
151
 
use vars qw ($AUTOLOAD @RES %OK_FIELD @ISA $ID);
 
148
use vars qw ($AUTOLOAD @RES %OK_FIELD $ID);
152
149
 
153
 
@ISA = qw(Bio::Root::Root Bio::SeqFeature::Generic);
 
150
use base qw(Bio::Root::Root Bio::SeqFeature::Generic);
154
151
 
155
152
BEGIN {
156
 
 @RES=qw(); # nothing here yet, not sure what we want!
157
 
 
 
153
 @RES = qw(); # nothing here yet, not sure what we want!
158
154
 foreach my $attr (@RES) {$OK_FIELD{$attr}++}
159
155
}
160
156
 
161
157
$ID = 'Bio::Tools::Analysis::Nucleotide::PrimedSeq';
162
158
 
163
 
 
164
159
sub AUTOLOAD {
165
160
 my $self = shift;
166
161
 my $attr = $AUTOLOAD;
171
166
}
172
167
 
173
168
 
174
 
 
175
 
 
176
169
=head2 new
177
170
 
178
171
 Title   : new()
179
 
 Usage   : $primed_sequence = new Bio::SeqFeature::Primer( -seq => $sequence,
180
 
                                                            -left_primer => $left_primer,
181
 
                                                            -right_primer => $right_primer);
 
172
 Usage   : $primed_sequence = new Bio::SeqFeature::Primer( 
 
173
                                     -seq => $sequence,
 
174
                                     -left_primer => $left_primer,
 
175
                                     -right_primer => $right_primer);
182
176
 Function: A constructor for an object representing a primed sequence 
183
177
 Returns : A Bio::Seq::PrimedSeq object
184
 
 Args    : 
185
 
     -target_sequence => a Bio::Seq object (required)
186
 
     -left_primer => a Bio::SeqFeature::Primer object (required)
187
 
     -right_primer => a Bio::SeqFeature::Primer object (required)
 
178
 Args    :  -seq => a Bio::Seq object (required)
 
179
            -left_primer => a Bio::SeqFeature::Primer object (required)
 
180
            -right_primer => a Bio::SeqFeature::Primer object (required)
188
181
 
189
 
     Many other parameters can be included including all of the output
190
 
     parameters from the primer3 program. At the moment some of these
191
 
     parameters (most...) will not do anything.
 
182
           Many other parameters can be included including all of the output
 
183
           parameters from the primer3 program. At the moment most of these
 
184
           parameters will not do anything.
192
185
 
193
186
=cut
194
187
 
195
188
sub new {
196
189
 
197
 
 # note, I have cleaned up a lot of the script that Chad had written here,
198
 
 # and I have removed the part where he removed the - before the tags.
199
 
 # Very confusing.
 
190
        # note, I have cleaned up a lot of the script that Chad had written here,
 
191
        # and I have removed the part where he removed the - before the tags.
 
192
        # Very confusing.
200
193
 
201
 
 my($class,%args) = @_;
202
 
 my $self = $class->SUPER::new(%args);
 
194
        my($class,%args) = @_;
 
195
        my $self = $class->SUPER::new(%args);
203
196
   # these are the absolute minimum components required to make
204
197
   # a primedseq
205
198
 
206
 
 foreach my $key (keys %args) {
207
 
     if ($key eq "-seq" || $key eq "-SEQ") {
208
 
         $self->{target_sequence} = $args{$key};
209
 
         next;
210
 
     } else {
211
 
         my $okey;
212
 
         ($okey=$key)=~s/^-//;;
213
 
         if (($okey eq "left_primer" || $okey eq "right_primer") && 
214
 
             ref($args{$key}) && $args{$key}->isa('Bio::SeqI') ) {
215
 
             # we have been parsed a bio seq object. 
216
 
             # Make it a Bio::SeqFeature::Primer object
217
 
             $self->{$okey} = Bio::SeqFeature::Primer->new(-seq=>$args{$key});
218
 
             push @{$self->{'arguments'}},$okey;
219
 
             next;
220
 
         }
 
199
        foreach my $key (keys %args) {
 
200
                if ($key =~  /^-seq/i) {
 
201
                        $self->{target_sequence} = $args{$key};
 
202
                        next;
 
203
                } else {
 
204
                        my $okey;
 
205
                        ($okey = $key) =~ s/^-//;
 
206
                        if (($okey eq "left_primer" || $okey eq "right_primer") && 
 
207
                                 ref($args{$key}) && $args{$key}->isa('Bio::SeqI') ) {
 
208
                                # We have been given a Bio::Seq object, 
 
209
                                # make it a Bio::SeqFeature::Primer object
 
210
                                $self->{$okey} = Bio::SeqFeature::Primer->new(-seq => $args{$key});
 
211
                                push @{$self->{'arguments'}},$okey;
 
212
                                next;
 
213
                        }
221
214
 
222
 
         $self->{$okey} = $args{$key};
223
 
         push @{$self->{'arguments'}},$okey;
224
 
     }
225
 
 }
226
 
 # and now the insurance- make sure that things are ok
227
 
 if (!$self->{target_sequence} || !$self->{left_primer} || !$self->{right_primer} ) {
228
 
   $self->throw("You must provide a -target_sequence, -left_primer, and -right_primer to create this object.");
229
 
 }
 
215
                        $self->{$okey} = $args{$key};
 
216
                        push @{$self->{'arguments'}},$okey;
 
217
                }
 
218
        }
 
219
        # and now the insurance - make sure that things are ok
 
220
        if (!$self->{target_sequence} || !$self->{left_primer} || 
 
221
                 !$self->{right_primer} ) {
 
222
                $self->throw("You must provide a -seq, -left_primer, and -right_primer to create this object.");
 
223
        }
230
224
  
231
 
 if (! ref($self->{target_sequence}) ||
232
 
     ! $self->{target_sequence}->isa('Bio::SeqI') ) {
233
 
     $self->throw("The target_sequence must be a Bio::Seq to create this object.");
234
 
 }
235
 
 if (! ref($self->{left_primer}) ||
236
 
     ! $self->{left_primer}->isa("Bio::SeqFeature::Primer") || 
237
 
     ! ref($self->{right_primer}) ||
238
 
     ! $self->{right_primer}->isa("Bio::SeqFeature::Primer")) {
239
 
     $self->throw("You must provide a left_primer and right_primer, both as Bio::SeqFeature::Primer to create this object.");
240
 
 }
 
225
        if (! ref($self->{target_sequence}) ||
 
226
                 ! $self->{target_sequence}->isa('Bio::SeqI') ) {
 
227
                $self->throw("The target_sequence must be a Bio::Seq to create this object.");
 
228
 }
 
229
        if (! ref($self->{left_primer}) ||
 
230
                 ! $self->{left_primer}->isa("Bio::SeqFeature::Primer") || 
 
231
                 ! ref($self->{right_primer}) ||
 
232
                 ! $self->{right_primer}->isa("Bio::SeqFeature::Primer")) {
 
233
                $self->throw("You must provide a left_primer and right_primer, both as Bio::SeqFeature::Primer to create this object.");
 
234
        }
241
235
 
242
 
 # now we have the sequences, lets find out where they are
243
 
 $self->_place_seqs();
244
 
 return $self;
 
236
        # now we have the sequences, lets find out where they are
 
237
        $self->_place_seqs();
 
238
        return $self;
245
239
}
246
240
 
247
241
 
248
242
=head2 get_primer
249
243
 
250
244
 Title   : get_primer();
251
 
 Usage   : $primer = $primedseq->get_primer(l, left, left_primer, -left_primer); to return the left primer
252
 
           or 
253
 
           $primer = $primedseq->get_primer(r, right, right_primer, -right_primer); to return the right primer
254
 
           or
255
 
           $primer = $primedseq->get_primer(b, both, both_primers, -both_primers); to return the left primer, right primer array
 
245
 Usage   : $primer = $primedseq->get_primer(l, left, left_primer, 
 
246
           -left_primer) to return the left primer or 
 
247
                $primer = $primedseq->get_primer(r, right, right_primer, 
 
248
           -right_primer) to return the right primer or
 
249
                $primer = $primedseq->get_primer(b, both, both_primers, 
 
250
           -both_primers)
 
251
           to return the left primer, right primer array
256
252
 Function: A getter for the left primer in thie PrimedSeq object.
257
253
 Returns : A Bio::SeqFeature::Primer object
258
 
 Args    : Either of (l, left, left_primer, -left_primer) to get left primer.
259
 
           Either of (r, right, right_primer, -right_primer) to get right primer
260
 
           Either of (b, both, both_primers, -both_primers) to get both primers. Note that this is plural.
 
254
 Args    : Either of (l, left, left_primer, -left_primer) to get left 
 
255
           primer.
 
256
           Either of (r, right, right_primer, -right_primer) to get 
 
257
           right primer
 
258
                Either of (b, both, both_primers, -both_primers) to get 
 
259
           both primers. 
 
260
           Note that this is plural. [default]
261
261
 
262
262
=cut
263
263
 
264
264
sub get_primer() {
265
265
 my ($self, $arg) = @_;
266
 
 if ($arg =~ /^l/ || $arg =~ /^-l/) { 
 
266
 if (! defined $arg ) {
 
267
  return ($self->{'left_primer'}, $self->{'right_primer'});
 
268
 } elsif( $arg =~ /^l/ || $arg =~ /^-l/) { 
267
269
  # what a cheat, I couldn't be bothered to write all those or statements!
268
270
  # Hah, now you can write leprechaun to get the left primer.
269
271
  return $self->{'left_primer'}
278
280
 
279
281
 Title   : annotated_sequence
280
282
 Usage   : $annotated_sequence_object = $primedseq->annotated_sequence()
281
 
 Function: Get an annotated sequence object containg the left and right primers
 
283
 Function: Get an annotated sequence object containg the left and right 
 
284
           primers
282
285
 Returns : An annotated sequence object or 0 if not defined.
283
286
 Args    : 
284
 
 Note    : Use this method to return a sequence object that you can write out (e.g. in GenBank format)
285
 
           See the example above.
 
287
 Note    : Use this method to return a sequence object that you can write
 
288
           out (e.g. in GenBank format). See the example above.
286
289
 
287
290
=cut
288
291
 
306
309
sub amplicon {
307
310
 my ($self,@args) = @_;
308
311
 my $id = $self->{'-seq'}->{'id'};
309
 
 unless ($id) {$id=""}# this just prevents a warning when $self->{'-seq'}->{'id'} is not defined
 
312
 unless ($id) {$id=""}
 
313
 # this just prevents a warning when $self->{'-seq'}->{'id'} is not defined
310
314
 $id = "Amplicon from ".$id;
311
315
 
312
316
 my $seqobj=Bio::Seq->new(-id=>$id, seq=>$self->{'amplicon_sequence'});
314
318
}
315
319
 
316
320
 
317
 
 
318
321
=head2 seq
319
322
 
320
323
 Title   : seq
335
338
 
336
339
 Title   : _place_seqs
337
340
 Usage   : $self->_place_seqs()
338
 
 Function: An internal method to place the primers on the sequence, and set up the ranges of the sequences
 
341
 Function: An internal method to place the primers on the sequence and 
 
342
           set up the ranges of the sequences
339
343
 Returns : Nothing
340
344
 Args    : None
341
345
 Note    : Internal use only
343
347
=cut
344
348
 
345
349
sub _place_seqs {
346
 
 my $self = shift;
347
 
 
348
 
 # we are going to pull out the target sequence, and then the primer sequences
349
 
 my $target_sequence = $self->{'target_sequence'}->seq();
350
 
 
351
 
 # left primer
352
 
 my $left_seq = $self->{'left_primer'}->seq()->seq();
353
 
 
354
 
 my $rprc = $self->{'right_primer'}->seq()->revcom();
355
 
 
356
 
 my $right_seq=$rprc->seq();
357
 
  
358
 
  
359
 
 # now just change the case, because we keep getting screwed on this
360
 
 $target_sequence=uc($target_sequence);
361
 
 $left_seq=uc($left_seq);
362
 
 $right_seq=uc($right_seq);
363
 
 
364
 
 unless ($target_sequence =~ /(.*)$left_seq(.*)$right_seq(.*)/) {
365
 
  unless ($target_sequence =~ /$left_seq/) {$self->throw("Can't place left sequence on target!")}
366
 
  unless ($target_sequence =~ /$right_seq/) {$self->throw("Can't place right sequence on target!")}
367
 
 }
368
 
 
369
 
 my ($before, $middle, $after) = ($1, $2, $3); # note didn't use $`, $', and $& because they are bad. Just use length instead.
370
 
 
371
 
 # cool now we can figure out lengths and what not.
372
 
 # we'll figure out the position and compare it to known positions (e.g. from primer3)
373
 
 
374
 
 my $left_location = length($before). ",". length($left_seq);
375
 
 my $right_location = (length($target_sequence)-length($after)-1).",".length($right_seq);
376
 
 my $amplicon_size = length($left_seq)+length($middle)+length($right_seq);
377
 
 
378
 
 if (exists $self->{'left_primer'}->{'PRIMER_LEFT'}) {
379
 
  # this is the left primer from primer3 input
380
 
  # just check to make sure it is right
381
 
  unless ($self->{'left_primer'}->{'PRIMER_LEFT'} eq $left_location) {
382
 
   $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.");
383
 
  }
384
 
 }
385
 
 else {
386
 
  $self->{'left_primer'}->{'PRIMER_LEFT'}=$left_location;
387
 
 }
388
 
 
389
 
 if (exists $self->{'right_primer'}->{'PRIMER_RIGHT'}) {
390
 
  # this is the right primer from primer3 input
391
 
  # just check to make sure it is right
392
 
  unless ($self->{'right_primer'}->{'PRIMER_RIGHT'} eq $right_location) {
393
 
   $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.");
394
 
  }
395
 
 }
396
 
 else {
397
 
  $self->{'right_primer'}->{'PRIMER_RIGHT'}=$right_location;
398
 
 }
399
 
 
400
 
 if (exists $self->{'PRIMER_PRODUCT_SIZE'}) {
401
 
  # this is the product size from primer3 input
402
 
  # just check to make sure it is right
403
 
  unless ($self->{'PRIMER_PRODUCT_SIZE'} eq $amplicon_size) {
404
 
   $self->warn("Note got |".$self->{'PRIMER_PRODUCT_SIZE'}."| from primer3 and |$amplicon_size| for the size. You should email redwards\@utmem.edu about this.");
405
 
  }
406
 
 }
407
 
 else {
408
 
  $self->{'PRIMER_PRODUCT_SIZE'} = $amplicon_size;
409
 
 }
410
 
 
411
 
 $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
412
 
 
413
 
 $self->_set_seqfeature;
 
350
        my $self = shift;
 
351
 
 
352
        # we are going to pull out the target sequence, and then the primer sequences
 
353
        my $target_sequence = $self->{'target_sequence'}->seq();
 
354
 
 
355
        # left primer
 
356
        my $left_seq = $self->{'left_primer'}->seq()->seq();
 
357
 
 
358
        my $rprc = $self->{'right_primer'}->seq()->revcom();
 
359
 
 
360
        my $right_seq=$rprc->seq();
 
361
  
 
362
        # now just change the case, because we keep getting screwed on this
 
363
        $target_sequence=uc($target_sequence);
 
364
        $left_seq=uc($left_seq);
 
365
        $right_seq=uc($right_seq);
 
366
 
 
367
        unless ($target_sequence =~ /(.*)$left_seq(.*)$right_seq(.*)/) {
 
368
                unless ($target_sequence =~ /$left_seq/) {$self->throw("Can't place left sequence on target!")}
 
369
                unless ($target_sequence =~ /$right_seq/) {$self->throw("Can't place right sequence on target!")}
 
370
 }
 
371
 
 
372
        my ($before, $middle, $after) = ($1, $2, $3); # note didn't use $`, $', and $& because they are bad. Just use length instead.
 
373
 
 
374
        # cool now we can figure out lengths and what not.
 
375
        # we'll figure out the position and compare it to known positions (e.g. from primer3)
 
376
 
 
377
        my $left_location = length($before). ",". length($left_seq);
 
378
        my $right_location = (length($target_sequence)-length($after)-1).",".length($right_seq);
 
379
        my $amplicon_size = length($left_seq)+length($middle)+length($right_seq);
 
380
 
 
381
        if (exists $self->{'left_primer'}->{'PRIMER_LEFT'}) {
 
382
                # this is the left primer from primer3 input
 
383
                # just check to make sure it is right
 
384
                unless ($self->{'left_primer'}->{'PRIMER_LEFT'} eq $left_location) {
 
385
                        $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.");
 
386
                }
 
387
        }
 
388
        else {
 
389
                $self->{'left_primer'}->{'PRIMER_LEFT'}=$left_location;
 
390
        }
 
391
 
 
392
        if (exists $self->{'right_primer'}->{'PRIMER_RIGHT'}) {
 
393
                # this is the right primer from primer3 input
 
394
                # just check to make sure it is right
 
395
                unless ($self->{'right_primer'}->{'PRIMER_RIGHT'} eq $right_location) {
 
396
                        $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.");
 
397
                }
 
398
        }
 
399
        else {
 
400
                $self->{'right_primer'}->{'PRIMER_RIGHT'}=$right_location;
 
401
        }
 
402
 
 
403
        if (exists $self->{'PRIMER_PRODUCT_SIZE'}) {
 
404
                # this is the product size from primer3 input
 
405
                # just check to make sure it is right
 
406
                unless ($self->{'PRIMER_PRODUCT_SIZE'} eq $amplicon_size) {
 
407
                        $self->warn("Note got |".$self->{'PRIMER_PRODUCT_SIZE'}."| from primer3 and |$amplicon_size| for the size. You should email redwards\@utmem.edu about this.");
 
408
                }
 
409
        }
 
410
        else {
 
411
                $self->{'PRIMER_PRODUCT_SIZE'} = $amplicon_size;
 
412
        }
 
413
 
 
414
        $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
 
415
        
 
416
        $self->_set_seqfeature;
414
417
}
415
418
 
416
419
=head2 _set_seqfeature
417
420
 
418
421
 Title   : _set_seqfeature
419
422
 Usage   : $self->_set_seqfeature()
420
 
 Function: An internal method to create Bio::SeqFeature::Generic objects for the primed seq
 
423
 Function: An internal method to create Bio::SeqFeature::Generic objects
 
424
           for the primed seq
421
425
 Returns : Nothing
422
426
 Args    : None
423
 
 Note    : Internal use only. Should only call this once left and right primers have been placed on the sequence
424
 
           This will then set them as sequence features so hopefully we can get a nice output with write_seq
 
427
 Note    : Internal use only. Should only call this once left and right 
 
428
           primers have been placed on the sequence. This will then set 
 
429
           them as sequence features so hopefully we can get a nice output 
 
430
           with write_seq.
425
431
 
426
432
=cut
427
433
 
428
434
 
429
435
sub _set_seqfeature {
430
 
 my $self = shift;
431
 
 unless ($self->{'left_primer'}->{'PRIMER_LEFT'} && $self->{'right_primer'}->{'PRIMER_RIGHT'}) {
432
 
  $self->warn("hmmm. Haven't placed primers, but trying to make annotated sequence");
433
 
  return 0;
434
 
 }
435
 
 my ($start, $length) = split /,/, $self->{'left_primer'}->{'PRIMER_LEFT'};
436
 
 my $tm=$self->{'left_primer'}->{'PRIMER_LEFT_TM'} || $self->{'left_primer'}->Tm || 0;
437
 
 
438
 
 my $seqfeatureL=new Bio::SeqFeature::Generic(
439
 
   -start => $start, -end => $start+$length, -strand => 1,
440
 
   -primary_tag => 'left_primer', -source => 'primer3',
441
 
   -tag    => {new => 1, author => 'Bio::Seq::PrimedSeq', Tm => $tm}
442
 
    );
443
 
 
444
 
 ($start, $length) = split /,/, $self->{'right_primer'}->{'PRIMER_RIGHT'};
445
 
 $tm=$self->{'right_primer'}->{'PRIMER_RIGHT_TM'} || $self->{'right_primer'}->Tm || 0;
446
 
 
447
 
 my $seqfeatureR=new Bio::SeqFeature::Generic(
448
 
   -start => $start-$length, -end => $start, -strand => -1,
 
436
        my $self = shift;
 
437
        unless ($self->{'left_primer'}->{'PRIMER_LEFT'} && 
 
438
                          $self->{'right_primer'}->{'PRIMER_RIGHT'}) {
 
439
                $self->warn("hmmm. Haven't placed primers, but trying to make annotated sequence");
 
440
                return 0;
 
441
        }
 
442
        my ($start, $length) = split /,/, $self->{'left_primer'}->{'PRIMER_LEFT'};
 
443
        my $tm=$self->{'left_primer'}->{'PRIMER_LEFT_TM'} || $self->{'left_primer'}->Tm || 0;
 
444
 
 
445
        my $seqfeatureL=new Bio::SeqFeature::Generic(
 
446
                                                  -start => $start+1, -end => $start+$length, -strand => 1,
 
447
                    -primary_tag => 'left_primer', -source => 'primer3',
 
448
                    -tag    => {new => 1, author => 'Bio::Seq::PrimedSeq', Tm => $tm}
 
449
                                                                                                                          );
 
450
 
 
451
        ($start, $length) = split /,/, $self->{'right_primer'}->{'PRIMER_RIGHT'};
 
452
        $tm=$self->{'right_primer'}->{'PRIMER_RIGHT_TM'} || $self->{'right_primer'}->Tm || 0;
 
453
 
 
454
        my $seqfeatureR=new Bio::SeqFeature::Generic(
 
455
   -start => $start-$length+2, -end => $start+1, -strand => -1,
449
456
   -primary_tag => 'right_primer', -source => 'primer3',
450
457
   -tag    => {new => 1, author => 'Bio::Seq::PrimedSeq', Tm => $tm}
451
 
    );
 
458
                                                                                                                          );
452
459
 
453
 
 # now add the sequences to a annotated sequence
454
 
 $self->{annotated_sequence} = $self->{target_sequence};
455
 
 $self->{annotated_sequence}->add_SeqFeature($seqfeatureL);
456
 
 $self->{annotated_sequence}->add_SeqFeature($seqfeatureR);
 
460
        # now add the sequences to a annotated sequence
 
461
        $self->{annotated_sequence} = $self->{target_sequence};
 
462
        $self->{annotated_sequence}->add_SeqFeature($seqfeatureL);
 
463
        $self->{annotated_sequence}->add_SeqFeature($seqfeatureR);
457
464
}
458
465
 
459
466
1;