scripts/convert_quals.pl

   1 #!/usr/bin/perl -w
   2
   3 #
   4 # convert_quals.pl
   5 #
   6 # Modify scale/encoding of quality values in a FASTQ file.
   7 #
   8 # Author: Ben Langmead
   9 #   Date: 5/5/2009
  10 #
  11 # p = probability that base is miscalled
  12 # Qphred = -10 * log10 (p)
  13 # Qsolexa = -10 * log10 (p / (1 - p))
  14 # See: http://en.wikipedia.org/wiki/FASTQ_format
  15 #
  16
  17 use strict;
  18 use warnings;
  19 use Getopt::Long;
  20
  21 my $inphred   = 33;
  22 my $insolexa  = 0;
  23 my $outphred  = 0;
  24 my $outsolexa = 64;
  25
  26 # Default: convert 33-based Phred quals into 64-based Solexa qualss
  27
  28 my $result =
  29         GetOptions ("inphred=i"   => \$inphred,
  30                     "insolexa=i"  => \$insolexa,
  31                     "outphred=i"  => \$outphred,
  32                     "outsolexa=i" => \$outsolexa);
  33 $result == 1 || die "One or more errors parsing script arguments";
  34
  35 if($inphred > 0) {
  36         $inphred >= 33 || die "Input base must be >= 33, was $inphred";
  37 } else {
  38         $insolexa >= 33 || die "Input base must be >= 33, was $insolexa";
  39 }
  40
  41 sub log10($) {
  42         return log(shift) / log(10.0);
  43 }
  44
  45 sub round {
  46     my($number) = shift;
  47     return int($number + .5 * ($number <=> 0));
  48 }
  49
  50 # Convert from phred qual to probability of miscall
  51 sub phredToP($) {
  52         my $phred = shift;
  53         my $p = (10.0 ** (($phred) / -10.0));
  54         ($p >= 0.0 && $p <= 1.0) || die "Bad prob: $p, from sol $phred";
  55         return $p;
  56 }
  57
  58 # Convert from solexa qual to probability of miscall
  59 sub solToP($) {
  60         my $sol = shift;
  61         my $x = (10.0 ** (($sol) / -10.0));
  62         my $p = $x / (1.0 + $x);
  63         ($p >= 0.0 && $p <= 1.0) || die "Bad prob: $p, from x $x, phred $sol";
  64         return $p;
  65 }
  66
  67 # Convert from probability of miscall to phred qual
  68 sub pToPhred($) {
  69         my $p = shift;
  70         ($p >= 0.0 && $p <= 1.0) || die "Bad prob: $p";
  71         return round(-10.0 * log10($p));
  72 }
  73
  74 # Convert from probability of miscall to solexa qual
  75 sub pToSol($) {
  76         my $p = shift;
  77         ($p >= 0.0 && $p <= 1.0) || die "Bad prob: $p";
  78         return 0.0 if $p == 1.0;
  79         return round(-10.0 * log10($p / (1.0 - $p)));
  80 }
  81
  82 while(<>) {
  83         my $name = $_;  print $name;
  84         my $seq = <>;   print $seq;
  85         my $name2 = <>; print $name2;
  86         my $quals = <>;
  87         chomp($quals);
  88         my @qual = split(//, $quals);
  89         for(my $i = 0; $i <= $#qual; $i++) {
  90                 my $co = ord($qual[$i]);
  91                 my $p;
  92                 # Convert input qual to p
  93                 if($inphred > 0) {
  94                         $co -= $inphred;
  95                         $co >= 0 || die "Bad Phred input quality: $co";
  96                         $p = phredToP($co);
  97                 } else {
  98                         $co -= $insolexa;
  99                         $p = solToP($co);
 100                 }
 101                 # Convert p to output qual
 102                 if($outphred > 0) {
 103                         $co = pToPhred($p);
 104                         $co >= 0 || die "Bad Phred output quality: $co";
 105                         $co += $outphred;
 106                 } else {
 107                         $co = pToSol($p);
 108                         $co += $outsolexa;
 109                 }
 110                 $co >= 33 || die "Error: Output qual " . $co . " char is less than 33.  Try a larger output base.";
 111                 print chr($co);
 112         }
 113         print "\n";
 114 }