SNVQ recalibration for flow based reads

src/main/java/org/broadinstitute/hellbender/tools/walkers/featuremapping/AddFlowSNVQualityArgumentCollection.java

@@ -1,5 +1,6 @@
 package org.broadinstitute.hellbender.tools.walkers.featuremapping;
 
+import org.broadinstitute.barclay.argparser.Advanced;
 import org.broadinstitute.barclay.argparser.Argument;
 import org.broadinstitute.barclay.argparser.Hidden;
 
@@ -11,6 +12,13 @@
  */
 public class AddFlowSNVQualityArgumentCollection implements Serializable{
     private static final long serialVersionUID = 1L;
+    public static final String MAX_PHRED_SCORE_FULL_NAME = "max-phred-score";
+    public static final String KEEP_SUPPLEMENTARY_ALIGNMENTS_FULL_NAME = "keep-supplementary-alignments";
+    public static final String INCLUDE_QC_FAILED_READ_FULL_NAME = "include-qc-failed-read";
+    public static final String SNVQ_MODE_FULL_NAME = "snvq-mode";
+    public static final String OUTPUT_QUALITY_ATTRIBUTE_FULL_NAME = "output-quality-attribute";
+    public static final String DEBUG_READ_NAME_FULL_NAME = "debug-read-name";
+    public static final String DEBUG_COLLECT_STATS_INTO_FULL_NAME = "debug-collect-stats-into";
 
     public enum SnvqModeEnum {
         Legacy,
@@ -23,25 +31,40 @@ public enum SnvqModeEnum {
      *  maximum value for
      *  delta in score
      **/
-    @Argument(fullName = "limit-phred-score", doc = "Limit value for phred scores", optional = true)
-    public double limitPhredScore = Double.NaN;
+    @Argument(fullName = MAX_PHRED_SCORE_FULL_NAME, doc = "Limit value for phred scores", optional = true)
+    public double maxPhredScore = Double.NaN;
 
     /**
      *  keep supplementary alignments?
      **/
-    @Argument(fullName = "keep-supplementary-alignments", doc = "keep supplementary alignments ?", optional = true)
+    @Argument(fullName = KEEP_SUPPLEMENTARY_ALIGNMENTS_FULL_NAME, doc = "keep supplementary alignments ?", optional = true)
     public boolean keepSupplementaryAlignments = true;
 
+    @Advanced
+    @Argument(fullName= INCLUDE_QC_FAILED_READ_FULL_NAME, doc = "include reads with QC failed flag", optional = true)
+    public boolean includeQcFailedReads = true;
+
     /**
      * snvq computation mode
      */
-    @Argument(fullName = "snvq-mode", doc = "ksnvq computation mode", optional = true)
+    @Argument(fullName = SNVQ_MODE_FULL_NAME, doc = "snvq calculation mode.", optional = true)
     public SnvqModeEnum snvMode = SnvqModeEnum.Geometric;
 
+    /**
+     * alternate quality attribute to set instead of the usual quality string
+     */
+    @Argument(fullName = OUTPUT_QUALITY_ATTRIBUTE_FULL_NAME, doc = "alternate quality attribute to set instead of the usual quality string.", optional = true)
+    public String outputQualityAttribute = null;
+
     /**
      *  debug read names?
      **/
     @Hidden
-    @Argument(fullName = "debug-read-name", doc = "debug specific reads?", optional = true)
+    @Argument(fullName = DEBUG_READ_NAME_FULL_NAME, doc = "Read names of reads to output details of as part of debugging. ", optional = true)
     public List<String> debugReadName = null;
+
+    @Advanced
+    @Hidden
+    @Argument(fullName= DEBUG_COLLECT_STATS_INTO_FULL_NAME, doc = "Statistics about the reads will be output to given filename.", optional = true)
+    public String debugCollectStatsInto = null;
 }

src/main/java/org/broadinstitute/hellbender/tools/walkers/groundtruth/AddFlowBaseQuality.java

@@ -239,6 +239,25 @@ private static double[][] extractErrorProbBands(final FlowBasedRead flowRead, fi
         return result;
     }
 
+    /**
+     * The following functions estimate the error probability for an hmer. specifically two error
+     * probability values are generated: one for the first base of the hmer and another for the
+     * rest of its bases.
+     *
+     * The computation itself is performed in a subsequent function: generateSidedHmerBaseErrorProbability
+     * It iterates over the possible valid combinations of errors and sums them up.
+     *
+     * @param key - key (hmer length) in flow space
+     * @param errorProbBands - for each flow (position in the key) three error probabilities are provided:
+     *                       [0] - for the hmer being one base shorter
+     *                       [1] - for the hmer to be at its length
+     *                       [2] - for the hmer to be one base longer
+     * @param flow - the flow (index) for which to generate the probabilities (0 <= flow < key.length)
+     * @param flowOrderLength - the cycle length of of the flow order (usually 4)
+     * @return an array of two probabilities:
+     * [0] - probability for the first base of the hmer
+     * [1] - probability for the rest of the bases of the hmer
+     */
     @VisibleForTesting
     protected static double[] generateHmerBaseErrorProbabilities(final int[] key, final double[][] errorProbBands, final int flow, final int flowOrderLength) {
 

src/main/java/org/broadinstitute/hellbender/tools/walkers/groundtruth/GroundTruthReadsBuilder.java

@@ -858,9 +858,9 @@ else if ( hasQ || hasZ ) {
         cols.put("ReadName", read.getName());
 
         // haplotypes and reference scores
-        cols.put("PaternalHaplotypeScore", paternal.score);
-        cols.put("MaternalHaplotypeScore", maternal.score);
-        cols.put("RefHaplotypeScore", refScore);
+        cols.put("PaternalHaplotypeScore", String.format("%.6f", paternal.score));
+        cols.put("MaternalHaplotypeScore", String.format("%.6f", maternal.score));
+        cols.put("RefHaplotypeScore", String.format("%.6f", refScore));
 
         // build haplotype keys
         final FlowBasedReadUtils.ReadGroupInfo rgInfo = FlowBasedReadUtils.getReadGroupInfo(getHeaderForReads(), read);

src/main/java/org/broadinstitute/hellbender/tools/walkers/groundtruth/SeriesStats.java

@@ -67,21 +67,6 @@ public void add(double v) {
         }
     }
 
-    public void aux(int v, int auxValue) {
-        aux((double)v, auxValue);
-    }
-
-    public void aux(double v, int auxValue) {
-
-        if ( auxBins.containsKey(v) ) {
-            auxBins.get(v).add(auxValue);
-        } else {
-            SeriesStats ss = new SeriesStats();
-            ss.add(auxValue);
-            auxBins.put(v, ss);
-        }
-    }
-
     public double getLast() {
         return last;
     }

src/main/java/org/broadinstitute/hellbender/utils/read/FlowBasedRead.java

@@ -254,33 +254,6 @@ public FlowBasedRead(final SAMRecord samRecord, final String flowOrder, final in
 
     }
 
-    public FlowBasedRead(final FlowBasedRead other, final boolean deepCopy) {
-        super(other.samRecord);
-        if ( deepCopy ) {
-            forwardSequence = other.forwardSequence.clone();
-            key = other.key.clone();
-            flow2base = other.flow2base.clone();
-            flowOrder = other.flowOrder.clone();
-            flowMatrix = other.flowMatrix.clone();
-        } else {
-            forwardSequence = other.forwardSequence;
-            key = other.key;
-            flow2base = other.flow2base;
-            flowOrder = other.flowOrder;
-            flowMatrix = other.flowMatrix;
-        }
-        maxHmer = other.maxHmer;
-        validKey = other.validKey;
-        direction = other.direction;
-        baseClipped = other.baseClipped;
-        trimLeftBase = other.trimLeftBase;
-        trimRightBase = other.trimRightBase;
-        fbargs = other.fbargs;
-        readInsQuals = other.readInsQuals;
-        readDelQuals = other.readDelQuals;
-        overallGCP = other.overallGCP;
-    }
-
     //since the last unclipped flow is uncertain (we give high probabilities to
     //also hmers higher than the called hmer)
     private void spreadFlowLengthProbsAcrossCountsAtFlow(final int flowToSpread) {
@@ -332,8 +305,6 @@ private void readFlowMatrix(final String _flowOrder) {
         // access qual, convert to flow representation
         final byte[]      quals = samRecord.getBaseQualities();
         final byte[]      tp = samRecord.getSignedByteArrayAttribute(FLOW_MATRIX_TAG_NAME);
-        // initialize matrix
-
 
         boolean specialTreatmentForZeroCalls = false;
         final byte[]      t0 = SAMUtils.fastqToPhred(samRecord.getStringAttribute(FLOW_MATRIX_T0_TAG_NAME));
@@ -864,12 +835,15 @@ private void applyFilteringFlowMatrix(){
     }
 
     /**
-     * clip probability values to fbargs.probabilityRatioThreshold
+     * clip probability values in a way that probability between perHmerMinErrorProb and 3*perHmerMinErrorProb
+     * is automatically clipped to perHmerMinErrorProb. This is done to avoid issues with rounding of
+     * small error probabilities in the basecalling
      */
     private void clipProbs() {
+        double probabilityThreshold = perHmerMinErrorProb*3;
         for ( int i = 0 ; i < getMaxHmer(); i++ ) {
             for ( int j =0; j < getNFlows(); j++) {
-                if ((flowMatrix[i][j] <= perHmerMinErrorProb*3) &&
+                if ((flowMatrix[i][j] <= probabilityThreshold) &&
                         (key[j]!=i)) {
                     flowMatrix[i][j] = perHmerMinErrorProb;
                 }

src/main/java/org/broadinstitute/hellbender/utils/read/FlowBasedReadUtils.java

@@ -480,4 +480,10 @@ public static CycleSkipStatus getCycleSkipStatus(FlowBasedRead read, ReferenceCo
         return status;
     }
 
+    public static int calcFlowOrderLength(String flowOrder) {
+
+        final int i = flowOrder.indexOf(flowOrder.charAt(0), 1);
+
+        return (i < 0) ? flowOrder.length() : i;
+    }
 }

src/test/java/org/broadinstitute/hellbender/tools/walkers/featuremapping/AddFlowSNVQualityIntegrationTest.java

@@ -32,7 +32,7 @@ public void testBasic() throws IOException {
                 "-O", outputFile.getAbsolutePath(),
                 "-I", publicTestDir + FlowTestConstants.ADD_FLOW_SNVQ_DATA_DIR + "/add_flow_snvq_input.bam",
                 "-L", "chr1:1-15000",
-                "--limit-phred-score", "50",
+                "--max-phred-score", "50",
                 "--debug-read-name", "30020185_2-UGAv3-182-1989782468",
                 "--verbosity", "INFO"
         };
@@ -45,7 +45,38 @@ public void testBasic() throws IOException {
 
         // walk the output and expected files, compare non-comment lines
         if ( !UPDATE_EXACT_MATCH_EXPECTED_OUTPUTS ) {
-            IntegrationTestSpec.assertEqualTextFiles(outputFile, expectedFile, "@PG");
+            IntegrationTestSpec.assertEqualTextFiles(outputFile, expectedFile, "@");
+        }
+    }
+
+    @Test
+    public void testBasicBQ() throws IOException {
+
+        final File outputDir = createTempDir("testAddFlowSNVQTest");
+        final String filename = "add_flow_snvq_output_bq.sam";
+        final File expectedFile = new File(testDir + "/" + filename);
+        final File outputFile = UPDATE_EXACT_MATCH_EXPECTED_OUTPUTS ? expectedFile : new File(outputDir + "/" + filename);
+
+        final String[] args = new String[] {
+                "-R", largeFileTestDir + "/Homo_sapiens_assembly38.fasta.gz",
+                "-O", outputFile.getAbsolutePath(),
+                "-I", publicTestDir + FlowTestConstants.ADD_FLOW_SNVQ_DATA_DIR + "/add_flow_snvq_input.bam",
+                "-L", "chr1:1-15000",
+                "--max-phred-score", "50",
+                "--debug-read-name", "30020185_2-UGAv3-182-1989782468",
+                "--verbosity", "INFO",
+                "--output-quality-attribute", "BQ"
+        };
+
+        // run the tool
+        runCommandLine(args);  // no assert, just make sure we don't throw
+
+        // make sure we've generated the output file
+        Assert.assertTrue(outputFile.exists());
+
+        // walk the output and expected files, compare non-comment lines
+        if ( !UPDATE_EXACT_MATCH_EXPECTED_OUTPUTS ) {
+            IntegrationTestSpec.assertEqualTextFiles(outputFile, expectedFile, "@");
         }
     }
 }

src/test/java/org/broadinstitute/hellbender/tools/walkers/featuremapping/AddFlowSNVQualityUnitTest.java

@@ -0,0 +1,39 @@
+package org.broadinstitute.hellbender.tools.walkers.featuremapping;
+
+import org.broadinstitute.hellbender.testutils.BaseTest;
+import org.testng.Assert;
+import org.testng.annotations.DataProvider;
+import org.testng.annotations.Test;
+
+public class AddFlowSNVQualityUnitTest extends BaseTest {
+
+    @DataProvider(name = "sliceIsValidDataProvider")
+    Object[][] sliceIsValidDataProvider() {
+        return new Object[][] {
+                { new int[] {0, 0, 0, 0, 0, 0}, 4, false },
+                { new int[] {0, 0, 1, 0, 0, 1}, 4, true },
+                { new int[] {1, 0, 0, 0, 1, 1}, 4, false },
+        };
+    }
+
+    @Test(dataProvider = "sliceIsValidDataProvider")
+    void testSliceIsValid(final int[] slice, final int flowOrderLength, final boolean isValid) {
+        Assert.assertEquals(AddFlowSNVQuality.sliceIsValidForConsideration(slice, flowOrderLength), isValid);
+    }
+
+    @DataProvider(name = "getSnvqDataProvider")
+    Object[][] getSnvqDataProvider() {
+        return new Object[][] {
+                { 0.5, 0.2, 0.3, AddFlowSNVQualityArgumentCollection.SnvqModeEnum.Legacy, 0.5 },
+                { 0.5, 0.2, 0.3, AddFlowSNVQualityArgumentCollection.SnvqModeEnum.Optimistic, 0.06 },
+                { 0.5, 0.2, 0.3, AddFlowSNVQualityArgumentCollection.SnvqModeEnum.Pessimistic, 0.44 },
+                { 0.5, 0.2, 0.3, AddFlowSNVQualityArgumentCollection.SnvqModeEnum.Geometric, Math.sqrt(0.06 * 0.44) },
+        };
+    }
+
+    @Test(dataProvider = "getSnvqDataProvider")
+    void testGetSnvq(final double sliceP, final double p1, final double p2, AddFlowSNVQualityArgumentCollection.SnvqModeEnum snvMode, final double expected) {
+        Assert.assertEquals(AddFlowSNVQuality.getSnvq(sliceP, p1, p2, snvMode), expected, 0.0001);
+    }
+
+}