diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/indels/PairHMMIndelErrorModel.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/indels/PairHMMIndelErrorModel.java index 0b0fa020e..318779cd2 100644 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/indels/PairHMMIndelErrorModel.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/indels/PairHMMIndelErrorModel.java @@ -205,10 +205,21 @@ public class PairHMMIndelErrorModel { } } - private LinkedHashMap trimHaplotypes(final LinkedHashMap haplotypeMap, - long startLocationInRefForHaplotypes, - long stopLocationInRefForHaplotypes, - final ReferenceContext ref){ + /** + * Trims the haplotypes in the given map to the provided start/stop. + * + * @param haplotypeMap the input map + * @param startLocationInRefForHaplotypes the start location of the trim + * @param stopLocationInRefForHaplotypes the stop location of the trim + * @param ref the reference context (used for debugging only, so can be null) + * @return a non-null mapping corresponding to the trimmed version of the original; + * some elements may be lost if trimming cannot be performed on them (e.g. they fall outside of the region to keep) + */ + protected static Map trimHaplotypes(final Map haplotypeMap, + long startLocationInRefForHaplotypes, + long stopLocationInRefForHaplotypes, + final ReferenceContext ref) { + if ( haplotypeMap == null ) throw new IllegalArgumentException("The input allele to haplotype map cannot be null"); final LinkedHashMap trimmedHaplotypeMap = new LinkedHashMap<>(); for (final Allele a: haplotypeMap.keySet()) { @@ -225,10 +236,13 @@ public class PairHMMIndelErrorModel { final long indStart = startLocationInRefForHaplotypes - haplotype.getStartPosition(); final long indStop = stopLocationInRefForHaplotypes - haplotype.getStartPosition(); + if ( indStart >= indStop ) + continue; - if (DEBUG) - System.out.format("indStart: %d indStop: %d WinStart:%d WinStop:%d start: %d stop: %d\n", - indStart, indStop, ref.getWindow().getStart(), ref.getWindow().getStop(), startLocationInRefForHaplotypes, stopLocationInRefForHaplotypes); + // commented out here because we need to make this method static for unit testing + //if (DEBUG) + // System.out.format("indStart: %d indStop: %d WinStart:%d WinStop:%d start: %d stop: %d\n", + // indStart, indStop, ref.getWindow().getStart(), ref.getWindow().getStop(), startLocationInRefForHaplotypes, stopLocationInRefForHaplotypes); // get the trimmed haplotype-bases array and create a new haplotype based on it. Pack this into the new map final byte[] trimmedHaplotypeBases = Arrays.copyOfRange(haplotype.getBases(), (int)indStart, (int)indStop); diff --git a/protected/java/test/org/broadinstitute/sting/gatk/walkers/indels/PairHMMIndelErrorModelUnitTest.java b/protected/java/test/org/broadinstitute/sting/gatk/walkers/indels/PairHMMIndelErrorModelUnitTest.java index bbbef43d3..3480b6775 100644 --- a/protected/java/test/org/broadinstitute/sting/gatk/walkers/indels/PairHMMIndelErrorModelUnitTest.java +++ b/protected/java/test/org/broadinstitute/sting/gatk/walkers/indels/PairHMMIndelErrorModelUnitTest.java @@ -50,9 +50,12 @@ package org.broadinstitute.sting.gatk.walkers.indels; import net.sf.picard.reference.IndexedFastaSequenceFile; import net.sf.samtools.SAMFileHeader; import org.broadinstitute.sting.BaseTest; +import org.broadinstitute.sting.utils.UnvalidatingGenomeLoc; import org.broadinstitute.sting.utils.fasta.CachingIndexedFastaSequenceFile; +import org.broadinstitute.sting.utils.haplotype.Haplotype; import org.broadinstitute.sting.utils.sam.ArtificialSAMUtils; import org.broadinstitute.sting.utils.sam.GATKSAMRecord; +import org.broadinstitute.variant.variantcontext.Allele; import org.testng.Assert; import org.testng.annotations.BeforeClass; import org.testng.annotations.DataProvider; @@ -138,4 +141,16 @@ public class PairHMMIndelErrorModelUnitTest extends BaseTest { Assert.assertEquals(PairHMMIndelErrorModel.mustClipDownstream(read, 13), true); Assert.assertEquals(PairHMMIndelErrorModel.mustClipDownstream(read, 14), false); } -} \ No newline at end of file + + @Test + public void trimHaplotypesToNullAlleleTest() { + // we need a case where start and stop > haplotype coordinates + final int start = 100, stop = 100; + final Haplotype h = new Haplotype(new byte[]{(byte)'A'}, new UnvalidatingGenomeLoc("1", 0, 10, 10)); + final Map input = new HashMap(1); + input.put(Allele.create("A"), h); + + final Map output = PairHMMIndelErrorModel.trimHaplotypes(input, start, stop, null); + Assert.assertTrue(output.isEmpty()); + } +}