Search in sources :

Example 1 with AllelicCountWithPhasePosteriorsCollection

use of org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountWithPhasePosteriorsCollection in project gatk-protected by broadinstitute.

the class CalculatePulldownPhasePosteriors method doWork.

@Override
public Object doWork() {
    if (!new HDF5Library().load(null)) {
        //Note: passing null means using the default temp dir.
        throw new UserException.HardwareFeatureException("Cannot load the required HDF5 library. " + "HDF5 is currently supported on x86-64 architecture and Linux or OSX systems.");
    }
    //read counts, segments, and parameters from files
    final AllelicCountCollection counts = new AllelicCountCollection(snpCountsFile);
    final List<ACNVModeledSegment> segments = SegmentUtils.readACNVModeledSegmentFile(segmentsFile);
    final AlleleFractionState state = reconstructState(segments, parametersFile);
    //load allelic-bias panel of normals if provided
    final AllelicPanelOfNormals allelicPoN = allelicPoNFile != null ? AllelicPanelOfNormals.read(allelicPoNFile) : AllelicPanelOfNormals.EMPTY_PON;
    //calculate phase posteriors
    final List<SimpleInterval> unmodeledSegments = segments.stream().map(ACNVModeledSegment::getInterval).collect(Collectors.toList());
    final AllelicCountWithPhasePosteriorsCollection countsWithPhasePosteriors = calculatePhasePosteriors(counts, unmodeledSegments, state, allelicPoN);
    //write phase posteriors to file with same verbosity as input file
    countsWithPhasePosteriors.write(outputFile, counts.getVerbosity());
    return "SUCCESS";
}
Also used : AllelicPanelOfNormals(org.broadinstitute.hellbender.tools.pon.allelic.AllelicPanelOfNormals) AllelicCountWithPhasePosteriorsCollection(org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountWithPhasePosteriorsCollection) HDF5Library(org.broadinstitute.hdf5.HDF5Library) AllelicCountCollection(org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountCollection) SimpleInterval(org.broadinstitute.hellbender.utils.SimpleInterval)

Example 2 with AllelicCountWithPhasePosteriorsCollection

use of org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountWithPhasePosteriorsCollection in project gatk-protected by broadinstitute.

the class CalculatePulldownPhasePosteriorsIntegrationTest method testCalculatePhasePosteriors.

/**
     * Uses {@link AlleleFractionSimulatedData} to test recovery of phase indicators.  Phase with highest posterior
     * probability is compared to the true phase; we require that
     * {@link CalculatePulldownPhasePosteriorsIntegrationTest#FRACTION_OF_INDICATORS_CORRECT_THRESHOLD} of the
     * indicators are recovered correctly.
     */
