perl/etc/kaks/get1986property.pl

   1 #!/bin/env perl
   2 use strict;
   3 use warnings;
   4 use Data::Dump qw(ddx);
   5
   6 my (%gen_code,%t,@t,%start_codes,$code,%AA,@AAs);
   7 $code=<<Ecode;
   8   AAs  = FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG
   9 Starts = ---M---------------M---------------M----------------------------
  10 Base1  = TTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGG
  11 Base2  = TTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGG
  12 Base3  = TCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAG
  13 Ecode
  14
  15 @t=split /\n/,$code;
  16 map {s/\s//g;@_=split /=/;$t{$_[0]}=[split //,$_[1]]} @t;
  17 for (@{$t{AAs}}) {
  18         ++$AA{$_};
  19         my $base=(shift @{$t{Base1}}).(shift @{$t{Base2}}).(shift @{$t{Base3}});
  20 #print "$base";
  21         my $start=shift @{$t{Starts}};
  22         ++$start_codes{$base} if $start eq 'M';
  23         $gen_code{$base}=$_;
  24 #print " -> [$_]\n";
  25 }
  26 @AAs = sort keys %AA;
  27
  28 ddx \%gen_code;
  29 ddx \%start_codes;
  30 ddx scalar @AAs,\@AAs;
  31
  32 my @Bases = sort qw( A T C G );
  33 my %mutBases;
  34 $mutBases{$_} = [] for @Bases;
  35 for my $mut (@Bases) {
  36         for (@Bases) {
  37                 push @{$mutBases{$_}},$mut if $mut ne $_;
  38         }
  39 }
  40 ddx \%mutBases;
  41
  42 for my $pos (0..2) {
  43         my ($cnsSyn,$cntAAc) = (0,0);
  44         for my $codon (sort keys %gen_code) {
  45                 my $ntApos = substr $codon,$pos,1;
  46                 my $oriAA = $gen_code{$codon};
  47                 next if $oriAA eq '*';
  48                 print "$pos\t$codon:$oriAA\t";
  49                 for my $mut ( @{$mutBases{$ntApos}} ) {
  50                         my $mutCodon = $codon;
  51                         substr $mutCodon,$pos,1,$mut;
  52                         my $newAA = $gen_code{$mutCodon};
  53                         next if $newAA eq '*';
  54                         print "$mutCodon:$newAA";
  55                         if ( $oriAA eq $newAA ) {
  56                                 ++$cnsSyn;
  57                                 print ".";
  58                         } else {
  59                                 ++$cntAAc;
  60                                 print "x";
  61                         }
  62                         print " ";
  63                 }
  64                 print "\n";
  65         }
  66         print ">>$pos: $cnsSyn,$cntAAc,",$cnsSyn/($cnsSyn+$cntAAc),"\n";
  67 }
  68
  69 __END__
  70 $ perl get1986property.pl |grep \>\>
  71 >>0: 8,166,0.0459770114942529
  72 >>1: 0,176,0
  73 >>2: 126,50,0.715909090909091
  74
  75