zzh
/
gatk-3.8
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

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.
This commit is contained in:
Ryan Poplin 2013-04-10 11:23:29 -04:00
parent 2fb8b61dd0
commit 850be5e9da
3 changed files with 18 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

32
protected/java/test/org/broadinstitute/sting/utils/SWPairwiseAlignmentUnitTest.java
View File

@ -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<Object[]> tests = new ArrayList<Object[]>();
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);
}

5
public/java/src/org/broadinstitute/sting/utils/SWPairwiseAlignment.java
View File

@ -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
}

6
public/java/test/org/broadinstitute/sting/utils/haplotypeBAMWriter/HaplotypeBAMWriterUnitTest.java
View File

@ -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";