Bio/Tools/Gel.pm

   1 # $Id$
   2 #
   3 # BioPerl module for Bio::Tools::Gel
   4 # Copyright Allen Day <allenday@ucla.edu>
   5 # You may distribute this module under the same terms as perl itself
   6
   7 # POD documentation - main docs before the code
   8
   9 =head1 NAME
  10
  11 Bio::Tools::Gel - Calculates relative electrophoretic migration distances
  12
  13 =head1 SYNOPSIS
  14
  15     use Bio::PrimarySeq;
  16     use Bio::Tools::RestrictionAnalysis;
  17     use Bio::Tools::Gel;
  18
  19     # get a sequence
  20     my $d = 'AAAAAAAAAGAATTCTTTTTTTTTTTTTTGAATTCGGGGGGGGGGGGGGGGGGGG';
  21     my $seq1 = Bio::Seq->new(-id=>'groundhog day',-seq=>$d);
  22
  23     # cut it with an enzyme
  24     my $ra=Bio::Restriction::Analysis->new(-seq=>$seq1);
  25     @cuts = $ra->fragments('EcoRI'), 3;
  26
  27     # analyse the fragments in a gel
  28     my $gel = Bio::Tools::Gel->new(-seq=>\@cuts,-dilate=>10);
  29     my %bands = $gel->bands;
  30     foreach my $band (sort {$b <=> $a} keys %bands){
  31       print $band,"\t", sprintf("%.1f", $bands{$band}),"\n";
  32     }
  33
  34     #prints:
  35     #20   27.0
  36     #25   26.0
  37     #10   30.0
  38
  39
  40 =head1 DESCRIPTION
  41
  42 This takes a set of sequences or Bio::Seq objects, and calculates their
  43 respective migration distances using:
  44     distance = dilation * (4 - log10(length(dna));
  45
  46 Source: Molecular Cloning, a Laboratory Manual. Sambrook, Fritsch, Maniatis.
  47 CSHL Press, 1989.
  48
  49 Bio::Tools::Gel currently calculates migration distances based solely on
  50 the length of the nucleotide sequence.  Secondary or tertiary structure,
  51 curvature, and other biophysical attributes of a sequence are currently
  52 not considered.  Polypeptide migration is currently not supported.
  53
  54 =head1 FEEDBACK
  55
  56 =head2 Mailing Lists
  57
  58 User feedback is an integral part of the evolution of this and other
  59 Bioperl modules. Send your comments and suggestions preferably to
  60 the Bioperl mailing list.  Your participation is much appreciated.
  61
  62   bioperl-l@bioperl.org                  - General discussion
  63   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
  64
  65 =head2 Reporting Bugs
  66
  67 Report bugs to the Bioperl bug tracking system to help us keep track
  68 of the bugs and their resolution. Bug reports can be submitted via the
  69 web:
  70
  71   http://bugzilla.open-bio.org/
  72
  73 =head1 AUTHOR - Allen Day
  74
  75 Email allenday@ucla.edu
  76
  77 =head1 APPENDIX
  78
  79 The rest of the documentation details each of the object methods.
  80 Internal methods are usually preceded with a _
  81
  82 =cut
  83
  84
  85 # Let the code begin...
  86
  87
  88 package Bio::Tools::Gel;
  89 use strict;
  90
  91 use Bio::PrimarySeq;
  92
  93 use base qw(Bio::Root::Root);
  94
  95 =head2 new
  96
  97  Title   : new
  98  Usage   : my $gel = Bio::Tools::Gel->new(-seq => $sequence,-dilate => 3);
  99  Function: Initializes a new Gel
 100  Returns : Bio::Tools::Gel
 101  Args    : -seq      => Bio::Seq(s), scalar(s) or list of either/both
 102                         (default: none)
 103            -dilate   => Expand band migration distances (default: 1)
 104
 105 =cut
 106
 107 sub new {
 108   my($class,@args) = @_;
 109
 110   my $self = $class->SUPER::new(@args);
 111   my ($seqs,$dilate) = $self->_rearrange([qw(SEQ DILATE)],
 112                                           @args);
 113   if( ! ref($seqs)  ) {
 114       $self->add_band([$seqs]);
 115   } elsif( ref($seqs) =~ /array/i ||
 116            $seqs->isa('Bio::PrimarySeqI') ) {
 117       $self->add_band($seqs);
 118   }
 119   $self->dilate($dilate || 1);
 120
 121   return $self;
 122 }
 123
 124
 125 =head2 add_band
 126
 127  Title   : add_band
 128  Usage   : $gel->add_band($seq);
 129  Function: Calls _add_band with a (possibly created) Bio::Seq object.
 130  Returns :
 131  Args    : Bio::Seq, scalar sequence, or list of either/both.
 132
 133 =cut
 134
 135 sub add_band {
 136   my($self,$args) = @_;
 137
 138   foreach my $arg (@$args){
 139       my $seq;
 140       if( ! ref($arg) ) {
 141           if( $arg =~ /^\d+/ ) {
 142               $seq= Bio::PrimarySeq->new(-seq=>"N"x$arg, -id => $arg);
 143           } else {
 144               $seq= Bio::PrimarySeq->new(-seq=>$arg,-id=>length($arg));
 145           }
 146       } elsif( $arg->isa('Bio::PrimarySeqI') ) {
 147           $seq = $arg;
 148       }
 149
 150     $seq->validate_seq or $seq->throw("invalid symbol in sequence".$seq->seq()."\n");
 151     $self->_add_band($seq);
 152   }
 153 }
 154
 155 =head2 _add_band
 156
 157  Title   : _add_band
 158  Usage   : $gel->_add_band($seq);
 159  Function: Adds a new band to the gel.
 160  Returns :
 161  Args    : Bio::Seq object
 162
 163 =cut
 164
 165 sub _add_band {
 166   my($self,$arg) = @_;
 167   if( defined $arg) {
 168       push (@{$self->{'bands'}},$arg);
 169   }
 170 }
 171
 172 =head2 dilate
 173
 174  Title   : dilate
 175  Usage   : $gel->dilate(1);
 176  Function: Sets/retrieves the dilation factor.
 177  Returns : dilation factor
 178  Args    : Float or none
 179
 180 =cut
 181
 182 sub dilate {
 183   my($self,$arg) = @_;
 184   return $self->{dilate} unless $arg;
 185   $self->throw("-dilate should be numeric") if defined $arg and $arg =~ /[^e\d\.]/;
 186   $self->{dilate} = $arg;
 187   return $self->{dilate};
 188 }
 189
 190 sub migrate {
 191   my ($self,$arg) = @_;
 192   $arg = $self unless $arg;
 193   if ( $arg ) {
 194       return 4 - log10($arg);
 195   } else { return 0; }
 196 }
 197
 198 =head2 bands
 199
 200  Title   : bands
 201  Usage   : $gel->bands;
 202  Function: Calculates migration distances of sequences.
 203  Returns : hash of (seq_id => distance)
 204  Args    :
 205
 206 =cut
 207
 208 sub bands {
 209   my $self = shift;
 210   $self->throw("bands() is read-only") if @_;
 211
 212   my %bands = ();
 213
 214   foreach my $band (@{$self->{bands}}){
 215       my $distance = $self->dilate * migrate($band->length);
 216       $bands{$band->id} = $distance;
 217   }
 218
 219   return %bands;
 220 }
 221
 222 =head2 log10
 223
 224  Title   : log10
 225  Usage   : log10($n);
 226  Function: returns base 10 log of $n.
 227  Returns : float
 228  Args    : float
 229
 230 =cut
 231
 232 #from programming perl
 233 sub log10 {
 234     my $n = shift;
 235     return log($n)/log(10);
 236 }
 237
 238 1;