@Test
public void testCalculatePhasePosteriors() {
    final double averageHetsPerSegment = 100;
    final int numSegments = 100;
    final int averageDepth = 100;
    final double biasMean = 1.1;
    final double biasVariance = 0.01;
    final double outlierProbability = 0.02;
    final AlleleFractionSimulatedData simulatedData = new AlleleFractionSimulatedData(averageHetsPerSegment, numSegments, averageDepth, biasMean, biasVariance, outlierProbability);
    final SegmentedGenome segmentedGenome = simulatedData.getSegmentedGenome();
    final AlleleFractionState trueState = simulatedData.getTrueState();
    final AlleleFractionSimulatedData.PhaseIndicators truePhases = simulatedData.getTruePhases();
    final AllelicCountCollection counts = new AllelicCountCollection();
    //note that chromosomes are in lexicographical order
    segmentedGenome.getGenome().getSNPs().targets().stream().forEach(counts::add);
    final AllelicCountWithPhasePosteriorsCollection countsWithPhasePosteriors = CalculatePulldownPhasePosteriors.calculatePhasePosteriors(counts, segmentedGenome.getSegments(), trueState, AllelicPanelOfNormals.EMPTY_PON);
    int numIndicatorsCorrect = 0;
    //order is ALT_MINOR, REF_MINOR, OUTLIER
    final Iterator<AlleleFractionIndicator> phaseIterator = truePhases.iterator();
    final Iterator<AllelicCountWithPhasePosteriors> countWithPhasePosteriorsIterator = countsWithPhasePosteriors.getCounts().iterator();
    while (phaseIterator.hasNext() && countWithPhasePosteriorsIterator.hasNext()) {
        final AlleleFractionIndicator truePhase = phaseIterator.next();
        final AllelicCountWithPhasePosteriors countWithPhasePosteriors = countWithPhasePosteriorsIterator.next();
        final List<Double> phaseProbabilities = Arrays.asList(countWithPhasePosteriors.getAltMinorProb(), countWithPhasePosteriors.getRefMinorProb(), countWithPhasePosteriors.getOutlierProb());
        final int indexOfMaxProbPhase = phaseProbabilities.indexOf(Collections.max(phaseProbabilities));
        final AlleleFractionIndicator maxProbPhase = AlleleFractionIndicator.values()[indexOfMaxProbPhase];
        if (maxProbPhase.equals(truePhase)) {
            numIndicatorsCorrect++;
        }
    }
    final double fractionOfIndicatorsCorrect = (double) numIndicatorsCorrect / countsWithPhasePosteriors.getCounts().size();
    Assert.assertTrue(fractionOfIndicatorsCorrect >= FRACTION_OF_INDICATORS_CORRECT_THRESHOLD);
}
Also used : AllelicCountWithPhasePosteriorsCollection(org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountWithPhasePosteriorsCollection) AlleleFractionIndicator(org.broadinstitute.hellbender.tools.exome.allelefraction.AlleleFractionIndicator) AllelicCountWithPhasePosteriors(org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountWithPhasePosteriors) AlleleFractionState(org.broadinstitute.hellbender.tools.exome.allelefraction.AlleleFractionState) AlleleFractionSimulatedData(org.broadinstitute.hellbender.tools.exome.allelefraction.AlleleFractionSimulatedData) AllelicCountCollection(org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountCollection) Test(org.testng.annotations.Test) CommandLineProgramTest(org.broadinstitute.hellbender.CommandLineProgramTest)

Example 3 with AllelicCountWithPhasePosteriorsCollection

use of org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountWithPhasePosteriorsCollection in project gatk by broadinstitute.

the class CalculatePulldownPhasePosteriorsIntegrationTest method testRun.

private void testRun(final String[] arguments, final String outputFileName) {
    runCommandLine(arguments);
    //only check that file is created and can be read as an AllelicCountWithPhasePosteriorsCollection,
    //do not check for correctness of results (which is tested by testCalculatePhasePosteriors)
    final File outputFile = new File(outputFileName);
    Assert.assertTrue(outputFile.isFile(), outputFile.getAbsolutePath() + " is not a file.");
    Assert.assertTrue(outputFile.length() > 0);
    try {
        //check that file has:
        //  - at least two lines and all either start with "#" or contain at least one "\t"
        //  - at least two lines with tab (column names + 1 SNP)
        final List<String> outputLines = FileUtils.readLines(outputFile);
        Assert.assertTrue(outputLines.size() >= 2);
        Assert.assertEquals(outputLines.stream().filter(l -> l.contains("\t") || l.startsWith("#")).count(), outputLines.size());
        Assert.assertTrue(outputLines.stream().filter(l -> l.split("\t").length > 2 && !l.startsWith("#")).count() > 2, "File: " + outputFile + " does not seem to have at least one SNP and a header.");
        final AllelicCountWithPhasePosteriorsCollection counts = new AllelicCountWithPhasePosteriorsCollection(outputFile);
    } catch (final Exception e) {
        Assert.fail("Could not create AllelicCountWithPhasePosteriorsCollection from file: " + outputFile, e);
    }
}
Also used : ExomeStandardArgumentDefinitions(org.broadinstitute.hellbender.cmdline.ExomeStandardArgumentDefinitions) Arrays(java.util.Arrays) AllelicCountWithPhasePosteriors(org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountWithPhasePosteriors) Iterator(java.util.Iterator) AlleleFractionIndicator(org.broadinstitute.hellbender.tools.exome.allelefraction.AlleleFractionIndicator) FileUtils(org.apache.commons.io.FileUtils) StandardArgumentDefinitions(org.broadinstitute.hellbender.cmdline.StandardArgumentDefinitions) Test(org.testng.annotations.Test) AllelicCountWithPhasePosteriorsCollection(org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountWithPhasePosteriorsCollection) CommandLineProgramTest(org.broadinstitute.hellbender.CommandLineProgramTest) File(java.io.File) List(java.util.List) AlleleFractionState(org.broadinstitute.hellbender.tools.exome.allelefraction.AlleleFractionState) Assert(org.testng.Assert) AlleleFractionSimulatedData(org.broadinstitute.hellbender.tools.exome.allelefraction.AlleleFractionSimulatedData) AllelicPanelOfNormals(org.broadinstitute.hellbender.tools.pon.allelic.AllelicPanelOfNormals) AllelicCountCollection(org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountCollection) Collections(java.util.Collections) AllelicCountWithPhasePosteriorsCollection(org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountWithPhasePosteriorsCollection) File(java.io.File)

