tools/etc/artalnmisrate.pl

   1 #!/bin/env perl
   2 use strict;
   3 use warnings;
   4
   5 die "Usage: $0 <art.aln file>" unless @ARGV == 1;
   6
   7 open I,'<',$ARGV[0] or die;
   8 my ($ref,$seq,@seq,@ref,%Mis);
   9 my ($BaseCount,$MisBaseCount) = (0,0);
  10 while(<I>) {
  11         chomp($ref=<I>);
  12         chomp($seq=<I>);
  13         #print "[$ref] [$seq]\n";
  14         @seq=split //,$seq;
  15         @ref=split //,$ref;
  16         my $mis=0;
  17         for (@seq) {
  18                 my $base=shift @ref;
  19                 ++$mis if $base ne $_;
  20         }
  21         ++$Mis{int(0.5+100*$mis/(length $seq))};
  22         $BaseCount += length $seq;
  23         $MisBaseCount += $mis;
  24         #print "[$ref] [$seq] $mis\n";
  25 }
  26
  27 print "mis%\tcount\n",'-'x75,"\n";
  28 for (sort {$a<=>$b} keys %Mis) {
  29         print "$_\t$Mis{$_}\n";
  30 }
  31 print "\nMisRate: $MisBaseCount/$BaseCount = ",$MisBaseCount/$BaseCount,"\n";