Examples with AlleleFractionGlobalParameters org.broadinstitute.hellbender.tools.exome.allelefraction.AlleleFractionGlobalParameters used on opensource projects

Example 1 with AlleleFractionGlobalParameters

use of org.broadinstitute.hellbender.tools.exome.allelefraction.AlleleFractionGlobalParameters in project gatk by broadinstitute.

the class AlleleFractionSegmenter method relearnAdditionalParameters.

@Override
protected void relearnAdditionalParameters(final ExpectationStep eStep) {
    final Function<AlleleFractionGlobalParameters, Double> emissionLogLikelihood = params -> {
        double logLikelihood = 0.0;
        for (int position = 0; position < numPositions(); position++) {
            for (int state = 0; state < numStates(); state++) {
                final double eStepPosterior = eStep.pStateAtPosition(state, position);
                logLikelihood += eStepPosterior < NEGLIGIBLE_POSTERIOR_FOR_M_STEP ? 0 : eStepPosterior * AlleleFractionHMM.logEmissionProbability(data.get(position), getState(state), params, allelicPoN);
            }
        }
        return logLikelihood;
    };
    final Function<Double, Double> meanBiasObjective = mean -> emissionLogLikelihood.apply(globalParameters.copyWithNewMeanBias(mean));
    final double newMeanBias = OptimizationUtils.argmax(meanBiasObjective, 0, AlleleFractionInitializer.MAX_REASONABLE_MEAN_BIAS, globalParameters.getMeanBias(), RELATIVE_TOLERANCE_FOR_OPTIMIZATION, ABSOLUTE_TOLERANCE_FOR_OPTIMIZATION, MAX_EVALUATIONS_FOR_OPTIMIZATION);
    final Function<Double, Double> biasVarianceObjective = variance -> emissionLogLikelihood.apply(globalParameters.copyWithNewBiasVariance(variance));
    final double newBiasVariance = OptimizationUtils.argmax(biasVarianceObjective, 0, AlleleFractionInitializer.MAX_REASONABLE_BIAS_VARIANCE, globalParameters.getBiasVariance(), RELATIVE_TOLERANCE_FOR_OPTIMIZATION, ABSOLUTE_TOLERANCE_FOR_OPTIMIZATION, MAX_EVALUATIONS_FOR_OPTIMIZATION);
    final Function<Double, Double> outlierProbabilityObjective = pOutlier -> emissionLogLikelihood.apply(globalParameters.copyWithNewOutlierProbability(pOutlier));
    final double newOutlierProbability = OptimizationUtils.argmax(outlierProbabilityObjective, 0, AlleleFractionInitializer.MAX_REASONABLE_OUTLIER_PROBABILITY, globalParameters.getOutlierProbability(), RELATIVE_TOLERANCE_FOR_OPTIMIZATION, ABSOLUTE_TOLERANCE_FOR_OPTIMIZATION, MAX_EVALUATIONS_FOR_OPTIMIZATION);
    globalParameters = new AlleleFractionGlobalParameters(newMeanBias, newBiasVariance, newOutlierProbability);
    logger.info(String.format("Global allelic bias parameters learned.  Mean allelic bias: %f, variance of allelic bias: %f, outlier probability: %f.", newMeanBias, newBiasVariance, newOutlierProbability));
}

Example 2 with AlleleFractionGlobalParameters

use of org.broadinstitute.hellbender.tools.exome.allelefraction.AlleleFractionGlobalParameters in project gatk-protected by broadinstitute.

the class AlleleFractionSegmenter method relearnAdditionalParameters.

