Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

mutect2 expose param #8447

Merged
merged 4 commits into from
Aug 14, 2023
Merged

mutect2 expose param #8447

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
fdbaf98
Allow users to turn off fix from 4.1.9.0 for consistency's sake
ldgauthier Aug 1, 2023
2f9bb34
David R says advanced args don't need tests
ldgauthier Aug 4, 2023
44231ca
Implement rename/refactor suggestion
ldgauthier Aug 10, 2023
1c3b062
Remove constant from engine and move documentation to arg colelction
ldgauthier Aug 10, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 9 additions & 1 deletion src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/M2ArgumentCollection.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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 USER_DEFINED_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;
Expand Down Expand Up @@ -249,6 +249,14 @@ 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.
*/
@Advanced
@Argument(fullName = USER_DEFINED_BASE_QUAL_CORRECTION, optional = true, doc = "Set to zero to turn off the error model changes included in GATK 4.1.9.0.")
public double userDefinedBaseQualCorrection = 1;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I lean toward making this private int activeRegionMultipleSubstitutionBaseQualCorrection = Math.round(10 * Math.log(3)) (maybe make the name a bit more graceful than that). That way 1) setting it to zero is inherently the desired behavior of turning off the change without requiring interpretation as a sentinel value, 2) it eliminates the indirection of the user-defined multiplier and goes straight to the desired variable, 3) it makes the variable an int from the beginning.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think I get what you're saying, but in that case it would be public since it's an arg right?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For sure, I'm not sure why I wrote private.


/**
* In tumor-only mode, we discard variants with population allele frequencies greater than this threshold.
*/
Expand Down
15 changes: 12 additions & 3 deletions src/main/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2Engine.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -122,8 +122,8 @@ public final class Mutect2Engine implements AssemblyRegionEvaluator, AutoCloseab

private final Optional<F1R2CountsCollector> 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,
Expand All @@ -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.userDefinedBaseQualCorrection);
normalPileupQualBuffer = new PileupQualBuffer(MTAC.userDefinedBaseQualCorrection);

// optimize set operations for the common cases of no normal and one normal
if (MTAC.normalSamples.isEmpty()) {
normalSamples = Collections.emptySet();
Expand Down Expand Up @@ -677,11 +680,17 @@ private static class PileupQualBuffer {
// 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 USER_DEFINED_QUAL_CORRECTION = 1;

// 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<ByteArrayList> buffers = IntStream.range(0,6).mapToObj(n -> new ByteArrayList()).collect(Collectors.toList());

public PileupQualBuffer() { }

public PileupQualBuffer(final double userDefinedQualCorrection) {
USER_DEFINED_QUAL_CORRECTION = userDefinedQualCorrection;
}

public void accumulateQuals(final ReadPileup pileup, final byte refBase, final int pcrErrorQual) {
clear();
final int position = pileup.getLocation().getStart();
Expand Down Expand Up @@ -731,7 +740,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 + ONE_THIRD_QUAL_CORRECTION * USER_DEFINED_QUAL_CORRECTION, QualityUtils.MAX_QUAL));
}
}

Expand Down
4 changes: 3 additions & 1 deletion src/test/java/org/broadinstitute/hellbender/tools/walkers/mutect/Mutect2EngineUnitTest.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,8 +119,10 @@ public void testSmallNumAltExact(final int numRef, final double errorRate) {

Assert.assertEquals(calculated, expected, precision);
}
}


@Test
public void testOldErrorModelBehavior() {
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Delete


}

Expand Down