6 # Modify scale/encoding of quality values in a FASTQ file.
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
26 # Default: convert 33-based Phred quals into 64-based Solexa qualss
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";
36 $inphred >= 33 || die "Input base must be >= 33, was $inphred";
38 $insolexa >= 33 || die "Input base must be >= 33, was $insolexa";
42 return log(shift) / log(10.0);
47 return int($number + .5 * ($number <=> 0));
50 # Convert from phred qual to probability of miscall
53 my $p = (10.0 ** (($phred) / -10.0));
54 ($p >= 0.0 && $p <= 1.0) || die "Bad prob: $p, from sol $phred";
58 # Convert from solexa qual to probability of miscall
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";
67 # Convert from probability of miscall to phred qual
70 ($p >= 0.0 && $p <= 1.0) || die "Bad prob: $p";
71 return round(-10.0 * log10($p));
74 # Convert from probability of miscall to solexa qual
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)));
83 my $name = $_; print $name;
84 my $seq = <>; print $seq;
85 my $name2 = <>; print $name2;
88 my @qual = split(//, $quals);
89 for(my $i = 0; $i <= $#qual; $i++) {
90 my $co = ord($qual[$i]);
92 # Convert input qual to p
95 $co >= 0 || die "Bad Phred input quality: $co";
101 # Convert p to output qual
104 $co >= 0 || die "Bad Phred output quality: $co";
110 $co >= 33 || die "Error: Output qual " . $co . " char is less than 33. Try a larger output base.";