diff --git a/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyWalker.java b/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyWalker.java index a722a334c..6e78f1d60 100755 --- a/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyWalker.java +++ b/java/src/org/broadinstitute/sting/playground/gatk/walkers/AlleleFrequencyWalker.java @@ -8,7 +8,6 @@ import org.broadinstitute.sting.playground.utils.AlleleFrequencyEstimate; import net.sf.samtools.SAMRecord; import java.util.List; -import java.util.Arrays; public class AlleleFrequencyWalker extends BasicLociWalker { @@ -88,19 +87,59 @@ public class AlleleFrequencyWalker extends BasicLociWalker 0); + Object SQ_field = reads.get(0).getAttribute("SQ"); + if (SQ_field == null) { + // Set all nonref qual scores to their share of the remaining probality not "used" by the reference base's qual + double nonref_quals = (1.0 - quals[i][callednum]) / 3; + for (int b=0; b<4; b++) + if (b != callednum) + quals[i][b] = nonref_quals; + }else{ + assert (SQ_field instanceof byte[]); + byte[] hex_quals = (byte[]) SQ_field; + System.out.printf("SQ field (hex): %s\n", bytesToHexString(hex_quals)); + System.out.printf("SAM record: %s\n", read.format()); + + int hex_qual = hex_quals[offset]; + int called2num = hex_qual & 0x3; + double qual2 = (double)(hex_qual >> 2) / 100.0; + System.out.printf("2ND %x %d %f\n", hex_qual, called2num, qual2); + quals[i][called2num] = qual2; + + // + double nonref_quals = (1.0 - quals[i][callednum] - quals[i][called2num]) / 3; + for (int b=0; b<4; b++) + if (b != callednum && b != called2num) + quals[i][b] = nonref_quals; + } } + //print_base_qual_matrix(quals); return quals; } + // Code pulled from SAMUtils for debugging + static String bytesToHexString(final byte[] data) { + final char[] chars = new char[2 * data.length]; + for (int i = 0; i < data.length; i++) { + final byte b = data[i]; + chars[2*i] = toHexDigit((b >> 4) & 0xF); + chars[2*i+1] = toHexDigit(b & 0xF); + } + return new String(chars); + } + + private static char toHexDigit(final int value) { + return (char) ((value < 10) ? ('0' + value) : ('A' + value - 10)); + } + public AlleleFrequencyEstimate AlleleFrequencyEstimator(String location, int N, byte[] bases, double[][] quals, int refnum, int altnum, int depth) {