Merge pull request #441 from broadinstitute/jt_gvcf_idx_user_error

Jt gvcf idx user error

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<List<VariantContext>, 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<List<VariantContext>, 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 ) {

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",
+        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);
+                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",
+        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");
+                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);
+    }
 }

public/java/src/org/broadinstitute/sting/gatk/arguments/GATKArgumentCollection.java

@@ -473,9 +473,11 @@ public class GATKArgumentCollection {
      */
 
     @Argument(fullName="variant_index_type",shortName = "variant_index_type",doc="which type of IndexCreator to use for VCF/BCF indices",required=false)
+    @Advanced
     public GATKVCFIndexType variant_index_type = GATKVCFUtils.DEFAULT_INDEX_TYPE;
 
     @Argument(fullName="variant_index_parameter",shortName = "variant_index_parameter",doc="the parameter (bin width or features per bin) to pass to the VCF/BCF IndexCreator",required=false)
+    @Advanced
     public int variant_index_parameter = GATKVCFUtils.DEFAULT_INDEX_PARAMETER;
 }
 

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

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;