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

Example 1 with CopyNumberTriState

use of org.broadinstitute.hellbender.tools.exome.germlinehmm.CopyNumberTriState in project gatk-protected by broadinstitute.

the class ConvertGSVariantsToSegments method apply.

@Override
public void apply(final VariantContext variant, final ReadsContext readsContext, final ReferenceContext referenceContext, final FeatureContext featureContext) {
    final SimpleInterval interval = new SimpleInterval(variant);
    final int targetCount = targets.indexRange(interval).size();
    final int[] callCounts = new int[CopyNumberTriState.values().length];
    for (final Genotype genotype : variant.getGenotypes().iterateInSampleNameOrder()) {
        final String sample = genotype.getSampleName();
        final double mean = doubleFrom(genotype.getExtendedAttribute(GS_COPY_NUMBER_FRACTION));
        final int copyNumber = intFrom(genotype.getExtendedAttribute(GS_COPY_NUMBER_FORMAT));
        final CopyNumberTriState call = copyNumber == neutralCopyNumber ? CopyNumberTriState.NEUTRAL : (copyNumber < neutralCopyNumber) ? CopyNumberTriState.DELETION : CopyNumberTriState.DUPLICATION;
        callCounts[call.ordinal()]++;
        final double[] probs = doubleArrayFrom(genotype.getExtendedAttribute(GS_COPY_NUMBER_POSTERIOR));
        final double log10PostQualCall = calculateLog10CallQuality(probs, call);
        final double log10PostQualNonRef = calculateLog10CallQualityNonRef(probs);
        final double phredProbCall = -10.0 * log10PostQualCall;
        final double phredProbNonRef = -10.0 * log10PostQualNonRef;
        final HiddenStateSegment<CopyNumberTriState, Target> segment = new HiddenStateSegment<>(interval, targetCount, mean, // GS VCF does not contain any stddev or var estimate for coverage fraction.
        0.0, call, // GS does not provide an EQ, we approximate it to be the 1 - sum of all call compatible CN corresponding posterior probs
        phredProbCall, // GS does not provide a SQ, we leave is a NaN.
        Double.NaN, // GS does not provide a START Q.
        Double.NaN, // GS does not provide a END Q.
        Double.NaN, phredProbNonRef);
        final HiddenStateSegmentRecord<CopyNumberTriState, Target> record = new HiddenStateSegmentRecord<>(sample, segment);
        try {
            outputWriter.writeRecord(record);
        } catch (final IOException ex) {
            throw new UserException.CouldNotCreateOutputFile(outputFile, ex);
        }
    }
}

Example 2 with CopyNumberTriState

use of org.broadinstitute.hellbender.tools.exome.germlinehmm.CopyNumberTriState in project gatk-protected by broadinstitute.

the class ConvertGSVariantsToSegmentsIntegrationTest method composeExpectedSegments.

private List<HiddenStateSegmentRecord<CopyNumberTriState, Target>> composeExpectedSegments(final File vcf, final TargetCollection<Target> targets) throws IOException {
    final VCFFileReader reader = new VCFFileReader(vcf, false);
    final List<HiddenStateSegmentRecord<CopyNumberTriState, Target>> result = new ArrayList<>();
    reader.iterator().forEachRemaining(vc -> {
        final int targetCount = targets.indexRange(vc).size();
        for (final Genotype genotype : vc.getGenotypes()) {
            final int cn = Integer.parseInt(genotype.getExtendedAttribute("CN").toString());
            final double[] cnp = Stream.of(genotype.getExtendedAttribute("CNP").toString().replaceAll("\\[\\]", "").split(",")).mapToDouble(Double::parseDouble).toArray();
            final double cnpSum = MathUtils.approximateLog10SumLog10(cnp);
            final CopyNumberTriState call = expectedCall(cn);
            final double exactLog10Prob = expectedExactLog10(call, cnp);
            final HiddenStateSegment<CopyNumberTriState, Target> expectedSegment = new HiddenStateSegment<>(new SimpleInterval(vc), targetCount, Double.parseDouble(genotype.getExtendedAttribute("CNF").toString()), 0.000, call, -10.0 * exactLog10Prob, Double.NaN, Double.NaN, Double.NaN, -10.0 * (cnp[ConvertGSVariantsToSegments.NEUTRAL_COPY_NUMBER_DEFAULT] - cnpSum));
            result.add(new HiddenStateSegmentRecord<>(genotype.getSampleName(), expectedSegment));
        }
    });
    return result;
}

