Search in sources :

Example 1 with Transition

use of org.broadinstitute.hellbender.tools.picard.analysis.artifacts.Transition in project gatk by broadinstitute.

the class FilterByOrientationBias method onTraversalSuccess.

@Override
public Object onTraversalSuccess() {
    logger.info("Tagging whether genotypes are in one of the artifact modes.");
    // Calculate how many artifacts need to be cut
    double fdrThreshold = 0.01;
    final List<VariantContext> finalVariants = OrientationBiasFilterer.annotateVariantContextsWithFilterResults(fdrThreshold, relevantTransitions, firstPassVariants, transitionToPreAdapterScoreMap);
    logger.info("Writing variants to VCF...");
    finalVariants.forEach(vcfWriter::add);
    logger.info("Writing a simple summary table...");
    List<String> sampleNames = new ArrayList<>();
    if (finalVariants.size() != 0) {
        sampleNames = finalVariants.get(0).getSampleNamesOrderedByName();
    }
    final List<Pair<String, Transition>> sampleTransitionCombinations = new ArrayList<>();
    for (Transition relevantTransition : relevantTransitions) {
        for (String sampleName : sampleNames) {
            sampleTransitionCombinations.add(Pair.of(sampleName, relevantTransition));
        }
    }
    OrientationBiasUtils.writeOrientationBiasSummaryTable(sampleTransitionCombinations, finalVariants, transitionToPreAdapterScoreMap, new File(outputFile.getAbsolutePath() + SUMMARY_FILE_SUFFIX));
    return null;
}
Also used : Transition(org.broadinstitute.hellbender.tools.picard.analysis.artifacts.Transition) VariantContext(htsjdk.variant.variantcontext.VariantContext) MetricsFile(htsjdk.samtools.metrics.MetricsFile) File(java.io.File) Pair(org.apache.commons.lang3.tuple.Pair)

Example 2 with Transition

use of org.broadinstitute.hellbender.tools.picard.analysis.artifacts.Transition in project gatk-protected by broadinstitute.

the class OrientationBiasFilterer method createTransitionToNumCutPrePreAdapterQ.

private static Map<Transition, Long> createTransitionToNumCutPrePreAdapterQ(double fdrThresh, String sampleName, long unfilteredGenotypeCount, final SortedMap<Genotype, VariantContext> genotypesToConsiderForFiltering, final Map<Transition, Long> transitionCount) {
    final long allTransitionCount = transitionCount.values().stream().mapToLong(Long::longValue).sum();
    final int totalNumToCut = calculateTotalNumToCut(fdrThresh, unfilteredGenotypeCount, genotypesToConsiderForFiltering);
    logger.info(sampleName + ": Cutting (total) pre-preAdapterQ: " + String.valueOf(totalNumToCut));
    // Adjust the number to cut based on artifact mode
    final Map<Transition, Long> transitionNumToCut = new HashMap<>();
    transitionCount.keySet().stream().forEach(transition -> transitionNumToCut.put(transition, 0L));
    for (final Transition transition : transitionNumToCut.keySet()) {
        transitionNumToCut.put(transition, Long.valueOf(Math.round(totalNumToCut * transitionCount.get(transition) / allTransitionCount)));
        logger.info(sampleName + ": Cutting (" + transition + ") pre-preAdapterQ: " + transitionNumToCut.get(transition));
    }
    return transitionNumToCut;
}
Also used : Transition(org.broadinstitute.hellbender.tools.picard.analysis.artifacts.Transition)

Example 3 with Transition

use of org.broadinstitute.hellbender.tools.picard.analysis.artifacts.Transition in project gatk-protected by broadinstitute.

the class OrientationBiasFiltererUnitTest method testAnnotateVariantContextWithPreprocessingValuesMultiArtifact.

@Test
public void testAnnotateVariantContextWithPreprocessingValuesMultiArtifact() {
    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);
    for (final VariantContext vc : featureDataSource) {
        final VariantContext updatedVariantContext = OrientationBiasFilterer.annotateVariantContextWithPreprocessingValues(vc, relevantTransitions, preAdapterQFakeScoreMap);
        final Genotype genotypeTumor = updatedVariantContext.getGenotype("TUMOR");
        final Genotype genotypeNormal = updatedVariantContext.getGenotype("NORMAL");
        // This is mostly just to make sure that nobody breaks the test itself.  I.e. that this test will test all tumor genotype paths be artifact or non-artifact.
        boolean wasGenotypeTumorTested = false;
        // Check whether this genotype is reverse complement or actual artifact mode
        wasGenotypeTumorTested |= assertArtifact(amGTPreAdapterQ, genotypeTumor, relevantTransitions.first());
        wasGenotypeTumorTested |= assertArtifact(amCTPreAdapterQ, genotypeTumor, relevantTransitions.last());
        // Check any variants that are not an artifact mode but are SNP
        if (!OrientationBiasUtils.isGenotypeInTransitionsWithComplement(genotypeTumor, relevantTransitions)) {
            assertNotTransition(genotypeTumor);
            wasGenotypeTumorTested = true;
        } else {
            // Check attributes common to all variants in artifact mode
            Assert.assertNotEquals(genotypeTumor.getExtendedAttribute(OrientationBiasFilterConstants.FOB, VCFConstants.EMPTY_ALLELE), VCFConstants.EMPTY_ALLELE);
            Assert.assertNotEquals(genotypeTumor.getExtendedAttribute(OrientationBiasFilterConstants.P_ARTIFACT_FIELD_NAME, VCFConstants.EMPTY_ALLELE), VCFConstants.EMPTY_ALLELE);
        }
        // The NORMAL is always ref/ref in the example file.
        assertNormal(genotypeNormal);
        Assert.assertTrue(wasGenotypeTumorTested, "The test seems to be broken...  A variant context was tested, but it had no tumor genotype.");
    }
}
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)

