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

UnifiedGenotyper bugfix: don't create haplotypes with 0 bases 

-- The PairHMM no longer allows us to create haplotypes with 0 bases.  The UG indel caller used to create such haplotypes.  Now we assign -Double.MAX_VALUE likelihoods to such haplotypes.
-- Add integration test to cover this case, along with private/testdata BAM
-- [Fixes #47523579]
This commit is contained in:
Mark DePristo 2013-04-13 14:53:30 -04:00
parent 6c9360b020
commit 3144eae51c
2 changed files with 24 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
protected/java/src/org/broadinstitute/sting/gatk/walkers/indels/PairHMMIndelErrorModel.java
View File

@ -388,20 +388,22 @@ public class PairHMMIndelErrorModel {
System.out.format("indStart: %d indStop: %d WinStart:%d WinStop:%d start: %d stop: %d readLength: %d C:%s\n",
indStart, indStop, ref.getWindow().getStart(), ref.getWindow().getStop(), startLocationInRefForHaplotypes, stopLocationInRefForHaplotypes, read.getReadLength(), read.getCigar().toString());
final byte[] haplotypeBases = Arrays.copyOfRange(haplotype.getBases(),
(int)indStart, (int)indStop);
final byte[] haplotypeBases = Arrays.copyOfRange(haplotype.getBases(), (int)indStart, (int)indStop);
if (firstHap) {
//no need to reallocate arrays for each new haplotype, as length won't change
pairHMM.initialize(readBases.length, haplotypeBases.length);
firstHap = false;
// it's possible that the indel starts at the last base of the haplotypes
if ( haplotypeBases.length == 0 ) {
readLikelihood = -Double.MAX_VALUE;
} else {
if (firstHap) {
//no need to reallocate arrays for each new haplotype, as length won't change
pairHMM.initialize(readBases.length, haplotypeBases.length);
firstHap = false;
}
readLikelihood = pairHMM.computeReadLikelihoodGivenHaplotypeLog10(haplotypeBases, readBases, readQuals,
baseInsertionQualities, baseDeletionQualities, contextLogGapContinuationProbabilities, firstHap);
}
readLikelihood = pairHMM.computeReadLikelihoodGivenHaplotypeLog10(haplotypeBases, readBases, readQuals,
baseInsertionQualities, baseDeletionQualities, contextLogGapContinuationProbabilities, firstHap);
if (DEBUG) {
System.out.println("H:"+new String(haplotypeBases));
System.out.println("R:"+new String(readBases));

11
protected/java/test/org/broadinstitute/sting/gatk/walkers/genotyper/UnifiedGenotyperIndelCallingIntegrationTest.java
View File

@ -51,6 +51,7 @@ import org.testng.annotations.Test;
import java.io.File;
import java.util.Arrays;
import java.util.Collections;
import java.util.List;
public class UnifiedGenotyperIndelCallingIntegrationTest extends WalkerTest {
@ -194,4 +195,14 @@ public class UnifiedGenotyperIndelCallingIntegrationTest extends WalkerTest {
Arrays.asList("3f07efb768e08650a7ce333edd4f9a52"));
executeTest("test minIndelFraction 1.0", spec);
}
// No testing of MD5 here, we previously blew up due to a 0 length haplotypes, so we just need to pass
@Test
public void testHaplotype0Length() {
WalkerTest.WalkerTestSpec spec = new WalkerTest.WalkerTestSpec(
"-T UnifiedGenotyper -I " + privateTestDir + "haplotype0.bam -L 20:47507681 -R " + b37KGReference + " -baq CALCULATE_AS_NECESSARY -glm BOTH -o /dev/null",
0,
Collections.<String>emptyList());
executeTest("testHaplotype0Length", spec);
}
}