From 850be5e9da3131eb056758d4d380924d1d000fbe Mon Sep 17 00:00:00 2001 From: Ryan Poplin Date: Wed, 10 Apr 2013 11:23:29 -0400 Subject: [PATCH] Bug fix in SWPairwiseAlignment. -- When the alignments are sufficiently apart from each other all the scores in the sw matrix could be negative which screwed up the max score calculation since it started at zero. --- .../utils/SWPairwiseAlignmentUnitTest.java | 32 +++++++++---------- .../sting/utils/SWPairwiseAlignment.java | 5 ++- .../HaplotypeBAMWriterUnitTest.java | 6 ---- 3 files changed, 18 insertions(+), 25 deletions(-) diff --git a/protected/java/test/org/broadinstitute/sting/utils/SWPairwiseAlignmentUnitTest.java b/protected/java/test/org/broadinstitute/sting/utils/SWPairwiseAlignmentUnitTest.java index 6d3c310b7..c55b4147d 100644 --- a/protected/java/test/org/broadinstitute/sting/utils/SWPairwiseAlignmentUnitTest.java +++ b/protected/java/test/org/broadinstitute/sting/utils/SWPairwiseAlignmentUnitTest.java @@ -72,22 +72,22 @@ public class SWPairwiseAlignmentUnitTest extends BaseTest { Assert.assertEquals(sw.getCigar().toString(), expectedCigar); } - // TODO - // TODO - // TODO this example demonstrates some kind of failure mode of SW that results in the read not being aligned - // TODO to the reference at all. It has something to do with the specific parameters provided to the - // TODO SW code. With the default parameters the result is the one expected. With the specified parameters - // TODO the code fails - // TODO - // TODO - @Test(enabled = false) - public void testOddNoAlignment() { - final String reference = "AAAGACTACTG"; - final String read = "AACGGACACTG"; - final int expectedStart = 0; - final String expectedCigar = "11M"; - final SWPairwiseAlignment sw = new SWPairwiseAlignment(reference.getBytes(), read.getBytes(), 5.0, -10.0, -22.0, -1.2); - sw.printAlignment(reference.getBytes(), read.getBytes()); + @DataProvider(name = "OddNoAlignment") + public Object[][] makeOddNoAlignment() { + List tests = new ArrayList(); + + final String ref1 = "AAAGACTACTG"; + final String read1 = "AACGGACACTG"; + tests.add(new Object[]{ref1, read1, 5.0, -10.0, -22.0, -1.2, 1, "2M2I3M1D4M"}); + tests.add(new Object[]{ref1, read1, 20.0, -5.0, -30.0, -2.2, 0, "11M"}); + + return tests.toArray(new Object[][]{}); + } + + @Test(dataProvider = "OddNoAlignment", enabled = true) + public void testOddNoAlignment(final String reference, final String read, final double match, final double mismatch, final double gap, final double gap_extend, + final int expectedStart, final String expectedCigar) { + final SWPairwiseAlignment sw = new SWPairwiseAlignment(reference.getBytes(), read.getBytes(), match, mismatch, gap, gap_extend); Assert.assertEquals(sw.getAlignmentStart2wrt1(), expectedStart); Assert.assertEquals(sw.getCigar().toString(), expectedCigar); } diff --git a/public/java/src/org/broadinstitute/sting/utils/SWPairwiseAlignment.java b/public/java/src/org/broadinstitute/sting/utils/SWPairwiseAlignment.java index e501cf40a..6c8beb32d 100644 --- a/public/java/src/org/broadinstitute/sting/utils/SWPairwiseAlignment.java +++ b/public/java/src/org/broadinstitute/sting/utils/SWPairwiseAlignment.java @@ -205,10 +205,9 @@ public final class SWPairwiseAlignment { private void calculateCigar(int n, int m, double [] sw, int [] btrack) { // p holds the position we start backtracking from; we will be assembling a cigar in the backwards order - //PrimitivePair.Int p = new PrimitivePair.Int(); int p1 = 0, p2 = 0; - double maxscore = 0.0; + double maxscore = Double.NEGATIVE_INFINITY; // sw scores are allowed to be negative int segment_length = 0; // length of the segment (continuous matches, insertions or deletions) // look for largest score. we use >= combined with the traversal direction @@ -259,7 +258,7 @@ public final class SWPairwiseAlignment { // move to next best location in the sw matrix: switch( new_state ) { - case MSTATE: data_offset -= (m+2); p1--; p2--; break; // move back along the diag in th esw matrix + case MSTATE: data_offset -= (m+2); p1--; p2--; break; // move back along the diag in the sw matrix case ISTATE: data_offset -= step_length; p2 -= step_length; break; // move left case DSTATE: data_offset -= (m+1)*step_length; p1 -= step_length; break; // move up } diff --git a/public/java/test/org/broadinstitute/sting/utils/haplotypeBAMWriter/HaplotypeBAMWriterUnitTest.java b/public/java/test/org/broadinstitute/sting/utils/haplotypeBAMWriter/HaplotypeBAMWriterUnitTest.java index 89d87a3c3..db16582b8 100644 --- a/public/java/test/org/broadinstitute/sting/utils/haplotypeBAMWriter/HaplotypeBAMWriterUnitTest.java +++ b/public/java/test/org/broadinstitute/sting/utils/haplotypeBAMWriter/HaplotypeBAMWriterUnitTest.java @@ -127,12 +127,6 @@ public class HaplotypeBAMWriterUnitTest extends BaseTest { } } - // test that reads without a good alignment to hap get excluded - { - final GATKSAMRecord read = makeRead("NNNNN"); - tests.add(new Object[]{read, allM, 10, -1, null}); - } - // example case of bad alignment because SW doesn't necessarily left-align indels { final String hap = "ACTGTGGGTTCCTCTTATTTTATTTCTACATCAATGTTCATATTTAACTTATTATTTTATCTTATTTTTAAATTTCTTTTATGTTGAGCCTTGATGAAAGCCATAGGTTCTCTCATATAATTGTATGTGTATGTATGTATATGTACATAATATATACATATATGTATATGTATGTGTATGTACATAATATATACGTATATGTATGTGTATGTACATAATATATACGTATATGTATGTGTATGTACATAATATATACGTATATGTATGTGTATGTACATAATATATACGTATATGTATGTGTATGTACATAATATATACGTATATGTATGTGTATGTGTATTACATAATATATACATATATGTATATATTATGTATATGTACATAATATATACATATATG";