New HaplotypeCaller priors for variants sites and homRef blocks #4793

ldgauthier · 2018-05-21T20:10:45Z

HC priors will let us reduce the WEx GVCF footprint. As a consequence, CalculateGenotypePosteriors now supports indels.

I fixed a bug and changed the args for CGP. We didn't have great tests, but CGP results will be fixed/improved in some cases.

droazen · 2018-05-21T20:55:51Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCaller.java

@@ -190,10 +193,14 @@ public AssemblyRegionEvaluator assemblyRegionEvaluator() {

    @Override
    public void onTraversalStart() {
+        if(hcArgs.genotypeArgs.applyPriors && hcArgs.genotypeArgs.samplePloidy != HomoSapiensConstants.DEFAULT_PLOIDY) {
+            throw new UserException("Priors can only be applied for diploid samples.");


This check should be moved to HaplotypeCallerEngine.validateAndInitializeArgs()

Is this moveable? It's gross having it in the tool itself rather than the HaplotypeCallerEngine

droazen · 2018-05-21T20:59:25Z

.../broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java

+        runCommandLine(args);
+
+        // Test for an exact match against past results
+        IntegrationTestSpec.assertEqualTextFiles(output, expected);


Can you think of a way to construct a good test that doesn't require an exact match on the output file? I've been trying to move HaplotypeCallerIntegrationTest away from that style of test. Eg., could you call VariantContextTestUtils.readEntireVCFIntoMemory() and inspect the VCs themselves?

Sure. I can refactor the integration test in the same way as in PosteriorProbabilitiesUtilsUnitTest.

codecov-io · 2018-05-23T20:06:04Z

Codecov Report

Merging #4793 into master will increase coverage by 0.027%.
The diff coverage is 88.889%.

@@               Coverage Diff               @@
##              master     #4793       +/-   ##
===============================================
+ Coverage     80.332%   80.359%   +0.027%     
- Complexity     17625     17670       +45     
===============================================
  Files           1088      1088               
  Lines          63849     63958      +109     
  Branches       10276     10307       +31     
===============================================
+ Hits           51291     51396      +105     
- Misses          8558      8559        +1     
- Partials        4000      4003        +3

Impacted Files	Coverage Δ	Complexity Δ
...lkers/ReferenceConfidenceVariantContextMerger.java	`95.257% <ø> (ø)`	`71 <0> (ø)`	⬇️
...tools/walkers/haplotypecaller/HaplotypeCaller.java	`93.939% <ø> (ø)`	`18 <0> (ø)`	⬇️
...lbender/utils/variant/GATKVariantContextUtils.java	`84.32% <ø> (ø)`	`216 <0> (ø)`	⬇️
...plotypecaller/HaplotypeCallerGenotypingEngine.java	`78.014% <100%> (+0.316%)`	`32 <0> (ø)`	⬇️
...hellbender/utils/test/VariantContextTestUtils.java	`83.544% <100%> (+0.283%)`	`66 <2> (+2)`	⬆️
...notyper/GenotypeCalculationArgumentCollection.java	`100% <100%> (ø)`	`2 <0> (ø)`	⬇️
...kers/haplotypecaller/ReferenceConfidenceModel.java	`93.532% <100%> (+0.339%)`	`68 <6> (+7)`	⬆️
...tools/walkers/genotyper/AlleleSubsettingUtils.java	`91.736% <100%> (+0.507%)`	`45 <0> (ø)`	⬇️
...kers/variantutils/PosteriorProbabilitiesUtils.java	`82.051% <84.211%> (+4.085%)`	`58 <3> (+15)`	⬆️
...walkers/haplotypecaller/HaplotypeCallerEngine.java	`79.643% <85.714%> (+0.676%)`	`71 <0> (+4)`	⬆️
... and 5 more

droazen · 2018-05-23T20:40:48Z

Adding @davidbenjamin as a second reviewer

davidbenjamin

Just some blue-collar coding suggestions.

davidbenjamin · 2018-05-24T01:55:20Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/AlleleSubsettingUtils.java

-            final GenotypeBuilder gb = useNewLikelihoods ? new GenotypeBuilder(g).PL(newLikelihoods).log10PError(newLog10GQ) : new GenotypeBuilder(g).noPL().noGQ();
-
+            final GenotypeBuilder gb;
+            final Map<String, Object> attributes = new HashMap<>(g.getExtendedAttributes());


attributes should be local to the if block.

davidbenjamin · 2018-05-24T01:57:18Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/genotyper/AlleleSubsettingUtils.java

@@ -79,8 +79,16 @@ public static GenotypesContext subsetAlleles(final GenotypesContext originalGs,
            }

            final boolean useNewLikelihoods = newLikelihoods != null && (depth != 0 || GATKVariantContextUtils.isInformative(newLikelihoods));
-            final GenotypeBuilder gb = useNewLikelihoods ? new GenotypeBuilder(g).PL(newLikelihoods).log10PError(newLog10GQ) : new GenotypeBuilder(g).noPL().noGQ();
-
+            final GenotypeBuilder gb;


Why not absorb the new GenotypeBuilder(g) into this declaration?

davidbenjamin · 2018-05-24T02:01:28Z

...broadinstitute/hellbender/tools/walkers/genotyper/GenotypeCalculationArgumentCollection.java

+    /**
+     * Apply population priors to genotype likelihoods.  If specified, a population resource VCF must be supplied via the -supporting argument.
+     */
+    @Argument(fullName = "apply-priors", optional = true)


It seems that this can be inferred from whether supportVariants is null.

I wanted to use the same arguments as CalculateGenotypePosteriors, in which case supporting-callsets is not particularly informative.

But why not eliminate this argument entirely, because it looks like it is intended to be true if and only if supportVariants != null?

davidbenjamin · 2018-05-24T02:05:05Z

...broadinstitute/hellbender/tools/walkers/genotyper/GenotypeCalculationArgumentCollection.java

+    /**
+     * When a variant is not seen in a panel, this argument controls whether to infer (and with what effective strength)
+     * that only reference alleles were observed at that site. E.g. "If not seen in 1000Genomes, treat it as AC=0,
+     * AN=2000". This is applied across all external panels, so if numRefIsMissing = 10, and the variant is absent in


This behavior for multiple panels doesn't seem like the most natural one. IMO a single number for the total imputed ref count is more intuitive. Maybe @droazen can break the tie.

I don't love any of the solutions here, but I do think if a variant is absent from all the resources then it should get a higher prior than if it's absent from just one. I suppose we make the assumption that the sets of samples in the resource files are disjoint and it's certainly not great that we have to apply the same number of each callset. But we can't determine the number of samples from the header because it might be a sites-only file and we can't use AN without traversing the whole VCF. It would greatly simplify things to limit CGP to one supporting callset as well, but I'm always hesitant to reduce functionality.

The thing is, if there are multiple resources the user can take that into account when choosing this value. What worries me most is the case of a single big resource, like gnomAD, along with a few boutique add-ons for whatever reason. Then the prior comes almost entirely from the big one and should hardly change when the little ones are added.

davidbenjamin · 2018-05-24T02:06:24Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCaller.java

        final ReferenceSequenceFile referenceReader = getReferenceReader(referenceArguments);
        hcEngine = new HaplotypeCallerEngine(hcArgs, createOutputBamIndex, createOutputBamMD5, getHeaderForReads(), referenceReader);

-        // The HC engine will make the right kind (VCF or GVCF) of writer for us
+            // The HC engine will make the right kind (VCF or GVCF) of writer for us


I agree with the old indentation here.

davidbenjamin · 2018-05-24T03:13:34Z

src/main/java/org/broadinstitute/hellbender/utils/test/VariantContextTestUtils.java

@@ -429,6 +429,12 @@ public static Boolean alleleSpecificAnnotationEquals(VariantContext actual, Vari
        return true;
    }

+    public static void assertGenotypeAttributeWasRemoved(final VariantContext actual, final VariantContext expected) {


The method's name doesn't suggest that it is hard-coded to check one particular attribute.

Oops. I initially had an attribute parameter before I realized the test expected a BiConsumer. Suggestions?

assertGenotypePosteriorsAttributeWasRemoved?

davidbenjamin · 2018-05-24T03:14:25Z

src/main/java/org/broadinstitute/hellbender/utils/variant/writers/HomRefBlock.java

+        if( genotype.hasExtendedAttribute(GATKVCFConstants.PHRED_SCALED_POSTERIORS_KEY)) {
+            if (minPPs == null ) {
+            minPPs = PosteriorProbabilitiesUtils.parsePosteriorsIntoPhredSpace(genotype);
+        } else { // otherwise take the min with the provided genotype's PLs


Indentation issues here.

davidbenjamin · 2018-05-24T03:17:45Z

...titute/hellbender/tools/walkers/variantutils/CalculateGenotypePosteriorsIntegrationTest.java

+                " -O %s" +
+                        " -R " + b37_reference_20_21 +    //NOTE: we need a reference for -L
+                        " -L 20:10,000,000-10,010,000" +
+                        " -supporting " + largeDir + "1000G.phase3.broad.withGenotypes.chr20.10100000.vcf" +


Use a variable for -supporting.

davidbenjamin · 2018-05-24T03:18:01Z

...titute/hellbender/tools/walkers/variantutils/CalculateGenotypePosteriorsIntegrationTest.java

+                " -O %s" +
+                        " -R " + b37_reference_20_21 +    //NOTE: we need a reference for -L
+                        " -L 20:10,000,000-10,010,000" +
+                        " -supporting " + largeDir + "1000G.phase3.broad.withGenotypes.chr20.10100000.vcf" +


Use variable for supporting.

davidbenjamin · 2018-05-24T03:18:13Z

...titute/hellbender/tools/walkers/variantutils/CalculateGenotypePosteriorsIntegrationTest.java

+                        " -supporting " + largeDir + "1000G.phase3.broad.withGenotypes.chr20.10100000.vcf" +
+                        " --" + StandardArgumentDefinitions.ADD_OUTPUT_VCF_COMMANDLINE +" false" +
+                        " -V " + dir + "NA12878.Jan2013.haplotypeCaller.subset.indels.vcf" +
+                        " --num-reference-samples-if-no-call 2500",


Use variable.

CalculateGenotypePosteriors now supports indels

…alid

ldgauthier · 2018-05-29T17:13:05Z

@davidbenjamin back to you

davidbenjamin · 2018-05-29T19:17:23Z

@ldgauthier Looks good to me. 👍

droazen reviewed May 21, 2018

View reviewed changes

droazen requested a review from davidbenjamin May 23, 2018 20:40

droazen assigned droazen and davidbenjamin May 23, 2018

davidbenjamin requested changes May 24, 2018

View reviewed changes

ldgauthier added 7 commits May 29, 2018 13:11

New HaplotypeCaller priors for variants sites and homRef sites

d8c9d1b

CalculateGenotypePosteriors now supports indels

Change tests to not bitwise equivalent

421ace5

Have GenotypeGVCFs remove PPs since removal of NON_REF makes them inv…

4474994

…alid

Turn off cast warning

7176386

Oh boy, I broke everything. Here's the fix.

a49e182

Oops

0e8772c

Address review comments

722afbb

ldgauthier force-pushed the ldg_HCwithPopPriors branch from 3547c12 to 722afbb Compare May 29, 2018 17:12

Merge branch 'master' into ldg_HCwithPopPriors

eb91616

ldgauthier merged commit ca2a4e3 into master May 30, 2018

ldgauthier deleted the ldg_HCwithPopPriors branch May 30, 2018 18:20

davidbenjamin mentioned this pull request May 31, 2018

PosteriorProbabilitiesUtils bug? #4760

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

New HaplotypeCaller priors for variants sites and homRef blocks #4793

New HaplotypeCaller priors for variants sites and homRef blocks #4793

ldgauthier commented May 21, 2018

droazen May 21, 2018

droazen May 24, 2018

droazen May 21, 2018

ldgauthier May 22, 2018

codecov-io commented May 23, 2018 •

edited

Loading

droazen commented May 23, 2018

davidbenjamin left a comment

davidbenjamin May 24, 2018

davidbenjamin May 24, 2018

davidbenjamin May 24, 2018

ldgauthier May 25, 2018

davidbenjamin May 27, 2018

davidbenjamin May 24, 2018

ldgauthier May 25, 2018

davidbenjamin May 27, 2018

davidbenjamin May 24, 2018

davidbenjamin May 24, 2018

ldgauthier May 25, 2018

davidbenjamin May 27, 2018

davidbenjamin May 24, 2018

davidbenjamin May 24, 2018

davidbenjamin May 24, 2018

davidbenjamin May 24, 2018

ldgauthier commented May 29, 2018

davidbenjamin commented May 29, 2018

New HaplotypeCaller priors for variants sites and homRef blocks #4793

New HaplotypeCaller priors for variants sites and homRef blocks #4793

Conversation

ldgauthier commented May 21, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

codecov-io commented May 23, 2018 • edited Loading

Codecov Report

droazen commented May 23, 2018

davidbenjamin left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ldgauthier commented May 29, 2018

davidbenjamin commented May 29, 2018

codecov-io commented May 23, 2018 •

edited

Loading