implement heuristics that align haplotypes to reads without calling Smith-Waterman

src/main/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerEngine.java

@@ -39,6 +39,7 @@
 import org.broadinstitute.hellbender.utils.read.GATKRead;
 import org.broadinstitute.hellbender.utils.read.ReadUtils;
 import org.broadinstitute.hellbender.utils.smithwaterman.SmithWatermanAligner;
+import org.broadinstitute.hellbender.utils.smithwaterman.SmithWatermanJavaAligner;
 import org.broadinstitute.hellbender.utils.variant.GATKVCFConstants;
 import org.broadinstitute.hellbender.utils.variant.GATKVCFHeaderLines;
 import org.broadinstitute.hellbender.utils.variant.GATKVariantContextUtils;
@@ -699,7 +700,9 @@ private List<VariantContext> referenceModelForNoVariation(final AssemblyRegion r
      */
     public void shutdown() {
         likelihoodCalculationEngine.close();
+
         aligner.close();
+
         if ( haplotypeBAMWriter.isPresent() ) {
             haplotypeBAMWriter.get().close();
         }
@@ -710,8 +713,6 @@ public void shutdown() {
                 throw new RuntimeIOException(e);
             }
         }
-
-
     }
 
     private Set<GATKRead> filterNonPassingReads( final AssemblyRegion activeRegion ) {

src/main/java/org/broadinstitute/hellbender/utils/Utils.java

@@ -1057,6 +1057,25 @@ public static int lastIndexOf(final byte[] reference, final byte[] query) {
         return -1;
     }
 
+    public static int lastIndexOfAtMostOneMismatch(final byte[] reference, final byte[] query) {
+        int queryLength = query.length;
+
+        // start search from the last possible matching position and search to the left
+        for (int refIndex = reference.length - queryLength; refIndex >= 0; refIndex--) {
+            int mismatchCount = 0;
+            for (int queryIndex = 0; queryIndex < queryLength && mismatchCount < 2; queryIndex++) {
+                if (reference[refIndex+queryIndex] != query[queryIndex]) {
+                    mismatchCount++;
+                }
+            }
+
+            if (mismatchCount < 2) {
+                return refIndex;
+            }
+        }
+        return -1;
+    }
+
     /**
      * Simple wrapper for sticking elements of a int[] array into a List<Integer>
      * @param ar - the array whose elements should be listified

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

@@ -74,6 +74,7 @@ public static GATKRead createReadAlignedToRef(final GATKRead originalRead,
 
         // compute the smith-waterman alignment of read -> haplotype
         final SmithWatermanAlignment swPairwiseAlignment = aligner.align(haplotype.getBases(), originalRead.getBases(), CigarUtils.ALIGNMENT_TO_BEST_HAPLOTYPE_SW_PARAMETERS, SWOverhangStrategy.SOFTCLIP);
+
         if ( swPairwiseAlignment.getAlignmentOffset() == -1 ) {
             // sw can fail (reasons not clear) so if it happens just don't realign the read
             return originalRead;

src/main/java/org/broadinstitute/hellbender/utils/smithwaterman/SWNativeAlignerWrapper.java

@@ -25,6 +25,24 @@ public SWNativeAlignerWrapper(final SWAlignerNativeBinding aligner) {
         this.aligner = aligner;
     }
 
+
+    //methods for printing out number of SW/non-SW alignments. not needed
+    public int getNumOfAlignments()
+    {
+        System.out.println("SW native wrapper being used");
+        return 0;
+    }
+    public int noSW()
+    {
+        System.out.println("native SW being used");
+        return 0;
+    }
+    public int yesSW()
+    {
+        System.out.println("native SW being used");
+        return 0;
+    }
+
     @Override
     public SmithWatermanAlignment align(final byte[] reference, final byte[] alternate, final SWParameters parameters, final SWOverhangStrategy overhangStrategy){
         long startTime = System.nanoTime();

src/main/java/org/broadinstitute/hellbender/utils/smithwaterman/SmithWatermanJavaAligner.java

@@ -24,6 +24,8 @@
 public final class SmithWatermanJavaAligner implements SmithWatermanAligner {
     private static final SmithWatermanJavaAligner ALIGNER = new SmithWatermanJavaAligner();
     private long totalComputeTime = 0;
+    private long totalHeuristicTime = 0;
+    private long totalMismatchHeuristicTime = 0;
 
     /**
      * return the stateless singleton instance of SmithWatermanJavaAligner
@@ -74,24 +76,38 @@ public SmithWatermanAlignment align(final byte[] reference, final byte[] alterna
 
         final SmithWatermanAlignment alignmentResult;
 
+        //if exact match
         if (matchIndex != -1) {
             // generate the alignment result when the substring search was successful
-            final List<CigarElement> lce = new ArrayList<>(alternate.length);
-            lce.add(makeElement(State.MATCH, alternate.length));
+            final List<CigarElement> lce = Collections.singletonList(makeElement(State.MATCH, alternate.length));
             alignmentResult = new SWPairwiseAlignmentResult(AlignmentUtils.consolidateCigar(new Cigar(lce)), matchIndex);
         }
         else {
-            // run full Smith-Waterman
-            final int n = reference.length+1;
-            final int m = alternate.length+1;
-            final int[][] sw = new int[n][m];
-            final int[][] btrack=new int[n][m];
-
-            calculateMatrix(reference, alternate, sw, btrack, overhangStrategy, parameters);
-            alignmentResult = calculateCigar(sw, btrack, overhangStrategy); // length of the segment (continuous matches, insertions or deletions)
+            int matchIndex2 = -1;
+            if (overhangStrategy == SWOverhangStrategy.SOFTCLIP || overhangStrategy == SWOverhangStrategy.IGNORE) {
+                //look for one mismatch
+                matchIndex2 = Utils.lastIndexOfAtMostOneMismatch(reference, alternate);
+            }
+
+            if (matchIndex2 != -1) {
+                // generate the alignment result when the substring search was successful
+                final List<CigarElement> lce = Collections.singletonList(makeElement(State.MATCH, alternate.length));
+                alignmentResult = new SWPairwiseAlignmentResult(AlignmentUtils.consolidateCigar(new Cigar(lce)), matchIndex2);
+            }
+            else{
+                // run full Smith-Waterman
+                final int n = reference.length+1;
+                final int m = alternate.length+1;
+                final int[][] sw = new int[n][m];
+                final int[][] btrack=new int[n][m];
+
+                calculateMatrix(reference, alternate, sw, btrack, overhangStrategy, parameters);
+                alignmentResult = calculateCigar(sw, btrack, overhangStrategy); // length of the segment (continuous matches, insertions or deletions)
+            }
         }
 
         totalComputeTime += System.nanoTime() - startTime;
+
         return alignmentResult;
     }
 

src/test/java/org/broadinstitute/hellbender/utils/UtilsUnitTest.java

@@ -613,6 +613,109 @@ public void testLastIndexOfLastBoundaries() {
         Assert.assertEquals(result, expected);
     }
 
+    @Test
+    public void testLastIndexOfAtMostOneMismatchLastBoundaries() {
+        final String reference = "AAAACCCCTTTTGGGG";
+
+        // match right boundary of reference
+        final String query1 = "TGAGG";
+        int result = Utils.lastIndexOfAtMostOneMismatch(reference.getBytes(), query1.getBytes());
+        int expected = reference.length() - query1.length();
+        Assert.assertEquals(result, expected);
+
+        // match left boundary of reference
+        final String query2 = "AAGAC";
+        result = Utils.lastIndexOfAtMostOneMismatch(reference.getBytes(), query2.getBytes());
+        Assert.assertEquals(result, 0);
+    }
+
+    @Test
+    public void testLastIndexOfAtMostOneMismatchRandomExactMatch() {
+        final int numTests = 100;
+        final int referenceLength = 100;
+        final int queryLength = 10;
+
+        byte[] reference = new byte[referenceLength];
+        byte[] query = new byte[queryLength];
+
+        final Random rng = Utils.getRandomGenerator();
+
+        for (int i = 0; i < numTests; i++) {
+            randomByteString(rng, reference);
+            randomByteString(rng, query);
+
+            int index = -1;
+
+            //add query to reference at a random lo cation for 75% of the tests
+            if (i % 4 > 0) {
+                index = rng.nextInt(referenceLength - queryLength);
+                System.arraycopy(query,0, reference, index, queryLength);
+            }
+
+            final int result = Utils.lastIndexOfAtMostOneMismatch(reference, query);
+            final int expected = index;
+            Assert.assertEquals(result, expected);
+        }
+    }
+
+    @Test
+    public void testLastIndexOfAtMostOneMismatchTwoMismatches() {
+        final String reference = "AAAACCCCTTTTGGGG";
+
+        // match right boundary of reference
+        final String query1 = "AGAGG";
+        int result = Utils.lastIndexOfAtMostOneMismatch(reference.getBytes(), query1.getBytes());
+        Assert.assertEquals(result, -1);
+
+        // match right boundary of reference
+        final String query2 = "GGAAC";
+        int result2 = Utils.lastIndexOfAtMostOneMismatch(reference.getBytes(), query2.getBytes());
+        Assert.assertEquals(result2, -1);
+    }
+
+     @Test
+    public void testLastIndexOfAtMostOneMismatchOneMismatchTwice() {
+        final String reference = "AGGCCCCCTTTTGGGG";
+
+        //matches both boundaries of reference
+        final String query1 = "AGGG";
+        int result = Utils.lastIndexOfAtMostOneMismatch(reference.getBytes(), query1.getBytes());
+        //will catch first match it encounters from the right
+        int expected = reference.length() - query1.length();
+        Assert.assertEquals(result, expected);
+    }
+
+    @Test
+    public void testLastIndexOfAtMostOneMismatchRandom() {
+        final int numTests = 100;
+        final int referenceLength = 100;
+        final int queryLength = 10;
+
+        byte[] reference = new byte[referenceLength];
+        byte[] query = new byte[queryLength];
+
+        final Random rng = Utils.getRandomGenerator();
+
+        for (int i = 0; i < numTests; i++) {
+            randomByteString(rng, reference);
+            randomByteString(rng, query);
+
+            int index = -1;
+
+            // add one-off query to reference at a random location for 75% of the tests
+            if (i % 4 > 0) {
+                index = rng.nextInt(referenceLength - queryLength);
+                System.arraycopy(query,0, reference, index, queryLength);
+                int mismatchPosition = rng.nextInt(queryLength);
+                reference[index + mismatchPosition] = BaseUtils.simpleComplement(query[mismatchPosition]);
+            }
+
+            final int result = Utils.lastIndexOfAtMostOneMismatch(reference, query);
+            final int expected = index;
+            Assert.assertEquals(result, expected);
+        }
+    }
+
     private void randomByteString(Random rng, byte[] bytes) {
         for (int i = 0; i < bytes.length; i++) {
             bytes[i] = (byte)(rng.nextInt(94) + 32);