private double getLog10GQ(List<Allele> genotypeAlleles,List<Allele> contextAlleles) { int allele1Index = contextAlleles.indexOf(genotypeAlleles.get(0)); int allele2Index = contextAlleles.indexOf(genotypeAlleles.get(1)); int plIndex = calculatePLindex(allele1Index,allele2Index); return getGQLog10FromLikelihoods(plIndex,getAsVector()); }
/** * 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; }
/** * This genotype has this PL value, converted from double[]. SLOW * @return */ public GenotypeBuilder PL(final double[] GLs) { this.PL = GenotypeLikelihoods.fromLog10Likelihoods(GLs).getAsPLs(); return this; }
public int[] getAsPLs() { final double[] GLs = getAsVector(); return GLs == null ? null : GLsToPLs(GLs); }
validateTrue(g.getPloidy() == 2, "only implemented for ploidy 2 for now."); final int expectedNumLikelihoods = GenotypeLikelihoods.numLikelihoods(allelesToKeep.size(), 2); newLog10GQ = GenotypeLikelihoods.getGQLog10FromLikelihoods(indexOfMostLikely, GenotypeLikelihoods.fromPLs(newPLs).getAsVector()); } else { newPLs = null; } else { gb = new GenotypeBuilder(g).PL(newPLs).log10PError(newLog10GQ); final List<Integer> originalDiploidAlleles = GenotypeLikelihoods.getAlleles(MathUtil.indexOfMin(newPLs), 2); gb.alleles(originalDiploidAlleles.stream().map(allelesToKeep::get).collect(Collectors.toList()));
/** * This is really dangerous and returns completely wrong results for genotypes from a multi-allelic context. * Use getLog10GQ(Genotype,VariantContext) or getLog10GQ(Genotype,List<Allele>) in place of it. * * If you **know** you're biallelic, use getGQLog10FromLikelihoods directly. * @param genotype - actually a genotype type (no call, hom ref, het, hom var) * @return an unsafe quantity that could be negative. In the bi-allelic case, the GQ resulting from best minus next best (if the type is the best). */ @Deprecated public double getLog10GQ(GenotypeType genotype){ return getGQLog10FromLikelihoods(genotype.ordinal() - 1 /* NO_CALL IS FIRST */, getAsVector()); }
} else { final int expectedNumLikelihoods = GenotypeLikelihoods.numLikelihoods(originalVC.getNAlleles(), g.getPloidy()); final double[] originalLikelihoods = g.getLikelihoods().getAsVector(); if ( likelihoodIndexesToUse == null ) { newLikelihoods = originalLikelihoods; gb.PL(newLikelihoods); final int PLindex = MathUtils.maxElementIndex(newLikelihoods); gb.log10PError(GenotypeLikelihoods.getGQLog10FromLikelihoods(PLindex, newLikelihoods));
@Test public void testGetLog10GQ(){ GenotypeLikelihoods gl = GenotypeLikelihoods.fromPLField(vPLString); //GQ for the best guess genotype Assert.assertEquals(gl.getLog10GQ(GenotypeType.HET),-3.9); double[] test = GeneralUtils.normalizeFromLog10(gl.getAsVector()); //GQ for the other genotypes Assert.assertEquals(gl.getLog10GQ(GenotypeType.HOM_REF), Math.log10(1.0 - test[GenotypeType.HOM_REF.ordinal()-1])); Assert.assertEquals(gl.getLog10GQ(GenotypeType.HOM_VAR), Math.log10(1.0 - test[GenotypeType.HOM_VAR.ordinal()-1])); //Test missing likelihoods gl = GenotypeLikelihoods.fromPLField("."); Assert.assertEquals(gl.getLog10GQ(GenotypeType.HOM_REF),Double.NEGATIVE_INFINITY); Assert.assertEquals(gl.getLog10GQ(GenotypeType.HET),Double.NEGATIVE_INFINITY); Assert.assertEquals(gl.getLog10GQ(GenotypeType.HOM_VAR),Double.NEGATIVE_INFINITY); }
public EnumMap<GenotypeType,Double> getAsMap(boolean normalizeFromLog10){ //Make sure that the log10likelihoods are set double[] likelihoods = normalizeFromLog10 ? GeneralUtils.normalizeFromLog10(getAsVector()) : getAsVector(); if(likelihoods == null) return null; EnumMap<GenotypeType,Double> likelihoodsMap = new EnumMap<GenotypeType, Double>(GenotypeType.class); likelihoodsMap.put(GenotypeType.HOM_REF,likelihoods[GenotypeType.HOM_REF.ordinal()-1]); likelihoodsMap.put(GenotypeType.HET,likelihoods[GenotypeType.HET.ordinal()-1]); likelihoodsMap.put(GenotypeType.HOM_VAR, likelihoods[GenotypeType.HOM_VAR.ordinal() - 1]); return likelihoodsMap; }
@Test public void testFromString2() { GenotypeLikelihoods gl = GenotypeLikelihoods.fromGLField(vGLString); assertDoubleArraysAreEqual(gl.getAsVector(), v); Assert.assertEquals(gl.getAsString(), vPLString); }
@Test public void testFromString1() { GenotypeLikelihoods gl = GenotypeLikelihoods.fromPLField(vPLString); assertDoubleArraysAreEqual(gl.getAsVector(), new double[]{-9.3, 0, -3.9}); Assert.assertEquals(gl.getAsString(), vPLString); }
private int numPLs() { return GenotypeLikelihoods.numLikelihoods(numAltAlleles+1, 2); }
private static GenotypeLikelihoodsAllelePair[] calculatePLcache(final int altAlleles) { final int numLikelihoods = numLikelihoods(1 + altAlleles, 2); final GenotypeLikelihoodsAllelePair[] cache = new GenotypeLikelihoodsAllelePair[numLikelihoods]; // for all possible combinations of 2 alleles for ( int allele1 = 0; allele1 <= altAlleles; allele1++ ) { for ( int allele2 = allele1; allele2 <= altAlleles; allele2++ ) { cache[calculatePLindex(allele1, allele2)] = new GenotypeLikelihoodsAllelePair(allele1, allele2); } } // a bit of sanity checking for ( int i = 0; i < cache.length; i++ ) { if ( cache[i] == null ) throw new IllegalStateException("BUG: cache entry " + i + " is unexpected null"); } return cache; }
@Test public void testFromVector2() { GenotypeLikelihoods gl = GenotypeLikelihoods.fromLog10Likelihoods(v); assertDoubleArraysAreEqual(gl.getAsVector(), v); Assert.assertEquals(gl.getAsString(), vPLString); }
gb.PL(decodeInts(genotypeValues.get(i))); } else if (gtKey.equals(VCFConstants.GENOTYPE_LIKELIHOODS_KEY)) { gb.PL(GenotypeLikelihoods.fromGLField(genotypeValues.get(i)).getAsPLs()); GenotypeLikelihoods glField = GenotypeLikelihoods.fromGLField(genotypeValues.get(i)); gb.PL(glField.getAsPLs()); gb.attribute(gtKey, genotypeValues.get(i)); } else if (gtKey.equals(VCFConstants.DEPTH_KEY)) {
for ( final Genotype genotype : vc1.getGenotypes() ) { if (!genotype.hasExtendedAttribute(GATKVCFConstants.PHRED_SCALED_POSTERIORS_KEY)){ likelihoods.add(genotype.hasLikelihoods() ? genotype.getLikelihoods().getAsVector() : null ); final String PPstring = (String)PPfromVCF; if (PPstring.charAt(0)=='.') //samples not in trios will have PP tag like ".,.,." if family priors are applied likelihoods.add(genotype.hasLikelihoods() ? genotype.getLikelihoods().getAsVector() : null ); else { final String[] likelihoodsAsStringVector = PPstring.split(","); GATKVariantContextUtils.GenotypeAssignmentMethod.USE_PLS_TO_ASSIGN, posteriors.get(genoIdx), vc1.getAlleles()); builder.attribute(GATKVCFConstants.PHRED_SCALED_POSTERIORS_KEY, Utils.listFromPrimitives(GenotypeLikelihoods.fromLog10Likelihoods(posteriors.get(genoIdx)).getAsPLs())); GenotypeLikelihoods.fromLog10Likelihoods(getDirichletPrior(alleleCounts, vc1.getMaxPloidy(2),useFlatPriors)).getAsPLs());
@Test (expectedExceptions = TribbleException.class) public void testErrorBadFormat() { GenotypeLikelihoods gl = GenotypeLikelihoods.fromPLField("adf,b,c"); gl.getAsVector(); }
final double[] normalizedLikelihoodsUnrounded = MathUtils.normalizeFromLog10(g.getLikelihoods().getAsVector()); double[] normalizedLikelihoods = new double[normalizedLikelihoodsUnrounded.length]; if (returnRounded) { alleleCounts.put(a, currentAlleleCounts + 2*normalizedLikelihoods[GenotypeLikelihoods.calculatePLindex(altIndex,altIndex)]); continue; idxAB = GenotypeLikelihoods.calculatePLindex(Math.min(i,altIndex),Math.max(i,altIndex)); final double aHetCounts = hetCounts.get(a); hetCounts.put(a, aHetCounts + normalizedLikelihoods[idxAB]);
/** * 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]); }
GenotypeLikelihoods glField = GenotypeLikelihoods.fromGLField((String) value); double[] glv = glField.getAsVector(); lhList = new ArrayList<>(glv.length); for (double d : glv) {