← Back to branch summary

~ubuntu-branches/ubuntu/lucid/bioperl/lucid

« back to all changes in this revision

Viewing changes to Bio/CodonUsage/IO.pm

  • Committer: Bazaar Package Importer
  • Author(s): Matt Hope
  • Date: 2004-04-18 14:24:11 UTC
  • mfrom: (1.2.1 upstream) (2.1.1 warty)
  • Revision ID: james.westby@ubuntu.com-20040418142411-gr92uexquw4w8liq
Tags: 1.4-1
* New upstream release
* Examples and working code are installed by default to usr/bin,
  this has been moved to usr/share/doc/bioperl/bin

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/CodonUsage/IO.pm
 
1
# $Id: IO.pm,v 1.4 2003/09/08 13:05:00 radams Exp $
 
2
#
 
3
# BioPerl module for Bio::CodonUsage::IO
 
4
#
 
5
# Cared for by Richard Adams (richard.adams@ed.ac.uk)
 
6
#
 
7
# Copyright Richard Adams
 
8
#
 
9
# You may distribute this module under the same terms as perl itself
 
10
# POD documentation - main docs before the code
 
11
 
 
12
=head1 NAME
 
13
 
 
14
Bio::CodonUsage::IO - for reading and writing codon usage tables to file
 
15
 
 
16
=head1 SYNOPSIS
 
17
 
 
18
        use Bio::CodonUsage::IO;
 
19
 
 
20
        ##read in a codon usage file
 
21
        my $io = Bio::CodonUsage::IO->new(-file => "in");
 
22
        my $cut = $io->next_data();
 
23
 
 
24
        ##write it out again
 
25
        my $out = Bio::CodonUsage::IO->new(-file => ">out");
 
26
        $out->write_data($cut);
 
27
 
 
28
=head1 DESCRIPTION
 
29
 
 
30
This class provides standard IO methods for reading and writing text files
 
31
of codon usage tables. These tables can initially be retrieved using
 
32
Bio::DB::CUTG. At present only this format is supported for read/write. 
 
33
 Reading a CUTG will return a Bio::CodonUsage::Table object. 
 
34
 
 
35
=head1 SEE ALSO
 
36
 
 
37
L<Bio::Tools::CodonTable>, 
 
38
L<Bio::WebAgent>,
 
39
L<Bio::CodonUsage::Table>,
 
40
L<Bio::CodonUsage::IO>
 
41
 
 
42
=head1 FEEDBACK
 
43
 
 
44
=head2 Mailing Lists
 
45
 
 
46
 
 
47
User feedback is an integral part of the evolution of this and other
 
48
Bioperl modules. Send your comments and suggestions preferably to one
 
49
of the Bioperl mailing lists.  Your participation is much appreciated.
 
50
 
 
51
  bioperl-l@bioperl.org                       - General discussion
 
52
  http://bio.perl.org/MailList.html           - About the mailing lists
 
53
 
 
54
=head2 Reporting Bugs
 
55
 
 
56
Report bugs to the Bioperl bug tracking system to help us keep track
 
57
the bugs and their resolution.  Bug reports can be submitted via email
 
58
or the web:
 
59
 
 
60
  bioperl-bugs@bio.perl.org
 
61
  http://bugzilla.bioperl.org/
 
62
 
 
63
=head1 AUTHORS
 
64
 
 
65
Richard Adams, Richard.Adams@ed.ac.uk
 
66
 
 
67
=head1 APPENDIX
 
68
 
 
69
The rest of the documentation details each of the object
 
70
methods. Internal methods are usually preceded with a _
 
71
 
 
72
=cut
 
73
 
 
74
 
 
75
# Let the code begin
 
76
 
 
77
package Bio::CodonUsage::IO;
 
78
use Bio::Root::IO;
 
79
use vars qw(@ISA);
 
80
 
 
81
@ISA = qw(Bio::Root::IO);
 
82
 
 
83
=head2  new
 
84
 
 
85
 Title  : new
 
86
 Usage  : my $io = Bio::CodonUsage::IO->new(-file => "CUTfile");
 
87
 Purpose: To  read/write a Bio:CodonUsage::Table object  
 
88
 Returns: A  Bio:CodonUsage::IO object
 
89
 Args   : a file or file handle 
 
90
 
 
91
=cut
 
92
 
 
93
sub new  {
 
94
        my ($class , @args) = @_;
 
95
        my $self = $class->SUPER::new(@args);
 
96
}
 
97
 
 
98
 
 
99
=head2  next_data
 
100
 
 
101
 Title  : next_data
 
102
 Usage  : my $cut = $io->next_data();
 
103
 Purpose: To  obtain a Bio:CodonUsage::Table object 
 
104
 Returns: A  Bio:CodonUsage::Table object
 
105
 Args   : none
 
106
 
 
107
=cut
 
108
 
 
109
sub next_data {
 
110
        my $self = shift;
 
111
        my $cut = $self->_parse;
 
112
        return $cut;
 
113
        }
 
114
 
 
115
=head2  write_data
 
116
 
 
117
 Title  : write_data
 
118
 Usage  : $io->write_data($cut);
 
119
 Purpose: To  write a CUT to file
 
120
 Returns: void
 
121
 Args   : a Bio::CodonUsage::Table object reference 
 
122
 
 
123
=cut
 
