diff --git a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java index d795f5b90..6b7373273 100644 --- a/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java +++ b/protected/java/src/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCaller.java @@ -92,6 +92,7 @@ import org.broadinstitute.sting.utils.help.HelpConstants; import org.broadinstitute.sting.utils.pairhmm.PairHMM; import org.broadinstitute.sting.utils.sam.GATKSAMRecord; import org.broadinstitute.sting.utils.sam.ReadUtils; +import org.broadinstitute.sting.utils.variant.GATKVCFIndexType; import org.broadinstitute.sting.utils.variant.GATKVariantContextUtils; import org.broadinstitute.variant.variantcontext.*; import org.broadinstitute.variant.variantcontext.writer.VariantContextWriter; @@ -536,6 +537,10 @@ public class HaplotypeCaller extends ActiveRegionWalker, In ReferenceConfidenceModel referenceConfidenceModel = null; + // as determined experimentally Nov-Dec 2013 + protected final static GATKVCFIndexType OPTIMAL_GVCF_INDEX_TYPE = GATKVCFIndexType.LINEAR; + protected final static int OPTIMAL_GVCF_INDEX_PARAMETER = 128000; + //--------------------------------------------------------------------------------------------------------------- // // initialize @@ -601,6 +606,12 @@ public class HaplotypeCaller extends ActiveRegionWalker, In if ( samples.size() != 1 ) throw new UserException.BadArgumentValue("emitRefConfidence", "Can only be used in single sample mode currently"); headerInfo.addAll(referenceConfidenceModel.getVCFHeaderLines()); if ( emitReferenceConfidence == ReferenceConfidenceMode.GVCF ) { + // a kluge to enforce the use of this indexing strategy + if (getToolkit().getArguments().variant_index_type != OPTIMAL_GVCF_INDEX_TYPE || + getToolkit().getArguments().variant_index_parameter != OPTIMAL_GVCF_INDEX_PARAMETER) { + throw new UserException.GVCFIndexException(OPTIMAL_GVCF_INDEX_TYPE, OPTIMAL_GVCF_INDEX_PARAMETER); + } + try { vcfWriter = new GVCFWriter(vcfWriter, GVCFGQBands); } catch ( IllegalArgumentException e ) { diff --git a/protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerGVCFIntegrationTest.java b/protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerGVCFIntegrationTest.java index e430cd8d1..97744f126 100644 --- a/protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerGVCFIntegrationTest.java +++ b/protected/java/test/org/broadinstitute/sting/gatk/walkers/haplotypecaller/HaplotypeCallerGVCFIntegrationTest.java @@ -47,6 +47,8 @@ package org.broadinstitute.sting.gatk.walkers.haplotypecaller; import org.broadinstitute.sting.WalkerTest; +import org.broadinstitute.sting.utils.exceptions.UserException; +import org.broadinstitute.sting.utils.variant.GATKVCFIndexType; import org.testng.annotations.DataProvider; import org.testng.annotations.Test; @@ -79,8 +81,8 @@ public class HaplotypeCallerGVCFIntegrationTest extends WalkerTest { */ @Test(dataProvider = "MyDataProvider") public void testHCWithGVCF(String bam, HaplotypeCaller.ReferenceConfidenceMode mode, String intervals, String md5) { - final String commandLine = String.format("-T HaplotypeCaller --disableDithering --pcr_indel_model NONE -R %s -I %s %s -ERC %s --no_cmdline_in_header", - b37KGReference, bam, intervals, mode); + final String commandLine = String.format("-T HaplotypeCaller --disableDithering --pcr_indel_model NONE -R %s -I %s %s -ERC %s --no_cmdline_in_header -variant_index_type %s -variant_index_parameter %d", + b37KGReference, bam, intervals, mode, HaplotypeCaller.OPTIMAL_GVCF_INDEX_TYPE, HaplotypeCaller.OPTIMAL_GVCF_INDEX_PARAMETER); final String name = "testHCWithGVCF bam=" + bam + " intervals= " + intervals + " gvcf= " + mode; final WalkerTestSpec spec = new WalkerTestSpec(commandLine + " -o %s", Arrays.asList(md5)); executeTest(name, spec); @@ -88,10 +90,42 @@ public class HaplotypeCallerGVCFIntegrationTest extends WalkerTest { @Test public void testERCRegionWithNoCalledHaplotypes() { - final String commandLine = String.format("-T HaplotypeCaller --pcr_indel_model NONE -R %s -I %s -L %s -ERC GVCF", - b37KGReference, privateTestDir + "noCallRefModel.bam", "20:17000001-18000001"); + final String commandLine = String.format("-T HaplotypeCaller --pcr_indel_model NONE -R %s -I %s -L %s -ERC GVCF -variant_index_type %s -variant_index_parameter %d", + b37KGReference, privateTestDir + "noCallRefModel.bam", "20:17000001-18000001", HaplotypeCaller.OPTIMAL_GVCF_INDEX_TYPE, HaplotypeCaller.OPTIMAL_GVCF_INDEX_PARAMETER); final WalkerTestSpec spec = new WalkerTestSpec(commandLine + " -o %s", Arrays.asList("")); spec.disableShadowBCF(); executeTest("testERCRegionWithNoCalledHaplotypes", spec); } + + @Test() + public void testMissingGVCFIndexException() { + final String commandLine = String.format("-T HaplotypeCaller --pcr_indel_model NONE -R %s -I %s -L %s -ERC GVCF", + b37KGReference, privateTestDir + "noCallRefModel.bam", "20:17000001-18000001"); + final WalkerTestSpec spec = new WalkerTestSpec(commandLine + " -o %s", 1, UserException.GVCFIndexException.class); + spec.disableShadowBCF(); + executeTest("testMissingGVCFIndexingStrategyException", spec); + } + + @Test() + public void testWrongParameterGVCFIndexException() { + final String commandLine = String.format("-T HaplotypeCaller --pcr_indel_model NONE -R %s -I %s -L %s -ERC GVCF -variant_index_type %s -variant_index_parameter %d", + b37KGReference, privateTestDir + "noCallRefModel.bam", "20:17000001-18000001", HaplotypeCaller.OPTIMAL_GVCF_INDEX_TYPE, HaplotypeCaller.OPTIMAL_GVCF_INDEX_PARAMETER + 1); + final WalkerTestSpec spec = new WalkerTestSpec(commandLine + " -o %s", 1, UserException.GVCFIndexException.class); + spec.disableShadowBCF(); + executeTest("testMissingGVCFIndexingStrategyException", spec); + } + + @Test() + public void testWrongTypeGVCFIndexException() { + // ensure non-optimal, if optimal changes + GATKVCFIndexType type = GATKVCFIndexType.DYNAMIC_SEEK; + if (HaplotypeCaller.OPTIMAL_GVCF_INDEX_TYPE == GATKVCFIndexType.DYNAMIC_SEEK) + type = GATKVCFIndexType.DYNAMIC_SIZE; + + final String commandLine = String.format("-T HaplotypeCaller --pcr_indel_model NONE -R %s -I %s -L %s -ERC GVCF -variant_index_type %s -variant_index_parameter %d", + b37KGReference, privateTestDir + "noCallRefModel.bam", "20:17000001-18000001", type, HaplotypeCaller.OPTIMAL_GVCF_INDEX_PARAMETER); + final WalkerTestSpec spec = new WalkerTestSpec(commandLine + " -o %s", 1, UserException.GVCFIndexException.class); + spec.disableShadowBCF(); + executeTest("testMissingGVCFIndexingStrategyException", spec); + } } diff --git a/public/java/src/org/broadinstitute/sting/utils/exceptions/UserException.java b/public/java/src/org/broadinstitute/sting/utils/exceptions/UserException.java index 6126116c2..40a730029 100644 --- a/public/java/src/org/broadinstitute/sting/utils/exceptions/UserException.java +++ b/public/java/src/org/broadinstitute/sting/utils/exceptions/UserException.java @@ -33,6 +33,7 @@ import org.broadinstitute.sting.utils.GenomeLoc; import org.broadinstitute.sting.utils.help.DocumentedGATKFeature; import org.broadinstitute.sting.utils.help.HelpConstants; import org.broadinstitute.sting.utils.sam.ReadUtils; +import org.broadinstitute.sting.utils.variant.GATKVCFIndexType; import org.broadinstitute.variant.variantcontext.VariantContext; import java.io.File; @@ -455,6 +456,14 @@ public class UserException extends ReviewedStingException { } } + public static class GVCFIndexException extends UserException { + public GVCFIndexException (GATKVCFIndexType indexType, int indexParameter) { + super(String.format("GVCF output requires a specific indexing strategy. Please re-run including the arguments " + + "-variant_index_type %s -variant_index_parameter %d.", + indexType, indexParameter)); + } + } + /** * A special exception that happens only in the case where * the filesystem, by design or configuration, is completely unable diff --git a/public/java/test/org/broadinstitute/sting/WalkerTest.java b/public/java/test/org/broadinstitute/sting/WalkerTest.java index e8e791a8a..994a2419c 100644 --- a/public/java/test/org/broadinstitute/sting/WalkerTest.java +++ b/public/java/test/org/broadinstitute/sting/WalkerTest.java @@ -201,6 +201,7 @@ public class WalkerTest extends BaseTest { this.testClass = getCallingTestClass(); } + // @Test(expectedExceptions) doesn't work in integration tests, so use this instead public WalkerTestSpec(String args, int nOutputFiles, Class expectedException) { this.args = args; this.nOutputFiles = nOutputFiles;