Examples with Target org.broadinstitute.hellbender.tools.exome.Target used on opensource projects

Example 81 with Target

use of org.broadinstitute.hellbender.tools.exome.Target in project gatk by broadinstitute.

the class EvaluateCopyNumberTriStateCallsIntegrationTest method checkOutputTargetNumbers.

private void checkOutputTargetNumbers(final File targetsFile, final File vcfOutput) {
    final List<VariantContext> outputVariants = readVCFFile(vcfOutput);
    final TargetCollection<Target> targets = TargetArgumentCollection.readTargetCollection(targetsFile);
    for (final VariantContext context : outputVariants) {
        final List<Target> overlappingTargets = targets.targets(context);
        Assert.assertEquals(context.getAttributeAsInt(XHMMSegmentGenotyper.NUMBER_OF_TARGETS_KEY, -1), overlappingTargets.size());
    }
}

Example 82 with Target

use of org.broadinstitute.hellbender.tools.exome.Target in project gatk by broadinstitute.

the class EvaluateCopyNumberTriStateCallsIntegrationTest method checkOutputTruthConcordance.

private void checkOutputTruthConcordance(final File truthFile, final File targetsFile, final File vcfOutput, final EvaluationFiltersArgumentCollection filteringOptions) {
    final List<VariantContext> truthVariants = readVCFFile(truthFile);
    final List<VariantContext> outputVariants = readVCFFile(vcfOutput);
    final Set<String> outputSamples = outputVariants.get(0).getSampleNames();
    final TargetCollection<Target> targets = TargetArgumentCollection.readTargetCollection(targetsFile);
    for (final VariantContext truth : truthVariants) {
        final List<Target> overlappingTargets = targets.targets(truth);
        final List<VariantContext> overlappingOutput = outputVariants.stream().filter(vc -> new SimpleInterval(vc).overlaps(truth)).collect(Collectors.toList());
        if (overlappingTargets.isEmpty()) {
            Assert.assertTrue(overlappingOutput.isEmpty());
            continue;
        }
        Assert.assertFalse(overlappingOutput.isEmpty());
        Assert.assertEquals(overlappingOutput.stream().filter(vc -> new SimpleInterval(truth).equals(new SimpleInterval(vc))).count(), 1);
        @SuppressWarnings("all") final Optional<VariantContext> prospectiveMatchingOutput = overlappingOutput.stream().filter(vc -> new SimpleInterval(truth).equals(new SimpleInterval(vc))).findFirst();
        Assert.assertTrue(prospectiveMatchingOutput.isPresent());
        final VariantContext matchingOutput = prospectiveMatchingOutput.get();
        final int[] truthAC = calculateACFromTruth(truth);
        final long truthAN = MathUtils.sum(truthAC);
        final double[] truthAF = IntStream.of(Arrays.copyOfRange(truthAC, 1, truthAC.length)).mapToDouble(d -> d / (double) truthAN).toArray();
        Assert.assertEquals(matchingOutput.getAttributeAsInt(VariantEvaluationContext.TRUTH_ALLELE_NUMBER_KEY, -1), truthAN);
        assertEquals(GATKProtectedVariantContextUtils.getAttributeAsDoubleArray(matchingOutput, VariantEvaluationContext.TRUTH_ALLELE_FREQUENCY_KEY, () -> new double[2], 0.0), truthAF, 0.001);
        assertOutputVariantFilters(filteringOptions, overlappingTargets, matchingOutput, truthAF);
        for (final String sample : outputSamples) {
            final Genotype outputGenotype = matchingOutput.getGenotype(sample);
            final Genotype truthGenotype = truth.getGenotype(sample);
            final int truthCN = GATKProtectedVariantContextUtils.getAttributeAsInt(truthGenotype, ConvertGSVariantsToSegments.GS_COPY_NUMBER_FORMAT, -1);
            final int truthGT = GATKProtectedVariantContextUtils.getAttributeAsInt(outputGenotype, VariantEvaluationContext.TRUTH_GENOTYPE_KEY, -1);
            final Object truthQualObject = outputGenotype.getAnyAttribute(VariantEvaluationContext.TRUTH_QUALITY_KEY);
            Assert.assertNotNull(truthQualObject, "" + truthGenotype);
            final double truthQual = Double.parseDouble(String.valueOf(truthQualObject));
            if (truthQual < filteringOptions.minimumTruthSegmentQuality) {
                Assert.assertEquals(outputGenotype.getFilters(), EvaluationFilter.LowQuality.acronym);
            } else {
                Assert.assertTrue(outputGenotype.getFilters() == null || outputGenotype.getFilters().equals(VCFConstants.PASSES_FILTERS_v4));
            }
            final double[] truthPosteriors = GATKProtectedVariantContextUtils.getAttributeAsDoubleArray(truthGenotype, ConvertGSVariantsToSegments.GS_COPY_NUMBER_POSTERIOR, () -> new double[0], Double.NEGATIVE_INFINITY);
            final double truthDelPosterior = MathUtils.log10SumLog10(truthPosteriors, 0, EvaluateCopyNumberTriStateCalls.REFERENCE_COPY_NUMBER_DEFAULT);
            final double truthRefPosterior = truthPosteriors[EvaluateCopyNumberTriStateCalls.REFERENCE_COPY_NUMBER_DEFAULT];
            final double truthDupPosterior = truthPosteriors.length < EvaluateCopyNumberTriStateCalls.REFERENCE_COPY_NUMBER_DEFAULT ? Double.NEGATIVE_INFINITY : MathUtils.log10SumLog10(truthPosteriors, EvaluateCopyNumberTriStateCalls.REFERENCE_COPY_NUMBER_DEFAULT + 1, truthPosteriors.length);
            final CopyNumberTriStateAllele truthAllele = truthGT == -1 ? null : CopyNumberTriStateAllele.ALL_ALLELES.get(truthGT);
            if (truthCN < EvaluateCopyNumberTriStateCalls.REFERENCE_COPY_NUMBER_DEFAULT) {
                Assert.assertEquals(truthAllele, CopyNumberTriStateAllele.DEL);
                Assert.assertEquals(truthQual * -.1, MathUtils.log10SumLog10(new double[] { truthRefPosterior, truthDupPosterior }), 0.01);
            } else if (truthCN > EvaluateCopyNumberTriStateCalls.REFERENCE_COPY_NUMBER_DEFAULT) {
                Assert.assertEquals(truthAllele, CopyNumberTriStateAllele.DUP);
                Assert.assertEquals(truthQual * -.1, MathUtils.log10SumLog10(new double[] { truthDelPosterior, truthRefPosterior }), 0.01, "" + truthGenotype + " " + outputGenotype);
            } else {
                Assert.assertEquals(truthAllele, CopyNumberTriStateAllele.REF);
                Assert.assertEquals(truthQual * -.1, MathUtils.log10SumLog10(new double[] { truthDelPosterior, truthDupPosterior }), 0.01);
            }
            final double outputTruthFraction = GATKProtectedVariantContextUtils.getAttributeAsDouble(outputGenotype, VariantEvaluationContext.TRUTH_COPY_FRACTION_KEY, -1);
            final double inputTruthFraction = GATKProtectedVariantContextUtils.getAttributeAsDouble(truthGenotype, ConvertGSVariantsToSegments.GS_COPY_NUMBER_FRACTION, -1);
            Assert.assertEquals(outputTruthFraction, inputTruthFraction, 0.01);
        }
    }
}