124
 
 
125
 
 
126
sub write_data {
 
127
        my ($self, $cut) = @_;
 
128
        if (!$cut || ! $cut->isa(Bio::CodonUsage::Table)) {
 
129
                $self->throw("must supply a Bio::CodonUsage::Table object for writing\n");
 
130
                        }
 
131
        my $outstring = "Codon usage table\n\n";
 
132
 
 
133
        my $sp_string = $cut->species . "[" . $cut->_gb_db . "]  " .
 
134
                                        $cut->cds_count . "  CDS's\n\n";
 
135
        $outstring .= $sp_string;
 
136
        my $colhead = sprintf("%-9s%-9s%15s%12s%12s\n\n", "AmAcid",
 
137
                                                         "Codon", "Number", "/1000", "Fraction");
 
138
        $outstring .= $colhead;
 
139
 
 
140
        ### now write bulk of codon data  ##
 
141
        my $ctable =  Bio::Tools::CodonTable->new;
 
142
 
 
143
        for my $f (qw(G A T C)) {
 
144
                for my $s (qw(G A T C)) {
 
145
                        for my $t (qw(G A T C)) {
 
146
                                $cod = $f . $s . $t;
 
147
                                my $aa =$Bio::SeqUtils::THREECODE {$ctable->translate($cod)};
 
148
                                my $codstr = sprintf("%-9s%-9s%15.2f%12.2f%12.2f\n",            
 
149
 
 
150
                                                $aa, $cod, my $tt = $cut->codon_count($cod)|| 0.00, 
 
151
                                                my $ll =$cut->{'_table'}{$aa}{$cod}{'per1000'}|| 0.00,
 
152
                                                my $ss = $cut->codon_rel_frequency($cod) || 0.00);
 
153
                                $outstring .= $codstr;
 
154
                        }
 
155
                $outstring .= "\n";
 
156
                }
 
157
        }
 
158
        $outstring .= "\n\n";
 
159
 
 
160
        ## now append GC data
 
161
        $outstring .= "Coding GC ". $cut->get_coding_gc('all'). "%\n";
 
162
        $outstring .= "1st letter GC ". $cut->get_coding_gc(1). "%\n";
 
163
        $outstring .= "2nd letter GC ". $cut->get_coding_gc(2). "%\n";
 
164
        $outstring .= "3rd letter GC ". $cut->get_coding_gc(3). "%\n";
 
165
        $outstring .= "Genetic code " . $cut->genetic_code() ."\n\n\n";
 
166
 
 
167
        
 
168
                        
 
169
$self->_print ($outstring);
 
170
$self->flush();
 
171
 
 
172
}
 
173
sub _parse {
 
174
        my $self = shift;
 
175
        my $cdtableobj = Bio::CodonUsage::Table->new();
 
176
        while (my $line = $self->_readline() ) {
 
177
                next if $line =~ /^$/ ;
 
178
                $line =~ s/End/Ter/;
 
179
                ## now parse in species name, cds number
 
180
 
 
181
                if ($line =~ /^(.+?)\[(\w+)\].+?(\d+)/) {
 
182
                        $cdtableobj->species($1);
 
183
                        $cdtableobj->{'_gb_db'} = $2;
 
184
                        $cdtableobj->cds_count($3);
 
185
                        }
 
186
 
 
187
                ## now parse in bulk of codon usage table
 
188
                elsif ( $line =~ /^(\w\w\w)\s+(\w+)\s+(\d+\.\d+)
 
189
                                        \s+(\d+\.\d+)\s+(\d+\.\d+)/x){
 
190
                        if (defined ($1)) {
 
191
                                $cdtableobj->{'_table'}{$1}{$2} = {
 
192
                                                'abs_count'=>$3,
 
193
                                                 'per1000'=> $4,
 
194
                                                 'rel_freq'=> $5,
 
195
                                                };
 
196
                                }
 
197
                        }
 
198
 
 
199
                ## now parse in gc data  ####
 
200
                if($line =~ /^Cod.+?(\d\d\.\d\d)/ ){
 
201
                                $cdtableobj->{'_coding_gc'}{'all'} = $1;
 
202
                                        }
 
203
                elsif ($line =~ /^1st.+?(\d\d\.\d\d)/){ 
 
204
                                $cdtableobj->{'_coding_gc'}{'1'} = $1;
 
205
                        }
 
206
                elsif($line =~ /^2nd.+?(\d\d\.\d\d)/){ 
 
207
                                $cdtableobj->{'_coding_gc'}{'2'} = $1;
 
208
                                }
 
209
                elsif ($line =~ /^3rd.+?(\d\d\.\d\d)/){ 
 
210
                                $cdtableobj->{'_coding_gc'}{'3'} = $1;
 
211
                                }
 
212
 
 
213
                elsif   ($line =~ /^Gen.+?(\d+)/){ 
 
214
                                $cdtableobj->{'_genetic_code'} = $1;
 
215
                        }
 
216
        }
 
217
                ## check has been parsed ok ##
 
218
                if (scalar keys %{$cdtableobj->{'_table'}} != 21) {
 
219
                        $cdtableobj->warn("probable parsing error - should be 21 entries for 20aa + stop codon");
 
220
                }
 
221
                return $cdtableobj;
 
222
                
 
223
}
 
224
 
 
225
 
 
226