@Override
protected void relearnAdditionalParameters(final ExpectationStep eStep) {
    final Function<AlleleFractionGlobalParameters, Double> emissionLogLikelihood = params -> {
        double logLikelihood = 0.0;
        for (int position = 0; position < numPositions(); position++) {
            for (int state = 0; state < numStates(); state++) {
                final double eStepPosterior = eStep.pStateAtPosition(state, position);
                logLikelihood += eStepPosterior < NEGLIGIBLE_POSTERIOR_FOR_M_STEP ? 0 : eStepPosterior * AlleleFractionHMM.logEmissionProbability(data.get(position), getState(state), params, allelicPoN);
            }
        }
        return logLikelihood;
    };
    final Function<Double, Double> meanBiasObjective = mean -> emissionLogLikelihood.apply(globalParameters.copyWithNewMeanBias(mean));
    final double newMeanBias = OptimizationUtils.argmax(meanBiasObjective, 0, AlleleFractionInitializer.MAX_REASONABLE_MEAN_BIAS, globalParameters.getMeanBias(), RELATIVE_TOLERANCE_FOR_OPTIMIZATION, ABSOLUTE_TOLERANCE_FOR_OPTIMIZATION, MAX_EVALUATIONS_FOR_OPTIMIZATION);
    final Function<Double, Double> biasVarianceObjective = variance -> emissionLogLikelihood.apply(globalParameters.copyWithNewBiasVariance(variance));
    final double newBiasVariance = OptimizationUtils.argmax(biasVarianceObjective, 0, AlleleFractionInitializer.MAX_REASONABLE_BIAS_VARIANCE, globalParameters.getBiasVariance(), RELATIVE_TOLERANCE_FOR_OPTIMIZATION, ABSOLUTE_TOLERANCE_FOR_OPTIMIZATION, MAX_EVALUATIONS_FOR_OPTIMIZATION);
    final Function<Double, Double> outlierProbabilityObjective = pOutlier -> emissionLogLikelihood.apply(globalParameters.copyWithNewOutlierProbability(pOutlier));
    final double newOutlierProbability = OptimizationUtils.argmax(outlierProbabilityObjective, 0, AlleleFractionInitializer.MAX_REASONABLE_OUTLIER_PROBABILITY, globalParameters.getOutlierProbability(), RELATIVE_TOLERANCE_FOR_OPTIMIZATION, ABSOLUTE_TOLERANCE_FOR_OPTIMIZATION, MAX_EVALUATIONS_FOR_OPTIMIZATION);
    globalParameters = new AlleleFractionGlobalParameters(newMeanBias, newBiasVariance, newOutlierProbability);
    logger.info(String.format("Global allelic bias parameters learned.  Mean allelic bias: %f, variance of allelic bias: %f, outlier probability: %f.", newMeanBias, newBiasVariance, newOutlierProbability));
}

Example 3 with AlleleFractionGlobalParameters

use of org.broadinstitute.hellbender.tools.exome.allelefraction.AlleleFractionGlobalParameters in project gatk-protected by broadinstitute.

the class AlleleFractionHMMUnitTest method constructorTest.

@Test
public void constructorTest() {
    final double memoryLength = 1e6;
    final List<Double> minorAlleleFractions = Arrays.asList(0.1, 0.5, 0.23);
    final List<Double> weights = Arrays.asList(0.2, 0.2, 0.6);
    final AlleleFractionGlobalParameters params = new AlleleFractionGlobalParameters(0.1, 0.01, 0.03);
    final AlleleFractionHMM model = new AlleleFractionHMM(minorAlleleFractions, weights, memoryLength, AllelicPanelOfNormals.EMPTY_PON, params);
    Assert.assertEquals(memoryLength, model.getMemoryLength());
    for (int n = 0; n < weights.size(); n++) {
        Assert.assertEquals(weights.get(n), model.getWeight(n));
        Assert.assertEquals(minorAlleleFractions.get(n), model.getMinorAlleleFraction(n));
    }
    Assert.assertEquals(model.getParameters().getMeanBias(), params.getMeanBias());
    Assert.assertEquals(model.getParameters().getBiasVariance(), params.getBiasVariance());
    Assert.assertEquals(model.getParameters().getOutlierProbability(), params.getOutlierProbability());
}

Example 4 with AlleleFractionGlobalParameters

use of org.broadinstitute.hellbender.tools.exome.allelefraction.AlleleFractionGlobalParameters in project gatk-protected by broadinstitute.

