<commit_before>package org.broadinstitute.hellbender.tools.walkers.contamination;

import htsjdk.samtools.util.OverlapDetector;
import org.apache.commons.lang.mutable.MutableDouble;
import org.apache.commons.lang3.tuple.ImmutablePair;
import org.apache.commons.lang3.tuple.Pair;
import org.apache.commons.math3.optim.univariate.UnivariatePointValuePair;
import org.apache.commons.math3.stat.descriptive.rank.Percentile;
import org.apache.commons.math3.util.FastMath;
import org.broadinstitute.hellbender.exceptions.GATKException;
import org.broadinstitute.hellbender.tools.walkers.qc.Pileup;
import org.broadinstitute.hellbender.utils.*;

import java.util.*;
import java.util.function.DoubleUnaryOperator;
import java.util.function.ToDoubleFunction;
import java.util.function.ToIntFunction;
import java.util.stream.Collectors;
import java.util.stream.IntStream;

/**
 * This is the probabilistic contamination model that we use to distinguish homs from hets
 * The model is similar to that of ContEst, in that it assumes that each contaminant read is independently drawn from
 * the population.  However, it also accounts for allelic CNVs and doesn't just apply to hom alt sites.
 *
 * Because this model makes assumptions about the form of contamination and doesn't model CNVs precisely (although it
 * is robust to them in that they don't lead to bad solutions), we don't actually use it for the final contamination
 * estimate.  Instead, we learn the model, then use it to predict which sites are hom alts.  It's important to do this right
 * because contamination makes hom alts look more like hets, and allelic CNVs make hets look more like hom alts, hence
 * distinguishing correctly is not trivial.
 */
public class ContaminationModel {
    public static final double INITIAL_MAF_THRESHOLD = 0.40;
    public static final double MAF_TO_SWITCH_TO_HOM_REF = 0.25;
    public static final double MAF_TO_SWITCH_TO_UNSCRUPULOUS_HOM_REF = 0.20;
    public static final double UNSCRUPULOUS_HOM_REF_ALLELE_FRACTION = 0.15;
    public static final double UNSCRUPULOUS_HOM_REF_FRACTION_TO_REMOVE_FOR_POSSIBLE_LOH = 0.1;
    public static final double UNSCRUPULOUS_HOM_REF_PERCENTILE = 100 * ( 1 - UNSCRUPULOUS_HOM_REF_FRACTION_TO_REMOVE_FOR_POSSIBLE_LOH);
    public static final double MAF_STEP_SIZE = 0.04;
    private final double contamination;
    private final double errorRate;
    private final List<Double> minorAlleleFractions;
    private final List<List<PileupSummary>> segments;

    public static final int HOM_REF = 0;
    public static final int HOM_ALT = 3;
    private static final int NUM_ITERATIONS = 3;
    private static final double MIN_FRACTION_OF_SITES_TO_USE = 0.25;
    private static final double MIN_RELATIVE_ERROR = 0.2;
    private static final List<Double> CONTAMINATION_INITIAL_GUESSES = Arrays.asList(0.02, 0.05, 0.1, 0.2);

    public ContaminationModel(List<PileupSummary> sites) {
        errorRate = calculateErrorRate(sites);

        // partition genome into minor allele fraction (MAF) segments to better distinguish hom alts from LoH hets.
        segments = ContaminationSegmenter.findSegments(sites);
        final int numSegments = segments.size();

        final List<Double> minorAlleleFractionsGuess = new ArrayList<>(Collections.nCopies(segments.size(), 0.5));
        final MutableDouble contaminationGuess = new MutableDouble(0);
        for (int n = 0; n < NUM_ITERATIONS; n++) {
            IntStream.range(0, numSegments).forEach(s -> minorAlleleFractionsGuess.set(s, calculateMinorAlleleFraction(contaminationGuess.doubleValue(), errorRate, segments.get(s))));
            final Pair<List<List<PileupSummary>>, List<Double>> nonLOHSegmentsAndMafs = getNonLOHSegments(segments, minorAlleleFractionsGuess);
            contaminationGuess.setValue(calculateContamination(errorRate, nonLOHSegmentsAndMafs.getLeft(), nonLOHSegmentsAndMafs.getRight()));
        }

        minorAlleleFractions = minorAlleleFractionsGuess;
        contamination = contaminationGuess.doubleValue();
    }

