use of htsjdk.variant.variantcontext.Genotype in project gatk-protected by broadinstitute.
the class Mutect2FilteringEngine method applyStrandArtifactFilter.
private void applyStrandArtifactFilter(final M2FiltersArgumentCollection MTFAC, final VariantContext vc, final Collection<String> filters) {
Genotype tumorGenotype = vc.getGenotype(tumorSample);
final double[] posteriorProbabilities = GATKProtectedVariantContextUtils.getAttributeAsDoubleArray(tumorGenotype, (StrandArtifact.POSTERIOR_PROBABILITIES_KEY), () -> null, -1);
final double[] mapAlleleFractionEstimates = GATKProtectedVariantContextUtils.getAttributeAsDoubleArray(tumorGenotype, (StrandArtifact.MAP_ALLELE_FRACTIONS_KEY), () -> null, -1);
if (posteriorProbabilities == null || mapAlleleFractionEstimates == null) {
return;
}
final int maxZIndex = MathUtils.maxElementIndex(posteriorProbabilities);
if (maxZIndex == StrandArtifact.StrandArtifactZ.NO_ARTIFACT.ordinal()) {
return;
}
if (posteriorProbabilities[maxZIndex] > MTFAC.STRAND_ARTIFACT_POSTERIOR_PROB_THRESHOLD && mapAlleleFractionEstimates[maxZIndex] < MTFAC.STRAND_ARTIFACT_ALLELE_FRACTION_THRESHOLD) {
filters.add(GATKVCFConstants.STRAND_ARTIFACT_FILTER_NAME);
}
}
use of htsjdk.variant.variantcontext.Genotype in project gatk-protected by broadinstitute.
the class XHMMSegmentGenotyperIntegrationTest method assertVariantsInfoFieldsAreConsistent.
private void assertVariantsInfoFieldsAreConsistent(final List<VariantContext> variants) {
for (final VariantContext variant : variants) {
final int expectedAN = variant.getGenotypes().size();
final int[] expectedAC = new int[CopyNumberTriState.values().length];
final List<Allele> alleles = variant.getAlleles();
for (final Genotype genotype : variant.getGenotypes()) {
Assert.assertEquals(genotype.getAlleles().size(), 1);
final int alleleIndex = alleles.indexOf(genotype.getAllele(0));
Assert.assertTrue(alleleIndex >= 0);
expectedAC[alleleIndex]++;
}
Assert.assertEquals(variant.getAlleles(), CopyNumberTriStateAllele.ALL_ALLELES);
Assert.assertTrue(variant.hasAttribute(XHMMSegmentGenotyper.NUMBER_OF_TARGETS_KEY));
Assert.assertTrue(variant.hasAttribute(VCFConstants.ALLELE_COUNT_KEY));
Assert.assertTrue(variant.hasAttribute(VCFConstants.END_KEY));
Assert.assertTrue(variant.hasAttribute(VCFConstants.ALLELE_FREQUENCY_KEY));
Assert.assertTrue(variant.hasAttribute(VCFConstants.ALLELE_NUMBER_KEY));
final int expectedANAnotherWay = IntStream.of(expectedAC).sum();
Assert.assertEquals(expectedANAnotherWay, expectedAN);
Assert.assertEquals(variant.getAttributeAsInt(VCFConstants.ALLELE_NUMBER_KEY, -1), expectedAN);
final double[] expectedAF = IntStream.of(expectedAC).mapToDouble(c -> c / (double) expectedAN).toArray();
final double[] observedAFWithoutRef = variant.getAttributeAsList(VCFConstants.ALLELE_FREQUENCY_KEY).stream().mapToDouble(o -> Double.parseDouble(String.valueOf(o))).toArray();
Assert.assertEquals(observedAFWithoutRef.length, expectedAF.length - 1);
for (int i = 0; i < observedAFWithoutRef.length; i++) {
Assert.assertEquals(observedAFWithoutRef[i], expectedAF[i + 1], 0.001);
}
final int[] observedACWithoutRef = variant.getAttributeAsList(VCFConstants.ALLELE_COUNT_KEY).stream().mapToInt(o -> Integer.parseInt(String.valueOf(o))).toArray();
Assert.assertEquals(observedACWithoutRef.length, expectedAC.length - 1);
for (int i = 0; i < observedACWithoutRef.length; i++) {
Assert.assertEquals(observedACWithoutRef[i], expectedAC[i + 1]);
}
Assert.assertEquals(variant.getAttributeAsInt(XHMMSegmentGenotyper.NUMBER_OF_TARGETS_KEY, -1), XHMMSegmentCallerBaseIntegrationTest.REALISTIC_TARGETS.targetCount(variant));
}
}
use of htsjdk.variant.variantcontext.Genotype in project gatk-protected by broadinstitute.
the class XHMMSegmentGenotyperIntegrationTest method assertVariantSegmentsAreCovered.
private void assertVariantSegmentsAreCovered(final List<VariantContext> variants, final List<HiddenStateSegmentRecord<CopyNumberTriState, Target>> variantSegments) {
for (final HiddenStateSegmentRecord<CopyNumberTriState, Target> variantSegment : variantSegments) {
final Optional<VariantContext> match = variants.stream().filter(vc -> new SimpleInterval(vc).equals(variantSegment.getSegment().getInterval())).findFirst();
Assert.assertTrue(match.isPresent());
final VariantContext matchedVariant = match.get();
final Genotype genotype = matchedVariant.getGenotype(variantSegment.getSampleName());
final String discovery = genotype.getAnyAttribute(XHMMSegmentGenotyper.DISCOVERY_KEY).toString();
Assert.assertTrue(discovery.equals(XHMMSegmentGenotyper.DISCOVERY_TRUE));
final CopyNumberTriState call = variantSegment.getSegment().getCall();
final List<Allele> gt = genotype.getAlleles();
Assert.assertEquals(gt.size(), 1);
// The call may not be the same for case where the event-quality is relatively low.
if (variantSegment.getSegment().getEventQuality() > 10) {
Assert.assertEquals(CopyNumberTriStateAllele.valueOf(gt.get(0)).state, call, genotype.toString());
}
final String[] SQ = genotype.getAnyAttribute(XHMMSegmentGenotyper.SOME_QUALITY_KEY).toString().split(VCFConstants.INFO_FIELD_ARRAY_SEPARATOR);
final double someQual = variantSegment.getSegment().getSomeQuality();
Assert.assertEquals(Double.parseDouble(SQ[call == CopyNumberTriState.DELETION ? 0 : 1]), someQual, XHMMSegmentGenotyper.PHRED_SCORE_PRECISION, variantSegment.getSampleName() + " => " + genotype.toString());
final String[] LQ = genotype.getAnyAttribute(XHMMSegmentGenotyper.START_QUALITY_KEY).toString().split(VCFConstants.INFO_FIELD_ARRAY_SEPARATOR);
final double startQuality = variantSegment.getSegment().getStartQuality();
Assert.assertEquals(Double.parseDouble(LQ[call == CopyNumberTriState.DELETION ? 0 : 1]), startQuality, XHMMSegmentGenotyper.PHRED_SCORE_PRECISION, variantSegment.getSampleName() + " => " + genotype.toString());
final String[] RQ = genotype.getAnyAttribute(XHMMSegmentGenotyper.END_QUALITY_KEY).toString().split(VCFConstants.INFO_FIELD_ARRAY_SEPARATOR);
final double endQuality = variantSegment.getSegment().getEndQuality();
Assert.assertEquals(Double.parseDouble(RQ[call == CopyNumberTriState.DELETION ? 0 : 1]), endQuality, XHMMSegmentGenotyper.PHRED_SCORE_PRECISION, variantSegment.getSampleName() + " => " + genotype.toString());
// Check the PL.
final int[] PL = genotype.getPL();
final int observedGQFromPL = Math.min(XHMMSegmentGenotyper.MAX_GQ, PL[CopyNumberTriStateAllele.REF.index()] - PL[CopyNumberTriStateAllele.valueOf(call).index()]);
final double expectedCallPL = GATKProtectedMathUtils.roundPhred(QualityUtils.phredScaleErrorRate(QualityUtils.qualToProb(variantSegment.getSegment().getExactQuality())), HMMPostProcessor.PHRED_SCORE_PRECISION);
final double expectedRefPL = GATKProtectedMathUtils.roundPhred(QualityUtils.phredScaleCorrectRate(QualityUtils.qualToProb(variantSegment.getSegment().getEventQuality())), HMMPostProcessor.PHRED_SCORE_PRECISION);
final int expectedGQFromPL = Math.min(XHMMSegmentGenotyper.MAX_GQ, (int) Math.round(expectedRefPL - expectedCallPL));
Assert.assertTrue(Math.abs(observedGQFromPL - expectedGQFromPL) <= 1, genotype.toString() + " " + variantSegment.getSegment().getEventQuality());
}
}
use of htsjdk.variant.variantcontext.Genotype in project gatk-protected by broadinstitute.
the class XHMMSegmentGenotyperIntegrationTest method assertVariantsAreCoveredBySegments.
private void assertVariantsAreCoveredBySegments(final List<VariantContext> variants, final List<HiddenStateSegmentRecord<CopyNumberTriState, Target>> variantSegments) {
for (final VariantContext variant : variants) {
final List<HiddenStateSegmentRecord<CopyNumberTriState, Target>> matches = variantSegments.stream().filter(s -> new SimpleInterval(variant).equals(s.getSegment().getInterval())).collect(Collectors.toList());
Assert.assertFalse(matches.isEmpty());
for (final Genotype genotype : variant.getGenotypes()) {
final boolean discovery = genotype.getExtendedAttribute(XHMMSegmentGenotyper.DISCOVERY_KEY).toString().equals(XHMMSegmentGenotyper.DISCOVERY_TRUE);
if (discovery) {
Assert.assertTrue(matches.stream().anyMatch(s -> s.getSampleName().equals(genotype.getSampleName())));
} else {
Assert.assertTrue(matches.stream().noneMatch(s -> s.getSampleName().equals(genotype.getSampleName())));
}
}
}
}
use of htsjdk.variant.variantcontext.Genotype in project gatk-protected by broadinstitute.
the class XHMMSegmentGenotyperIntegrationTest method assertVariantsPlGtAndGQAreConsistent.
private void assertVariantsPlGtAndGQAreConsistent(final List<VariantContext> variants) {
for (final VariantContext vc : variants) {
for (final Genotype gt : vc.getGenotypes()) {
final int[] PL = gt.getPL();
Assert.assertNotNull(PL);
final int[] twoLowestPLIndices = IntStream.range(0, PL.length).boxed().sorted((a, b) -> Integer.compare(PL[a], PL[b])).limit(2).mapToInt(n -> n).toArray();
final int minPLIndex = twoLowestPLIndices[0];
final int secondPLIndex = twoLowestPLIndices[1];
Assert.assertEquals(vc.getAlleles().indexOf(gt.getAlleles().get(0)), minPLIndex);
final int expectedGQ = Math.min(XHMMSegmentGenotyper.MAX_GQ, PL[secondPLIndex] - PL[minPLIndex]);
Assert.assertEquals(gt.getGQ(), expectedGQ);
}
}
}
Aggregations