Example 83 with Target

use of org.broadinstitute.hellbender.tools.exome.Target in project gatk 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());
    }
}

Example 84 with Target

use of org.broadinstitute.hellbender.tools.exome.Target in project gatk by broadinstitute.

the class XHMMModelUnitTest method testLogTransitionProbabilityOverlappingTargets.

@Test(dataProvider = "testData", dependsOnMethods = "testInstantiation")
public void testLogTransitionProbabilityOverlappingTargets(final double eventRate, final double meanEventSize, final double deletionMean, final double duplicationMean) {
    final XHMMModel model = new XHMMModel(eventRate, meanEventSize, deletionMean, duplicationMean);
    final Target fromTarget = new Target("target1", new SimpleInterval("1", 1, 100));
    final Target toTarget = new Target("target2", new SimpleInterval("1", 50, 150));
    assertTransitionProbabilities(eventRate, meanEventSize, model, fromTarget, toTarget);
}

Example 85 with Target

use of org.broadinstitute.hellbender.tools.exome.Target in project gatk by broadinstitute.

the class XHMMModelUnitTest method testLogTransitionProbabilityWithNoDistance.

@Test(dataProvider = "testData", dependsOnMethods = "testInstantiation", expectedExceptions = IllegalArgumentException.class)
public void testLogTransitionProbabilityWithNoDistance(final double eventStartProbability, final double meanEventSize, final double deletionMean, final double duplicationMean) {
    final XHMMModel model = new XHMMModel(eventStartProbability, meanEventSize, deletionMean, duplicationMean);
    final Target fromTarget = new Target("target1");
    final Target toTarget = new Target("target2");
    assertTransitionProbabilities(eventStartProbability, meanEventSize, model, fromTarget, toTarget);
}