Search in sources :

Example 36 with Genotype

use of htsjdk.variant.variantcontext.Genotype in project gatk by broadinstitute.

the class HomRefBlock method toVariantContext.

/**
     * Convert a HomRefBlock into a VariantContext
     *
     * @param sampleName sample name to give this variant context
     * @return a VariantContext representing the gVCF encoding for this block.
     * It will return {@code null} if input {@code block} is {@code null}, indicating that there
     * is no variant-context to be output into the VCF.
     */
public VariantContext toVariantContext(String sampleName) {
    final VariantContextBuilder vcb = new VariantContextBuilder(getStartingVC());
    // clear the attributes
    vcb.attributes(new LinkedHashMap<>(2));
    vcb.stop(getEnd());
    vcb.attribute(VCFConstants.END_KEY, getEnd());
    final Genotype genotype = createHomRefGenotype(sampleName);
    return vcb.genotypes(genotype).make();
}
Also used : VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) Genotype(htsjdk.variant.variantcontext.Genotype)

Example 37 with Genotype

use of htsjdk.variant.variantcontext.Genotype in project gatk 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 38 with Genotype

use of htsjdk.variant.variantcontext.Genotype in project gatk 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)

Example 39 with Genotype

use of htsjdk.variant.variantcontext.Genotype in project gatk 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 40 with Genotype

use of htsjdk.variant.variantcontext.Genotype in project gatk by broadinstitute.

the class OrientationBiasFiltererUnitTest method testCreateSampleToGenotypeVCMap.

@Test
public void testCreateSampleToGenotypeVCMap() {
    // Setup the test
    final FeatureDataSource<VariantContext> featureDataSource = new FeatureDataSource<>(new File(smallM2VcfMore));
    SortedSet<Transition> relevantTransitions = new TreeSet<>();
    relevantTransitions.add(Transition.transitionOf('G', 'T'));
    relevantTransitions.add(Transition.transitionOf('C', 'T'));
    final Map<Transition, Double> preAdapterQFakeScoreMap = new HashMap<>();
    final double amGTPreAdapterQ = 20.0;
    final double amCTPreAdapterQ = 25.0;
    // preAdapterQ suppression will do nothing.
    preAdapterQFakeScoreMap.put(relevantTransitions.first(), amGTPreAdapterQ);
    // preAdapterQ suppression will do nothing.
    preAdapterQFakeScoreMap.put(relevantTransitions.last(), amCTPreAdapterQ);
    final List<VariantContext> updatedVariants = new ArrayList<>();
    for (final VariantContext vc : featureDataSource) {
        final VariantContext updatedVariantContext = OrientationBiasFilterer.annotateVariantContextWithPreprocessingValues(vc, relevantTransitions, preAdapterQFakeScoreMap);
        updatedVariants.add(updatedVariantContext);
    }
    final List<String> sampleNames = updatedVariants.get(0).getSampleNamesOrderedByName();
    // Do the test
    // Create a mapping from sample name to a genotype->variant context map with the second map sorted by p_artifact (i.e. OrientationBiasFilterConstants.P_ARTIFACT_FIELD_NAME)
    final Map<String, SortedMap<Genotype, VariantContext>> sampleNameToVariants = OrientationBiasFilterer.createSampleToGenotypeVariantContextSortedMap(sampleNames, updatedVariants);
    Assert.assertEquals(sampleNameToVariants.keySet().size(), 2);
    Assert.assertTrue(sampleNameToVariants.keySet().contains("TUMOR"));
    Assert.assertTrue(sampleNameToVariants.keySet().contains("NORMAL"));
    Assert.assertEquals(sampleNameToVariants.get("TUMOR").keySet().size(), 8);
    // None of the normal genotypes should have a pvalue, so cannot/shouldn't be added to the sorted map
    Assert.assertEquals(sampleNameToVariants.get("NORMAL").keySet().size(), 0);
    // Check that the sorted map is getting smaller (or same) values of p_artifact and not staying put.
    double previousPArtifact = Double.POSITIVE_INFINITY;
    for (final Genotype genotypeTumor : sampleNameToVariants.get("TUMOR").keySet()) {
        final Double pArtifact = OrientationBiasUtils.getGenotypeDouble(genotypeTumor, OrientationBiasFilterConstants.P_ARTIFACT_FIELD_NAME, Double.POSITIVE_INFINITY);
        Assert.assertNotNull(pArtifact);
        Assert.assertTrue(pArtifact <= previousPArtifact);
        Assert.assertNotEquals(pArtifact, Double.POSITIVE_INFINITY);
    }
}
Also used : VariantContext(htsjdk.variant.variantcontext.VariantContext) Genotype(htsjdk.variant.variantcontext.Genotype) Transition(org.broadinstitute.hellbender.tools.picard.analysis.artifacts.Transition) File(java.io.File) FeatureDataSource(org.broadinstitute.hellbender.engine.FeatureDataSource) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

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