File Coverage

blib/lib/Bio/GeneDesign/CodonJuggle.pm
Criterion Covered Total %
statement 117 118 99.1
branch 9 12 75.0
condition 3 3 100.0
subroutine 11 11 100.0
pod n/a
total 140 144 97.2


line stmt bran cond sub pod time code
1             =head1 NAME
2              
3             Bio::GeneDesign::CodonJuggle
4              
5             =head1 VERSION
6              
7             Version 5.52
8              
9             =head1 DESCRIPTION
10              
11             Codon Juggle a sequence using rscu data to determine replacement likelihood
12              
13             =head1 AUTHOR
14              
15             Sarah Richardson <SMRichardson@lbl.gov>.
16              
17             =cut
18              
19             package Bio::GeneDesign::CodonJuggle;
20             require Exporter;
21              
22 11     11   6169 use Bio::GeneDesign::Random qw(_random_index _weighted_rand);
  11         36  
  11         819  
23 11     11   62 use Bio::GeneDesign::Basic qw(_compare_sequences);
  11         24  
  11         479  
24 11     11   62 use Carp;
  11         17  
  11         591  
25              
26 11     11   54 use strict;
  11         21  
  11         340  
27 11     11   64 use warnings;
  11         19  
  11         454  
28              
29             our $VERSION = 5.52;
30              
31 11     11   56 use base qw(Exporter);
  11         29  
  11         18369  
