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
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
36
# For example, start with a fasta file
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:
38
# For example, start with a fasta file
39
40
use Bio::Tools::Run::Primer3;
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;
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
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');
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);
56
#(ii) to create a genbank file for a sequence and its cognate primers
59
Or, (ii), to create a genbank file for a sequence and its cognate primers:
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
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);
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
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);
82
# This should output a genbank file with the two primers labeled.
83
This should output a genbank file with the two primers labeled.
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.
91
92
The primers are represented as Bio::SeqFeature::Primer objects, and should
94
The simplest way to initiate a PrimedSeq object is as follows:
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.
95
A simple way to create a PrimedSeq object is as follows:
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,
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.
106
107
This is based on the PrimedSeq.pm module started by Chad Matsalla, with
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(
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
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)
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.
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.
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.
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
206
foreach my $key (keys %args) {
207
if ($key eq "-seq" || $key eq "-SEQ") {
208
$self->{target_sequence} = $args{$key};
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;
199
foreach my $key (keys %args) {
200
if ($key =~ /^-seq/i) {
201
$self->{target_sequence} = $args{$key};
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;
222
$self->{$okey} = $args{$key};
223
push @{$self->{'arguments'}},$okey;
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.");
215
$self->{$okey} = $args{$key};
216
push @{$self->{'arguments'}},$okey;
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.");
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.");
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.");
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.");
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.");
242
# now we have the sequences, lets find out where they are
243
$self->_place_seqs();
236
# now we have the sequences, lets find out where they are
237
$self->_place_seqs();
248
242
=head2 get_primer
250
244
Title : get_primer();
251
Usage : $primer = $primedseq->get_primer(l, left, left_primer, -left_primer); to return the left primer
253
$primer = $primedseq->get_primer(r, right, right_primer, -right_primer); to return the right primer
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,
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
256
Either of (r, right, right_primer, -right_primer) to get
258
Either of (b, both, both_primers, -both_primers) to get
260
Note that this is plural. [default]
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'}
345
349
sub _place_seqs {
348
# we are going to pull out the target sequence, and then the primer sequences
349
my $target_sequence = $self->{'target_sequence'}->seq();
352
my $left_seq = $self->{'left_primer'}->seq()->seq();
354
my $rprc = $self->{'right_primer'}->seq()->revcom();
356
my $right_seq=$rprc->seq();
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);
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!")}
369
my ($before, $middle, $after) = ($1, $2, $3); # note didn't use $`, $', and $& because they are bad. Just use length instead.
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)
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);
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.");
386
$self->{'left_primer'}->{'PRIMER_LEFT'}=$left_location;
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.");
397
$self->{'right_primer'}->{'PRIMER_RIGHT'}=$right_location;
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.");
408
$self->{'PRIMER_PRODUCT_SIZE'} = $amplicon_size;
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
413
$self->_set_seqfeature;
352
# we are going to pull out the target sequence, and then the primer sequences
353
my $target_sequence = $self->{'target_sequence'}->seq();
356
my $left_seq = $self->{'left_primer'}->seq()->seq();
358
my $rprc = $self->{'right_primer'}->seq()->revcom();
360
my $right_seq=$rprc->seq();
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);
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!")}
372
my ($before, $middle, $after) = ($1, $2, $3); # note didn't use $`, $', and $& because they are bad. Just use length instead.
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)
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);
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.");
389
$self->{'left_primer'}->{'PRIMER_LEFT'}=$left_location;
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.");
400
$self->{'right_primer'}->{'PRIMER_RIGHT'}=$right_location;
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.");
411
$self->{'PRIMER_PRODUCT_SIZE'} = $amplicon_size;
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
416
$self->_set_seqfeature;
416
419
=head2 _set_seqfeature
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
421
425
Returns : Nothing
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
429
435
sub _set_seqfeature {
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");
435
my ($start, $length) = split /,/, $self->{'left_primer'}->{'PRIMER_LEFT'};
436
my $tm=$self->{'left_primer'}->{'PRIMER_LEFT_TM'} || $self->{'left_primer'}->Tm || 0;
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}
444
($start, $length) = split /,/, $self->{'right_primer'}->{'PRIMER_RIGHT'};
445
$tm=$self->{'right_primer'}->{'PRIMER_RIGHT_TM'} || $self->{'right_primer'}->Tm || 0;
447
my $seqfeatureR=new Bio::SeqFeature::Generic(
448
-start => $start-$length, -end => $start, -strand => -1,
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");
442
my ($start, $length) = split /,/, $self->{'left_primer'}->{'PRIMER_LEFT'};
443
my $tm=$self->{'left_primer'}->{'PRIMER_LEFT_TM'} || $self->{'left_primer'}->Tm || 0;
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}
451
($start, $length) = split /,/, $self->{'right_primer'}->{'PRIMER_RIGHT'};
452
$tm=$self->{'right_primer'}->{'PRIMER_RIGHT_TM'} || $self->{'right_primer'}->Tm || 0;
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}
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);