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.

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

Fix Mutect2 IndexOutOfBoundException with germline resource #7979

Merged
merged 4 commits into from
Aug 9, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
Original file line number Diff line number Diff line change
Expand Up @@ -490,10 +490,9 @@ public ActivityProfileState isActive(final AlignmentContext context, final Refer

for (final VariantContext germlineVC : germline) {
final List<Double> germlineAlleleFrequencies = getAttributeAsDoubleList(germlineVC, VCFConstants.ALLELE_FREQUENCY_KEY, 0.0);
final List<Allele> germlineAlts = germlineVC.getAlternateAlleles();
final Allele germlineRef = germlineVC.getReference();

for (int germlineAltIdx = 0; germlineAltIdx < germlineAlts.size(); germlineAltIdx++) {
for (int germlineAltIdx = 0; germlineAltIdx < germlineAlleleFrequencies.size(); germlineAltIdx++) {
if (germlineAlleleFrequencies.get(germlineAltIdx) < MTAC.maxPopulationAlleleFrequency) {
continue;
}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,8 @@
import org.apache.commons.math3.linear.DefaultRealMatrixChangingVisitor;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.util.FastMath;
import org.apache.logging.log4j.LogManager;
import org.apache.logging.log4j.Logger;
import org.broadinstitute.hellbender.engine.FeatureContext;
import org.broadinstitute.hellbender.engine.ReferenceContext;
import org.broadinstitute.hellbender.tools.walkers.annotator.VariantAnnotatorEngine;
Expand All @@ -34,6 +36,7 @@
import java.util.stream.IntStream;

public class SomaticGenotypingEngine implements AutoCloseable {
protected static final Logger logger = LogManager.getLogger(SomaticGenotypingEngine.class);

private final M2ArgumentCollection MTAC;
private final Set<String> normalSamples;
Expand Down Expand Up @@ -372,6 +375,10 @@ private static Map<String, Object> getNegativeLogPopulationAFAnnotation(List<Var
static double[] getGermlineAltAlleleFrequencies(final List<Allele> allAlleles, final Optional<VariantContext> germlineVC, final double afOfAllelesNotInGermlineResource) {
Utils.validateArg(!allAlleles.isEmpty(), "allAlleles are empty -- there is not even a reference allele.");
if (germlineVC.isPresent()) {
if (! germlineVC.get().hasAttribute(VCFConstants.ALLELE_FREQUENCY_KEY)) {
logger.warn("Germline resource variant at " + germlineVC.get().getContig() + ":" + germlineVC.get().getStart() +" missing AF attribute");
return Doubles.toArray(Collections.nCopies(allAlleles.size() - 1, afOfAllelesNotInGermlineResource));
}
List<OptionalInt> germlineIndices = GATKVariantContextUtils.alleleIndices(allAlleles, germlineVC.get().getAlleles());
final List<Double> germlineAltAFs = Mutect2Engine.getAttributeAsDoubleList(germlineVC.get(), VCFConstants.ALLELE_FREQUENCY_KEY, afOfAllelesNotInGermlineResource);

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ Object[][] germlineAltAlleleFrequenciesData() {

//biallelic, same alt allele in different representations
{ Arrays.asList(Allele.REF_A, Allele.ALT_G), Arrays.asList(ATref, GT), new double[] {0.1}, new double[] {0.1}},
{ Arrays.asList(ATref, GT), Arrays.asList(Allele.REF_A, Allele.ALT_G), new double[] {0.1}, new double[] {0.1}},
{ Arrays.asList(ATref, GT), Arrays.asList(Allele.REF_A, Allele.ALT_G), new double[] {0.1}, new double[] {0.1}}
};
}

Expand All @@ -56,4 +56,20 @@ public void testGetGermlineAltAlleleFrequencies(final List<Allele> calledAlleles
final double[] result = SomaticGenotypingEngine.getGermlineAltAlleleFrequencies(calledAlleles, Optional.of(vc1), DEFAULT_AF);
Assert.assertEquals(result, expected, 1.0e-10);
}

@DataProvider(name = "missingAFData")
Object[][] missingAFData() {
return new Object[][]{
{new VariantContextBuilder("SOURCE", "1", 1, 1, Arrays.asList(Allele.REF_A, Allele.ALT_C))
.attribute(VCFConstants.ALLELE_FREQUENCY_KEY, VCFConstants.MISSING_VALUE_v4).make()},
{new VariantContextBuilder("SOURCE", "1", 1, 1, Arrays.asList(Allele.REF_A, Allele.ALT_C))
.make()}
};
}

@Test(dataProvider = "missingAFData")
public void testGetGermlineAltAlleleFrequenciesWithMissingAF(final VariantContext vc) {
final double[] result = SomaticGenotypingEngine.getGermlineAltAlleleFrequencies(vc.getAlleles(), Optional.of(vc), DEFAULT_AF);
Assert.assertEquals(result, new double[] {DEFAULT_AF}, 1.0e-10);
}
}