From 3f085296f16ec71c9a1066b7c98b8703dfb76a99 Mon Sep 17 00:00:00 2001
From: David Benjamin <davidben@broadinstitute.org>
Date: Wed, 1 Aug 2018 12:05:54 -0400
Subject: [PATCH] handle overlapping mates in M2 isActive

---
 .../tools/walkers/mutect/M2ArgumentCollection.java   |  7 +++++++
 .../tools/walkers/mutect/Mutect2Engine.java          | 12 ++++++++----
 2 files changed, 15 insertions(+), 4 deletions(-)

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 5abadc28ace..f12c9a23172 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
@@ -23,6 +23,7 @@ public class M2ArgumentCollection extends AssemblyBasedCallerArgumentCollection
     public static final String EMISSION_LOG_SHORT_NAME = "emit-lod";
     public static final String INITIAL_TUMOR_LOD_LONG_NAME = "initial-tumor-lod";
     public static final String INITIAL_TUMOR_LOD_SHORT_NAME = "init-lod";
+    public static final String INITIAL_PCR_ERROR_QUAL = "initial-pcr-qual";
     public static final String MAX_POPULATION_AF_LONG_NAME = "max-population-af";
     public static final String MAX_POPULATION_AF_SHORT_NAME = "max-af";
     public static final String DOWNSAMPLING_STRIDE_LONG_NAME = "downsampling-stride";
@@ -108,6 +109,12 @@ public double getDefaultAlleleFrequency() {
     @Argument(fullName = INITIAL_TUMOR_LOD_LONG_NAME, shortName = INITIAL_TUMOR_LOD_SHORT_NAME, optional = true, doc = "LOD threshold to consider pileup active.")
     public double initialTumorLod = 2.0;
 
+    /**
+     * PCR error rate for overlapping fragments in isActive()
+     */
+    @Argument(fullName = INITIAL_PCR_ERROR_QUAL, optional = true, doc = "PCR error rate for overlapping fragments in isActive()")
+    public int initialPCRErrorQual = 40;
+
     /**
      * 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 ddf2737af6f..5013b298d86 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
@@ -259,14 +259,14 @@ public ActivityProfileState isActive(final AlignmentContext context, final Refer
 
         final ReadPileup pileup = context.getBasePileup();
         final ReadPileup tumorPileup = pileup.getPileupForSample(tumorSample, header);
-        final List<Byte> tumorAltQuals = altQuals(tumorPileup, refBase);
+        final List<Byte> tumorAltQuals = altQuals(tumorPileup, refBase, MTAC.initialPCRErrorQual);
         final double tumorLog10Odds = MathUtils.logToLog10(lnLikelihoodRatio(tumorPileup.size()-tumorAltQuals.size(), tumorAltQuals));
 
         if (tumorLog10Odds < MTAC.initialTumorLod) {
             return new ActivityProfileState(refInterval, 0.0);
         } else if (hasNormal() && !MTAC.genotypeGermlineSites) {
             final ReadPileup normalPileup = pileup.getPileupForSample(normalSample, header);
-            final List<Byte> normalAltQuals = altQuals(normalPileup, refBase);
+            final List<Byte> normalAltQuals = altQuals(normalPileup, refBase, MTAC.initialPCRErrorQual);
             final int normalAltCount = normalAltQuals.size();
             final double normalQualSum = normalAltQuals.stream().mapToDouble(Byte::doubleValue).sum();
             if (normalAltCount > normalPileup.size() * MAX_ALT_FRACTION_IN_NORMAL && normalQualSum > MAX_NORMAL_QUAL_SUM) {
@@ -297,8 +297,9 @@ private static byte indelQual(final int indelLength) {
         return (byte) Math.min(INDEL_START_QUAL + (indelLength - 1) * INDEL_CONTINUATION_QUAL, Byte.MAX_VALUE);
     }
 
-    private static List<Byte> altQuals(final ReadPileup pileup, final byte refBase) {
+    private static List<Byte> altQuals(final ReadPileup pileup, final byte refBase, final int pcrErrorQual) {
         final List<Byte> result = new ArrayList<>();
+        final int position = pileup.getLocation().getStart();
 
         for (final PileupElement pe : pileup) {
             final int indelLength = getCurrentOrFollowingIndelLength(pe);
@@ -307,7 +308,10 @@ private static List<Byte> altQuals(final ReadPileup pileup, final byte refBase)
             } else if (isNextToUsefulSoftClip(pe)) {
                 result.add(indelQual(1));
             } else if (pe.getBase() != refBase && pe.getQual() > MINIMUM_BASE_QUALITY) {
-                result.add(pe.getQual());
+                final GATKRead read = pe.getRead();
+                final int mateStart = read.mateIsUnmapped() ? Integer.MAX_VALUE : read.getMateStart();
+                final boolean overlapsMate = mateStart <= position && position < mateStart + read.getLength();
+                result.add(overlapsMate ? (byte) FastMath.min(pe.getQual(), pcrErrorQual/2) : pe.getQual());
             }
         }