the class AlleleFractionSegmenterUnitTest method testSegmentation.

@Test
public void testSegmentation() {
    final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(563));
    final List<Double> trueWeights = Arrays.asList(0.2, 0.5, 0.3);
    final List<Double> trueMinorAlleleFractions = Arrays.asList(0.12, 0.32, 0.5);
    final double trueMemoryLength = 1e5;
    final AlleleFractionGlobalParameters trueParams = new AlleleFractionGlobalParameters(1.0, 0.01, 0.01);
    final AlleleFractionHMM trueModel = new AlleleFractionHMM(trueMinorAlleleFractions, trueWeights, trueMemoryLength, AllelicPanelOfNormals.EMPTY_PON, trueParams);
    // randomly set positions
    final int chainLength = 10000;
    final List<SimpleInterval> positions = CopyRatioSegmenterUnitTest.randomPositions("chr1", chainLength, rng, trueMemoryLength / 4);
    final List<Integer> trueStates = trueModel.generateHiddenStateChain(positions);
    final List<Double> truthMinorFractions = trueStates.stream().map(trueModel::getMinorAlleleFraction).collect(Collectors.toList());
    final AllelicCountCollection counts = generateCounts(truthMinorFractions, positions, rng, trueParams);
    final AlleleFractionSegmenter segmenter = new AlleleFractionSegmenter(10, counts, AllelicPanelOfNormals.EMPTY_PON);
    final List<ModeledSegment> segments = segmenter.getModeledSegments();
    final double[] segmentMinorFractions = segments.stream().flatMap(s -> Collections.nCopies((int) s.getTargetCount(), s.getSegmentMean()).stream()).mapToDouble(x -> x).toArray();
    final double averageMinorFractionError = IntStream.range(0, truthMinorFractions.size()).mapToDouble(n -> Math.abs(segmentMinorFractions[n] - truthMinorFractions.get(n))).average().getAsDouble();
    Assert.assertEquals(averageMinorFractionError, 0, 0.01);
}

Example 5 with AlleleFractionGlobalParameters

use of org.broadinstitute.hellbender.tools.exome.allelefraction.AlleleFractionGlobalParameters in project gatk by broadinstitute.

the class AlleleFractionSegmenterUnitTest method testChromosomesOnDifferentSegments.

@Test
public void testChromosomesOnDifferentSegments() {
    final RandomGenerator rng = RandomGeneratorFactory.createRandomGenerator(new Random(563));
    final double[] trueMinorAlleleFractions = new double[] { 0.12, 0.32, 0.5 };
    final double trueMemoryLength = 1e5;
    final AlleleFractionGlobalParameters trueParams = new AlleleFractionGlobalParameters(1.0, 0.01, 0.01);
    // randomly set positions
    final int chainLength = 100;
    final List<SimpleInterval> positions = CopyRatioSegmenterUnitTest.randomPositions("chr1", chainLength, rng, trueMemoryLength / 4);
    positions.addAll(CopyRatioSegmenterUnitTest.randomPositions("chr2", chainLength, rng, trueMemoryLength / 4));
    positions.addAll(CopyRatioSegmenterUnitTest.randomPositions("chr3", chainLength, rng, trueMemoryLength / 4));
    //fix everything to the same state 2
    final int trueState = 2;
    final List<Double> minorAlleleFractionSequence = Collections.nCopies(positions.size(), trueMinorAlleleFractions[trueState]);
    final AllelicCountCollection counts = generateCounts(minorAlleleFractionSequence, positions, rng, trueParams);
    final AlleleFractionSegmenter segmenter = new AlleleFractionSegmenter(10, counts, AllelicPanelOfNormals.EMPTY_PON);
    final List<ModeledSegment> segments = segmenter.getModeledSegments();
    //check that each chromosome has at least one segment
    final int numDifferentContigsInSegments = (int) segments.stream().map(ModeledSegment::getContig).distinct().count();
    Assert.assertEquals(numDifferentContigsInSegments, 3);
}