/** * get the allele index pair for the given PL using the deprecated PL ordering: * AA,AB,AC,AD,BB,BC,BD,CC,CD,DD instead of AA,AB,BB,AC,BC,CC,AD,BD,CD,DD. * Although it's painful to keep this conversion around, our DiploidSNPGenotypeLikelihoods class uses the deprecated * ordering and I know with certainty that external users have built code on top of it; changing it now would * cause a whole lot of heartache for our collaborators, so for now at least there's a standard conversion method. * This method assumes at most 3 alternate alleles. * * @param PLindex the PL index * @return the allele index pair */ @Deprecated public static GenotypeLikelihoodsAllelePair getAllelePairUsingDeprecatedOrdering(final int PLindex) { return getAllelePair(PLindexConversion[PLindex]); }
/** * get the allele index pair for the given PL using the deprecated PL ordering: * AA,AB,AC,AD,BB,BC,BD,CC,CD,DD instead of AA,AB,BB,AC,BC,CC,AD,BD,CD,DD. * Although it's painful to keep this conversion around, our DiploidSNPGenotypeLikelihoods class uses the deprecated * ordering and I know with certainty that external users have built code on top of it; changing it now would * cause a whole lot of heartache for our collaborators, so for now at least there's a standard conversion method. * This method assumes at most 3 alternate alleles. * * @param PLindex the PL index * @return the allele index pair * * @deprecated since 2/5/13 */ @Deprecated public static GenotypeLikelihoodsAllelePair getAllelePairUsingDeprecatedOrdering(final int PLindex) { return getAllelePair(PLindexConversion[PLindex]); }
/** * get the allele index pair for the given PL using the deprecated PL ordering: * AA,AB,AC,AD,BB,BC,BD,CC,CD,DD instead of AA,AB,BB,AC,BC,CC,AD,BD,CD,DD. * Although it's painful to keep this conversion around, our DiploidSNPGenotypeLikelihoods class uses the deprecated * ordering and I know with certainty that external users have built code on top of it; changing it now would * cause a whole lot of heartache for our collaborators, so for now at least there's a standard conversion method. * This method assumes at most 3 alternate alleles. * * @param PLindex the PL index * @return the allele index pair * * @deprecated since 2/5/13 */ @Deprecated public static GenotypeLikelihoodsAllelePair getAllelePairUsingDeprecatedOrdering(final int PLindex) { return getAllelePair(PLindexConversion[PLindex]); }
/** * Get the allele ploidy indices for the given PL index * Must use the same ploidy as @see #initializeAnyploidPLIndexToAlleleIndices * * @param PLindex the PL index * @param ploidy number of chromosomes * @return the ploidy allele indices * @throws IllegalStateException if @see #anyploidPloidyToPLIndexToAlleleIndices does not contain the requested ploidy or PL index */ public static synchronized List<Integer> getAlleles(final int PLindex, final int ploidy) { if ( ploidy == 2 ) { // diploid final GenotypeLikelihoodsAllelePair pair = getAllelePair(PLindex); return Arrays.asList(pair.alleleIndex1, pair.alleleIndex2); } else { // non-diploid if (!anyploidPloidyToPLIndexToAlleleIndices.containsKey(ploidy)) throw new IllegalStateException("Must initialize the cache of allele anyploid indices for ploidy " + ploidy); if (PLindex < 0 || PLindex >= anyploidPloidyToPLIndexToAlleleIndices.get(ploidy).size()) { final String msg = "The PL index " + PLindex + " does not exist for " + ploidy + " ploidy, " + (PLindex < 0 ? "cannot have a negative value." : "initialized the cache of allele anyploid indices with the incorrect number of alternate alleles."); throw new IllegalStateException(msg); } return anyploidPloidyToPLIndexToAlleleIndices.get(ploidy).get(PLindex); } }
/** * Get the allele ploidy indices for the given PL index * Must use the same ploidy as @see #initializeAnyploidPLIndexToAlleleIndices * * @param PLindex the PL index * @param ploidy number of chromosomes * @return the ploidy allele indices * @throws IllegalStateException if @see #anyploidPloidyToPLIndexToAlleleIndices does not contain the requested ploidy or PL index */ public static synchronized List<Integer> getAlleles(final int PLindex, final int ploidy) { if ( ploidy == 2 ) { // diploid final GenotypeLikelihoodsAllelePair pair = getAllelePair(PLindex); return Arrays.asList(pair.alleleIndex1, pair.alleleIndex2); } else { // non-diploid if (!anyploidPloidyToPLIndexToAlleleIndices.containsKey(ploidy)) throw new IllegalStateException("Must initialize the cache of allele anyploid indices for ploidy " + ploidy); if (PLindex < 0 || PLindex >= anyploidPloidyToPLIndexToAlleleIndices.get(ploidy).size()) { final String msg = "The PL index " + PLindex + " does not exist for " + ploidy + " ploidy, " + (PLindex < 0 ? "cannot have a negative value." : "initialized the cache of allele anyploid indices with the incorrect number of alternate alleles."); throw new IllegalStateException(msg); } return anyploidPloidyToPLIndexToAlleleIndices.get(ploidy).get(PLindex); } }
@Test(dataProvider = "testGetAllelePairData") public void testGetAllelePair(final int PLindex, final int allele1, final int allele2) { Assert.assertEquals(GenotypeLikelihoods.getAllelePair(PLindex).alleleIndex1, allele1, "allele index " + allele1 + " from PL index " + PLindex + " was not calculated correctly"); Assert.assertEquals(GenotypeLikelihoods.getAllelePair(PLindex).alleleIndex2, allele2, "allele index " + allele2 + " from PL index " + PLindex + " was not calculated correctly"); }
/** * Given a list of original alleles and a subset of new alleles to retain, find the array of old PL indices that correspond * to new PL indices i.e. result[7] = old PL index of genotype containing same alleles as the new genotype with PL index 7. * <p> * This method is written in terms f indices rather than subsetting PLs directly in order to produce output that can be * recycled from sample to sample, provided that the ploidy is the same. * * @param originalAlleles List of original alleles * @param newAlleles New alleles -- must be a subset of {@code originalAlleles} * @return old PL indices of new genotypes */ public static int[] subsettedPLIndices(final List<Allele> originalAlleles, final List<Allele> newAlleles) { final int[] result = new int[GenotypeLikelihoods.numLikelihoods(newAlleles.size(), 2)]; for (int oldPLIndex = 0; oldPLIndex < GenotypeLikelihoods.numLikelihoods(originalAlleles.size(), 2); oldPLIndex++) { final GenotypeLikelihoods.GenotypeLikelihoodsAllelePair allelePairFromPLIndex = GenotypeLikelihoods.getAllelePair(oldPLIndex); final Allele allele1 = originalAlleles.get(allelePairFromPLIndex.alleleIndex1); final Allele allele2 = originalAlleles.get(allelePairFromPLIndex.alleleIndex2); final boolean containsOnlyNewAlleles = newAlleles.contains(allele1) && newAlleles.contains(allele2); if (containsOnlyNewAlleles) { // thus we want this PL...but we need to figure out where its new index is.... final int newPLIndex = GenotypeLikelihoods.calculatePLindex(newAlleles.indexOf(allele1), newAlleles.indexOf(allele2)); result[newPLIndex] = oldPLIndex; } } return result; }
/** * Given a list of original alleles and a subset of new alleles to retain, find the array of old PL indices that correspond * to new PL indices i.e. result[7] = old PL index of genotype containing same alleles as the new genotype with PL index 7. * <p> * This method is written in terms f indices rather than subsetting PLs directly in order to produce output that can be * recycled from sample to sample, provided that the ploidy is the same. * * @param originalAlleles List of original alleles * @param newAlleles New alleles -- must be a subset of {@code originalAlleles} * @return old PL indices of new genotypes */ public static int[] subsettedPLIndices(final List<Allele> originalAlleles, final List<Allele> newAlleles) { final int[] result = new int[GenotypeLikelihoods.numLikelihoods(newAlleles.size(), 2)]; for (int oldPLIndex = 0; oldPLIndex < GenotypeLikelihoods.numLikelihoods(originalAlleles.size(), 2); oldPLIndex++) { final GenotypeLikelihoods.GenotypeLikelihoodsAllelePair allelePairFromPLIndex = GenotypeLikelihoods.getAllelePair(oldPLIndex); final Allele allele1 = originalAlleles.get(allelePairFromPLIndex.alleleIndex1); final Allele allele2 = originalAlleles.get(allelePairFromPLIndex.alleleIndex2); final boolean containsOnlyNewAlleles = newAlleles.contains(allele1) && newAlleles.contains(allele2); if (containsOnlyNewAlleles) { // thus we want this PL...but we need to figure out where its new index is.... final int newPLIndex = GenotypeLikelihoods.calculatePLindex(newAlleles.indexOf(allele1), newAlleles.indexOf(allele2)); result[newPLIndex] = oldPLIndex; } } return result; }
@Override protected void reduceScopeCalculateLikelihoodSums(final VariantContext vc, final int defaultPloidy, final LikelihoodSum[] likelihoodSums) { final ArrayList<double[]> GLs = getGLs(vc.getGenotypes(), true); for ( final double[] likelihoods : GLs ) { final int PLindexOfBestGL = MathUtils.maxElementIndex(likelihoods); if ( PLindexOfBestGL != PL_INDEX_OF_HOM_REF ) { final GenotypeLikelihoods.GenotypeLikelihoodsAllelePair alleles = GenotypeLikelihoods.getAllelePair(PLindexOfBestGL); final int alleleLikelihoodIndex1 = alleles.alleleIndex1 - 1; final int alleleLikelihoodIndex2 = alleles.alleleIndex2 - 1; if ( alleles.alleleIndex1 != 0 ) likelihoodSums[alleleLikelihoodIndex1].sum += likelihoods[PLindexOfBestGL] - likelihoods[PL_INDEX_OF_HOM_REF]; // don't double-count it if ( alleles.alleleIndex2 != 0 && alleles.alleleIndex2 != alleles.alleleIndex1 ) likelihoodSums[alleleLikelihoodIndex2].sum += likelihoods[PLindexOfBestGL] - likelihoods[PL_INDEX_OF_HOM_REF]; } } }
private boolean combineAltAlleleLikelihoods(final Genotype g, final int plMaxIndex, final double[] dest, final double[] hetLikelihoods, final double[] homAltLikelihoods) { final int[] pls = g.getPL(); if (pls == null) return false; int hetNextIndex = 0; int homAltNextIndex = 0; for (int plIndex = 1; plIndex < plMaxIndex; plIndex++) { final GenotypeLikelihoods.GenotypeLikelihoodsAllelePair alleles = GenotypeLikelihoods.getAllelePair(plIndex); if (alleles.alleleIndex1 == 0 || alleles.alleleIndex2 == 0) hetLikelihoods[hetNextIndex++] = pls[plIndex] * PHRED_2_LOG10_COEFF; else homAltLikelihoods[homAltNextIndex++] = pls[plIndex] * PHRED_2_LOG10_COEFF; } dest[0] = pls[0] * PHRED_2_LOG10_COEFF; dest[1] = MathUtils.approximateLog10SumLog10(hetLikelihoods); dest[2] = MathUtils.approximateLog10SumLog10(homAltLikelihoods); return true; } }
protected List<Allele> determineAlternateAlleles(final byte ref, final List<SampleGenotypeData> sampleDataList) { final int baseIndexOfRef = BaseUtils.simpleBaseToBaseIndex(ref); final int PLindexOfRef = DiploidGenotype.createDiploidGenotype(ref, ref).ordinal(); for ( int i = 0; i < 4; i++ ) likelihoodSums[i] = 0.0; // based on the GLs, find the alternate alleles with enough probability for ( SampleGenotypeData sampleData : sampleDataList ) { final double[] likelihoods = sampleData.GL.getLikelihoods(); final int PLindexOfBestGL = MathUtils.maxElementIndex(likelihoods); if ( PLindexOfBestGL != PLindexOfRef ) { GenotypeLikelihoods.GenotypeLikelihoodsAllelePair alleles = GenotypeLikelihoods.getAllelePair(PLindexOfBestGL); if ( alleles.alleleIndex1 != baseIndexOfRef ) likelihoodSums[alleles.alleleIndex1] += likelihoods[PLindexOfBestGL] - likelihoods[PLindexOfRef]; // don't double-count it if ( alleles.alleleIndex2 != baseIndexOfRef && alleles.alleleIndex2 != alleles.alleleIndex1 ) likelihoodSums[alleles.alleleIndex2] += likelihoods[PLindexOfBestGL] - likelihoods[PLindexOfRef]; } } final List<Allele> allelesToUse = new ArrayList<Allele>(3); for ( int i = 0; i < 4; i++ ) { if ( likelihoodSums[i] > 0.0 ) allelesToUse.add(Allele.create(BaseUtils.baseIndexToSimpleBase(i), false)); } return allelesToUse; }
private void pushData(final ExactACset targetSet, final ExactACset dependentSet, final int PLsetIndex, final ArrayList<double[]> genotypeLikelihoods) { final int totalK = targetSet.getACsum(); final double[] targetLog10Likelihoods = targetSet.getLog10Likelihoods(); final double[] dependentLog10Likelihoods = dependentSet.getLog10Likelihoods(); final int[] targetACcounts = targetSet.getACcounts().getCounts(); // skip impossible conformations, namely those for which there aren't enough possible chromosomes (2*j in the if loop below) // to fill up the totalK alternate alleles needed; we do want to ensure that there's at least one sample included (hence the Math.max) final int firstIndex = Math.max(1, (totalK + 1) / 2); // find the 2 alleles that are represented by this PL index final GenotypeLikelihoods.GenotypeLikelihoodsAllelePair alleles = GenotypeLikelihoods.getAllelePair(PLsetIndex); // if neither allele is reference then the coefficient is constant throughout this invocation final Double constantCoefficient = (alleles.alleleIndex1 == 0) ? null : determineCoefficient(alleles, 0, targetACcounts, totalK); for ( int j = firstIndex; j < targetLog10Likelihoods.length; j++ ) { final double[] gl = genotypeLikelihoods.get(j); final double coefficient = (constantCoefficient != null) ? constantCoefficient : determineCoefficient(alleles, j, targetACcounts, totalK); final double conformationValue = coefficient + dependentLog10Likelihoods[j-1] + gl[PLsetIndex]; targetLog10Likelihoods[j] = MathUtils.approximateLog10SumLog10(targetLog10Likelihoods[j], conformationValue); } }
int xxOffset = 0; for ( int index = 0; index < plCount; index++ ) { final GenotypeLikelihoods.GenotypeLikelihoodsAllelePair pair = GenotypeLikelihoods.getAllelePair(index); int i = pair.alleleIndex1; int j = pair.alleleIndex2;