Example 3 with CopyNumberTriState

use of org.broadinstitute.hellbender.tools.exome.germlinehmm.CopyNumberTriState in project gatk-protected by broadinstitute.

the class XHMMSegmentCallerIntegrationTest method assertOutputIsInOrder.

private void assertOutputIsInOrder(final List<HiddenStateSegmentRecord<CopyNumberTriState, Target>> outputRecords, final TargetCollection<Target> targets) {
    for (int i = 1; i < outputRecords.size(); i++) {
        final HiddenStateSegmentRecord<CopyNumberTriState, Target> nextRecord = outputRecords.get(i);
        final HiddenStateSegmentRecord<CopyNumberTriState, Target> previousRecord = outputRecords.get(i - 1);
        final IndexRange nextRange = targets.indexRange(nextRecord.getSegment());
        final IndexRange previousRange = targets.indexRange(previousRecord.getSegment());
        Assert.assertTrue(nextRange.from >= previousRange.from);
    }
}

Example 4 with CopyNumberTriState

use of org.broadinstitute.hellbender.tools.exome.germlinehmm.CopyNumberTriState in project gatk-protected by broadinstitute.

the class XHMMModelUnitTest method assertTransitionProbabilities.

// CopyNumberTriStateTransitionProbabilityCache is thoroughly tested.
// Here we simply test that it is invoked  properly
private void assertTransitionProbabilities(double eventStartProbability, double meanEventSize, XHMMModel model, Target fromTarget, Target toTarget) {
    final CopyNumberTriStateTransitionProbabilityCache logTransitionProbabilityCache = new CopyNumberTriStateTransitionProbabilityCache(meanEventSize, eventStartProbability);
    final double distance = XHMMModel.calculateDistance(fromTarget, toTarget);
    for (final CopyNumberTriState from : CopyNumberTriState.values()) {
        for (final CopyNumberTriState to : CopyNumberTriState.values()) {
            Assert.assertEquals(model.logTransitionProbability(from, fromTarget, to, toTarget), logTransitionProbabilityCache.logProbability((int) distance, to, from), EPSILON);
        }
    }
}

Example 5 with CopyNumberTriState

use of org.broadinstitute.hellbender.tools.exome.germlinehmm.CopyNumberTriState in project gatk-protected by broadinstitute.

the class XHMMSegmentCallerIntegrationTest method testRunCommandLine.

//TODO: this test used to contain a test of concordance with XHMM.  It no longer does that because our model has
//TODO: diverged from XHMM's.  Eventually the right thing to do is use the simulateChain() method to generate
//TODO: simulated data for some artificial set of CNV segments and to test concordance with those segments.
//TODO: however we still use XHMM's emission model, which is both not generative and quite silly.  Once we
//TODO: have a generative model of coverage we can modify simulateChain() accordingly and then write a concordance
//TODO: test here.  Until then, we do not have an integration test but we do have our ongoing evaluations, which
//TODO: show the superiority of our modifications versus the original XHMM model.
@Test(dataProvider = "simulatedChainData")
public File testRunCommandLine(final XHMMData chain) {
    final File inputFile = writeChainInTempFile(chain);
    final File outputFile = createTempFile("output", ".tab");
    runCommandLine(chain, inputFile, outputFile);
    Assert.assertTrue(outputFile.exists());
    final TargetCollection<Target> targets = TargetArgumentCollection.readTargetCollection(REALISTIC_TARGETS_FILE);
    final List<HiddenStateSegmentRecord<CopyNumberTriState, Target>> outputRecords = readOutputRecords(outputFile);
    assertOutputIsInOrder(outputRecords, targets);
    assertOutputHasConsistentNumberOfTargets(outputRecords, targets);
    final Map<String, List<HiddenStateSegmentRecord<CopyNumberTriState, Target>>> outputRecordsBySample = splitOutputRecordBySample(outputRecords);
    assertSampleNames(outputRecordsBySample.keySet(), chain);
    for (final List<HiddenStateSegmentRecord<CopyNumberTriState, Target>> sampleRecords : outputRecordsBySample.values()) {
        assertSampleSegmentsCoverAllTargets(sampleRecords, targets);
        assertSampleSegmentsCoordinates(sampleRecords, targets);
    }
    return outputFile;
}