Examples with HiddenStateSegmentRecord org.broadinstitute.hellbender.utils.hmm.segmentation.HiddenStateSegmentRecord used on opensource projects

Example 1 with HiddenStateSegmentRecord

use of org.broadinstitute.hellbender.utils.hmm.segmentation.HiddenStateSegmentRecord 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 HiddenStateSegmentRecord

use of org.broadinstitute.hellbender.utils.hmm.segmentation.HiddenStateSegmentRecord 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 HiddenStateSegmentRecord

use of org.broadinstitute.hellbender.utils.hmm.segmentation.HiddenStateSegmentRecord 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;
}

Example 4 with HiddenStateSegmentRecord

use of org.broadinstitute.hellbender.utils.hmm.segmentation.HiddenStateSegmentRecord in project gatk 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;
}

Example 5 with HiddenStateSegmentRecord

use of org.broadinstitute.hellbender.utils.hmm.segmentation.HiddenStateSegmentRecord 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())));
            }
        }
    }
}