This is in fact a shallow copy if {@link GenotypeLikelihoodCalculators#ploidyLog10}
and is not meant to be modified by * this class. */ - private final double[] log10; + private final double[] ploidyLog10; /** * Buffer field use as a temporal container for sorted allele counts when calculating the likelihood of a @@ -202,24 +202,22 @@ public class GenotypeLikelihoodCalculator { * Creates a new calculator providing its ploidy and number of genotyping alleles. */ protected GenotypeLikelihoodCalculator(final int ploidy, final int alleleCount, - final int[][] alleleFirstGenotypeOffsetByPloidy, - final GenotypeAlleleCounts[][] genotypeTableByPloidy, - final double[] ploidyLog10) { + final int[][] alleleFirstGenotypeOffsetByPloidy, + final GenotypeAlleleCounts[][] genotypeTableByPloidy, + final double[] ploidyLog10) { + this.alleleFirstGenotypeOffsetByPloidy = alleleFirstGenotypeOffsetByPloidy; genotypeAlleleCounts = genotypeTableByPloidy[ploidy]; this.alleleCount = alleleCount; this.ploidy = ploidy; genotypeCount = this.alleleFirstGenotypeOffsetByPloidy[ploidy][alleleCount]; - if (genotypeCount == GenotypeLikelihoodCalculators.GENOTYPE_COUNT_OVERFLOW) - throw new IllegalArgumentException( - String.format("the combination of ploidy (%s) and number of alleles (%s) results in a very large number of genotypes (> %s). You need to limit ploidy or the number of alternative alleles to analyze this locus", - ploidy,alleleCount,Integer.MAX_VALUE)); + alleleHeap = new IntMaxHeap(ploidy); readLikelihoodsByGenotypeIndex = new double[genotypeCount][]; - log10 = ploidyLog10; + this.ploidyLog10 = ploidyLog10; // The number of possible components is limited by distinct allele count and ploidy. maximumDistinctAllelesInGenotype = Math.min(ploidy, alleleCount); - genotypeAllelesAndCounts = new int[maximumDistinctAllelesInGenotype << 1]; + genotypeAllelesAndCounts = new int[maximumDistinctAllelesInGenotype*2]; } /** @@ -349,7 +347,7 @@ public class GenotypeLikelihoodCalculator { */ private double[] genotypeLikelihoods(final double[][] readLikelihoodsByGenotypeIndex, final int readCount) { final double[] result = new double[genotypeCount]; - final double denominator = readCount * log10[ploidy]; // instead of dividing each read likelihood by ploidy + final double denominator = readCount * ploidyLog10[ploidy]; // instead of dividing each read likelihood by ploidy // ( so subtract log10(ploidy) ) we multiply them all and the divide by ploidy^readCount (so substract readCount * log10(ploidy) ) for (int g = 0; g < genotypeCount; g++) { final double[] likelihoodsByRead = readLikelihoodsByGenotypeIndex[g]; @@ -464,7 +462,9 @@ public class GenotypeLikelihoodCalculator { * exactly one allele present in the genotype. */ private void singleComponentGenotypeLikelihoodByRead(final GenotypeAlleleCounts genotypeAlleleCounts, - final double[] likelihoodByRead, final double[] readLikelihoodComponentsByAlleleCount, final int readCount) { + final double[] likelihoodByRead, + final double[] readLikelihoodComponentsByAlleleCount, + final int readCount) { final int allele = genotypeAlleleCounts.alleleIndexAt(0); // the count of the only component must be = ploidy. int offset = (allele * (ploidy + 1) + ploidy) * readCount; @@ -493,7 +493,7 @@ public class GenotypeLikelihoodCalculator { // p = 2 because the frequency == 1 we already have it. for (int frequency = 2, destinationOffset = frequency1Offset + readCount; frequency <= ploidy; frequency++) { - final double log10frequency = log10[frequency]; + final double log10frequency = ploidyLog10[frequency]; for (int r = 0, sourceOffset = frequency1Offset; r < readCount; r++) readAlleleLikelihoodByAlleleCount[destinationOffset++] = readAlleleLikelihoodByAlleleCount[sourceOffset++] + log10frequency; @@ -620,7 +620,11 @@ public class GenotypeLikelihoodCalculator { * @param destination where to store the new genotype index mapping to old. * @param sortedAlleleCountsBuffer a buffer to re-use to get the genotype-allele-count's sorted allele counts. */ - private void genotypeIndexMapPerGenotypeIndex(final int newGenotypeIndex, final GenotypeAlleleCounts alleleCounts, final int[] oldToNewAlleleIndexMap, final int[] destination, final int[] sortedAlleleCountsBuffer) { + private void genotypeIndexMapPerGenotypeIndex(final int newGenotypeIndex, + final GenotypeAlleleCounts alleleCounts, + final int[] oldToNewAlleleIndexMap, + final int[] destination, + final int[] sortedAlleleCountsBuffer) { final int distinctAlleleCount = alleleCounts.distinctAlleleCount(); alleleCounts.copyAlleleCounts(sortedAlleleCountsBuffer,0); for (int j = 0, jj = 0; j < distinctAlleleCount; j++) { diff --git a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/GenotypeLikelihoodCalculators.java b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/GenotypeLikelihoodCalculators.java index 83b280bff..da0ff0975 100644 --- a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/GenotypeLikelihoodCalculators.java +++ b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/genotyper/GenotypeLikelihoodCalculators.java @@ -51,7 +51,11 @@ package org.broadinstitute.gatk.tools.walkers.genotyper; +import org.broadinstitute.gatk.utils.MathUtils; +import org.broadinstitute.gatk.utils.exceptions.GATKException; + import java.util.Arrays; +import java.util.stream.IntStream; /** * Genotype likelihood calculator utility. @@ -116,7 +120,10 @@ public class GenotypeLikelihoodCalculators { private volatile static GenotypeAlleleCounts[][] genotypeTableByPloidy = buildGenotypeAlleleCountsTable(maximumPloidy,maximumAllele,alleleFirstGenotypeOffsetByPloidy); - + /** + * Cached log10 values for the first integer up to the maximum ploidy requested thus far. + */ + private volatile static double[] ploidyLog10 = IntStream.range(0, maximumPloidy + 1).mapToDouble(Math::log10).toArray(); /** * Build the table with the genotype offsets based on ploidy and the maximum allele index with representation @@ -291,40 +298,29 @@ public class GenotypeLikelihoodCalculators { return result; } - /** - * Cached log10 values for the first integer up to the maximum ploidy requested thus far. - */ - private volatile static double[] ploidyLog10; - - // Initialize {@link #ploidyLog10}. - static { - ploidyLog10 = new double[maximumPloidy + 1]; - for (int i = 0; i <= maximumPloidy; i++) - ploidyLog10[i] = Math.log10(i); - } - /** * Returns an instance given its ploidy and the number of alleles. * * @param alleleCount the required allele-count. * @param ploidy the required ploidy-count. * - * @throws IllegalArgumentException if either {@code ploidy} or {@code alleleCount} is {@code null}, or - * the resulting number of genotypes is too large. + * @throws IllegalArgumentException if either {@code ploidy} or {@code alleleCount} is negative, or the resulting number of genotypes is too large. * * @return never {@code null}. */ - public static GenotypeLikelihoodCalculator getInstance(final int ploidy, - final int alleleCount) { + public static GenotypeLikelihoodCalculator getInstance(final int ploidy, final int alleleCount) { checkPloidyAndMaximumAllele(ploidy, alleleCount); // Non-thread safe (fast) check on tables capacities, - // if not enough capacity we expand the tables in a thread-safe manner: - if (alleleCount > maximumAllele || ploidy > maximumPloidy) - ensureCapacity(alleleCount, ploidy); + // if not enough capacity we expand the tables in a thread-safe manner + // also checks if the requested ploidy and allele count result in a genotype count too large to deal with + if(calculateGenotypeCountUsingTables(ploidy, alleleCount) == GENOTYPE_COUNT_OVERFLOW){ + final double largeGenotypeCount = MathUtils.binomialCoefficient(ploidy + alleleCount - 1, alleleCount - 1); + throw new IllegalArgumentException(String.format("the number of genotypes is too large for ploidy %d and allele %d: approx. %.0f", ploidy, alleleCount, largeGenotypeCount)); + } // At this point the tables must have at least the requested capacity, likely to be much more. - return new GenotypeLikelihoodCalculator(ploidy,alleleCount,alleleFirstGenotypeOffsetByPloidy,genotypeTableByPloidy,ploidyLog10); + return new GenotypeLikelihoodCalculator(ploidy, alleleCount, alleleFirstGenotypeOffsetByPloidy, genotypeTableByPloidy, ploidyLog10); } /** @@ -413,14 +409,59 @@ public class GenotypeLikelihoodCalculators { * @param ploidy the requested ploidy. * @param alleleCount the requested number of alleles. * - * @throws IllegalArgumentException if {@code ploidy} or {@code alleleCount} is negative. + * @throws IllegalArgumentException if {@code ploidy} or {@code alleleCount} is negative or + * the number of genotypes is too large (more than {@link Integer#MAX_VALUE}). * - * @return 0 or greater. + * @return the number of genotypes given ploidy and allele count (0 or greater). */ public final static int genotypeCount(final int ploidy, final int alleleCount) { + + final int result = calculateGenotypeCountUsingTables(ploidy, alleleCount); + if (result == GENOTYPE_COUNT_OVERFLOW) { + final double largeGenotypeCount = MathUtils.binomialCoefficient(ploidy + alleleCount - 1, alleleCount - 1); + throw new IllegalArgumentException(String.format("the number of genotypes is too large for ploidy %d and allele %d: approx. %.0f", ploidy, alleleCount, largeGenotypeCount)); + } + return result; + } + + /** + * Compute the maximally acceptable allele count (ref allele included) given the maximally acceptable genotype count. + * @param ploidy sample ploidy + * @param maxGenotypeCount maximum number of genotype count used to calculate upper bound on number of alleles given ploidy + * @throws IllegalArgumentException if {@code ploidy} or {@code alleleCount} is negative. + * @return the maximally acceptable allele count given ploidy and maximum number of genotypes acceptable + */ + public static int computeMaxAcceptableAlleleCount(final int ploidy, final int maxGenotypeCount){ + + checkPloidyAndMaximumAllele(ploidy, ploidy); // a hack to check ploidy makes sense (could duplicate code but choice must be made) + + final double log10MaxGenotypeCount = Math.log10(maxGenotypeCount); + + // Math explanation: genotype count is determined by ${P+A-1 \choose A-1}$, this leads to constraint + // $\log(\frac{(P+A-1)!}{(A-1)!}) \le \log(P!G)$, + // where $P$ is ploidy, $A$ is allele count, and $G$ is maxGenotypeCount + // The upper and lower bounds of the left hand side of the constraint are $P \log(A-1+P)$ and $P \log(A)$ + // which require $A$ to be searched in interval $[10^{\log(P!G)/P} - (P-1), 10^{\log(P!G)/P}]$ + // Denote $[10^{\log(P!G)/P}$ as $x$ in the code. + + final double x = Math.pow(10, (MathUtils.log10Factorial(ploidy) + log10MaxGenotypeCount)/ploidy ); + final int lower = (int)Math.floor(x) - ploidy - 1; + final int upper = (int)Math.ceil(x); + for(int a=upper; a>=lower; --a){// check one by one + + final double log10GTCnt = MathUtils.log10BinomialCoefficient(ploidy+a-1, a-1); + if(log10MaxGenotypeCount >= log10GTCnt) { + return a; + } + } + throw new GATKException("Code should never reach here."); + } + + private static int calculateGenotypeCountUsingTables(int ploidy, int alleleCount) { checkPloidyAndMaximumAllele(ploidy, alleleCount); - if (ploidy > maximumPloidy || alleleCount > maximumAllele) - ensureCapacity(alleleCount,ploidy); + if (ploidy > maximumPloidy || alleleCount > maximumAllele) { + ensureCapacity(alleleCount, ploidy); + } return alleleFirstGenotypeOffsetByPloidy[ploidy][alleleCount]; } } \ No newline at end of file diff --git a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/haplotypecaller/HaplotypeCallerGenotypingEngine.java b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/haplotypecaller/HaplotypeCallerGenotypingEngine.java index 30e7efe22..ad0f2b773 100644 --- a/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/haplotypecaller/HaplotypeCallerGenotypingEngine.java +++ b/protected/gatk-tools-protected/src/main/java/org/broadinstitute/gatk/tools/walkers/haplotypecaller/HaplotypeCallerGenotypingEngine.java @@ -56,6 +56,8 @@ import com.google.java.contract.Ensures; import com.google.java.contract.Requires; import htsjdk.samtools.util.StringUtil; import htsjdk.variant.variantcontext.*; +import org.apache.commons.lang.ArrayUtils; +import org.apache.commons.math3.stat.StatUtils; import org.broadinstitute.gatk.engine.arguments.GenotypeCalculationArgumentCollection; import org.broadinstitute.gatk.utils.*; import org.broadinstitute.gatk.utils.contexts.ReferenceContext; @@ -73,6 +75,7 @@ import org.broadinstitute.gatk.utils.variant.GATKVCFConstants; import org.broadinstitute.gatk.utils.variant.GATKVariantContextUtils; import java.util.*; +import java.util.stream.Collectors; /** * {@link HaplotypeCaller}'s genotyping strategy implementation. @@ -84,6 +87,9 @@ public class HaplotypeCallerGenotypingEngine extends GenotypingEngine
@@ -462,7 +579,7 @@ public class HaplotypeCallerGenotypingEngine extends GenotypingEngine