AH - implement changes for mitochondrial pipeline #6399

ahaessly · 2020-01-21T18:12:08Z

This PR converts the Mutect2Filtering engine to be allele specific. This required changes to SomaticClusteringModel and ThresholdCalculator as well as ErrorProbabilities and of course the filters themselves. There are some filters which have not yet been converted, but I am prioritizing the ones in this PR for Sarah Calvo and the mitochondria pipeline. This provides the implementation for dsp-spec-ops tickets 166, 168, 169

gatk-bot · 2020-01-21T21:01:09Z

Travis reported job failures from build 28724
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
unit	openjdk11	28724.13	logs
unit	openjdk8	28724.3	logs

gatk-bot · 2020-01-21T21:36:43Z

Travis reported job failures from build 28733
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
unit	openjdk11	28733.13	logs
unit	openjdk8	28733.3	logs

gatk-bot · 2020-01-22T15:37:39Z

Travis reported job failures from build 28741
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
unit	openjdk11	28741.13	logs
unit	openjdk8	28741.3	logs

gatk-bot · 2020-01-22T15:47:34Z

Travis reported job failures from build 28743
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
unit	openjdk11	28743.13	logs
unit	openjdk8	28743.3	logs

gatk-bot · 2020-01-22T20:48:25Z

Travis reported job failures from build 28760
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
unit	openjdk11	28760.13	logs
unit	openjdk8	28760.3	logs

davidbenjamin

Finished reviewing first half of code.

...va/org/broadinstitute/hellbender/tools/walkers/annotator/allelespecific/StrandBiasUtils.java

...va/org/broadinstitute/hellbender/tools/walkers/mutect/clustering/SomaticClusteringModel.java

...n/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/ContaminationFilter.java

davidbenjamin · 2020-01-22T20:42:07Z

...va/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/DuplicatedAltReadFilter.java

