Genotyping code for the Gnarly Pipeline (gnomAD v3) #4947

ldgauthier · 2018-06-26T15:18:22Z

This along with the GVCF reblocking branch constitute the new code I'm using for gnomAD v3 on the Gnarly Pipeline.

Some of the GDB hacks are gross, but I can't clean it up until after the protobuf update.

codecov-io · 2018-06-26T16:44:48Z

Codecov Report

❗ No coverage uploaded for pull request base (master@39a9d13). Click here to learn what that means.
The diff coverage is 6.58%.

@@            Coverage Diff             @@
##             master     #4947   +/-   ##
==========================================
  Coverage          ?   13.338%           
  Complexity        ?      6396           
==========================================
  Files             ?      2016           
  Lines             ?    151745           
  Branches          ?     16269           
==========================================
  Hits              ?     20240           
  Misses            ?    129234           
  Partials          ?      2271

Impacted Files	Coverage Δ	Complexity Δ
...ute/hellbender/utils/variant/GATKVCFConstants.java	`50% <ø> (ø)`	`2 <0> (?)`
...iantutils/PosteriorProbabilitiesUtilsUnitTest.java	`4.386% <ø> (ø)`	`1 <0> (?)`
...te/hellbender/utils/variant/writers/TLODBlock.java	`0% <ø> (ø)`	`0 <0> (?)`
...notyper/GenotypeCalculationArgumentCollection.java	`100% <ø> (ø)`	`2 <0> (?)`
...ender/utils/variant/writers/GVCFBlockCombiner.java	`0% <0%> (ø)`	`0 <0> (?)`
...titute/hellbender/tools/walkers/GenotypeGVCFs.java	`0% <0%> (ø)`	`0 <0> (?)`
...rs/annotator/allelespecific/AS_StrandBiasTest.java	`5.085% <0%> (ø)`	`2 <0> (?)`
...ender/utils/variant/writers/TLODBlockUnitTest.java	`4.348% <0%> (ø)`	`2 <0> (?)`
...lkers/annotator/allelespecific/AS_QualByDepth.java	`1.351% <0%> (ø)`	`1 <0> (?)`
...ls/variant/writers/GVCFBlockCombiningIterator.java	`0% <0%> (ø)`	`0 <0> (?)`
... and 39 more

ldgauthier · 2018-07-07T15:53:41Z

TODO: refactor duplicated VC generation code from PosteriorProbabilitiesUtilsUnitTest in ReblockGVCFUnitTest by extracting to VariantContextTestUtils

droazen · 2018-08-24T19:04:42Z

@ldgauthier Should this branch be reviewed before the GVCF reblocking PR goes in?

ldgauthier · 2018-08-24T19:26:29Z

They may actually be more independent than two PRs both modifying SNP/indel data usually are. I think there will be VariantContextTestUtils conflicts either way, but the VCF data will probably be the same.

ldgauthier · 2018-11-05T18:00:40Z

@droazen can someone take a look soon? Eric tells me we need to start the CCDG 60K callset next week and I'd rather run that off a release.

droazen · 2018-11-05T18:34:47Z

@ldgauthier Looks like this branch is in conflict -- can you rebase before we review?

ldgauthier · 2018-11-09T21:50:32Z

@droazen After a gruesome rebase, conflicts are resolved.

droazen · 2018-11-14T20:14:22Z

@ldgauthier Since I'm about to leave for ~2 weeks, I've recruited @lbergelson to review this PR in my place.

src/main/java/org/broadinstitute/hellbender/engine/FeatureDataSource.java

ldgauthier · 2018-12-20T15:24:31Z

Conflicts resolved and David's comments addressed. @lbergelson can you take a look?

ldgauthier · 2019-01-10T16:22:21Z

@lbergelson I fixed the PL wiggle in the failing test. Think you could review before 4.1? Or maybe @droazen will look again? (Getting it in would be nice, but not critical.)

lbergelson

@ldgauthier I'm not quite done looking through this. I just want to checkpoint what I've reviewed so far.

src/main/java/org/broadinstitute/hellbender/tools/genomicsdb/GenomicsDBImport.java

lbergelson · 2019-01-24T15:52:24Z

src/main/java/org/broadinstitute/hellbender/engine/FeatureDataSource.java