Example 4 with Transition

use of org.broadinstitute.hellbender.tools.picard.analysis.artifacts.Transition in project gatk-protected by broadinstitute.

the class OrientationBiasFiltererUnitTest method testAnnotateVariantContextWithPreprocessingValues.

@Test
public void testAnnotateVariantContextWithPreprocessingValues() {
    final FeatureDataSource<VariantContext> featureDataSource = new FeatureDataSource<>(new File(smallM2Vcf));
    SortedSet<Transition> relevantTransitions = new TreeSet<>();
    relevantTransitions.add(Transition.transitionOf('G', 'T'));
    final Map<Transition, Double> preAdapterQFakeScoreMap = new HashMap<>();
    final double amGTPreAdapterQ = 20.0;
    // preAdapterQ suppression will do nothing.
    preAdapterQFakeScoreMap.put(Transition.transitionOf('G', 'T'), amGTPreAdapterQ);
    for (final VariantContext vc : featureDataSource) {
        final VariantContext updatedVariantContext = OrientationBiasFilterer.annotateVariantContextWithPreprocessingValues(vc, relevantTransitions, preAdapterQFakeScoreMap);
        final Genotype genotypeTumor = updatedVariantContext.getGenotype("TUMOR");
        final Genotype genotypeNormal = updatedVariantContext.getGenotype("NORMAL");
        Assert.assertNotEquals(genotypeTumor.getExtendedAttribute(OrientationBiasFilterConstants.FOB, VCFConstants.EMPTY_ALLELE), VCFConstants.EMPTY_ALLELE);
        Assert.assertNotEquals(genotypeTumor.getExtendedAttribute(OrientationBiasFilterConstants.P_ARTIFACT_FIELD_NAME, VCFConstants.EMPTY_ALLELE), VCFConstants.EMPTY_ALLELE);
        assertArtifact(amGTPreAdapterQ, genotypeTumor, Transition.transitionOf('G', 'T'));
        // The NORMAL is always ref/ref in the example file.
        assertNormal(genotypeNormal);
    }
}
Also used : Transition(org.broadinstitute.hellbender.tools.picard.analysis.artifacts.Transition) VariantContext(htsjdk.variant.variantcontext.VariantContext) Genotype(htsjdk.variant.variantcontext.Genotype) File(java.io.File) FeatureDataSource(org.broadinstitute.hellbender.engine.FeatureDataSource) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Example 5 with Transition

use of org.broadinstitute.hellbender.tools.picard.analysis.artifacts.Transition 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

Transition (org.broadinstitute.hellbender.tools.picard.analysis.artifacts.Transition)32 BaseTest (org.broadinstitute.hellbender.utils.test.BaseTest)16 Test (org.testng.annotations.Test)16 VariantContext (htsjdk.variant.variantcontext.VariantContext)14 File (java.io.File)14 FeatureDataSource (org.broadinstitute.hellbender.engine.FeatureDataSource)12 VisibleForTesting (com.google.cloud.dataflow.sdk.repackaged.com.google.common.annotations.VisibleForTesting)6 MetricsFile (htsjdk.samtools.metrics.MetricsFile)6 Genotype (htsjdk.variant.variantcontext.Genotype)6 RealMatrix (org.apache.commons.math3.linear.RealMatrix)6 htsjdk.variant.variantcontext (htsjdk.variant.variantcontext)4 htsjdk.variant.vcf (htsjdk.variant.vcf)4 FileReader (java.io.FileReader)4 java.util (java.util)4 HashMap (java.util.HashMap)4 Collectors (java.util.stream.Collectors)4 IntStream (java.util.stream.IntStream)4 Array2DRowRealMatrix (org.apache.commons.math3.linear.Array2DRowRealMatrix)4 LogManager (org.apache.logging.log4j.LogManager)4 Logger (org.apache.logging.log4j.Logger)4