@@ -20,11 +21,11 @@ public DuplicatedAltReadFilter(final int uniqueAltReadCount) {
    public ErrorType errorType() { return ErrorType.ARTIFACT; }

    @Override
-    public boolean isArtifact(final VariantContext vc, final Mutect2FilteringEngine filteringEngine) {
-        return vc.getAttributeAsInt(UniqueAltReadCount.KEY, 1) <= uniqueAltReadCount;
+    public List<Boolean> areAllelesArtifacts(final VariantContext vc, final Mutect2FilteringEngine filteringEngine, ReferenceContext referenceContext) {


This always returns a singleton list, even if there is more than one alt allele. . .

good catch. So, UniqueAltReadCount is actually assuming a single alt allele. I guess I should change that class to write an allele specific annotation that can be used in an allele specific way in the filter.

davidbenjamin · 2020-01-22T20:57:59Z

...in/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/ErrorProbabilities.java

-            trueProbability *= (1 - errorProb);
+        // if vc has symbolic allele, remove it
+        if (vc.hasSymbolicAlleles()) {
+            // can we assume it's the last allele?


The VCF spec doesn't require it, so let's be cautious and not assume.

...in/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/ErrorProbabilities.java

davidbenjamin

Back to you, @ahaessly. It's fundamentally sound.

...va/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/MinAlleleFractionFilter.java

...ava/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/Mutect2FilteringEngine.java

...n/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/TumorEvidenceFilter.java

src/main/java/org/broadinstitute/hellbender/utils/variant/GATKVCFConstants.java

davidbenjamin · 2020-01-23T03:13:37Z

src/test/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2IntegrationTest.java

                args -> {
                    intervals.stream().map(SimpleInterval::new).forEach(args::addInterval);
                    return args;
                });

+        // add tests for DUPLICATE


What's the status of this comment?

...est/resources/org/broadinstitute/hellbender/tools/mutect/mito/expected_LowHetNone_output.txt

...nstitute/hellbender/tools/mutect/mito/expected_LowHetVariantWalkerIntegrationTest_output.txt

gatk-bot · 2020-01-24T15:47:44Z

Travis reported job failures from build 28795
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
variantcalling	openjdk8	28795.4	logs
integration	oraclejdk8	28795.11	logs
unit	openjdk8	28795.3	logs
integration	openjdk8	28795.2	logs

gatk-bot · 2020-01-24T19:06:25Z

Travis reported job failures from build 28802
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
unit	openjdk8	28802.3	logs

ldgauthier

I reviewed the annotation aspects and left the Mutect2-specific parts to David. I pushed a commit with some changes I should have made a long time ago and a fix for the AFs in the multi-allelic splitting. Let me know what you think. (There are probably a lot more methods we could move from RefConfMerger to AlleleSubsettingUtils, but I just moved the ones I used here.)

ldgauthier · 2020-01-23T14:26:12Z

...ava/org/broadinstitute/hellbender/tools/walkers/annotator/AS_StrandBiasMutectAnnotation.java

@@ -0,0 +1,44 @@
+package org.broadinstitute.hellbender.tools.walkers.annotator;


I'm on the fence about putting this in the allelespecific subpackage. On the one hand, it is allele-specific, but on the other hand, all those classes have interfaces and/or parent classes that this one doesn't.

hmm, now that I look at that class again, I think I'll add the AlleleSpecificAnnotation interface to it. but, yeah, I wasn't really sure which package it belonged in.

ldgauthier · 2020-01-23T14:36:37Z

...va/org/broadinstitute/hellbender/tools/walkers/annotator/allelespecific/StrandBiasUtils.java

+
+    public static Map<String, Object> computeSBAnnotation(VariantContext vc, AlleleLikelihoods<GATKRead, Allele> likelihoods, String key) {
+        // calculate the annotation from the likelihoods
+        // likelihoods can come from HaplotypeCaller call to VariantAnnotatorEngine


Well, sometimes likelihoods come from Mutect2, right?

ldgauthier · 2020-01-31T15:32:30Z

...ava/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/Mutect2FilteringEngine.java

        return vcb.make();
    }

+    /**
+     * Creates a comma separated string of all the filters that apply to the allele.


Does it really create a comma separated string? Also, it looks like the only duplicate filter you're checking for is PASS, so I would call the method "removeExtraPassFilters" or something less general than what's there now.

thanks - didn't update the doc after a refactor.
the call to distinct removes any duplicate filters (of which the only possibility should be PASS). the first if statement removes the PASS if there are filters that did not pass.

ldgauthier · 2020-01-31T15:37:46Z

src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/NuMTFilter.java

+    }
+
+    @Override
+    protected List<String> requiredAnnotations() { return Collections.emptyList(); }


This does have required annotations, they're just not info annotations. How important is this check?

well the code that uses this sees if the Variant Context has the annotations - so this wouldn't work here because it's on the genotype. If there are no required annotations the check passes and the error probabilities are calculated. I could rename the method to requiredInfoAnnotations or requiredVariantAnnotations?

If the upstream method has a check for INFO annotations, then I'm happy.

ok - renaming the method to requiredInfoAnnotations so it's more explicit

...ava/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/StrictStrandBiasFilter.java

...roadinstitute/hellbender/tools/walkers/mutect/filtering/MTLowHeteroplasmyFilterToolTest.java

gatk-bot · 2020-02-03T22:51:07Z

Travis reported job failures from build 28960
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
unit	openjdk11	28960.13	logs
unit	openjdk8	28960.3	logs

davidbenjamin · 2020-02-12T20:59:16Z

@ahaessly Here's a representative line in the Mutect2IntegrationTest output:

20      38245275        .       TAC     T,TACAC .       germline;multiallelic;normal_artifact;slippage  AS_FilterStatus=.|weak_evidence|PASS

Comments:

Could we make AS_FilterStatus an alt-only field, since the ref isn't filtered?
I don't feel good about an AS_FilterStatus of PASS when the site fails the variant filters, especially because the next step (which I intend to do soon after your PR) is to make every filter allele-specific so that a filtered variant just means that every allele had a filter in common. I think of the allele filters as additive on top of the variant filters, so that a better representation would be to add "." If you want to leave this change for later that's fine with me.
I think some users have complained about INFO fields with underscores violating the VCF spec. Let's rename this to someone like ASFILTER ("Status" can be omitted).

Otherwise, looks like you have addressed everything.

ahaessly · 2020-02-12T21:49:55Z

@davidbenjamin AS_FilterStatus was already a field used in the INFO section of a VCF (@ldgauthier added it I believe). Also, it seems like there are several other INFO field attributes that use the _ if you take a look at GATKVCFConstants file. I'm fine changing it just wanted to hear from laura as well.

davidbenjamin · 2020-02-13T03:19:28Z

I double-checked the VCF spec and it turns out that while FORMAT keys must be alphanumeric, INFO keys may have underscores. I still prefer ASFILTER for the sake of brevity and hope the existing constant is not set in stone yet, though there may be compelling reasons not to change it. Let's let Laura weigh in.

ldgauthier · 2020-02-19T20:23:31Z

@davidbenjamin I'm willing to concede number 3, and for number 2 I grudgingly agree. A dot in the FILTER field means filters haven't been applied, so by extension a dot in the ASFILTER field would mean that filters haven't been applied to that allele, which isn't true. But in the autosomal pipeline, a variant that fails a site filter before VQSR wouldn't be seen or filtered by VQSR, which I suppose does mean that filters have not been applied. I guess my objection boils down to the fact that by outputting . instead of PASS you're throwing away information.

gatk-bot · 2020-02-20T15:39:46Z

Travis reported job failures from build 29219
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
variantcalling	openjdk8	29219.4	logs

gatk-bot · 2020-02-20T20:36:00Z

Travis reported job failures from build 29224
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
integration	oraclejdk8	29224.12	logs
unit	openjdk11	29224.14	logs
variantcalling	openjdk8	29224.4	logs
integration	openjdk11	29224.13	logs
integration	openjdk8	29224.2	logs

gatk-bot · 2020-02-21T20:33:49Z

Travis reported job failures from build 29232
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
variantcalling	openjdk8	29232.4	logs
unit	openjdk8	29232.3	logs

gatk-bot · 2020-02-21T22:36:57Z

Travis reported job failures from build 29236
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
variantcalling	openjdk8	29236.4	logs
unit	openjdk8	29236.3	logs

davidbenjamin · 2020-02-23T05:41:27Z

I guess my objection boils down to the fact that by outputting . Instead of PASS you're throwing away information.

I agree that we need to address that. What if an allele whose failing filters were all passed to the site got a concise tag like site to indicate that it has no additional filters, is not a passing allele, and that allele-specific filters were applied?

Besides the aesthetics of the word "PASS" being associated in any way with something failing, I'm worried about users shooting themselves in the foot with pipelines that only check the allele filters. You could imagine cases where failing to check the site filters wouldn't be so egregious as to be obvious.

ldgauthier · 2020-02-24T14:14:43Z

I could be onboard with site. It's not entirely intuitive, but it will get a header line, right?

davidbenjamin · 2020-02-25T04:36:58Z

It will get a header line.

ldgauthier · 2020-02-25T15:09:34Z

@ahaessly after reviewing the allele-specific filtering discussion in hts-specs (samtools/hts-specs#115), I think the delimiters for the ASFILTER values should be switched, i.e. data for each allele should be separated by commas (ALLELE_SPECIFIC_PRINT_DELIM?) and the filters for each allele should be separated by pipes (ALLELE_SPECIFIC_RAW_DELIM).

meganshand

The WDLs look good to me (that's all I looked at). I learned some new tricks! Only have a couple of minor comments/questions. Thanks!

meganshand · 2020-03-25T20:05:55Z

scripts/mitochondria_m2_wdl/AlignAndCall.wdl

@@ -1,6 +1,6 @@
 version 1.0

-import "AlignmentPipeline.wdl" as AlignAndMarkDuplicates
+import "https://api.firecloud.org/ga4gh/v1/tools/mitochondria:AlignmentPipeline/versions/1/plain-WDL/descriptor" as AlignAndMarkDuplicates


I think for the tests to work this needs to be the file name "AlignmentPipeline.wdl". If there's a better solution for moving things back and forth between Terra and the gatk repo though I'm all for it.

scripts/mitochondria_m2_wdl/AlignAndCall.wdl

meganshand · 2020-03-25T20:21:46Z

scripts/mitochondria_m2_wdl/AlignAndCall.wdl

+  awk '{print $6}' output-data > major_hg.txt
+  awk '{print $8}' output-data > minor_hg.txt
+  awk '{print $14}' output-data > mean_het_major.txt
+  awk '{print $15}' output-data > mean_het_minor.txt


Do you think it's worth adding any sanity checks here to make sure this doesn't break in the future since it's relying on column order to extract the values?

scripts/mitochondria_m2_wdl/AlignAndCall.wdl

gatk-bot · 2020-04-01T18:39:42Z

Travis reported job failures from build 29837
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
integration	oraclejdk8	29837.12	logs
integration	openjdk11	29837.13	logs
integration	openjdk8	29837.2	logs

gatk-bot · 2020-04-01T18:48:00Z

Travis reported job failures from build 29839
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
integration	oraclejdk8	29839.12	logs
integration	openjdk11	29839.13	logs
integration	openjdk8	29839.2	logs

gatk-bot · 2020-04-01T21:58:08Z

Travis reported job failures from build 29842
Failures in the following jobs:

Test Type	JDK	Job ID	Logs
integration	oraclejdk8	29842.12	logs
integration	openjdk11	29842.13	logs
integration	openjdk8	29842.2	logs

meganshand · 2020-04-08T14:42:35Z

scripts/mitochondria_m2_wdl/AlignAndCall.wdl

@@ -470,7 +498,7 @@ task M2 {
        --max-mnp-distance 0
  >>>
  runtime {
-      docker: "us.gcr.io/broad-gatk/gatk:4.1.1.0"
+      docker: select_first([gatk_docker_override, "us.gcr.io/broad-gatk/gatk:4.1.1.0"])


I'm not sure what the solution here is, but does keeping this as the default version mean that we'll have to keep updating it when there are new releases? That's the way it already was, so it's not necessarily a problem, I just want to remember that once this PR is merged and released we should go back and update this default.

meganshand

I had one small question/comment but otherwise the WDL and WDL tests look good to me! Thanks for adding them.

davidbenjamin

@ahaessly Looks like you have resolved everything. Don't bother sending it back to me on account my one trivial comment, which you can take or leave at your discretion.

davidbenjamin · 2020-04-12T06:16:12Z

...va/org/broadinstitute/hellbender/tools/walkers/mutect/clustering/SomaticClusteringModel.java

+        // set tumorAD to 0 for symbolic alleles so it won't contribute to overall AD
+        List<Allele> symbolicAlleles = vc.getAlternateAlleles().stream().filter(allele -> allele.isSymbolic()).collect(Collectors.toList());
+        // convert allele index to alt allele index
+        List<Integer> symIndexes = vc.getAlleleIndices(symbolicAlleles).stream().map(i -> i-1).collect(Collectors.toList());


Can you go directly compute symbolic indices via new IndexRange(0, vc.getNAlleles()).filter(n -> vc.getAllele(n).isSymbolic());?

ldgauthier · 2020-04-21T18:31:53Z

...n/java/org/broadinstitute/hellbender/tools/walkers/mutect/filtering/Mutect2AlleleFilter.java

@@ -73,7 +73,7 @@


    /**
-     *
+     * Returns a list for each alternate allele which is the probability that the allele should be filtered out.


"a list with one element for each allele"?

new refactor adding allele specific filters in mutect2 fix headers * Fixed AF and SB splitting; also some javadoc I should have done in the past and some refactoring I should have done in the past * make unique alt read count allele specific * fix genotypes not included in vcf * changed getRequiredAnnotations to getRequiredInfoAnnotations to be more explicit * fix split multiallelics to work for all info fields * change count type for RPA * update example inputs to use broad public bucket for refs Co-authored-by: Laura Gauthier <[email protected]>

ahaessly requested review from davidbenjamin and ldgauthier January 22, 2020 14:55

ahaessly force-pushed the ah_SO166_morefilters branch from f3ba57b to 959a4a2 Compare January 22, 2020 15:07

davidbenjamin requested changes Jan 22, 2020

View reviewed changes

davidbenjamin requested changes Jan 23, 2020

View reviewed changes

davidbenjamin self-assigned this Jan 23, 2020

ahaessly force-pushed the ah_SO166_morefilters branch from ecdbbe7 to 1946785 Compare January 24, 2020 15:12

ldgauthier requested changes Jan 31, 2020

View reviewed changes

davidbenjamin mentioned this pull request Feb 6, 2020

Ah so166 asfilter #6290

Closed

ahaessly force-pushed the ah_SO166_morefilters branch from 3f2937e to f07b809 Compare February 19, 2020 22:37

ahaessly requested a review from meganshand March 25, 2020 19:48

meganshand reviewed Mar 25, 2020

View reviewed changes

ahaessly added 4 commits March 30, 2020 15:16

cleanup

272f6b8

updating wdl

35ab284

fix bug that removes existing filters

0dfab1d

add one more test

9b8f211

fix overwrite of filter

3b43c10

ahaessly added 7 commits April 2, 2020 14:14

add missing test files

8191f3c

should fix wdl tests

54b7fc1

fix wdl tests

37c8c81

fix path

4dcd088

fix path again

d6b1229

remove blacklisted shifted

e60ed86

add error checking for contam file format

717d87f

meganshand reviewed Apr 8, 2020

View reviewed changes

meganshand approved these changes Apr 8, 2020

View reviewed changes

davidbenjamin approved these changes Apr 12, 2020

View reviewed changes

ahaessly added 3 commits April 14, 2020 10:47

minor refactor

69cca57

update example inputs to use broad public bucket for refs

f51ba10

doc updates

59f3dc9

ldgauthier reviewed Apr 21, 2020

View reviewed changes

ldgauthier approved these changes Apr 21, 2020

View reviewed changes

update documentation

33c2942

ahaessly merged commit 3021e69 into master Apr 22, 2020

ahaessly deleted the ah_SO166_morefilters branch April 22, 2020 17:41

yl-h mentioned this pull request Aug 6, 2020

Default Mutect2 INFO field discrepancies between version 4.1.5.0 and 4.1.8.0 #6721

Open

2 tasks

		@@ -0,0 +1,44 @@
		package org.broadinstitute.hellbender.tools.walkers.annotator;

AH - implement changes for mitochondrial pipeline #6399

AH - implement changes for mitochondrial pipeline #6399

Conversation

ahaessly commented Jan 21, 2020 • edited Loading

gatk-bot commented Jan 21, 2020 • edited Loading

gatk-bot commented Jan 21, 2020 • edited Loading

gatk-bot commented Jan 22, 2020 • edited Loading

gatk-bot commented Jan 22, 2020 • edited Loading

gatk-bot commented Jan 22, 2020 • edited Loading

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

davidbenjamin left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

gatk-bot commented Jan 24, 2020 • edited Loading

gatk-bot commented Jan 24, 2020

ldgauthier 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

gatk-bot commented Feb 3, 2020 • edited Loading

davidbenjamin commented Feb 12, 2020

ahaessly commented Feb 12, 2020

davidbenjamin commented Feb 13, 2020

ldgauthier commented Feb 19, 2020

gatk-bot commented Feb 20, 2020

gatk-bot commented Feb 20, 2020 • edited Loading

gatk-bot commented Feb 21, 2020 • edited Loading

gatk-bot commented Feb 21, 2020 • edited Loading

davidbenjamin commented Feb 23, 2020 • edited Loading

ldgauthier commented Feb 24, 2020

davidbenjamin commented Feb 25, 2020

ldgauthier commented Feb 25, 2020

meganshand left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

gatk-bot commented Apr 1, 2020 • edited Loading

gatk-bot commented Apr 1, 2020 • edited Loading

gatk-bot commented Apr 1, 2020 • edited Loading

Choose a reason for hiding this comment

meganshand left a comment

Choose a reason for hiding this comment

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

ahaessly commented Jan 21, 2020 •

edited

Loading

gatk-bot commented Jan 21, 2020 •

edited

Loading

gatk-bot commented Jan 21, 2020 •

edited

Loading

gatk-bot commented Jan 22, 2020 •

edited

Loading

gatk-bot commented Jan 22, 2020 •

edited

Loading

gatk-bot commented Jan 22, 2020 •

edited

Loading

gatk-bot commented Jan 24, 2020 •

edited

Loading

gatk-bot commented Feb 3, 2020 •

edited

Loading

gatk-bot commented Feb 20, 2020 •

edited

Loading

gatk-bot commented Feb 21, 2020 •

edited

Loading

gatk-bot commented Feb 21, 2020 •

edited

Loading

davidbenjamin commented Feb 23, 2020 •

edited

Loading

gatk-bot commented Apr 1, 2020 •

edited

Loading

gatk-bot commented Apr 1, 2020 •

edited

Loading

gatk-bot commented Apr 1, 2020 •

edited

Loading