diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/M2ArgumentCollection.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/M2ArgumentCollection.java index 7de109846aa..d25c7a33fc5 100644 --- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/M2ArgumentCollection.java +++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/M2ArgumentCollection.java @@ -9,7 +9,6 @@ import org.broadinstitute.hellbender.cmdline.ReadFilterArgumentDefinitions; import org.broadinstitute.hellbender.engine.FeatureInput; import org.broadinstitute.hellbender.tools.walkers.haplotypecaller.FlowBasedAlignmentArgumentCollection; -import org.broadinstitute.hellbender.tools.FlowBasedArgumentCollection; import org.broadinstitute.hellbender.tools.walkers.haplotypecaller.*; import org.broadinstitute.hellbender.tools.walkers.haplotypecaller.readthreading.ReadThreadingAssembler; import org.broadinstitute.hellbender.tools.walkers.mutect.filtering.FilterMutectCalls; @@ -50,6 +49,7 @@ public class M2ArgumentCollection extends AssemblyBasedCallerArgumentCollection public static final String CALLABLE_DEPTH_LONG_NAME = "callable-depth"; public static final String PCR_SNV_QUAL_LONG_NAME = "pcr-snv-qual"; public static final String PCR_INDEL_QUAL_LONG_NAME = "pcr-indel-qual"; + public static final String MULTIPLE_SUBSTITUTION_BASE_QUAL_CORRECTION = "base-qual-correction-factor"; public static final String F1R2_TAR_GZ_NAME = "f1r2-tar-gz"; public static final double DEFAULT_AF_FOR_TUMOR_ONLY_CALLING = 5e-8; @@ -249,6 +249,18 @@ public double getInitialLogOdds() { @Argument(fullName = PCR_INDEL_QUAL_LONG_NAME, optional = true, doc = "Phred-scaled PCR indel qual for overlapping fragments") public int pcrIndelQual = 40; + /** + * A scale factor to modify the base qualities reported by the sequencer and used in the Mutect2 substitution error model. + * Set to zero to turn off the error model changes included in GATK 4.1.9.0. + * Our pileup likelihoods models assume that the base quality (qual) corresponds to the probability that a ref base is misread + * as the *particular* alt base, whereas the qual actually means the probability of *any* substitution error. + * Since there are three possible substitutions for each ref base we must divide the error probability by three + * which corresponds to adding 10*log10(3) = 4.77 ~ 5 to the qual. + */ + @Advanced + @Argument(fullName = MULTIPLE_SUBSTITUTION_BASE_QUAL_CORRECTION, optional = true, doc = "Set to zero to turn off the error model changes included in GATK 4.1.9.0.") + public int activeRegionMultipleSubstitutionBaseQualCorrection = (int)Math.round(10 * Math.log10(3)); + /** * In tumor-only mode, we discard variants with population allele frequencies greater than this threshold. */ diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2Engine.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2Engine.java index 9f0bb90345c..1d2096f3c2b 100644 --- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2Engine.java +++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2Engine.java @@ -122,8 +122,8 @@ public final class Mutect2Engine implements AssemblyRegionEvaluator, AutoCloseab private final Optional f1R2CountsCollector; - private PileupQualBuffer tumorPileupQualBuffer = new PileupQualBuffer(); - private PileupQualBuffer normalPileupQualBuffer = new PileupQualBuffer(); + private PileupQualBuffer tumorPileupQualBuffer; + private PileupQualBuffer normalPileupQualBuffer; /** * Create and initialize a new HaplotypeCallerEngine given a collection of HaplotypeCaller arguments, a reads header, @@ -144,6 +144,9 @@ public Mutect2Engine(final M2ArgumentCollection MTAC, AssemblyRegionArgumentColl aligner = SmithWatermanAligner.getAligner(MTAC.smithWatermanImplementation); samplesList = new IndexedSampleList(new ArrayList<>(ReadUtils.getSamplesFromHeader(header))); + tumorPileupQualBuffer = new PileupQualBuffer(MTAC.activeRegionMultipleSubstitutionBaseQualCorrection); + normalPileupQualBuffer = new PileupQualBuffer(MTAC.activeRegionMultipleSubstitutionBaseQualCorrection); + // optimize set operations for the common cases of no normal and one normal if (MTAC.normalSamples.isEmpty()) { normalSamples = Collections.emptySet(); @@ -671,16 +674,14 @@ private static class PileupQualBuffer { private static final int OTHER_SUBSTITUTION = 4; private static final int INDEL = 5; - // our pileup likelihoods models assume that the qual corresponds to the probability that a ref base is misread - // as the *particular* alt base, whereas the qual actually means the probability of *any* substitution error. - // since there are three possible substitutions for each ref base we must divide the error probability by three - // which corresponds to adding 10*log10(3) = 4.77 ~ 5 to the qual. - private static final int ONE_THIRD_QUAL_CORRECTION = 5; + private static double MULTIPLE_SUBSTITUTION_BASE_QUAL_CORRECTION; // indices 0-3 are A,C,G,T; 4 is other substitution (just in case it's some exotic protocol); 5 is indel private List buffers = IntStream.range(0,6).mapToObj(n -> new ByteArrayList()).collect(Collectors.toList()); - public PileupQualBuffer() { } + public PileupQualBuffer(final double multipleSubstitutionQualCorrection) { + MULTIPLE_SUBSTITUTION_BASE_QUAL_CORRECTION = multipleSubstitutionQualCorrection; + } public void accumulateQuals(final ReadPileup pileup, final byte refBase, final int pcrErrorQual) { clear(); @@ -731,7 +732,7 @@ private void accumulateSubstitution(final byte base, final byte qual) { if (index == -1) { // -1 is the hard-coded value for non-simple bases in BaseUtils buffers.get(OTHER_SUBSTITUTION).add(qual); } else { - buffers.get(index).add((byte) FastMath.min(qual + ONE_THIRD_QUAL_CORRECTION, QualityUtils.MAX_QUAL)); + buffers.get(index).add((byte) FastMath.min(qual + MULTIPLE_SUBSTITUTION_BASE_QUAL_CORRECTION, QualityUtils.MAX_QUAL)); } } diff --git a/src/test/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2EngineUnitTest.java b/src/test/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2EngineUnitTest.java index 18e2915897f..cb24556e4b6 100644 --- a/src/test/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2EngineUnitTest.java +++ b/src/test/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2EngineUnitTest.java @@ -119,9 +119,6 @@ public void testSmallNumAltExact(final int numRef, final double errorRate) { Assert.assertEquals(calculated, expected, precision); } - - - } @DataProvider(name = "fewAltData")