Java Genotype.getPL方法代码示例

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
/**
 * Add information from this Genotype to this band
 * @param g a non-null Genotype with GQ and DP attributes
 */
public void add(final int pos, final Genotype g) {
    if ( g == null ) throw new IllegalArgumentException("g cannot be null");
    if ( ! g.hasGQ() ) throw new IllegalArgumentException("g must have GQ field");
    if ( ! g.hasPL() ) throw new IllegalArgumentException("g must have PL field");
    if ( pos != stop + 1 ) throw new IllegalArgumentException("adding genotype at pos " + pos + " isn't contiguous with previous stop " + stop);
    if ( g.getPloidy() != ploidy)
        throw new IllegalArgumentException("cannot add a genotype with a different ploidy: " + g.getPloidy() + " != " + ploidy);

    if( minPLs == null )
        minPLs = g.getPL();
    else { // otherwise take the min with the provided genotype's PLs
        final int[] PL = g.getPL();
        if (PL.length != minPLs.length)
            throw new IllegalStateException("trying to merge different PL array sizes: " + PL.length + " != " + minPLs.length);
        for (int i = 0; i < PL.length; i++)
            if (minPLs[i] > PL[i])
                minPLs[i] = PL[i];
    }
    stop = pos;
    GQs.add(Math.min(g.getGQ(), 99)); // cap the GQs by the max. of 99 emission
    DPs.add(Math.max(g.getDP(),0));
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
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;
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
private void setGenotypeData(Variation variation, VCFHeader header, Genotype genotype) {
    Map<String, Object> genotypeInfo = new HashMap<>();
    if (genotype.getAD() != null) {
        genotypeInfo.put(header.getFormatHeaderLine("AD").getDescription(), genotype.getAD());
    }
    genotypeInfo.put(header.getFormatHeaderLine("GQ") != null ? header.getFormatHeaderLine("GQ")
            .getDescription() : "GQ", genotype.getGQ());
    genotypeInfo.put(header.getFormatHeaderLine("DP") != null ? header.getFormatHeaderLine("DP")
            .getDescription() : "DP", genotype.getDP());
    if (genotype.getPL() != null) {
        genotypeInfo.put(header.getFormatHeaderLine("PL") != null ? header.getFormatHeaderLine("PL")
                .getDescription() : "PL", genotype.getPL());
    }
    genotypeInfo.put(header.getFormatHeaderLine("GT") != null ? header.getFormatHeaderLine("GT")
            .getDescription() : "GT", genotype.getGenotypeString());

    variation.getGenotypeData().setInfo(genotypeInfo);
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
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] = GvcfMathUtils.approximateLog10SumLog10(hetLikelihoods);
    dest[2] = GvcfMathUtils.approximateLog10SumLog10(homAltLikelihoods);
    return true;
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
/**
 * Returns a new Genotype with the PLs of the multi-allelic original reduced to a bi-allelic case.
 *
 * <p>Uses the log-sum-exp trick in order to work well with very low PLs</p>
 *
 * <p>This is handled in the following way:</p>
 *
 * <p>Suppose we have for a A/B/C site the following GLs:</p>
 *
 * <p>AA AB BB AC BC CC</p>
 *
 * <p>and we want to get the bi-allelic GLs for X/B, where X is everything not B</p>
 *
 * <p>XX = AA + AC + CC (since X = A or C)<br/>
 * XB = AB + BC                           <br/>
 * BB = BB                                     <br/>
 * </p>
 * <p>
 *     This implementation use the log sum trick in order to avoid numeric inestability.
 * </p>
 *
 * @param original the original multi-allelic genotype
 * @param altIndex the index of the alt allele we wish to keep in the bialleic case -- with ref == 0
 * @param nAlts the total number of alt alleles
 * @return a new biallelic genotype with appropriate PLs
 */
protected Genotype combineGLsPrecise(final Genotype original, final int altIndex, final int nAlts ) {

    if ( original.isNonInformative() )
        return new GenotypeBuilder(original).PL(BIALLELIC_NON_INFORMATIVE_PLS).alleles(BIALLELIC_NOCALL).make();

    if ( altIndex < 1 || altIndex > nAlts ) throw new IllegalStateException("altIndex must be between 1 and nAlts " + nAlts);

    final int[] pls = original.getPL();

    final int nAlleles = nAlts + 1;

    final int plCount = pls.length;

    double BB = 0;
    final double[] XBvalues = new double[nAlleles - 1];
    final double[] XXvalues = new double[plCount - nAlleles];

    int xbOffset = 0;
    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;
        if (i == j) {
            if (i == altIndex) BB = PHRED_2_LOG10_COEFF * pls[index]; else XXvalues[xxOffset++] = PHRED_2_LOG10_COEFF * pls[index];
        } else if (i == altIndex || j == altIndex)
            XBvalues[xbOffset++] = PHRED_2_LOG10_COEFF * pls[index];
        else
            XXvalues[xxOffset++] = PHRED_2_LOG10_COEFF * pls[index];
    }

    final double XB = MathUtils.log10sumLog10(XBvalues);
    final double XX = MathUtils.log10sumLog10(XXvalues);

    final double[] GLs = new double[] { XX, XB, BB};
    return new GenotypeBuilder(original).PL(GLs).alleles(BIALLELIC_NOCALL).make();
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
private void parseInfo(final Variation variation, final Allele alt, final Long refId, final Genotype genotype) {
    if (genotype != null) {
        Map<String, Object> genotypeInfo = new HashMap<>();
        if (genotype.getAD() != null) {
            genotypeInfo.put("AD", genotype.getAD());
        }
        genotypeInfo.put("GQ", genotype.getGQ());
        genotypeInfo.put("DP", genotype.getDP());
        if (genotype.getPL() != null) {
            genotypeInfo.put("PL", genotype.getPL());
        }
        genotypeInfo.put("GT", genotype.getGenotypeString());
    }

    if (variation.getType() == VariationType.BND && variation.getBindInfo() == null) {
        variation.setBindInfo(new HashMap<>());
        for (Pattern pattern : BIND_PATTERNS) {
            Matcher matcher = pattern.matcher(alt.getDisplayString());
            if (matcher.matches()) {
                String chrName = matcher.group(1);
                Optional<Chromosome> chromosome = referenceGenomeManager.loadChromosomes(refId)
                        .stream()
                        .filter(c -> c.getName().equals(chrName) ||
                                c.getName().equals(Utils.changeChromosomeName(chrName)))
                        .findAny();

                variation.getBindInfo().put(BIND_CHR_ATTRIBUTE, chromosome.isPresent() ?
                        chromosome.get().getId() : chrName);
                variation.getBindInfo().put(BIND_POS_ATTRIBUTE, matcher.group(2));
                break;
            }
        }
    }
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
/**
 * Returns a new Genotype with the PLs of the multi-allelic original reduced to a bi-allelic case.
 *
 * <p>Uses the log-sum-exp trick in order to work well with very low PLs</p>
 *
 * <p>This is handled in the following way:</p>
 *
 * <p>Suppose we have for a A/B/C site the following GLs:</p>
 *
 * <p>AA AB BB AC BC CC</p>
 *
 * <p>and we want to get the bi-allelic GLs for X/B, where X is everything not B</p>
 *
 * <p>XX = AA + AC + CC (since X = A or C)<br/>
 * XB = AB + BC                           <br/>
 * BB = BB                                     <br/>
 * </p>
 * <p>
 *     This implementation use the log sum trick in order to avoid numeric inestability.
 * </p>
 *
 * @param original the original multi-allelic genotype
 * @param altIndex the index of the alt allele we wish to keep in the bialleic case -- with ref == 0
 * @param nAlts the total number of alt alleles
 * @return a new biallelic genotype with appropriate PLs
 */
@Requires({"original.hasLikelihoods()"})
@Ensures({"result.hasLikelihoods()", "result.getPL().length == 3"})
protected Genotype combineGLsPrecise(final Genotype original, final int altIndex, final int nAlts ) {

    if ( original.isNonInformative() )
        return new GenotypeBuilder(original).PL(BIALLELIC_NON_INFORMATIVE_PLS).alleles(BIALLELIC_NOCALL).make();

    if ( altIndex < 1 || altIndex > nAlts ) throw new IllegalStateException("altIndex must be between 1 and nAlts " + nAlts);

    final int[] pls = original.getPL();

    final int nAlleles = nAlts + 1;

    final int plCount = pls.length;

    double BB = 0;
    final double[] XBvalues = new double[nAlleles - 1];
    final double[] XXvalues = new double[plCount - nAlleles];

    int xbOffset = 0;
    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;
        if (i == j) {
          if (i == altIndex) BB = PHRED_2_LOG10_COEFF * pls[index]; else XXvalues[xxOffset++] = PHRED_2_LOG10_COEFF * pls[index];
        } else if (i == altIndex || j == altIndex)
          XBvalues[xbOffset++] = PHRED_2_LOG10_COEFF * pls[index];
        else
          XXvalues[xxOffset++] = PHRED_2_LOG10_COEFF * pls[index];
    }

    final double XB = GvcfMathUtils.log10sumLog10(XBvalues);
    final double XX = GvcfMathUtils.log10sumLog10(XXvalues);

    final double[] GLs = new double[] { XX, XB, BB};
    return new GenotypeBuilder(original).PL(GLs).alleles(BIALLELIC_NOCALL).make();
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
private static int[] generatePL(final Genotype g, final int[] genotypeIndexMapByPloidy) {
	final int[] PLs = new int[genotypeIndexMapByPloidy.length];
	final int[] oldPLs = g.getPL();
	for (int i = 0; i < PLs.length; i++)
		PLs[i] = oldPLs[genotypeIndexMapByPloidy[i]];
	return PLs;
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
private boolean genotypeCanBeMergedInCurrentBlock(final Genotype g) {
    return currentBlock != null && currentBlock.withinBounds(g.getGQ()) && currentBlock.getPloidy() == g.getPloidy()
            && (currentBlock.getMinPLs() == null || !g.hasPL() || (currentBlock.getMinPLs().length == g.getPL().length));
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
/**
 * Returns a new Genotype with the PLs of the multi-allelic original reduced to a bi-allelic case.
 *
 * <p>Uses the log-sum-exp trick in order to work well with very low PLs</p>
 *
 * <p>This is handled in the following way:</p>
 *
 * <p>Suppose we have for a A/B/C site the following GLs:</p>
 *
 * <p>AA AB BB AC BC CC</p>
 *
 * <p>and we want to get the bi-allelic GLs for X/B, where X is everything not B</p>
 *
 * <p>XX = AA + AC + CC (since X = A or C)<br/>
 * XB = AB + BC                           <br/>
 * BB = BB                                <br/>
 * </p>
 * <p>
 *     This implementation uses the log-sum-exp trick in order to avoid numeric instability (underflow).
 * </p>
 *
 * @param original the original multi-allelic genotype
 * @param alleleIndex the index of the alt allele we wish to keep in the bialleic case -- with ref == 0
 * @param numberOfAlleles the total number of alleles (alternatives + the reference).
 * @return a new biallelic genotype with appropriate PLs
 */
@Requires({"original.hasLikelihoods() && alleleIndex >= 0"})
@Ensures({"result.hasLikelihoods()"})
private Genotype combineGLs(final Genotype original, final int defaultPloidy, final int alleleIndex, final int numberOfAlleles ) {

    final int declaredPloidy = original.getPloidy();
    final int ploidy = declaredPloidy <= 0 ? defaultPloidy : declaredPloidy;
    if ( original.isNonInformative() )
        return new GenotypeBuilder(original).PL(biallelicNonInformativePls(ploidy)).alleles(biallelicNoCall(ploidy)).make();

    final int[] pls = original.getPL();

    final GenotypeLikelihoodCalculator calculator = GenotypeLikelihoodCalculators.getInstance(ploidy, numberOfAlleles);
    final double[] newPLs = new double[ploidy + 1];
    Arrays.fill(newPLs, Double.NEGATIVE_INFINITY);
    for (int i = 0; i < pls.length; i++) {
        final GenotypeAlleleCounts alleleCounts = calculator.genotypeAlleleCountsAt(i);
        final int alleleCount = alleleCounts.alleleCountFor(alleleIndex);
        final int newPLIndex = alleleIndex == 0 ? ploidy - alleleCount : alleleCount;
        newPLs[newPLIndex] = GvcfMathUtils.approximateLog10SumLog10(newPLs[newPLIndex], -.1 * pls[i]);
    }

    return new GenotypeBuilder(original).PL(newPLs).alleles(biallelicNoCall(ploidy)).make();
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
/**
 * Composes a new likelihood array given the original genotype and the genotype index map.
 *
 * @param g the original genotype.
 * @param genotypeIndexMapByPloidy genotype index map. The ith element indicates what genotype in {@code g} corresponds
 *                                 to the ith genotype in the return likelihoods array.
 *
 * @throws NullPointerException if {@code g} or {@code genotypeIndexMapByPloidy} is {@code null}, or if {@code g}
 *    does not contain likelihoods.
 * @throws IndexOutOfBoundsException if {@code genotypeIndexMapByPloidy} contain non valid
 *  genotype indices given the likelihood array in {@code g}.
 *
 * @return never {@code null} but an array of exactly {@code genotypeIndexMapByPloidy.length} positions.
 */
private static int[] generatePL(final Genotype g, final int[] genotypeIndexMapByPloidy) {
    final int[] PLs = new int[genotypeIndexMapByPloidy.length];
    final int[] oldPLs = g.getPL();
    for (int i = 0; i < PLs.length; i++)
        PLs[i] = oldPLs[genotypeIndexMapByPloidy[i]];
    return PLs;
}