    private static Pair<List<List<PileupSummary>>, List<Double>> getNonLOHSegments(final List<List<PileupSummary>> segments, final List<Double> mafs) {
        final int numSites = segments.stream().mapToInt(List::size).sum();
        for (double minMaf = INITIAL_MAF_THRESHOLD; minMaf > 0; minMaf -= MAF_STEP_SIZE) {
            final double threshold = minMaf;
            final int[] nonLOHIndices = IntStream.range(0, segments.size()).filter(s -> mafs.get(s) > threshold).toArray();
            final List<List<PileupSummary>> nonLOHSegments = Arrays.stream(nonLOHIndices).mapToObj(segments::get).collect(Collectors.toList());
            final List<Double> nonLOHMafs = Arrays.stream(nonLOHIndices).mapToObj(mafs::get).collect(Collectors.toList());

            final int numNonLOHSites = nonLOHSegments.stream().mapToInt(List::size).sum();
            if ((double) numNonLOHSites / numSites > MIN_FRACTION_OF_SITES_TO_USE) {
                return ImmutablePair.of(nonLOHSegments, nonLOHMafs);
            }
        }
        return ImmutablePair.of(segments, mafs);
    }

    /**
     * Calculate the contamination in a tumor sample using the hom alt and hom ref sites inferred from this model, which
     * could be derived from the tumor itself or a matched normal.
     * @return
     */
    public Pair<Double, Double> calculateContaminationFromHoms(final List<PileupSummary> tumorSites) {
        for (double minMaf = INITIAL_MAF_THRESHOLD; minMaf >= 0; minMaf -= MAF_STEP_SIZE) {

            final Pair<Double, Double> result;
            if (minMaf > MAF_TO_SWITCH_TO_HOM_REF) {
                result = calculateContamination(Strategy.HOM_ALT, tumorSites, minMaf);
            } else if (minMaf > MAF_TO_SWITCH_TO_UNSCRUPULOUS_HOM_REF) {
                result = calculateContamination(Strategy.HOM_REF, tumorSites, minMaf);
            } else {
                result = calculateContamination(Strategy.UNSCRUPULOUS_HOM_REF, tumorSites, minMaf);
            }
            if (result.getRight() < result.getLeft() * MIN_RELATIVE_ERROR) {
                return result;
            }
        }

        return calculateContamination(Strategy.UNSCRUPULOUS_HOM_REF, tumorSites, 0);
    }

    /**
     * Depending on how many sites we have found for use in our estimate, we are more or less strict as to which calculation<commit_msg>fixed a bug in CalculateContamination when contamination is indistinguishable from zero (#5971)<commit_after>package org.broadinstitute.hellbender.tools.walkers.contamination;

import htsjdk.samtools.util.OverlapDetector;
import org.apache.commons.lang.mutable.MutableDouble;
import org.apache.commons.lang3.tuple.ImmutablePair;
import org.apache.commons.lang3.tuple.Pair;
import org.apache.commons.math3.optim.univariate.UnivariatePointValuePair;
import org.apache.commons.math3.stat.descriptive.rank.Percentile;
import org.apache.commons.math3.util.FastMath;
import org.broadinstitute.hellbender.exceptions.GATKException;
import org.broadinstitute.hellbender.tools.walkers.qc.Pileup;
import org.broadinstitute.hellbender.utils.*;

import java.util.*;
import java.util.function.DoubleUnaryOperator;
import java.util.function.ToDoubleFunction;
import java.util.function.ToIntFunction;
import java.util.stream.Collectors;
import java.util.stream.IntStream;

/**
 * This is the probabilistic contamination model that we use to distinguish homs from hets
 * The model is similar to that of ContEst, in that it assumes that each contaminant read is independently drawn from
 * the population.  However, it also accounts for allelic CNVs and doesn't just apply to hom alt sites.
 *
 * Because this model makes assumptions about the form of contamination and doesn't model CNVs precisely (although it
 * is robust to them in that they don't lead to bad solutions), we don't actually use it for the final contamination
 * estimate.  Instead, we learn the model, then use it to predict which sites are hom alts.  It's important to do this right
 * because contamination makes hom alts look more like hets, and allelic CNVs make hets look more like hom alts, hence
 * distinguishing correctly is not trivial.
 */
public class ContaminationModel {
    public static final double INITIAL_MAF_THRESHOLD = 0.40;
    public static final double MAF_TO_SWITCH_TO_HOM_REF = 0.25;
    public static final double MAF_TO_SWITCH_TO_UNSCRUPULOUS_HOM_REF = 0.20;
    public static final double UNSCRUPULOUS_HOM_REF_ALLELE_FRACTION = 0.15;
    public static final double UNSCRUPULOUS_HOM_REF_FRACTION_TO_REMOVE_FOR_POSSIBLE_LOH = 0.1;
    public static final double UNSCRUPULOUS_HOM_REF_PERCENTILE = 100 * ( 1 - UNSCRUPULOUS_HOM_REF_FRACTION_TO_REMOVE_FOR_POSSIBLE_LOH);
    public static final double MAF_STEP_SIZE = 0.04;
    private final double contamination;
    private final double errorRate;
    private final List<Double> minorAlleleFractions;
    private final List<List<PileupSummary>> segments;

