Search in sources :

Example 86 with Genotype

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);
    }
}
Also used : Genotype(htsjdk.variant.variantcontext.Genotype)

Example 87 with Genotype

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));
    }
}
Also used : Genotype(htsjdk.variant.variantcontext.Genotype) HiddenStateSegmentRecord(org.broadinstitute.hellbender.utils.hmm.segmentation.HiddenStateSegmentRecord) IntStream(java.util.stream.IntStream) Allele(htsjdk.variant.variantcontext.Allele) htsjdk.variant.vcf(htsjdk.variant.vcf) java.util(java.util) GATKProtectedMathUtils(org.broadinstitute.hellbender.utils.GATKProtectedMathUtils) DataProvider(org.testng.annotations.DataProvider) StandardArgumentDefinitions(org.broadinstitute.hellbender.cmdline.StandardArgumentDefinitions) QualityUtils(org.broadinstitute.hellbender.utils.QualityUtils) Test(org.testng.annotations.Test) IOException(java.io.IOException) TargetArgumentCollection(org.broadinstitute.hellbender.tools.exome.TargetArgumentCollection) CopyNumberTriStateAllele(org.broadinstitute.hellbender.tools.exome.germlinehmm.CopyNumberTriStateAllele) SimpleInterval(org.broadinstitute.hellbender.utils.SimpleInterval) Collectors(java.util.stream.Collectors) File(java.io.File) HMMPostProcessor(org.broadinstitute.hellbender.utils.hmm.segmentation.HMMPostProcessor) Assert(org.testng.Assert) Target(org.broadinstitute.hellbender.tools.exome.Target) CopyNumberTriState(org.broadinstitute.hellbender.tools.exome.germlinehmm.CopyNumberTriState) VariantContext(htsjdk.variant.variantcontext.VariantContext) HiddenStateSegmentRecordReader(org.broadinstitute.hellbender.utils.hmm.segmentation.HiddenStateSegmentRecordReader) Allele(htsjdk.variant.variantcontext.Allele) CopyNumberTriStateAllele(org.broadinstitute.hellbender.tools.exome.germlinehmm.CopyNumberTriStateAllele) VariantContext(htsjdk.variant.variantcontext.VariantContext) Genotype(htsjdk.variant.variantcontext.Genotype)

Example 88 with Genotype

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());
    }
}
Also used : Genotype(htsjdk.variant.variantcontext.Genotype) HiddenStateSegmentRecord(org.broadinstitute.hellbender.utils.hmm.segmentation.HiddenStateSegmentRecord) IntStream(java.util.stream.IntStream) Allele(htsjdk.variant.variantcontext.Allele) htsjdk.variant.vcf(htsjdk.variant.vcf) java.util(java.util) GATKProtectedMathUtils(org.broadinstitute.hellbender.utils.GATKProtectedMathUtils) DataProvider(org.testng.annotations.DataProvider) StandardArgumentDefinitions(org.broadinstitute.hellbender.cmdline.StandardArgumentDefinitions) QualityUtils(org.broadinstitute.hellbender.utils.QualityUtils) Test(org.testng.annotations.Test) IOException(java.io.IOException) TargetArgumentCollection(org.broadinstitute.hellbender.tools.exome.TargetArgumentCollection) CopyNumberTriStateAllele(org.broadinstitute.hellbender.tools.exome.germlinehmm.CopyNumberTriStateAllele) SimpleInterval(org.broadinstitute.hellbender.utils.SimpleInterval) Collectors(java.util.stream.Collectors) File(java.io.File) HMMPostProcessor(org.broadinstitute.hellbender.utils.hmm.segmentation.HMMPostProcessor) Assert(org.testng.Assert) Target(org.broadinstitute.hellbender.tools.exome.Target) CopyNumberTriState(org.broadinstitute.hellbender.tools.exome.germlinehmm.CopyNumberTriState) VariantContext(htsjdk.variant.variantcontext.VariantContext) HiddenStateSegmentRecordReader(org.broadinstitute.hellbender.utils.hmm.segmentation.HiddenStateSegmentRecordReader) VariantContext(htsjdk.variant.variantcontext.VariantContext) Genotype(htsjdk.variant.variantcontext.Genotype) Target(org.broadinstitute.hellbender.tools.exome.Target) Allele(htsjdk.variant.variantcontext.Allele) CopyNumberTriStateAllele(org.broadinstitute.hellbender.tools.exome.germlinehmm.CopyNumberTriStateAllele) CopyNumberTriState(org.broadinstitute.hellbender.tools.exome.germlinehmm.CopyNumberTriState) SimpleInterval(org.broadinstitute.hellbender.utils.SimpleInterval)

