A test to ensure Bio::PrimarySeqI->trunc() doesn't use clone() for a Bio::Seq::RichSe...
[bioperl-live.git] / Bio / Tools / Gel.pm
blobed4d3eeec54e5a2a7a12e18c696aa2302c8dca9e
1 #
2 # BioPerl module for Bio::Tools::Gel
3 # Copyright Allen Day <allenday@ucla.edu>
4 # You may distribute this module under the same terms as perl itself
6 # POD documentation - main docs before the code
8 =head1 NAME
10 Bio::Tools::Gel - Calculates relative electrophoretic migration distances
12 =head1 SYNOPSIS
14 use Bio::PrimarySeq;
15 use Bio::Restriction::Analysis;
16 use Bio::Tools::Gel;
18 # get a sequence
19 my $d = 'AAAAAAAAAGAATTCTTTTTTTTTTTTTTGAATTCGGGGGGGGGGGGGGGGGGGG';
20 my $seq1 = Bio::Seq->new(-id=>'groundhog day',-seq=>$d);
22 # cut it with an enzyme
23 my $ra=Bio::Restriction::Analysis->new(-seq=>$seq1);
24 @cuts = $ra->fragments('EcoRI'), 3;
26 # analyse the fragments in a gel
27 my $gel = Bio::Tools::Gel->new(-seq=>\@cuts,-dilate=>10);
28 my %bands = $gel->bands;
29 foreach my $band (sort {$b <=> $a} keys %bands){
30 print $band,"\t", sprintf("%.1f", $bands{$band}),"\n";
33 #prints:
34 #20 27.0
35 #25 26.0
36 #10 30.0
38 =head1 DESCRIPTION
40 This takes a set of sequences or Bio::Seq objects, and calculates their
41 respective migration distances using:
42 distance = dilation * (4 - log10(length(dna));
44 Source: Molecular Cloning, a Laboratory Manual. Sambrook, Fritsch, Maniatis.
45 CSHL Press, 1989.
47 Bio::Tools::Gel currently calculates migration distances based solely on
48 the length of the nucleotide sequence. Secondary or tertiary structure,
49 curvature, and other biophysical attributes of a sequence are currently
50 not considered. Polypeptide migration is currently not supported.
52 =head1 FEEDBACK
54 =head2 Mailing Lists
56 User feedback is an integral part of the evolution of this and other
57 Bioperl modules. Send your comments and suggestions preferably to
58 the Bioperl mailing list. Your participation is much appreciated.
60 bioperl-l@bioperl.org - General discussion
61 http://bioperl.org/wiki/Mailing_lists - About the mailing lists
63 =head2 Support
65 Please direct usage questions or support issues to the mailing list:
67 I<bioperl-l@bioperl.org>
69 rather than to the module maintainer directly. Many experienced and
70 reponsive experts will be able look at the problem and quickly
71 address it. Please include a thorough description of the problem
72 with code and data examples if at all possible.
74 =head2 Reporting Bugs
76 Report bugs to the Bioperl bug tracking system to help us keep track
77 of the bugs and their resolution. Bug reports can be submitted via the
78 web:
80 https://github.com/bioperl/bioperl-live/issues
82 =head1 AUTHOR - Allen Day
84 Email allenday@ucla.edu
86 =head1 APPENDIX
88 The rest of the documentation details each of the object methods.
89 Internal methods are usually preceded with a _
91 =cut
94 package Bio::Tools::Gel;
95 use strict;
97 use Bio::PrimarySeq;
99 use base qw(Bio::Root::Root);
101 =head2 new
103 Title : new
104 Usage : my $gel = Bio::Tools::Gel->new(-seq => $sequence,-dilate => 3);
105 Function: Initializes a new Gel
106 Returns : Bio::Tools::Gel
107 Args : -seq => Bio::Seq(s), scalar(s) or list of either/both
108 (default: none)
109 -dilate => Expand band migration distances (default: 1)
111 =cut
113 sub new {
114 my ($class, @args) = @_;
116 my $self = $class->SUPER::new(@args);
117 my ($seqs, $dilate) = $self->_rearrange([qw(SEQ DILATE)], @args);
118 if( ! ref($seqs) ) {
119 $self->add_band([$seqs]);
120 } elsif( ref($seqs) =~ /array/i ||
121 $seqs->isa('Bio::PrimarySeqI') ) {
122 $self->add_band($seqs);
124 $self->dilate($dilate || 1);
126 return $self;
130 =head2 add_band
132 Title : add_band
133 Usage : $gel->add_band($seq);
134 Function: Calls _add_band with a (possibly created) Bio::Seq object.
135 Returns :
136 Args : Bio::Seq, scalar sequence, or list of either/both.
138 =cut
140 sub add_band {
141 my ($self, $args) = @_;
143 foreach my $arg (@$args){
144 my $seq;
145 if( ! ref $arg ) {
146 if( $arg =~ /^\d+/ ) {
147 # $arg is a number
148 $seq = Bio::PrimarySeq->new(-seq=>'N'x$arg, -id => $arg);
149 } else {
150 # $arg is a sequence string
151 $seq = Bio::PrimarySeq->new(-seq=>$arg, -id=>length $arg);
153 } elsif( $arg->isa('Bio::PrimarySeqI') ) {
154 # $arg is a sequence object
155 $seq = $arg;
158 $self->_add_band($seq);
160 return 1;
164 =head2 _add_band
166 Title : _add_band
167 Usage : $gel->_add_band($seq);
168 Function: Adds a new band to the gel.
169 Returns :
170 Args : Bio::Seq object
172 =cut
174 sub _add_band {
175 my ($self, $arg) = @_;
176 if ( defined $arg) {
177 push (@{$self->{'bands'}},$arg);
179 return 1;
183 =head2 dilate
185 Title : dilate
186 Usage : $gel->dilate(1);
187 Function: Sets/retrieves the dilation factor.
188 Returns : dilation factor
189 Args : Float or none
191 =cut
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};
202 sub migrate {
203 my ($self, $arg) = @_;
204 $arg = $self unless $arg;
205 if ( $arg ) {
206 return 4 - log10($arg);
207 } else {
208 return 0;
213 =head2 bands
215 Title : bands
216 Usage : $gel->bands;
217 Function: Calculates migration distances of sequences.
218 Returns : hash of (seq_id => distance)
219 Args :
221 =cut
223 sub bands {
224 my $self = shift;
225 $self->throw("bands() is read-only") if @_;
227 my %bands = ();
229 foreach my $band (@{$self->{bands}}){
230 my $distance = $self->dilate * migrate($band->length);
231 $bands{$band->id} = $distance;
234 return %bands;
238 =head2 log10
240 Title : log10
241 Usage : log10($n);
242 Function: returns base 10 log of $n.
243 Returns : float
244 Args : float
246 =cut
248 # from "Programming Perl"
249 sub log10 {
250 my $n = shift;
251 return log($n)/log(10);