branch-1-6/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::Restriction::Analysis;
  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 Support
  66
  67 Please direct usage questions or support issues to the mailing list:
  68
  69 I<bioperl-l@bioperl.org>
  70
  71 rather than to the module maintainer directly. Many experienced and
  72 reponsive experts will be able look at the problem and quickly
  73 address it. Please include a thorough description of the problem
  74 with code and data examples if at all possible.
  75
  76 =head2 Reporting Bugs
  77
  78 Report bugs to the Bioperl bug tracking system to help us keep track
  79 of the bugs and their resolution. Bug reports can be submitted via the
  80 web:
  81
  82   http://bugzilla.open-bio.org/
  83
  84 =head1 AUTHOR - Allen Day
  85
  86 Email allenday@ucla.edu
  87
  88 =head1 APPENDIX
  89
  90 The rest of the documentation details each of the object methods.
  91 Internal methods are usually preceded with a _
  92
  93 =cut
  94
  95
  96 # Let the code begin...
  97
  98
  99 package Bio::Tools::Gel;
 100 use strict;
 101
 102 use Bio::PrimarySeq;
 103
 104 use base qw(Bio::Root::Root);
 105
 106 =head2 new
 107
 108  Title   : new
 109  Usage   : my $gel = Bio::Tools::Gel->new(-seq => $sequence,-dilate => 3);
 110  Function: Initializes a new Gel
 111  Returns : Bio::Tools::Gel
 112  Args    : -seq      => Bio::Seq(s), scalar(s) or list of either/both
 113                         (default: none)
 114            -dilate   => Expand band migration distances (default: 1)
 115
 116 =cut
 117
 118 sub new {
 119   my($class,@args) = @_;
 120
 121   my $self = $class->SUPER::new(@args);
 122   my ($seqs,$dilate) = $self->_rearrange([qw(SEQ DILATE)],
 123                                           @args);
 124   if( ! ref($seqs)  ) {
 125       $self->add_band([$seqs]);
 126   } elsif( ref($seqs) =~ /array/i ||
 127            $seqs->isa('Bio::PrimarySeqI') ) {
 128       $self->add_band($seqs);
 129   }
 130   $self->dilate($dilate || 1);
 131
 132   return $self;
 133 }
 134
 135
 136 =head2 add_band
 137
 138  Title   : add_band
 139  Usage   : $gel->add_band($seq);
 140  Function: Calls _add_band with a (possibly created) Bio::Seq object.
 141  Returns :
 142  Args    : Bio::Seq, scalar sequence, or list of either/both.
 143
 144 =cut
 145
 146 sub add_band {
 147   my($self,$args) = @_;
 148
 149   foreach my $arg (@$args){
 150       my $seq;
 151       if( ! ref($arg) ) {
 152           if( $arg =~ /^\d+/ ) {
 153               $seq= Bio::PrimarySeq->new(-seq=>"N"x$arg, -id => $arg);
 154           } else {
 155               $seq= Bio::PrimarySeq->new(-seq=>$arg,-id=>length($arg));
 156           }
 157       } elsif( $arg->isa('Bio::PrimarySeqI') ) {
 158           $seq = $arg;
 159       }
 160
 161     $seq->validate_seq or $seq->throw("invalid symbol in sequence".$seq->seq()."\n");
 162     $self->_add_band($seq);
 163   }
 164 }
 165
 166 =head2 _add_band
 167
 168  Title   : _add_band
 169  Usage   : $gel->_add_band($seq);
 170  Function: Adds a new band to the gel.
 171  Returns :
 172  Args    : Bio::Seq object
 173
 174 =cut
 175
 176 sub _add_band {
 177   my($self,$arg) = @_;
 178   if( defined $arg) {
 179       push (@{$self->{'bands'}},$arg);
 180   }
 181 }
 182
 183 =head2 dilate
 184
 185  Title   : dilate
 186  Usage   : $gel->dilate(1);
 187  Function: Sets/retrieves the dilation factor.
 188  Returns : dilation factor
 189  Args    : Float or none
 190
 191 =cut
 192
 193 sub dilate {
 194   my($self,$arg) = @_;
 195   return $self->{dilate} unless $arg;
 196   $self->throw("-dilate should be numeric") if defined $arg and $arg =~ /[^e\d\.]/;
 197   $self->{dilate} = $arg;
 198   return $self->{dilate};
 199 }
 200
 201 sub migrate {
 202   my ($self,$arg) = @_;
 203   $arg = $self unless $arg;
 204   if ( $arg ) {
 205       return 4 - log10($arg);
 206   } else { return 0; }
 207 }
 208
 209 =head2 bands
 210
 211  Title   : bands
 212  Usage   : $gel->bands;
 213  Function: Calculates migration distances of sequences.
 214  Returns : hash of (seq_id => distance)
 215  Args    :
 216
 217 =cut
 218
 219 sub bands {
 220   my $self = shift;
 221   $self->throw("bands() is read-only") if @_;
 222
 223   my %bands = ();
 224
 225   foreach my $band (@{$self->{bands}}){
 226       my $distance = $self->dilate * migrate($band->length);
 227       $bands{$band->id} = $distance;
 228   }
 229
 230   return %bands;
 231 }
 232
 233 =head2 log10
 234
 235  Title   : log10
 236  Usage   : log10($n);
 237  Function: returns base 10 log of $n.
 238  Returns : float
 239  Args    : float
 240
 241 =cut
 242
 243 #from programming perl
 244 sub log10 {
 245     my $n = shift;
 246     return log($n)/log(10);
 247 }
 248
 249 1;