+             new GenomicsDBOptions());
+    }
+
+    public FeatureDataSource(final FeatureInput<T> featureInput, final int queryLookaheadBases, final Class<? extends Feature> targetFeatureType,


src/main/java/org/broadinstitute/hellbender/tools/genomicsdb/GenomicsDBOptions.java

src/main/java/org/broadinstitute/hellbender/tools/walkers/GnarlyGenotyper.java

src/test/java/org/broadinstitute/hellbender/tools/walkers/GenotypeGVCFsIntegrationTest.java

...roadinstitute/hellbender/tools/walkers/variantutils/PosteriorProbabilitiesUtilsUnitTest.java

src/test/java/org/broadinstitute/hellbender/utils/GenotypeUtilsUnitTest.java

src/testUtils/java/org/broadinstitute/hellbender/testutils/VariantContextTestUtils.java

ldgauthier · 2019-03-04T20:39:16Z

That should work for both my cases. It could be nice for SelectVariants to be able to specify whether genotypes should be called or not too. Other tools might want the sites-only option.

…

On Mon, Mar 4, 2019 at 12:40 PM droazen ***@***.***> wrote: ***@***.**** commented on this pull request. ------------------------------ In src/main/java/org/broadinstitute/hellbender/tools/genomicsdb/GenomicsDBUtils.java <#4947 (comment)>: > @@ -40,7 +40,7 @@ */ public static GenomicsDBExportConfiguration.ExportConfiguration createExportConfiguration(final File reference, final String workspace, final String callsetJson, final String vidmapJson, - final String vcfHeader) { + final String vcfHeader, final boolean doGnarlyGenotyping) { @lbergelson <https://github.com/lbergelson> @ldgauthier <https://github.com/ldgauthier> If tools had a way to inject custom GDB config (eg., via an overridable method in GATKTool), and the engine used this config when creating the Feature Manager on startup, would that solve the problem here? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#4947 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AGRhdOOjGpZBu39mqk7jekA7iOzWDTFrks5vTVqFgaJpZM4U4KK0> .

-- Laura Doyle Gauthier, Ph.D. Associate Director, Germline Methods Data Sciences Platform [email protected] Broad Institute of MIT & Harvard 320 Charles St. Cambridge MA 0214

droazen · 2019-03-28T15:49:22Z

What's the status of this one @ldgauthier ? Do you need a new reviewer with @lbergelson out?

ldgauthier · 2019-03-28T18:36:17Z

Yes please. Where we left it I think Louis was happy, but we wanted to ask Nalini if she had any suggestions to avoid threading the argument for genotypes all the way through the engine.

…

On Thu, Mar 28, 2019 at 11:49 AM droazen ***@***.***> wrote: What's the status of this one @ldgauthier <https://github.com/ldgauthier> ? Do you need a new reviewer with @lbergelson <https://github.com/lbergelson> out? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub <#4947 (comment)>, or mute the thread <https://github.com/notifications/unsubscribe-auth/AGRhdJBn0j6n9Dfp-sz763M3mP5b1oyDks5vbOSHgaJpZM4U4KK0> .

-- Laura Doyle Gauthier, Ph.D. Associate Director, Germline Methods Data Sciences Platform [email protected] Broad Institute of MIT & Harvard 320 Charles St. Cambridge MA 0214

src/main/java/org/broadinstitute/hellbender/tools/genomicsdb/GenomicsDBUtils.java

nalinigans · 2019-03-29T17:46:06Z

Yes please. Where we left it I think Louis was happy, but we wanted to ask Nalini if she had any suggestions to avoid threading the argument for genotypes all the way through the engine.

Not sure we can avoid threading the argument for genotypes, but would using GenomicsDBOptions instead work?

droazen · 2019-06-05T20:24:13Z

@ldgauthier Let me know once you've had the chance to do that GnarlyGenotyperEngine refactoring discussed above, and I'd be happy to give this a (hopefully) final look.

…rge WGS cohort "Gnarly Pipeline" to finalize annotation values for filtering Super fast on Hail-derived sites-only input, needs to calculate ACs for GenomicsDB input Does NOT subset alternate alleles, just drops genotypes for sites with more than 6 alts Does NOT perform joint discovery, i.e. aggregating small amounts of evidence across samples Fix ReblockGVCF NPE on no-calls Some annotation cleanup New scheme to preserve de novo calling -- more blocks, but drop PLs, MIN_DP, floor GQs for each block Add AS combine operation and capability for GenomicsDB

count

ldgauthier · 2019-07-24T15:51:48Z

@droazen a few minor tests I missed, but this should be good to go otherwise.

I did the refactor for Evoquer and I think everything else has been addressed.

droazen · 2019-07-26T15:29:36Z

@ldgauthier Can you move spanDel.exome.chr20.vcf (as well as any other test files over about 1-2 MB) into a directory under src/test/resources/large so that they'll be stored in git lfs rather than check in to the repo?

droazen

Back to @ldgauthier with a final round of mostly-minor comments. We should be able to merge once they're addressed and tests pass.

droazen · 2019-07-26T15:56:11Z

src/main/java/org/broadinstitute/hellbender/engine/FeatureManager.java

@@ -205,7 +207,8 @@ private void initializeFeatureSources( final int featureQueryLookahead, final Co
            // Only create a data source for Feature arguments that were actually specified
            if ( featureInput != null ) {
                final Class<? extends Feature> featureType = getFeatureTypeForFeatureInputField(featureArgument.getKey());
-                addToFeatureSources(featureQueryLookahead, featureInput, featureType, cloudPrefetchBuffer, cloudIndexPrefetchBuffer, reference);
+                addToFeatureSources(featureQueryLookahead, featureInput, featureType, cloudPrefetchBuffer, cloudIndexPrefetchBuffer,
+                        toolInstance instanceof VariantWalker ? ((VariantWalker) toolInstance).getGenomicsDBOptions() : null);


We could probably avoid this ugly instanceof check by making GATKTool.getGenomicsDBOptions() return null (or new GenomicsDBOptions(referenceArguments.getReferencePath())) instead of throwing an exception, and then having GATKTool.initializeFeatures() (and its overrides) pass the GenomicsDB options in to the FeatureManager constructor, which can then propagate them down here.

Could be done in a separate PR if you don't want to do it here.

droazen · 2019-07-26T16:15:32Z

src/main/java/org/broadinstitute/hellbender/engine/FeatureManager.java

@@ -217,6 +220,13 @@ public void dumpAllFeatureCacheStats() {
        }
    }

+    void addToFeatureSources(final int featureQueryLookahead, final FeatureInput<? extends Feature> featureInput,


Add javadoc for this overload (can copy from below)

droazen · 2019-07-26T16:19:30Z

src/main/java/org/broadinstitute/hellbender/engine/GATKTool.java

+     * @return By default, not every GATK tool can read from a GenomicsDB -- child classes can override
+     */
+    protected GenomicsDBOptions getGenomicsDBOptions() {
+        throw new IllegalArgumentException("This tool does not take a GenomicsDB as a feature input.");


If you have this return either null or a new GenomicsDBOptions(referenceArguments.getReferencePath()) instead of throwing an exception, you could likely get rid of the VariantWalker special-casing in FeatureManager.

src/main/java/org/broadinstitute/hellbender/engine/VariantWalker.java

droazen · 2019-07-26T16:49:00Z

src/main/java/org/broadinstitute/hellbender/tools/evoquer/GnarlyGenotyperEngine.java

+    }
+
+    @SuppressWarnings({"unchecked", "rawtypes"})
+    public VariantContext finalizeGenotype(final VariantContext variant, final VariantContextBuilder annotationDBBuilder) {


As discussed in person, modify the method to tolerate a null annotationDBBuilder, and add a one-arg overload of finalizeGenotype() that takes just a VariantContext arg.

droazen · 2019-07-26T17:44:22Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/StrandBiasTest.java

@@ -78,7 +78,7 @@
        boolean foundData = false;

        for( final Genotype g : genotypes ) {
-            if( g.isNoCall() || ! g.hasAnyAttribute(GATKVCFConstants.STRAND_BIAS_BY_SAMPLE_KEY) ) {
+            if( ! g.hasAnyAttribute(GATKVCFConstants.STRAND_BIAS_BY_SAMPLE_KEY) ) {


You don't want to skip no-calls here?

...ava/org/broadinstitute/hellbender/tools/walkers/annotator/allelespecific/AS_QualByDepth.java

droazen · 2019-07-26T17:47:45Z

...oadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerArgumentCollection.java

+     */
+    @Advanced
+    @Argument(fullName=HaplotypeCallerArgumentCollection.OUTPUT_BLOCK_LOWER_BOUNDS, doc = "Output the band lower bound for each GQ block regardless of the data it represents", optional = true)
+    public boolean floorBlocks = false;


Is there a test covering this new HC arg?

Just in ReblockGVCF. Do you want one in HaplotypeCaller?

src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/ReblockGVCF.java

src/main/java/org/broadinstitute/hellbender/utils/GenotypeUtils.java

matching local vs. Travis

droazen · 2019-08-01T15:26:50Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/annotator/AnnotationUtils.java

+     * @param annotation the annotation to be tested
+     * @return true if the annotation is expected to have values per-allele
+     */
+    public static boolean isAlleleSpecific(final InfoFieldAnnotation annotation) {


There is a danger that if someone adds a new allele-specific annotation one day, they will not realize that they have to update this method as well. Perhaps we could introduce an empty AlleleSpecificAnnotation marker interface, and have just these 4 classes implement it? That would be much less likely to be missed in the future...

droazen · 2019-08-01T15:43:26Z

...oadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerArgumentCollection.java

+     * Output the band lower bound for each GQ block instead of the min GQ -- for better compression
+     */
+    @Advanced
+    @Argument(fullName=HaplotypeCallerArgumentCollection.OUTPUT_BLOCK_LOWER_BOUNDS, doc = "Output the band lower bound for each GQ block regardless of the data it represents", optional = true)


It would be good to cover this new HC arg with a quick/simple direct test, since the ReblockGVCF arg is separate.

droazen

Latest version looks good overall and tests pass -- I still have a lingering maintenance concern about AnnotationUtils.isAlleleSpecific() that could be resolved by adding an empty marker interface for AS annotations, and I'd still like to see a direct test covering that new HC arg.

droazen · 2019-08-01T16:00:19Z

After discussing with @ldgauthier, I'm going to approve this PR as-is, and Laura will address the remaining TODOs in a separate PR. For the record, the three remaining issues that need addressing are:

Get rid of the instanceof VariantWalker check in FeatureManager by making GATKTool.getGenomicsDBOptions() return null (or new GenomicsDBOptions(referenceArguments.getReferencePath())) instead of throwing an exception, and then having GATKTool.initializeFeatures() (and its overrides) pass the GenomicsDB options in to the FeatureManager constructor, which can then propagate them down here.
Add a simple direct integration test for the new --floor-blocks HaplotypeCaller arg
Address my maintenance concerns about AnnotationUtils.isAlleleSpecific() by adding an empty marker interface for AS annotations (open to other ideas here if you don't like that one)

lbergelson mentioned this pull request Jul 6, 2018

GVCF "reblocking" code to compress homRef blocks #4940

Merged

ldgauthier force-pushed the ldg_getVQSRinput branch from 0d143f6 to 5e3bf34 Compare August 22, 2018 19:16

droazen self-requested a review August 24, 2018 19:28

droazen self-assigned this Aug 24, 2018

ldgauthier force-pushed the ldg_getVQSRinput branch 3 times, most recently from b909de3 to 6b7d045 Compare November 8, 2018 16:31

droazen assigned lbergelson and unassigned droazen Nov 14, 2018

droazen requested review from lbergelson and removed request for droazen November 14, 2018 20:14

droazen reviewed Nov 14, 2018

View reviewed changes

src/main/java/org/broadinstitute/hellbender/engine/FeatureDataSource.java Outdated Show resolved Hide resolved

ldgauthier force-pushed the ldg_getVQSRinput branch from 337c33a to 8da4f3d Compare December 20, 2018 15:21

lbergelson reviewed Jan 25, 2019

View reviewed changes

ldgauthier force-pushed the ldg_getVQSRinput branch from 94e68ef to d283f3e Compare February 22, 2019 17:41

nalinigans reviewed Mar 29, 2019

View reviewed changes

src/main/java/org/broadinstitute/hellbender/tools/genomicsdb/GenomicsDBUtils.java Outdated Show resolved Hide resolved

droazen unassigned lbergelson May 23, 2019

ldgauthier force-pushed the ldg_getVQSRinput branch from 37c7482 to f098d6c Compare May 28, 2019 15:42

ldgauthier force-pushed the ldg_getVQSRinput branch from 789cc37 to a5b245f Compare July 19, 2019 18:11

ldgauthier added 4 commits July 22, 2019 10:21

Fail on warnings gets me every time

76f1681

Fixes for tests

ae1cd9b

Add test for annotation DB, expect RAW_GT_COUNT needs corrected homRef

c6eb337

count

Better annotation DB test

b1e1ce1

Last tests fixed(?) and large files moved

44d2b35

droazen suggested changes Jul 26, 2019

View reviewed changes

ldgauthier added 7 commits July 29, 2019 14:21

Address comments, fix tests

6e97901

Maybe adding index fixes tests???

5f62f35

Output all the annotations for the failing test (will break Hail test)

d75a69b

Fix expected output to match Travis -- fixTooHighQD random number isn't

db83c09

matching local vs. Travis

Reset the random number generator.... of course

4ef5198

Botched the vcf hatchet job

3a991d9

Merge branch 'master' into ldg_getVQSRinput

fbd8916

droazen reviewed Aug 1, 2019

View reviewed changes

droazen suggested changes Aug 1, 2019

View reviewed changes

droazen approved these changes Aug 1, 2019

View reviewed changes

droazen merged commit d0c8739 into master Aug 1, 2019

droazen deleted the ldg_getVQSRinput branch August 1, 2019 16:03

droazen mentioned this pull request Aug 1, 2019

Fix ReblockGVCF NPE on no-calls #5863

Closed

This was referenced Aug 2, 2019

GnarlyGenotyper tech debt #6075

Merged

Add ability to specify AS-annotation combine operations #5250

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Genotyping code for the Gnarly Pipeline (gnomAD v3) #4947

Genotyping code for the Gnarly Pipeline (gnomAD v3) #4947

ldgauthier commented Jun 26, 2018

codecov-io commented Jun 26, 2018 •

edited by codecov bot

Loading

ldgauthier commented Jul 7, 2018

droazen commented Aug 24, 2018

ldgauthier commented Aug 24, 2018

ldgauthier commented Nov 5, 2018

droazen commented Nov 5, 2018

ldgauthier commented Nov 9, 2018

droazen commented Nov 14, 2018

ldgauthier commented Dec 20, 2018

ldgauthier commented Jan 10, 2019

lbergelson left a comment

lbergelson Jan 24, 2019

ldgauthier commented Mar 4, 2019 via email

droazen commented Mar 28, 2019

ldgauthier commented Mar 28, 2019 via email

nalinigans commented Mar 29, 2019 •

edited

Loading

droazen commented Jun 5, 2019

ldgauthier commented Jul 24, 2019

droazen commented Jul 26, 2019

droazen left a comment

droazen Jul 26, 2019

droazen Jul 26, 2019

droazen Jul 26, 2019

droazen Jul 26, 2019

droazen Jul 26, 2019

droazen Jul 26, 2019

ldgauthier Jul 26, 2019

droazen Aug 1, 2019

droazen Aug 1, 2019

droazen left a comment

droazen commented Aug 1, 2019

Genotyping code for the Gnarly Pipeline (gnomAD v3) #4947

Genotyping code for the Gnarly Pipeline (gnomAD v3) #4947

Conversation

ldgauthier commented Jun 26, 2018

codecov-io commented Jun 26, 2018 • edited by codecov bot Loading

Codecov Report

ldgauthier commented Jul 7, 2018

droazen commented Aug 24, 2018

ldgauthier commented Aug 24, 2018

ldgauthier commented Nov 5, 2018

droazen commented Nov 5, 2018

ldgauthier commented Nov 9, 2018

droazen commented Nov 14, 2018

ldgauthier commented Dec 20, 2018

ldgauthier commented Jan 10, 2019

lbergelson left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ldgauthier commented Mar 4, 2019 via email

droazen commented Mar 28, 2019

ldgauthier commented Mar 28, 2019 via email

nalinigans commented Mar 29, 2019 • edited Loading

droazen commented Jun 5, 2019

ldgauthier commented Jul 24, 2019

droazen commented Jul 26, 2019

droazen 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

droazen left a comment

Choose a reason for hiding this comment

droazen commented Aug 1, 2019

codecov-io commented Jun 26, 2018 •

edited by codecov bot

Loading

nalinigans commented Mar 29, 2019 •

edited

Loading