    public static final int HOM_REF = 0;
    public static final int HOM_ALT = 3;
    private static final int NUM_ITERATIONS = 3;
    private static final double MIN_FRACTION_OF_SITES_TO_USE = 0.25;
    private static final double MIN_RELATIVE_ERROR = 0.2;
    private static final double MIN_ABSOLUTE_ERROR = 0.001;
    private static final List<Double> CONTAMINATION_INITIAL_GUESSES = Arrays.asList(0.02, 0.05, 0.1, 0.2);

    public ContaminationModel(List<PileupSummary> sites) {
        errorRate = calculateErrorRate(sites);

        // partition genome into minor allele fraction (MAF) segments to better distinguish hom alts from LoH hets.
        segments = ContaminationSegmenter.findSegments(sites);
        final int numSegments = segments.size();

        final List<Double> minorAlleleFractionsGuess = new ArrayList<>(Collections.nCopies(segments.size(), 0.5));
        final MutableDouble contaminationGuess = new MutableDouble(0);
        for (int n = 0; n < NUM_ITERATIONS; n++) {
            IntStream.range(0, numSegments).forEach(s -> minorAlleleFractionsGuess.set(s, calculateMinorAlleleFraction(contaminationGuess.doubleValue(), errorRate, segments.get(s))));
            final Pair<List<List<PileupSummary>>, List<Double>> nonLOHSegmentsAndMafs = getNonLOHSegments(segments, minorAlleleFractionsGuess);
            contaminationGuess.setValue(calculateContamination(errorRate, nonLOHSegmentsAndMafs.getLeft(), nonLOHSegmentsAndMafs.getRight()));
        }

        minorAlleleFractions = minorAlleleFractionsGuess;
        contamination = contaminationGuess.doubleValue();
    }

    private static Pair<List<List<PileupSummary>>, List<Double>> getNonLOHSegments(final List<List<PileupSummary>> segments, final List<Double> mafs) {
        final int numSites = segments.stream().mapToInt(List::size).sum();
        for (double minMaf = INITIAL_MAF_THRESHOLD; minMaf > 0; minMaf -= MAF_STEP_SIZE) {
            final double threshold = minMaf;
            final int[] nonLOHIndices = IntStream.range(0, segments.size()).filter(s -> mafs.get(s) > threshold).toArray();
            final List<List<PileupSummary>> nonLOHSegments = Arrays.stream(nonLOHIndices).mapToObj(segments::get).collect(Collectors.toList());
            final List<Double> nonLOHMafs = Arrays.stream(nonLOHIndices).mapToObj(mafs::get).collect(Collectors.toList());

            final int numNonLOHSites = nonLOHSegments.stream().mapToInt(List::size).sum();
            if ((double) numNonLOHSites / numSites > MIN_FRACTION_OF_SITES_TO_USE) {
                return ImmutablePair.of(nonLOHSegments, nonLOHMafs);
            }
        }
        return ImmutablePair.of(segments, mafs);
    }

    /**
     * Calculate the contamination in a tumor sample using the hom alt and hom ref sites inferred from this model, which
     * could be derived from the tumor itself or a matched normal.
     * @return
     */
    public Pair<Double, Double> calculateContaminationFromHoms(final List<PileupSummary> tumorSites) {
        for (double minMaf = INITIAL_MAF_THRESHOLD; minMaf >= 0; minMaf -= MAF_STEP_SIZE) {

            final Pair<Double, Double> result;
            if (minMaf > MAF_TO_SWITCH_TO_HOM_REF) {
                result = calculateContamination(Strategy.HOM_ALT, tumorSites, minMaf);
            } else if (minMaf > MAF_TO_SWITCH_TO_UNSCRUPULOUS_HOM_REF) {
                result = calculateContamination(Strategy.HOM_REF, tumorSites, minMaf);
            } else {
                result = calculateContamination(Strategy.UNSCRUPULOUS_HOM_REF, tumorSites, minMaf);
            }
            if (result.getRight() < (result.getLeft() * MIN_RELATIVE_ERROR + MIN_ABSOLUTE_ERROR)) {
                return result;
            }
        }

        return calculateContamination(Strategy.UNSCRUPULOUS_HOM_REF, tumorSites, 0);
    }

    /**
     * Depending on how many sites we have found for use in our estimate, we are more or less strict as to which calculation
{'fraction_non_alphanumeric': 0.08286266924564797, 'fraction_numerical': 0.008820116054158608, 'mean_word_length': 3.9073652239939256, 'pattern_counts': {'":': 0, '<': 55, '<?xml version=': 0, '>': 67, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6652715', 'n_tokens_mistral': 4142, 'n_tokens_neox': 3904, 'n_words': 1007}