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   if(!defined($length)) { $length = length($read); }
  23   if($read =~ /\./) { next; }
  24   $count++;
  25   for(0..$length-1) {
  26     my $base = uc(substr($read,$_,1));
  27     $calls{$base}[$_] += 1;
  28   }
  29 }
  30
  31 for(0..$length-1) {
  32   print "$_\t",$calls{A}[$_]/$count,"\t",$calls{C}[$_]/$count,"\t",$calls{G}[$_]/$count,"\t",$calls{T}[$_]/$count,"\n";
  33 }