Bio/CodonUsage/IO.pm

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