scripts/gen_solqual_lookup.pl

   1 #!/usr/bin/perl -w
   2
   3 use warnings;
   4 use strict;
   5
   6 sub log10($) {
   7         return log(shift) / log(10.0);
   8 }
   9
  10 sub round {
  11     my($number) = shift;
  12     return int($number + .5 * ($number <=> 0));
  13 }
  14
  15 # Convert from solexa qual to probability of miscall
  16 sub phredToP($) {
  17         my $sol = shift;
  18         my $p = (10.0 ** (($sol) / -10.0));
  19         ($p >= 0.0 && $p <= 1.0) || die "Bad prob: $p, from sol $sol";
  20         return $p;
  21 }
  22
  23 # Convert from phred qual to probability of miscall
  24 sub solToP($) {
  25         my $phred = shift;
  26         my $x = (10.0 ** (($phred) / -10.0));
  27         my $p = $x / (1.0 + $x);
  28         ($p >= 0.0 && $p <= 1.0) || die "Bad prob: $p, from x $x, phred $phred";
  29         return $p;
  30 }
  31
  32 # Convert from probability of miscall to phred qual
  33 sub pToPhred($) {
  34         my $p = shift;
  35         ($p >= 0.0 && $p <= 1.0) || die "Bad prob: $p";
  36         return round(-10.0 * log10($p));
  37 }
  38
  39 # Convert from probability of miscall to solexa qual
  40 sub pToSol($) {
  41         my $p = shift;
  42         ($p >= 0.0 && $p <= 1.0) || die "Bad prob: $p";
  43         return 0 if($p == 1.0);
  44         return round(-10.0 * log10($p / (1.0 - $p)));
  45 }
  46
  47 # Print conversion table from Phred to Solexa
  48 print "uint8_t solToPhred[] = {";
  49 my $cols = 10;
  50 my $cnt = 0;
  51 for(my $i = -10; $i < 256; $i++) {
  52         # Solexa qual = $i
  53         my $p = solToP($i);
  54         my $ph = pToPhred($p);
  55         print "\n\t/* $i */ " if($cnt == 0);
  56         $cnt++;
  57         $cnt = 0 if($cnt == 10);
  58         print "$ph";
  59         print ", " if($i < 255);
  60 }
  61 print "\n};\n";