Example 89 with Genotype

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())));
            }
        }
    }
}
Also used : Genotype(htsjdk.variant.variantcontext.Genotype) HiddenStateSegmentRecord(org.broadinstitute.hellbender.utils.hmm.segmentation.HiddenStateSegmentRecord) IntStream(java.util.stream.IntStream) Allele(htsjdk.variant.variantcontext.Allele) htsjdk.variant.vcf(htsjdk.variant.vcf) java.util(java.util) GATKProtectedMathUtils(org.broadinstitute.hellbender.utils.GATKProtectedMathUtils) DataProvider(org.testng.annotations.DataProvider) StandardArgumentDefinitions(org.broadinstitute.hellbender.cmdline.StandardArgumentDefinitions) QualityUtils(org.broadinstitute.hellbender.utils.QualityUtils) Test(org.testng.annotations.Test) IOException(java.io.IOException) TargetArgumentCollection(org.broadinstitute.hellbender.tools.exome.TargetArgumentCollection) CopyNumberTriStateAllele(org.broadinstitute.hellbender.tools.exome.germlinehmm.CopyNumberTriStateAllele) SimpleInterval(org.broadinstitute.hellbender.utils.SimpleInterval) Collectors(java.util.stream.Collectors) File(java.io.File) HMMPostProcessor(org.broadinstitute.hellbender.utils.hmm.segmentation.HMMPostProcessor) Assert(org.testng.Assert) Target(org.broadinstitute.hellbender.tools.exome.Target) CopyNumberTriState(org.broadinstitute.hellbender.tools.exome.germlinehmm.CopyNumberTriState) VariantContext(htsjdk.variant.variantcontext.VariantContext) HiddenStateSegmentRecordReader(org.broadinstitute.hellbender.utils.hmm.segmentation.HiddenStateSegmentRecordReader) VariantContext(htsjdk.variant.variantcontext.VariantContext) HiddenStateSegmentRecord(org.broadinstitute.hellbender.utils.hmm.segmentation.HiddenStateSegmentRecord) Genotype(htsjdk.variant.variantcontext.Genotype) SimpleInterval(org.broadinstitute.hellbender.utils.SimpleInterval)

Example 90 with Genotype

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);
        }
    }
}
Also used : Genotype(htsjdk.variant.variantcontext.Genotype) HiddenStateSegmentRecord(org.broadinstitute.hellbender.utils.hmm.segmentation.HiddenStateSegmentRecord) IntStream(java.util.stream.IntStream) Allele(htsjdk.variant.variantcontext.Allele) htsjdk.variant.vcf(htsjdk.variant.vcf) java.util(java.util) GATKProtectedMathUtils(org.broadinstitute.hellbender.utils.GATKProtectedMathUtils) DataProvider(org.testng.annotations.DataProvider) StandardArgumentDefinitions(org.broadinstitute.hellbender.cmdline.StandardArgumentDefinitions) QualityUtils(org.broadinstitute.hellbender.utils.QualityUtils) Test(org.testng.annotations.Test) IOException(java.io.IOException) TargetArgumentCollection(org.broadinstitute.hellbender.tools.exome.TargetArgumentCollection) CopyNumberTriStateAllele(org.broadinstitute.hellbender.tools.exome.germlinehmm.CopyNumberTriStateAllele) SimpleInterval(org.broadinstitute.hellbender.utils.SimpleInterval) Collectors(java.util.stream.Collectors) File(java.io.File) HMMPostProcessor(org.broadinstitute.hellbender.utils.hmm.segmentation.HMMPostProcessor) Assert(org.testng.Assert) Target(org.broadinstitute.hellbender.tools.exome.Target) CopyNumberTriState(org.broadinstitute.hellbender.tools.exome.germlinehmm.CopyNumberTriState) VariantContext(htsjdk.variant.variantcontext.VariantContext) HiddenStateSegmentRecordReader(org.broadinstitute.hellbender.utils.hmm.segmentation.HiddenStateSegmentRecordReader) VariantContext(htsjdk.variant.variantcontext.VariantContext) Genotype(htsjdk.variant.variantcontext.Genotype)

Aggregations

Genotype (htsjdk.variant.variantcontext.Genotype)150 VariantContext (htsjdk.variant.variantcontext.VariantContext)97 Allele (htsjdk.variant.variantcontext.Allele)82 ArrayList (java.util.ArrayList)54 VariantContextBuilder (htsjdk.variant.variantcontext.VariantContextBuilder)52 GenotypeBuilder (htsjdk.variant.variantcontext.GenotypeBuilder)51 File (java.io.File)48 VCFHeader (htsjdk.variant.vcf.VCFHeader)46 IOException (java.io.IOException)45 Collectors (java.util.stream.Collectors)42 Test (org.testng.annotations.Test)37 HashSet (java.util.HashSet)35 VariantContextWriter (htsjdk.variant.variantcontext.writer.VariantContextWriter)29 List (java.util.List)29 VCFHeaderLine (htsjdk.variant.vcf.VCFHeaderLine)27 VCFFormatHeaderLine (htsjdk.variant.vcf.VCFFormatHeaderLine)25 SAMSequenceDictionaryProgress (com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress)23 HashMap (java.util.HashMap)23 java.util (java.util)22 Set (java.util.Set)22