32             our @EXPORT_OK = qw(
33             _codonJuggle_balanced
34             _codonJuggle_high
35             _codonJuggle_least_different_rscu
36             _codonJuggle_most_different_sequence
37             _codonJuggle_random
38             );
39             our %EXPORT_TAGS = (GD => \@EXPORT_OK);
40            
41             =head2 _codonJuggle_balanced
42              
43             =cut
44              
45             sub _codonJuggle_balanced
46             {
47 10000     10000   17331 my ($codon_table, $reverse_codon_table, $rscu_table, $nucseq) = @_;
48 10000         19585 $nucseq = uc $nucseq;
49 10000         13818 my %changehsh = ();
50 10000         43310 foreach my $aa (keys %$reverse_codon_table)
51             {
52 210000         326815 $changehsh{$aa} = {};
53 210000         236030 my @codons = @{$reverse_codon_table->{$aa}};
  210000         466410  
54 210000         260593 my $count = scalar(@codons);
55 210000         212556 my $checksum = 0;
56 210000         286613 foreach my $coda (@codons)
57             {
58 640000         826549 my $likely = ($rscu_table->{$coda}) / $count;
59 640000         970510 $changehsh{$aa}->{$coda} = $likely;
60 640000         1025776 $checksum += $likely;
61             }
62 210000 50       610127 if ($checksum == 0)
63             {
64 0         0 croak "This RSCU table has no positive values for $aa\n";
65             }
66             }
67 10000         25416 my $offset = 0;
68 10000         12254 my $newseq = q{};
69 10000         23955 while ($offset < length($nucseq))
70             {
71 50000         69663 my $curcod = substr($nucseq, $offset, 3);
72 50000         75922 my $aa = $codon_table->{$curcod};
73 50000         130352 my $newcod = _weighted_rand($changehsh{$aa});
74 50000         69385 $newseq .= $newcod;
75 50000         121776 $offset += 3;
76             }
77 10000         146062 return $newseq;
78             }
79              
80             =head2 _codonJuggle_high
81              
82            
83             =cut
84              
85             sub _codonJuggle_high
86             {
87 1     1   3 my ($codon_table, $reverse_codon_table, $rscu_table, $nucseq) = @_;
88 1         2 my $cod_highs = {};
89 1         3 foreach my $aa (keys %{$reverse_codon_table})
  1         6  
90             {
91 21         26 my $myrscu = -1;
92 21         23 foreach my $codon (@{$reverse_codon_table->{$aa}})
  21         43  
93             {
94 64 100       191 if ($rscu_table->{$codon} > $myrscu)
95             {
96 36         57 $cod_highs->{$aa} = $codon;
97 36         70 $myrscu = $rscu_table->{$codon};
98             }
99             }
100             }
101 1         4 my $offset = 0;
102 1         2 my $newseq = q{};
103 1         6 while ($offset < length($nucseq))
104             {
105 199         237 my $curcod = substr($nucseq, $offset, 3);
106 199         276 $newseq .= $cod_highs->{$codon_table->{$curcod}};
107 199         370 $offset += 3;
108             }
109 1         11 return $newseq;
110             }
111              
112             =head2 _codonJuggle_least_different_rscu
113              
114             =cut
115              
116             sub _codonJuggle_least_different_rscu
117             {
118 1     1   3 my ($codon_table, $reverse_codon_table, $rscu_table, $nucseq) = @_;
119 1         2 my %changehsh = ();
120 1         8 foreach my $aa (keys %$reverse_codon_table)
121             {
122 21         26 foreach my $coda (@{$reverse_codon_table->{$aa}})
  21         41  
123             {
124 64         190 $changehsh{$coda} = $coda;
125 56         144 my @posarr = sort {abs($rscu_table->{$a} - $rscu_table->{$coda})
  180         453  
126             <=> abs($rscu_table->{$b} - $rscu_table->{$coda})}
127 244         440 grep {abs($rscu_table->{$_} - $rscu_table->{$coda}) <= 1}
128 64         109 grep {$_ ne $coda}
129 64         70 @{$reverse_codon_table->{$aa}};
130 64 100       228 $changehsh{$coda} = $posarr[0] if (scalar @posarr);
131             }
132             }
133              
134 1         4 my $offset = 0;
135 1         2 my $newseq = q{};
136 1         5 while ($offset < length($nucseq))
137             {
138 199         343 my $curcod = substr($nucseq, $offset, 3);
139 199         268 $newseq .= $changehsh{$curcod};
140 199         435 $offset += 3;
141             }
142 1         14 return $newseq;
143             }
144              
145             =head2 _codonJuggle_most_different_sequence
146              
147             =cut
148              
149             sub _codonJuggle_most_different_sequence
150             {
151 1     1   4 my ($codon_table, $reverse_codon_table, $rscu_table, $nucseq) = @_;
152 1         3 my %changehsh = ();
153 1         7 foreach my $aa (keys %$reverse_codon_table)
154             {
155 21         35 foreach my $coda (@{$reverse_codon_table->{$aa}})
  21         63  
156             {
157 64         92 my %hsh;
158 64         92 foreach my $codb (@{$reverse_codon_table->{$aa}})
  64         179  
159             {
160 244         823 $hsh{$codb} = _compare_sequences($coda, $codb);
161             }
162 64 50 100     254 my @mdcod = sort {$hsh{$b}->{D} <=> $hsh{$a}->{D}
  308         1501  
163             || $hsh{$b}->{V} <=> $hsh{$a}->{V}
164             || abs($rscu_table->{$a} - $rscu_table->{$coda})
165             <=> abs($rscu_table->{$b} - $rscu_table->{$coda})}
166             keys %hsh;
167 64 100       184 if (scalar @mdcod > 1)
168             {
169 62 50       172 shift @mdcod if ($mdcod[0] eq $coda);
170             }
171 64         452 $changehsh{$coda} = $mdcod[0];
172             }
173             }
174              
175 1         7 my $offset = 0;
176 1         2 my $newseq = q{};
177 1         6 while ($offset < length($nucseq))
178             {
179 199         221 my $curcod = substr($nucseq, $offset, 3);
180 199         221 $newseq .= $changehsh{$curcod};
181 199         364 $offset += 3;
182             }
183 1         20 return $newseq;
184             }
185              
186              
187             =head2 _codonJuggle_random
188              
189             =cut
190              
191             sub _codonJuggle_random
192             {
193 10000     10000   22425 my ($codon_table, $reverse_codon_table, $rscu_table, $nucseq) = @_;
194 10000         16345 my $cod_highs = {};
195 10000         12365 foreach my $aa (keys %{$reverse_codon_table})
  10000         52691  
196             {
197 210000         398831 $cod_highs->{$aa} = [];
198 210000         230304 foreach my $codon (@{$reverse_codon_table->{$aa}})
  210000         391289  
199             {
200 640000         669147 push @{$cod_highs->{$aa}}, $codon;
  640000         1469108  
201             }
202             }
203 10000         27093 my $offset = 0;
204 10000         25158 my $newseq = q{};
205 10000         23696 while ($offset < length($nucseq))
206             {
207 50000         73933 my $curcod = substr($nucseq, $offset, 3);
208 50000         82950 my $aa = $codon_table->{$curcod};
209 50000         53046 my $index = _random_index(scalar @{$cod_highs->{$aa}});
  50000         164546  
210 50000         100672 $newseq .= $cod_highs->{$aa}->[$index];
211 50000         125045 $offset += 3;
212             }
213 10000         121213 return $newseq;
214             }
215              
216             1;
217              
218             __END__
219              
220             =head1 COPYRIGHT AND LICENSE
221              
222             Copyright (c) 2013, GeneDesign developers
223             All rights reserved.
224              
225             Redistribution and use in source and binary forms, with or without modification,
226             are permitted provided that the following conditions are met:
227              
228             * Redistributions of source code must retain the above copyright notice, this
229             list of conditions and the following disclaimer.
230              
231             * Redistributions in binary form must reproduce the above copyright notice, this
232             list of conditions and the following disclaimer in the documentation and/or
233             other materials provided with the distribution.
234              
235             * The names of Johns Hopkins, the Joint Genome Institute, the Lawrence Berkeley
236             National Laboratory, the Department of Energy, and the GeneDesign developers may
237             not be used to endorse or promote products derived from this software without
238             specific prior written permission.
239              
240             THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
241             ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
242             WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
243             DISCLAIMED. IN NO EVENT SHALL THE DEVELOPERS BE LIABLE FOR ANY DIRECT, INDIRECT,
244             INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
245             LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
246             PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
247             LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
248             OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
249             ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
250              
251             =cut