Bio/Seq/PrimedSeq.pm

   1 # $Id$
   2 #
   3 # This is the original copyright statement. I have relied on Chad's module
   4 # extensively for this module.
   5 #
   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.
   9 #
  10 # Copyright Chad Matsalla
  11 #
  12 # You may distribute this module under the same terms as perl itself
  13 # POD documentation - main docs before the code
  14 #
  15 # But I have modified lots of it, so I guess I should add:
  16 #
  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.
  20 #
  21 # Copyright Rob Edwards
  22 #
  23 # You may distribute this module under the same terms as perl itself
  24 # POD documentation - main docs before the code
  25
  26
  27 =head1 NAME
  28
  29 Bio::Seq::PrimedSeq - A representation of a sequence and two primers
  30 flanking a target region
  31
  32 =head1 SYNOPSIS
  33
  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:
  37
  38       # For example, start with a fasta file
  39       use Bio::SeqIO;
  40       use Bio::Tools::Run::Primer3;
  41
  42       my $file = shift || die "need a file to read";
  43       my $seqin = Bio::SeqIO->new(-file => $file);
  44       my $seq = $seqin->next_seq;
  45
  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
  49
  50       # create a file to write the results to
  51       my $seqout = Bio::SeqIO->new(-file => ">primed_sequence.gbk",
  52                                    -format => 'genbank');
  53
  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);
  57       }
  58
  59   # Or, (ii), to create a genbank file for a sequence and its cognate
  60   # primers:
  61
  62      use Bio::SeqIO;
  63      use Bio::Seq::PrimedSeq;
  64
  65      # have a sequence file ($file) with the template, and two primers
  66      # that match it, in fasta format
  67
  68      my $file = shift || die "$0 <file>";
  69      my $seqin = Bio::SeqIO->new(-file => $file);
  70
  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);
  83
  84   # This should output a genbank file with the two primers labeled.
  85
  86 =head1 DESCRIPTION
  87
  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.
  92
  93 The primers are represented as Bio::SeqFeature::Primer objects, and should
  94 be instantiated first.
  95
  96 A simple way to create a PrimedSeq object is as follows:
  97
  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,
 102   );
 103
 104 From the PrimedSeq object you should be able to retrieve
 105 information about melting temperatures and what not on each of the primers
 106 and the amplicon.
 107
 108 This is based on the PrimedSeq.pm module started by Chad Matsalla, with
 109 additions/improvements by Rob Edwards.
 110
 111 =head1 FEEDBACK
 112
 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.
 116
 117   bioperl-l@bioperl.org                  - General discussion
 118   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
 119
 120 =head2 Support
 121
 122 Please direct usage questions or support issues to the mailing list:
 123
 124 L<bioperl-l@bioperl.org>
 125
 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.
 130
 131 =head2 Reporting Bugs
 132
 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
 135 web:
 136
 137   http://bugzilla.open-bio.org/
 138
 139 =head1 AUTHOR
 140
 141 Rob Edwards, redwards@utmem.edu
 142
 143 Based on a module written by Chad Matsalla, bioinformatics1@dieselwurks.com
 144
 145 =head1 APPENDIX
 146
 147 The rest of the documentation details each of the object
 148 methods. Internal methods are usually preceded with a _
 149
 150 =cut
 151
 152
 153 # Let the code begin...
 154
 155 package Bio::Seq::PrimedSeq;
 156
 157 use strict;
 158 use Bio::SeqFeature::Primer;
 159
 160 use vars qw ($AUTOLOAD @RES %OK_FIELD $ID);
 161
 162 use base qw(Bio::Root::Root Bio::SeqFeature::Generic);
 163
 164 BEGIN {
 165  @RES = qw(); # nothing here yet, not sure what we want!
 166  foreach my $attr (@RES) {$OK_FIELD{$attr}++}
 167 }
 168
 169 $ID = 'Bio::Tools::Analysis::Nucleotide::PrimedSeq';
 170
 171 sub AUTOLOAD {
 172  my $self = shift;
 173  my $attr = $AUTOLOAD;
 174  $attr =~ s/.*:://;
 175  $self->throw("Unallowed parameter: $attr !") unless $OK_FIELD{$attr};
 176  $self->{$attr} = shift if @_;
 177  return $self->{$attr};
 178 }
 179
 180
 181 =head2 new
 182
 183  Title   : new()
 184  Usage   : $primed_sequence = Bio::SeqFeature::Primer->new(
 185                                      -seq => $sequence,
 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)
 193
 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.
 197
 198 =cut
 199
 200 sub new {
 201
 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.
 204         # Very confusing.
 205
 206         my($class,%args) = @_;
 207         my $self = $class->SUPER::new(%args);
 208    # these are the absolute minimum components required to make
 209    # a primedseq
 210
 211         foreach my $key (keys %args) {
 212                 if ($key =~  /^-seq/i) {
 213                         $self->{target_sequence} = $args{$key};
 214                         next;
 215                 } else {
 216                         my $okey;
 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;
 224                                 next;
 225                         }
 226
 227                         $self->{$okey} = $args{$key};
 228                         push @{$self->{'arguments'}},$okey;
 229                 }
 230         }
 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.");
 235         }
 236
 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.");
 240  }
 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.");
 246         }
 247
 248         # now we have the sequences, lets find out where they are
 249         $self->_place_seqs();
 250         return $self;
 251 }
 252
 253
 254 =head2 get_primer
 255
 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,
 262            -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
 267            primer.
 268            Either of (r, right, right_primer, -right_primer) to get
 269            right primer
 270                 Either of (b, both, both_primers, -both_primers) to get
 271            both primers.
 272            Note that this is plural. [default]
 273
 274 =cut
 275
 276 sub get_primer() {
 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'}
 284  }
 285  elsif ($arg =~ /^r/ || $arg =~ /^-r/) {return $self->{'right_primer'}}
 286  elsif ($arg =~ /^b/ || $arg =~ /^-b/) {return ($self->{'left_primer'}, $self->{'right_primer'})}
 287 }
 288
 289
 290
 291 =head2 annotated_sequence
 292
 293  Title   : annotated_sequence
 294  Usage   : $annotated_sequence_object = $primedseq->annotated_sequence()
 295  Function: Get an annotated sequence object containg the left and right
 296            primers
 297  Returns : An annotated sequence object or 0 if not defined.
 298  Args    :
 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.
 301
 302 =cut
 303
 304 sub annotated_sequence {
 305   my $self = shift;
 306   if (exists $self->{annotated_sequence}) {return $self->{annotated_sequence}}
 307   else {return 0}
 308 }
 309
 310 =head2 amplicon
 311
 312  Title   : amplicon
 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
 316  Args    : None
 317  Note    :
 318
 319 =cut
 320
 321 sub amplicon {
 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;
 327
 328  my $seqobj=Bio::Seq->new(-id=>$id, seq=>$self->{'amplicon_sequence'});
 329  return $seqobj;
 330 }
 331
 332
 333 =head2 seq
 334
 335  Title   : seq
 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
 339  Args    : None
 340  Note    :
 341
 342 =cut
 343
 344 sub seq {
 345  my $self = shift;
 346  return $self->{target_sequence};
 347 }
 348
 349 =head2 _place_seqs
 350
 351  Title   : _place_seqs
 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
 355  Returns : Nothing
 356  Args    : None
 357  Note    : Internal use only
 358
 359 =cut
 360
 361 sub _place_seqs {
 362         my $self = shift;
 363
 364         # we are going to pull out the target sequence, and then the primer sequences
 365         my $target_sequence = $self->{'target_sequence'}->seq();
 366
 367         # left primer
 368         my $left_seq = $self->{'left_primer'}->seq()->seq();
 369
 370         my $rprc = $self->{'right_primer'}->seq()->revcom();
 371
 372         my $right_seq=$rprc->seq();
 373
 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);
 378
 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!")}
 382  }
 383
 384         my ($before, $middle, $after) = ($1, $2, $3); # note didn't use $`, $', and $& because they are bad. Just use length instead.
 385
 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)
 388
 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);
 392
 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.");
 398                 }
 399         }
 400         else {
 401                 $self->{'left_primer'}->{'PRIMER_LEFT'}=$left_location;
 402         }
 403
 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.");
 409                 }
 410         }
 411         else {
 412                 $self->{'right_primer'}->{'PRIMER_RIGHT'}=$right_location;
 413         }
 414
 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.");
 420                 }
 421         }
 422         else {
 423                 $self->{'PRIMER_PRODUCT_SIZE'} = $amplicon_size;
 424         }
 425
 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
 427
 428         $self->_set_seqfeature;
 429 }
 430
 431 =head2 _set_seqfeature
 432
 433  Title   : _set_seqfeature
 434  Usage   : $self->_set_seqfeature()
 435  Function: An internal method to create Bio::SeqFeature::Generic objects
 436            for the primed seq
 437  Returns : Nothing
 438  Args    : None
 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
 442            with write_seq.
 443
 444 =cut
 445
 446
 447 sub _set_seqfeature {
 448         my $self = shift;
 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");
 452                 return 0;
 453         }
 454         my ($start, $length) = split /,/, $self->{'left_primer'}->{'PRIMER_LEFT'};
 455         my $tm=$self->{'left_primer'}->{'PRIMER_LEFT_TM'} || $self->{'left_primer'}->Tm || 0;
 456
 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}
 461                                                                                                                           );
 462
 463         ($start, $length) = split /,/, $self->{'right_primer'}->{'PRIMER_RIGHT'};
 464         $tm=$self->{'right_primer'}->{'PRIMER_RIGHT_TM'} || $self->{'right_primer'}->Tm || 0;
 465
 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}
 470                                                                                                                           );
 471
 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);
 476 }
 477
 478 1;