Example 4 with AllelicCountWithPhasePosteriorsCollection

use of org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountWithPhasePosteriorsCollection in project gatk by broadinstitute.

the class CalculatePulldownPhasePosteriorsIntegrationTest method testCalculatePhasePosteriors.

/**
     * Uses {@link AlleleFractionSimulatedData} to test recovery of phase indicators.  Phase with highest posterior
     * probability is compared to the true phase; we require that
     * {@link CalculatePulldownPhasePosteriorsIntegrationTest#FRACTION_OF_INDICATORS_CORRECT_THRESHOLD} of the
     * indicators are recovered correctly.
     */
@Test
public void testCalculatePhasePosteriors() {
    final double averageHetsPerSegment = 100;
    final int numSegments = 100;
    final int averageDepth = 100;
    final double biasMean = 1.1;
    final double biasVariance = 0.01;
    final double outlierProbability = 0.02;
    final AlleleFractionSimulatedData simulatedData = new AlleleFractionSimulatedData(averageHetsPerSegment, numSegments, averageDepth, biasMean, biasVariance, outlierProbability);
    final SegmentedGenome segmentedGenome = simulatedData.getSegmentedGenome();
    final AlleleFractionState trueState = simulatedData.getTrueState();
    final AlleleFractionSimulatedData.PhaseIndicators truePhases = simulatedData.getTruePhases();
    final AllelicCountCollection counts = new AllelicCountCollection();
    //note that chromosomes are in lexicographical order
    segmentedGenome.getGenome().getSNPs().targets().stream().forEach(counts::add);
    final AllelicCountWithPhasePosteriorsCollection countsWithPhasePosteriors = CalculatePulldownPhasePosteriors.calculatePhasePosteriors(counts, segmentedGenome.getSegments(), trueState, AllelicPanelOfNormals.EMPTY_PON);
    int numIndicatorsCorrect = 0;
    //order is ALT_MINOR, REF_MINOR, OUTLIER
    final Iterator<AlleleFractionIndicator> phaseIterator = truePhases.iterator();
    final Iterator<AllelicCountWithPhasePosteriors> countWithPhasePosteriorsIterator = countsWithPhasePosteriors.getCounts().iterator();
    while (phaseIterator.hasNext() && countWithPhasePosteriorsIterator.hasNext()) {
        final AlleleFractionIndicator truePhase = phaseIterator.next();
        final AllelicCountWithPhasePosteriors countWithPhasePosteriors = countWithPhasePosteriorsIterator.next();
        final List<Double> phaseProbabilities = Arrays.asList(countWithPhasePosteriors.getAltMinorProb(), countWithPhasePosteriors.getRefMinorProb(), countWithPhasePosteriors.getOutlierProb());
        final int indexOfMaxProbPhase = phaseProbabilities.indexOf(Collections.max(phaseProbabilities));
        final AlleleFractionIndicator maxProbPhase = AlleleFractionIndicator.values()[indexOfMaxProbPhase];
        if (maxProbPhase.equals(truePhase)) {
            numIndicatorsCorrect++;
        }
    }
    final double fractionOfIndicatorsCorrect = (double) numIndicatorsCorrect / countsWithPhasePosteriors.getCounts().size();
    Assert.assertTrue(fractionOfIndicatorsCorrect >= FRACTION_OF_INDICATORS_CORRECT_THRESHOLD);
}
Also used : AllelicCountWithPhasePosteriorsCollection(org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountWithPhasePosteriorsCollection) AlleleFractionIndicator(org.broadinstitute.hellbender.tools.exome.allelefraction.AlleleFractionIndicator) AllelicCountWithPhasePosteriors(org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountWithPhasePosteriors) AlleleFractionState(org.broadinstitute.hellbender.tools.exome.allelefraction.AlleleFractionState) AlleleFractionSimulatedData(org.broadinstitute.hellbender.tools.exome.allelefraction.AlleleFractionSimulatedData) AllelicCountCollection(org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountCollection) Test(org.testng.annotations.Test) CommandLineProgramTest(org.broadinstitute.hellbender.CommandLineProgramTest)

