28
Bio::Seq::PrimedSeq - A representation of a sequence and two primers
29
flanking a target region
28
Bio::Seq::PrimedSeq - A sequence and a pair of primers matching on it
33
# The easiest way to use this is probably either, (i), get the
34
# output from Bio::Tools::Run::Primer3, Bio::Tools::Primer3, or
35
# Bio::Tools::PCRSimulation:
37
# For example, start with a fasta file
39
use Bio::Tools::Run::Primer3;
41
my $file = shift || die "need a file to read";
42
my $seqin = Bio::SeqIO->new(-file => $file);
43
my $seq = $seqin->next_seq;
45
# use primer3 to design some primers
46
my $primer3run = Bio::Tools::Run::Primer3->new(-seq => $seq);
47
$primer3run -> run; # run it with the default parameters
49
# create a file to write the results to
50
my $seqout = Bio::SeqIO->new(-file => ">primed_sequence.gbk",
51
-format => 'genbank');
53
# now just get all the results and write them out.
54
while (my $results = $primer3run->next_primer) {
55
$seqout->write_seq($results->annotated_seq);
58
# Or, (ii), to create a genbank file for a sequence and its cognate
62
use Bio::Seq::PrimedSeq;
64
# have a sequence file ($file) with the template, and two primers
65
# that match it, in fasta format
67
my $file = shift || die "$0 <file>";
68
my $seqin = Bio::SeqIO->new(-file => $file);
70
# read three sequences
71
my ($template, $leftprimer, $rightprimer) =
72
($seqin->next_seq, $seqin->next_seq, $seqin->next_seq);
73
# set up the primed sequence object
74
my $primedseq = Bio::Seq::PrimedSeq->new(-seq => $template,
75
-left_primer => $leftprimer,
76
-right_primer => $rightprimer);
77
# open a file for output
78
my $seqout = Bio::SeqIO->new(-file => ">primed_sequence.gbk",
79
-format => 'genbank');
80
# print the sequence out
81
$seqout->write_seq($primedseq->annotated_sequence);
83
# This should output a genbank file with the two primers labeled.
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.
92
The primers are represented as Bio::SeqFeature::Primer objects, and should
93
be instantiated first.
95
A simple way to create a PrimedSeq object is as follows:
33
use Bio::Seq::PrimedSeq;
35
my $template = Bio::Seq->new( -seq => 'AGCTTTTCATTCTGACTGCAAC' );
36
my $left = Bio::Seq->new( -seq => 'AGCT' );
37
my $right = Bio::Seq->new( -seq => 'GTTGC' );
97
39
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,
40
-seq => $template, # sequence object
41
-left_primer => $left, # sequence or primer object
42
-right_primer => $right # sequence or primer object
103
From the PrimedSeq object you should be able to retrieve
104
information about melting temperatures and what not on each of the primers
107
This is based on the PrimedSeq.pm module started by Chad Matsalla, with
108
additions/improvements by Rob Edwards.
45
# Get the primers as Bio::SeqFeature::Primer objects
46
my @primer_objects = $primedseq->get_primer();
48
# Sequence object representing the PCR product, i.e. the section of the target
49
# sequence contained between the primers (primers included)
50
my $amplicon_seq = $primedseq->amplicon();
52
print 'Got amplicon sequence: '.$amplicon_seq->seq."\n";
53
# Amplicon should be: agctTTTCATTCTGACTgcaac
57
This module was created to address the fact that a primer is more than a
58
SeqFeature and that there is a need to represent the primer-sequence complex and
59
the attributes that are associated with the complex.
61
A PrimedSeq object is initialized with a target sequence and two primers. The
62
location of the primers on the target sequence is calculated if needed so that
63
one can then get the PCR product, or amplicon sequence. From the PrimedSeq object
64
you can also retrieve information about melting temperatures and what not on each
65
of the primers and the amplicon. The L<Bio::Tools::Primer3> module uses PrimedSeq
68
Note that this module does not simulate PCR. It assumes that the primers
69
will match in a single location on the target sequence and does not understand
76
Providing primers as sequence objects
78
If the primers are specified as sequence objects, e.g. L<Bio::PrimarySeq> or
79
L<Bio::Seq>, they are first converted to L<Bio::SeqFeature::Primer> objects.
80
Their location on the target sequence is then calculated and added to the
81
L<Bio::SeqFeature::Primer> objects, which you can retrieve using the get_primer()
86
Providing primers as primer objects
88
Because of the limitations of specifying primers as sequence objects, the
89
recommended method is to provide L<Bio::SeqFeature::Primer> objects. If you did
90
not record the location of the primers in the primer object, it will be
95
L<Bio::Seq::PrimedSeq> was initially started by Chad Matsalla, and later
96
improved on by Rob Edwards.
104
Run Primer3 to get PrimedSeq objects:
107
use Bio::Tools::Run::Primer3;
109
# Read a target sequences from a FASTA file
110
my $file = shift || die "Need a file to read";
111
my $seqin = Bio::SeqIO->new(-file => $file);
112
my $seq = $seqin->next_seq;
114
# Use Primer3 to design some primers
115
my $primer3 = Bio::Tools::Run::Primer3->new(-seq => $seq);
116
my $results = $primer3->run; # default parameters
118
# Write all the results in a Genbank file
119
my $seqout = Bio::SeqIO->new(-file => ">primed_sequence.gbk",
120
-format => 'genbank');
121
while (my $primedseq = $results->next_primer) {
122
$seqout->write_seq( $primedseq->annotated_seq );
127
Create a genbank file for a sequence and its cognate primers:
130
use Bio::Seq::PrimedSeq;
132
# Read a FASTA file that contains the target sequence and its two primers
133
my $file = shift || die "$0 <file>";
134
my $seqin = Bio::SeqIO->new(-file => $file);
135
my ($template, $leftprimer, $rightprimer) =
136
($seqin->next_seq, $seqin->next_seq, $seqin->next_seq);
138
# Set up a PrimedSeq object
139
my $primedseq = Bio::Seq::PrimedSeq->new(-seq => $template,
140
-left_primer => $leftprimer,
141
-right_primer => $rightprimer);
143
# Write the sequences in an output Genbank file
144
my $seqout = Bio::SeqIO->new(-file => ">primed_sequence.gbk",
145
-format => 'genbank');
146
$seqout->write_seq($primedseq->annotated_sequence);
157
197
use Bio::SeqFeature::Primer;
159
use vars qw ($AUTOLOAD @RES %OK_FIELD $ID);
161
199
use base qw(Bio::Root::Root Bio::SeqFeature::Generic);
164
@RES = qw(); # nothing here yet, not sure what we want!
165
foreach my $attr (@RES) {$OK_FIELD{$attr}++}
168
$ID = 'Bio::Tools::Analysis::Nucleotide::PrimedSeq';
172
my $attr = $AUTOLOAD;
174
$self->throw("Unallowed parameter: $attr !") unless $OK_FIELD{$attr};
175
$self->{$attr} = shift if @_;
176
return $self->{$attr};
200
# Since this module occupies the Bio::Seq::* namespace, it should probably
201
# inherit from Bio::Seq before it inherits from Bio::SeqFeature::Generic. But
202
# then, Bio::SeqI and Bio::SeqFeatureI both request a seq() method that return
203
# different things. So, being Bio::SeqI is incompatible with being Bio::SeqFeatureI
183
Usage : $primed_sequence = Bio::SeqFeature::Primer->new(
185
-left_primer => $left_primer,
186
-right_primer => $right_primer);
187
Function: A constructor for an object representing a primed sequence
209
Usage : my $primedseq = Bio::SeqFeature::Primer->new(
211
-left_primer => $left_primer,
212
-right_primer => $right_primer
214
Function: Construct a primed sequence.
188
215
Returns : A Bio::Seq::PrimedSeq object
189
216
Args : -seq => a Bio::Seq object (required)
190
-left_primer => a Bio::SeqFeature::Primer object (required)
191
-right_primer => a Bio::SeqFeature::Primer object (required)
217
-left_primer => a Bio::SeqFeature::Primer or sequence object (required)
218
-right_primer => a Bio::SeqFeature::Primer or sequence object (required)
220
If you pass a sequence object to specify a primer, it will be used to
221
construct a Bio::SeqFeature::Primer that you can retrieve with the
222
L<get_primer> method.
193
224
Many other parameters can be included including all of the output
194
parameters from the primer3 program. At the moment most of these
195
parameters will not do anything.
225
parameters from the primer3 program (see L<Bio::Tools::Primer3>). At
226
the moment these parameters will simply be stored and do anything.
201
# note, I have cleaned up a lot of the script that Chad had written here,
202
# and I have removed the part where he removed the - before the tags.
205
my($class,%args) = @_;
206
my $self = $class->SUPER::new(%args);
207
# these are the absolute minimum components required to make
210
foreach my $key (keys %args) {
211
if ($key =~ /^-seq/i) {
212
$self->{target_sequence} = $args{$key};
216
($okey = $key) =~ s/^-//;
217
if (($okey eq "left_primer" || $okey eq "right_primer") &&
218
ref($args{$key}) && $args{$key}->isa('Bio::SeqI') ) {
219
# We have been given a Bio::Seq object,
220
# make it a Bio::SeqFeature::Primer object
221
$self->{$okey} = Bio::SeqFeature::Primer->new(-seq => $args{$key});
222
push @{$self->{'arguments'}},$okey;
226
$self->{$okey} = $args{$key};
227
push @{$self->{'arguments'}},$okey;
230
# and now the insurance - make sure that things are ok
231
if (!$self->{target_sequence} || !$self->{left_primer} ||
232
!$self->{right_primer} ) {
233
$self->throw("You must provide a -seq, -left_primer, and -right_primer to create this object.");
236
if (! ref($self->{target_sequence}) ||
237
! $self->{target_sequence}->isa('Bio::SeqI') ) {
238
$self->throw("The target_sequence must be a Bio::Seq to create this object.");
240
if (! ref($self->{left_primer}) ||
241
! $self->{left_primer}->isa("Bio::SeqFeature::Primer") ||
242
! ref($self->{right_primer}) ||
243
! $self->{right_primer}->isa("Bio::SeqFeature::Primer")) {
244
$self->throw("You must provide a left_primer and right_primer, both as Bio::SeqFeature::Primer to create this object.");
247
# now we have the sequences, lets find out where they are
248
$self->_place_seqs();
231
my($class,%args) = @_;
232
my $self = $class->SUPER::new(%args);
234
# Need an amplicon sequence
235
$self->{seq} = delete $args{-seq} || delete $args{-target_sequence} ||
236
$self->throw("Need to provide a sequence during PrimedSeq object construction");
237
if (! ref($self->{seq}) || ! $self->{seq}->isa('Bio::SeqI') ) {
238
$self->throw("The target_sequence must be a Bio::Seq to create this object.");
241
# Need a left and right primers
242
for my $primer ( 'left', 'right' ) {
243
$self->{$primer} = delete $args{'-'.$primer.'_primer'} ||
244
$self->throw("Need to provide both primers during PrimedSeq object construction");
245
if ( ref $self->{$primer} && $self->{$primer}->isa('Bio::PrimarySeqI') ) {
246
# Convert Bio::Seq or Bio::PrimarySeq objects to Bio::SeqFeature::Primer
247
$self->{$primer} = Bio::SeqFeature::Primer->new(-seq => $self->{$primer});
249
if (not (ref $self->{$primer} && $self->{$primer}->isa("Bio::SeqFeature::Primer"))) {
250
$self->throw("Primers must be Bio::SeqFeature::Primer objects but got a ".ref($self->{$primer}));
254
# Save remaining arguments
255
while (my ($arg, $val) = each %args) {
256
$self->{$arg} = $val;
259
# Determine primer location on target if needed
260
if ( not( $self->{'left'}->start && $self->{'left'}->end &&
261
$self->{'right'}->start && $self->{'right'}->end ) ) {
262
$self->_place_primers();
253
269
=head2 get_primer
255
271
Title : get_primer();
256
Usage : $primer = $primedseq->get_primer(l, left, left_primer,
257
-left_primer) to return the left primer or
258
$primer = $primedseq->get_primer(r, right, right_primer,
259
-right_primer) to return the right primer or
260
$primer = $primedseq->get_primer(b, both, both_primers,
262
to return the left primer, right primer array
263
Function: A getter for the left primer in thie PrimedSeq object.
272
Usage : my @primers = $primedseq->get_primer();
274
my $primer = $primedseq->get_primer('-left_primer');
275
Function: Get the primers associated with the PrimedSeq object.
264
276
Returns : A Bio::SeqFeature::Primer object
265
Args : Either of (l, left, left_primer, -left_primer) to get left
267
Either of (r, right, right_primer, -right_primer) to get
269
Either of (b, both, both_primers, -both_primers) to get
271
Note that this is plural. [default]
277
Args : For the left primer, use: l, left, left_primer or -left_primer
278
For the right primer, use: r, right, right_primer or -right_primer
279
For both primers [default], use: b, both, both_primers or -both_primers
276
my ($self, $arg) = @_;
277
if (! defined $arg ) {
278
return ($self->{'left_primer'}, $self->{'right_primer'});
279
} elsif( $arg =~ /^l/ || $arg =~ /^-l/) {
280
# what a cheat, I couldn't be bothered to write all those or statements!
281
# Hah, now you can write leprechaun to get the left primer.
282
return $self->{'left_primer'}
284
elsif ($arg =~ /^r/ || $arg =~ /^-r/) {return $self->{'right_primer'}}
285
elsif ($arg =~ /^b/ || $arg =~ /^-b/) {return ($self->{'left_primer'}, $self->{'right_primer'})}
284
my ($self, $arg) = @_;
285
if (! defined $arg ) {
286
return ($self->{left}, $self->{right});
287
} elsif ( $arg =~ /^-?l/ ) {
288
# What a cheat, I couldn't be bothered to write all the exact statements!
289
# Hah, now you can write 'leprechaun' to get the left primer.
291
} elsif ( $arg =~ /^-?r/ ) {
292
return $self->{right};
293
} elsif ( $arg =~ /^-?b/ ) {
294
return ($self->{left}, $self->{right});
290
299
=head2 annotated_sequence
292
301
Title : annotated_sequence
293
Usage : $annotated_sequence_object = $primedseq->annotated_sequence()
302
Usage : my $annotated_sequence_object = $primedseq->annotate_sequence();
294
303
Function: Get an annotated sequence object containg the left and right
296
305
Returns : An annotated sequence object or 0 if not defined.
303
312
sub annotated_sequence {
305
if (exists $self->{annotated_sequence}) {return $self->{annotated_sequence}}
314
my $seq = $self->{'seq'}; ### clone??
315
$seq->add_SeqFeature($self->{'left'});
316
$seq->add_SeqFeature($self->{'right'});
312
Usage : my $amplicon = $primedseq->amplicon()
313
Function: Retrieve the amplicon as a sequence object
314
Returns : A seq object. To get the sequence use $amplicon->seq
324
Usage : my $amplicon = $primedseq->amplicon();
325
Function: Retrieve the amplicon as a sequence object. The amplicon sequnce is
326
only the section of the target sequence between the primer matches
328
Returns : A Bio::Seq object. To get the sequence use $amplicon->seq
321
my ($self,@args) = @_;
322
my $id = $self->{'-seq'}->{'id'};
323
unless ($id) {$id=""}
324
# this just prevents a warning when $self->{'-seq'}->{'id'} is not defined
325
$id = "Amplicon from ".$id;
327
my $seqobj=Bio::Seq->new(-id=>$id, seq=>$self->{'amplicon_sequence'});
336
my $left = $self->{left};
337
my $right = $self->{right};
338
my $target = $self->{seq};
339
return Bio::PrimarySeq->new(
340
-id => 'Amplicon_from_'.($target->id || 'unidentified'),
341
-seq => lc( $left->seq->seq ).
342
uc( $target->subseq($left->end+1, $right->start-1) ).
343
lc( $right->seq->revcom->seq ),
335
Usage : my $seqobj = $primedseq->seq()
351
Usage : my $seqobj = $primedseq->seq();
336
352
Function: Retrieve the target sequence as a sequence object
337
353
Returns : A seq object. To get the sequence use $seqobj->seq
363
# we are going to pull out the target sequence, and then the primer sequences
364
my $target_sequence = $self->{'target_sequence'}->seq();
367
my $left_seq = $self->{'left_primer'}->seq()->seq();
369
my $rprc = $self->{'right_primer'}->seq()->revcom();
371
my $right_seq=$rprc->seq();
373
# now just change the case, because we keep getting screwed on this
374
$target_sequence=uc($target_sequence);
375
$left_seq=uc($left_seq);
376
$right_seq=uc($right_seq);
378
unless ($target_sequence =~ /(.*)$left_seq(.*)$right_seq(.*)/) {
379
unless ($target_sequence =~ /$left_seq/) {$self->throw("Can't place left sequence on target!")}
380
unless ($target_sequence =~ /$right_seq/) {$self->throw("Can't place right sequence on target!")}
383
my ($before, $middle, $after) = ($1, $2, $3); # note didn't use $`, $', and $& because they are bad. Just use length instead.
385
# cool now we can figure out lengths and what not.
386
# we'll figure out the position and compare it to known positions (e.g. from primer3)
388
my $left_location = length($before). ",". length($left_seq);
389
my $right_location = (length($target_sequence)-length($after)-1).",".length($right_seq);
390
my $amplicon_size = length($left_seq)+length($middle)+length($right_seq);
392
if (exists $self->{'left_primer'}->{'PRIMER_LEFT'}) {
393
# this is the left primer from primer3 input
394
# just check to make sure it is right
395
unless ($self->{'left_primer'}->{'PRIMER_LEFT'} eq $left_location) {
396
$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.");
400
$self->{'left_primer'}->{'PRIMER_LEFT'}=$left_location;
403
if (exists $self->{'right_primer'}->{'PRIMER_RIGHT'}) {
404
# this is the right primer from primer3 input
405
# just check to make sure it is right
406
unless ($self->{'right_primer'}->{'PRIMER_RIGHT'} eq $right_location) {
407
$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.");
411
$self->{'right_primer'}->{'PRIMER_RIGHT'}=$right_location;
414
if (exists $self->{'PRIMER_PRODUCT_SIZE'}) {
415
# this is the product size from primer3 input
416
# just check to make sure it is right
417
unless ($self->{'PRIMER_PRODUCT_SIZE'} eq $amplicon_size) {
418
$self->warn("Note got |".$self->{'PRIMER_PRODUCT_SIZE'}."| from primer3 and |$amplicon_size| for the size. You should email redwards\@utmem.edu about this.");
422
$self->{'PRIMER_PRODUCT_SIZE'} = $amplicon_size;
425
$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
427
$self->_set_seqfeature;
430
=head2 _set_seqfeature
432
Title : _set_seqfeature
433
Usage : $self->_set_seqfeature()
434
Function: An internal method to create Bio::SeqFeature::Generic objects
438
Note : Internal use only. Should only call this once left and right
439
primers have been placed on the sequence. This will then set
440
them as sequence features so hopefully we can get a nice output
446
sub _set_seqfeature {
448
unless ($self->{'left_primer'}->{'PRIMER_LEFT'} &&
449
$self->{'right_primer'}->{'PRIMER_RIGHT'}) {
450
$self->warn("hmmm. Haven't placed primers, but trying to make annotated sequence");
453
my ($start, $length) = split /,/, $self->{'left_primer'}->{'PRIMER_LEFT'};
454
my $tm=$self->{'left_primer'}->{'PRIMER_LEFT_TM'} || $self->{'left_primer'}->Tm || 0;
456
my $seqfeatureL=Bio::SeqFeature::Generic->new(
457
-start => $start+1, -end => $start+$length, -strand => 1,
458
-primary_tag => 'left_primer', -source => 'primer3',
459
-tag => {new => 1, author => 'Bio::Seq::PrimedSeq', Tm => $tm}
462
($start, $length) = split /,/, $self->{'right_primer'}->{'PRIMER_RIGHT'};
463
$tm=$self->{'right_primer'}->{'PRIMER_RIGHT_TM'} || $self->{'right_primer'}->Tm || 0;
465
my $seqfeatureR=Bio::SeqFeature::Generic->new(
466
-start => $start-$length+2, -end => $start+1, -strand => -1,
467
-primary_tag => 'right_primer', -source => 'primer3',
468
-tag => {new => 1, author => 'Bio::Seq::PrimedSeq', Tm => $tm}
471
# now add the sequences to a annotated sequence
472
$self->{annotated_sequence} = $self->{target_sequence};
473
$self->{annotated_sequence}->add_SeqFeature($seqfeatureL);
474
$self->{annotated_sequence}->add_SeqFeature($seqfeatureR);
380
# Get the target and primer sequence strings, all in uppercase
381
my $left = $self->{left};
382
my $right = $self->{right};
383
my $target_seq = uc $self->{seq}->seq();
384
my $left_seq = uc $left->seq()->seq();
385
my $right_seq = uc $right->seq()->revcom()->seq();
387
# Locate primers on target sequence
388
my ($before, $middle, $after) = ($target_seq =~ /^(.*)$left_seq(.*)$right_seq(.*)$/);
389
if (not defined $before || not defined $middle || not defined $after) {
390
if ($target_seq !~ /$left_seq/) {
391
$self->throw("Could not place left primer on target");
393
if ($target_seq !~ /$right_seq/) {
394
$self->throw("Could not place right primer on target")
398
# Save location information in primer object
399
my $left_location = length($before) + 1;
400
my $right_location = length($target_seq) - length($after);
402
$left->start($left_location);
403
$left->end($left_location + $left->seq->length - 1);
405
$right->start($right_location - $right->seq->length + 1);
406
$right->end($right_location);
409
# If Primer3 information was recorded, compare it to what we calculated
410
if ( exists($left->{PRIMER_LEFT}) || exists($right->{PRIMER_RIGHT}) || exists($self->{PRIMER_PRODUCT_SIZE}) ) {
412
# Bio::Seq::PrimedSeq positions
413
my $amplicon_size = length($left_seq) + length($middle) + length($right_seq);
414
$left_location = $left_location.','.length($left_seq);
415
$right_location = $right_location.','.length($right_seq);
418
my $primer_product = $self->{PRIMER_PRODUCT_SIZE};
419
my $primer_left = $left->{PRIMER_LEFT};
420
my $primer_right = $right->{PRIMER_RIGHT};
422
if ( defined($primer_left) && (not $primer_left eq $left_location) ) {
423
$self->warn("Got |$primer_left| from primer3 but calculated ".
424
"|$left_location| for the left primer.");
426
if ( defined($primer_right) && (not $primer_right eq $right_location) ) {
427
$self->warn("Got |$primer_right| from primer3 but calculated ".
428
"|$right_location| for the right primer.");
430
if ( defined($primer_product) && (not $primer_product eq $amplicon_size) ) {
431
$self->warn("Got |$primer_product| from primer3 but calculated ".
432
"|$amplicon_size| for the size.");
added
removed
Bio/Seq/PrimedSeq.pm
