htswanalysis/scripts/count_bases.pm

   1 #!/usr/bin/perl -w
   2
   3 use strict;
   4 use warnings;
   5
   6 my %calls;
   7 my @a = (0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0);
   8 my @c = (0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0);
   9 my @g = (0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0);
  10 my @t = (0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0);
  11
  12 $calls{'A'} = \@a;
  13 $calls{'C'} = \@c;
  14 $calls{'G'} = \@g;
  15 $calls{'T'} = \@t;
  16
  17 my $count = 0;
  18 my $length;
  19
  20 while( my $read = <>) {
  21   chomp $read;
  22   $length = length($read);
  23   $count++;
  24   for(0..$length-1) {
  25     my $base = uc(substr($read,$_,1));
  26     $calls{$base}[$_] += 1;
  27   }
  28 }
  29
  30 for(0..$length-1) {
  31   print "$_\t",$calls{A}[$_]/$count,"\t",$calls{C}[$_]/$count,"\t",$calls{G}[$_]/$count,"\t",$calls{T}[$_]/$count,"\n";
  32 }