Example 5 with AllelicCountWithPhasePosteriorsCollection

use of org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountWithPhasePosteriorsCollection in project gatk-protected by broadinstitute.

the class CalculatePulldownPhasePosteriors method calculatePhasePosteriors.

@VisibleForTesting
protected static AllelicCountWithPhasePosteriorsCollection calculatePhasePosteriors(final AllelicCountCollection counts, final List<SimpleInterval> segments, final AlleleFractionState state, final AllelicPanelOfNormals allelicPoN) {
    final TargetCollection<SimpleInterval> segmentTargetCollection = new HashedListTargetCollection<>(segments);
    final AllelicCountWithPhasePosteriorsCollection countsWithPhasePosteriors = new AllelicCountWithPhasePosteriorsCollection();
    for (final AllelicCount count : counts.getCounts()) {
        final int segmentIndex = segmentTargetCollection.index(count.getInterval());
        if (segmentIndex < 0) {
            throw new UserException.EmptyIntersection(String.format("The AllelicCount at %s is not located within one of the input segments.", count.getInterval()));
        }
        final AlleleFractionGlobalParameters parameters = state.globalParameters();
        final double minorFraction = state.segmentMinorFraction(segmentIndex);
        final double refMinorLogProb = AlleleFractionLikelihoods.hetLogLikelihood(parameters, minorFraction, count, AlleleFractionIndicator.REF_MINOR, allelicPoN);
        final double altMinorLogProb = AlleleFractionLikelihoods.hetLogLikelihood(parameters, minorFraction, count, AlleleFractionIndicator.ALT_MINOR, allelicPoN);
        final double outlierLogProb = AlleleFractionLikelihoods.hetLogLikelihood(parameters, minorFraction, count, AlleleFractionIndicator.OUTLIER, allelicPoN);
        final AllelicCountWithPhasePosteriors countWithPhasePosteriors = new AllelicCountWithPhasePosteriors(count, refMinorLogProb, altMinorLogProb, outlierLogProb);
        countsWithPhasePosteriors.add(countWithPhasePosteriors);
    }
    return countsWithPhasePosteriors;
}
Also used : AllelicCountWithPhasePosteriorsCollection(org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountWithPhasePosteriorsCollection) AllelicCountWithPhasePosteriors(org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountWithPhasePosteriors) SimpleInterval(org.broadinstitute.hellbender.utils.SimpleInterval) AllelicCount(org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCount) VisibleForTesting(com.google.common.annotations.VisibleForTesting)

Aggregations

AllelicCountWithPhasePosteriorsCollection (org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountWithPhasePosteriorsCollection)8 AllelicCountCollection (org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountCollection)6 AllelicCountWithPhasePosteriors (org.broadinstitute.hellbender.tools.exome.alleliccount.AllelicCountWithPhasePosteriors)6 CommandLineProgramTest (org.broadinstitute.hellbender.CommandLineProgramTest)4 AlleleFractionIndicator (org.broadinstitute.hellbender.tools.exome.allelefraction.AlleleFractionIndicator)4 AlleleFractionSimulatedData (org.broadinstitute.hellbender.tools.exome.allelefraction.AlleleFractionSimulatedData)4 AlleleFractionState (org.broadinstitute.hellbender.tools.exome.allelefraction.AlleleFractionState)4 AllelicPanelOfNormals (org.broadinstitute.hellbender.tools.pon.allelic.AllelicPanelOfNormals)4 SimpleInterval (org.broadinstitute.hellbender.utils.SimpleInterval)4 Test (org.testng.annotations.Test)4 VisibleForTesting (com.google.common.annotations.VisibleForTesting)2 File (java.io.File)2 Arrays (java.util.Arrays)2 Collections (java.util.Collections)2 Iterator (java.util.Iterator)2 List (java.util.List)2 FileUtils (org.apache.commons.io.FileUtils)2 HDF5Library (org.broadinstitute.hdf5.HDF5Library)2 ExomeStandardArgumentDefinitions (org.broadinstitute.hellbender.cmdline.ExomeStandardArgumentDefinitions)2 StandardArgumentDefinitions (org.broadinstitute.hellbender.cmdline.StandardArgumentDefinitions)2