/** * 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; }
/** * 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; }
private synchronized GenotypeLikelihoods[] initializeIndelPLCache(final int ploidy) { if (indelPLCache.length <= ploidy) indelPLCache = Arrays.copyOf(indelPLCache,ploidy << 1); if (indelPLCache[ploidy] != null) return indelPLCache[ploidy]; final double denominator = - MathUtils.Log10Cache.get(ploidy); final GenotypeLikelihoods[] result = new GenotypeLikelihoods[MAX_N_INDEL_INFORMATIVE_READS + 1]; result[0] = GenotypeLikelihoods.fromLog10Likelihoods(new double[ploidy + 1]); for( int nInformativeReads = 1; nInformativeReads <= MAX_N_INDEL_INFORMATIVE_READS; nInformativeReads++ ) { double[] PLs = new double[ploidy + 1]; PLs[0] = nInformativeReads * NO_INDEL_LIKELIHOOD; for (int altCount = 1; altCount <= ploidy; altCount++) { final double refLikelihoodAccum = NO_INDEL_LIKELIHOOD + MathUtils.Log10Cache.get(ploidy - altCount); final double altLikelihoodAccum = INDEL_LIKELIHOOD + MathUtils.Log10Cache.get(altCount); PLs[altCount] = nInformativeReads * (MathUtils.approximateLog10SumLog10(refLikelihoodAccum ,altLikelihoodAccum) + denominator); } result[nInformativeReads] = GenotypeLikelihoods.fromLog10Likelihoods(PLs); } indelPLCache[ploidy] = result; return result; }
public void testGetQualFromLikelihoodsMultiAllelicBroken() { GenotypeLikelihoods gl = GenotypeLikelihoods.fromLog10Likelihoods(triAllelic); double actualGQ = gl.getLog10GQ(GenotypeType.HET); double expectedGQ = 1.6; Assert.assertEquals(actualGQ,expectedGQ); }
/** * This genotype has this PL value, converted from double[]. SLOW * @return */ @Requires("PL == null || PL.length > 0") @Ensures({"this.PL == PL"}) public GenotypeBuilder PL(final double[] GLs) { this.PL = GenotypeLikelihoods.fromLog10Likelihoods(GLs).getAsPLs(); return this; }
@Test public void testFromVector2() { GenotypeLikelihoods gl = GenotypeLikelihoods.fromLog10Likelihoods(v); assertDoubleArraysAreEqual(gl.getAsVector(), v); Assert.assertEquals(gl.getAsString(), vPLString); }
/** * Calculate the likelihoods given the list of alleles and the likelihood map. * * <p>This operation is <b>thread-unsafe</b>.</p> * * @param likelihoods the likelihood matrix all alleles vs all reads. * * @throws IllegalArgumentException if {@code alleleList} is {@code null} or {@code likelihoods} is {@code null} * or the alleleList size does not match the allele-count of this calculator, or there are missing allele vs * read combinations in {@code likelihoods}. * * @return never {@code null}. */ public <A extends Allele> GenotypeLikelihoods genotypeLikelihoods(final ReadLikelihoods.Matrix<A> likelihoods) { if (likelihoods == null) throw new IllegalArgumentException("the likelihood map cannot be null"); if (likelihoods.alleleCount() != alleleCount) throw new IllegalArgumentException("mismatch between allele list and alleleCount"); final int readCount = likelihoods.readCount(); ensureReadCapacity(readCount); /// [x][y][z] = z * LnLk(Read_x | Allele_y) final double[] readLikelihoodComponentsByAlleleCount = readLikelihoodComponentsByAlleleCount(likelihoods); final double[][] genotypeLikelihoodByRead = genotypeLikelihoodByRead(readLikelihoodComponentsByAlleleCount,readCount); final double[] readLikelihoodsByGenotypeIndex = genotypeLikelihoods(genotypeLikelihoodByRead, readCount); return GenotypeLikelihoods.fromLog10Likelihoods(readLikelihoodsByGenotypeIndex); }
private Genotype getUpdatedGenotype(final VariantContext vc, final Genotype genotype, final double jointLikelihood, final double jointPosteriorProb, final double[] log10Posteriors){ //Don't update null, missing or unavailable genotypes if(genotype == null || !hasCalledGT(genotype.getType())) return genotype; int phredScaledJL = -1; int phredScaledJP = -1; if(jointLikelihood != NO_JOINT_VALUE){ double dphredScaledJL = QualityUtils.phredScaleLog10ErrorRate(Math.log10(1-jointLikelihood)); phredScaledJL = dphredScaledJL < Byte.MAX_VALUE ? (byte)dphredScaledJL : Byte.MAX_VALUE; } if(jointPosteriorProb != NO_JOINT_VALUE){ double dphredScaledJP = QualityUtils.phredScaleLog10ErrorRate(Math.log10(1-jointPosteriorProb)); phredScaledJP = dphredScaledJP < Byte.MAX_VALUE ? (byte)dphredScaledJP : Byte.MAX_VALUE; } //Add the joint trio calculations final Map<String, Object> genotypeAttributes = new HashMap<>(); genotypeAttributes.putAll(genotype.getExtendedAttributes()); genotypeAttributes.put(GATKVCFConstants.JOINT_LIKELIHOOD_TAG_NAME, phredScaledJL); genotypeAttributes.put(GATKVCFConstants.JOINT_POSTERIOR_TAG_NAME, phredScaledJP); final GenotypeBuilder builder = new GenotypeBuilder(genotype); //final double[] log10Posteriors = MathUtils.toLog10(normalizedPosteriors); //update genotype types based on posteriors GATKVariantContextUtils.updateGenotypeAfterSubsetting(vc.getAlleles(), genotype.getPloidy(), builder, GATKVariantContextUtils.GenotypeAssignmentMethod.USE_PLS_TO_ASSIGN, log10Posteriors, vc.getAlleles()); builder.attribute(GATKVCFConstants.PHRED_SCALED_POSTERIORS_KEY, Utils.listFromPrimitives(GenotypeLikelihoods.fromLog10Likelihoods(log10Posteriors).getAsPLs())); builder.attributes(genotypeAttributes); return builder.make(); }
@Test public void testGetAsMap(){ GenotypeLikelihoods gl = GenotypeLikelihoods.fromLog10Likelihoods(v); //Log scale EnumMap<GenotypeType,Double> glMap = gl.getAsMap(false); Assert.assertEquals(v[GenotypeType.HOM_REF.ordinal()-1],glMap.get(GenotypeType.HOM_REF)); Assert.assertEquals(v[GenotypeType.HET.ordinal()-1],glMap.get(GenotypeType.HET)); Assert.assertEquals(v[GenotypeType.HOM_VAR.ordinal()-1],glMap.get(GenotypeType.HOM_VAR)); //Linear scale glMap = gl.getAsMap(true); double [] vl = GeneralUtils.normalizeFromLog10(v); Assert.assertEquals(vl[GenotypeType.HOM_REF.ordinal()-1],glMap.get(GenotypeType.HOM_REF)); Assert.assertEquals(vl[GenotypeType.HET.ordinal()-1],glMap.get(GenotypeType.HET)); Assert.assertEquals(vl[GenotypeType.HOM_VAR.ordinal()-1],glMap.get(GenotypeType.HOM_VAR)); //Test missing likelihoods gl = GenotypeLikelihoods.fromPLField("."); glMap = gl.getAsMap(false); Assert.assertNull(glMap); }
public void testGetQualFromLikelihoodsMultiAllelic() { GenotypeLikelihoods gl = GenotypeLikelihoods.fromLog10Likelihoods(triAllelic); Allele ref = Allele.create((byte)'A',true); Allele alt1 = Allele.create((byte)'C'); Allele alt2 = Allele.create((byte)'T'); List<Allele> allAlleles = Arrays.asList(ref,alt1,alt2); List<Allele> gtAlleles = Arrays.asList(alt1,alt2); GenotypeBuilder gtBuilder = new GenotypeBuilder(); gtBuilder.alleles(gtAlleles); double actualGQ = gl.getLog10GQ(gtBuilder.make(),allAlleles); double expectedGQ = 1.6; Assert.assertEquals(actualGQ,expectedGQ); }
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());
gb.alleles(noCall); if (UAC.annotateAllSitesWithPLs) gb.attribute(GATKVCFConstants.PL_FOR_ALL_SNP_ALLELES_KEY,GenotypeLikelihoods.fromLog10Likelihoods(MathUtils.normalizeFromLog10(allLikelihoods, false, true))); genotypes.add(gb.make());
final GenotypeLikelihoods snpGLs = GenotypeLikelihoods.fromLog10Likelihoods(homRefCalc.genotypeLikelihoods); final int nIndelInformativeReads = calcNIndelInformativeReads(pileup, refOffset, ref, indelInformativeDepthIndelSize); final GenotypeLikelihoods indelGLs = getIndelPLs(ploidy,nIndelInformativeReads);