Java Genotype.hasLikelihoods方法代码示例

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
/**
 * Unpack GenotypesContext into arraylist of doubel values
 * @param GLs            Input genotype context
 * @return               ArrayList of doubles corresponding to GL vectors
 */
protected static ArrayList<double[]> getGLs(final GenotypesContext GLs, final boolean includeDummy) {
    final ArrayList<double[]> genotypeLikelihoods = new ArrayList<>(GLs.size() + 1);

    if ( includeDummy ) genotypeLikelihoods.add(new double[]{0.0,0.0,0.0}); // dummy
    for ( Genotype sample : GLs.iterateInSampleNameOrder() ) {
        if ( sample.hasLikelihoods() ) {
            final double[] gls = sample.getLikelihoods().getAsVector();

            if ( MathUtils.sum(gls) < GATKVariantContextUtils.SUM_GL_THRESH_NOCALL )
                genotypeLikelihoods.add(gls);
        }
    }

    return genotypeLikelihoods;
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
/**
 * Unpack GenotypesContext into arraylist of double values
 * @param GLs            Input genotype context
 * @param keepUninformative Don't filter out uninformative genotype likelihoods (i.e. all log likelihoods near 0)
 *                          This is useful for VariantContexts with a NON_REF allele
 * @return               ArrayList of doubles corresponding to GL vectors
 */
protected static ArrayList<double[]> getGLs(final GenotypesContext GLs, final boolean includeDummy, final boolean keepUninformative) {
    final ArrayList<double[]> genotypeLikelihoods = new ArrayList<>(GLs.size() + 1);

    if ( includeDummy ) genotypeLikelihoods.add(new double[]{0.0,0.0,0.0}); // dummy
    for ( Genotype sample : GLs.iterateInSampleNameOrder() ) {
        if ( sample.hasLikelihoods() ) {
            final double[] gls = sample.getLikelihoods().getAsVector();
            
            if ( MathUtils.sum(gls) < GaeaGvcfVariantContextUtils.SUM_GL_THRESH_NOCALL || keepUninformative )
                genotypeLikelihoods.add(gls);
        }
    }

    return genotypeLikelihoods;
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
private double[] effectiveAlleleCounts(final VariantContext vc, final double[] log10AlleleFrequencies) {
	final int numAlleles = vc.getNAlleles();
	JointCallingUtils.validateArg(numAlleles == log10AlleleFrequencies.length, "number of alleles inconsistent");
	final double[] log10Result = new double[numAlleles];
	Arrays.fill(log10Result, Double.NEGATIVE_INFINITY);
	for (final Genotype g : vc.getGenotypes()) {
		if (!g.hasLikelihoods()) {
			continue;
		}
		final GenotypeLikelihoodCalculator glCalc = GL_CALCS.getInstance(g.getPloidy(), numAlleles);
		final double[] log10GenotypePosteriors = log10NormalizedGenotypePosteriors(g, glCalc,
				log10AlleleFrequencies);

		new IndexRange(0,
				glCalc.genotypeCount()).forEach(genotypeIndex -> glCalc.genotypeAlleleCountsAt(genotypeIndex)
						.forEachAlleleIndexAndCount((alleleIndex,
								count) -> log10Result[alleleIndex] = GvcfMathUtils.log10SumLog10(
										log10Result[alleleIndex], log10GenotypePosteriors[genotypeIndex]
												+ GvcfMathUtils.Log10Cache.get(count))));
	}
	return GvcfMathUtils.applyToArrayInPlace(log10Result, x -> Math.pow(10.0, x));
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
/**
 * Unpack GenotypesContext into arraylist of doubel values
 * @param GLs            Input genotype context
 * @return               ArrayList of doubles corresponding to GL vectors
 */
protected static ArrayList<double[]> getGLs(final GenotypesContext GLs, final boolean includeDummy) {
    ArrayList<double[]> genotypeLikelihoods = new ArrayList<double[]>(GLs.size() + 1);
    if ( includeDummy ) genotypeLikelihoods.add(new double[]{0.0,0.0,0.0}); // dummy
    for ( Genotype sample : GLs.iterateInSampleNameOrder() ) {
    	if ( sample.hasLikelihoods() ) {
            double[] gls = sample.getLikelihoods().getAsVector();
            if ( MathUtils.sum(gls) < GaeaVariantContextUtils.SUM_GL_THRESH_NOCALL )
                genotypeLikelihoods.add(gls);
        }
    }

    return genotypeLikelihoods;
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
/**
 * From a given genotype, extract a given subset of alleles and update genotype PLs and SACs.
 * @param g                                 genotype to subset
 * @param allelesToUse                      alleles to subset
 * @param vc                                variant context with alleles and genotypes
 * @param defaultPloidy                     ploidy to assume in case that {@code vc} does not contain that information for a sample.
 * @param assignGenotypes                   true: assign hard genotypes, false: leave as no-call
 * @return                                  Genotypes with new PLs and SACs
 */
private Genotype subsetGenotypeAlleles(final Genotype g, final List<Allele> allelesToUse, final VariantContext vc, final int defaultPloidy,
                                       boolean assignGenotypes) {
    final int ploidy = g.getPloidy() <= 0 ? defaultPloidy : g.getPloidy();
    if (!g.hasLikelihoods())
        return GenotypeBuilder.create(g.getSampleName(),GaeaGvcfVariantContextUtils.noCallAlleles(ploidy));
    else {
        // subset likelihood alleles
        final double[] newLikelihoods = subsetLikelihoodAlleles(g, allelesToUse, vc, ploidy);
        if (MathUtils.sum(newLikelihoods) > GaeaGvcfVariantContextUtils.SUM_GL_THRESH_NOCALL)
            return GenotypeBuilder.create(g.getSampleName(), GaeaGvcfVariantContextUtils.noCallAlleles(ploidy));
        else  // just now we would care about newSACs
            return subsetGenotypeAllelesWithLikelihoods(g, allelesToUse, vc, ploidy, assignGenotypes, newLikelihoods);
    }
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
/**
 * Create the new GenotypesContext with the subsetted PLs and ADs
 *
 * @param originalGs               the original GenotypesContext
 * @param vc                       the original VariantContext
 * @param allelesToUse             the actual alleles to use with the new Genotypes
 * @param likelihoodIndexesToUse   the indexes in the PL to use given the allelesToUse (@see #determineLikelihoodIndexesToUse())
 * @param assignGenotypes          assignment strategy for the (subsetted) PLs
 * @return a new non-null GenotypesContext
 */
private static GenotypesContext createGenotypesWithSubsettedLikelihoods(final GenotypesContext originalGs,
                                                                        final VariantContext vc,
                                                                        final List<Allele> allelesToUse,
                                                                        final List<Integer> likelihoodIndexesToUse,
                                                                        final GenotypeAssignmentMethod assignGenotypes) {
    // the new genotypes to create
    final GenotypesContext newGTs = GenotypesContext.create(originalGs.size());

    // make sure we are seeing the expected number of likelihoods per sample
    final int expectedNumLikelihoods = GenotypeLikelihoods.numLikelihoods(vc.getNAlleles(), 2);

    // the samples
    final List<String> sampleIndices = originalGs.getSampleNamesOrderedByName();

    // create the new genotypes
    for ( int k = 0; k < originalGs.size(); k++ ) {
        final Genotype g = originalGs.get(sampleIndices.get(k));
        final GenotypeBuilder gb = new GenotypeBuilder(g);

        // create the new likelihoods array from the alleles we are allowed to use
        double[] newLikelihoods;
        if ( !g.hasLikelihoods() ) {
            // we don't have any likelihoods, so we null out PLs and make G ./.
            newLikelihoods = null;
            gb.noPL();
        } else {
            final double[] originalLikelihoods = g.getLikelihoods().getAsVector();
            if ( likelihoodIndexesToUse == null ) {
                newLikelihoods = originalLikelihoods;
            } else if ( originalLikelihoods.length != expectedNumLikelihoods ) {
                newLikelihoods = null;
            } else {
                newLikelihoods = new double[likelihoodIndexesToUse.size()];
                int newIndex = 0;
                for ( final int oldIndex : likelihoodIndexesToUse )
                    newLikelihoods[newIndex++] = originalLikelihoods[oldIndex];

                // might need to re-normalize
                newLikelihoods = MathUtils.normalizeFromLog10(newLikelihoods, false, true);
            }

            if ( newLikelihoods == null || likelihoodsAreUninformative(newLikelihoods) )
                gb.noPL();
            else
                gb.PL(newLikelihoods);
        }

        updateGenotypeAfterSubsetting(g.getAlleles(), gb, assignGenotypes, newLikelihoods, allelesToUse);
        newGTs.add(gb.make());
    }

    return fixADFromSubsettedAlleles(newGTs, vc, allelesToUse);
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
public static Genotype removePLsAndAD(final Genotype g) {
    return ( g.hasLikelihoods() || g.hasAD() ) ? new GenotypeBuilder(g).noPL().noAD().make() : g;
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
/**
 * From a given variant context, extract a given subset of alleles, and update genotype context accordingly,
 * including updating the PL's, and assign genotypes accordingly
 *
 * @param vc              variant context with alleles and genotype likelihoods
 * @param defaultPloidy   ploidy to assume in case that {@code vc} does not contain that information
 *                        for a sample.
 * @param allelesToUse    alleles to subset
 * @param assignGenotypes true: assign hard genotypes, false: leave as no-call
 * @return GenotypesContext with new PLs
 */
public GenotypesContext subsetAlleles(final VariantContext vc, final int defaultPloidy,
                                      final List<Allele> allelesToUse,
                                      final boolean assignGenotypes) {
    // the genotypes with PLs
    final GenotypesContext oldGTs = vc.getGenotypes();

    // samples
    final List<String> sampleIndices = oldGTs.getSampleNamesOrderedByName();

    // the new genotypes to create
    final GenotypesContext newGTs = GenotypesContext.create();

    // we need to determine which of the alternate alleles (and hence the likelihoods) to use and carry forward
    final int numOriginalAltAlleles = vc.getAlternateAlleles().size();
    final int numNewAltAlleles = allelesToUse.size() - 1;


    // create the new genotypes
    for (int k = 0; k < oldGTs.size(); k++) {
        final Genotype g = oldGTs.get(sampleIndices.get(k));
        final int declaredPloidy = g.getPloidy();
        final int ploidy = declaredPloidy <= 0 ? defaultPloidy : declaredPloidy;
        if (!g.hasLikelihoods()) {
            newGTs.add(GenotypeBuilder.create(g.getSampleName(), GATKVariantContextUtils.noCallAlleles(ploidy)));
            continue;
        }

        // create the new likelihoods array from the alleles we are allowed to use
        final double[] originalLikelihoods = g.getLikelihoods().getAsVector();
        double[] newLikelihoods;

        // Optimization: if # of new alt alleles = 0 (pure ref call), keep original likelihoods so we skip normalization
        // and subsetting
        if (numOriginalAltAlleles == numNewAltAlleles || numNewAltAlleles == 0) {
            newLikelihoods = originalLikelihoods;
        } else {
            newLikelihoods = GeneralPloidyGenotypeLikelihoods.subsetToAlleles(originalLikelihoods, ploidy, vc.getAlleles(), allelesToUse);

            // might need to re-normalize
            newLikelihoods = MathUtils.normalizeFromLog10(newLikelihoods, false, true);
        }

        // if there is no mass on the (new) likelihoods, then just no-call the sample
        if (MathUtils.sum(newLikelihoods) > GATKVariantContextUtils.SUM_GL_THRESH_NOCALL) {
            newGTs.add(GenotypeBuilder.create(g.getSampleName(), GATKVariantContextUtils.noCallAlleles(ploidy)));
        } else {
            final GenotypeBuilder gb = new GenotypeBuilder(g);

            if (numNewAltAlleles == 0)
                gb.noPL();
            else
                gb.PL(newLikelihoods);

            // if we weren't asked to assign a genotype, then just no-call the sample
            if (!assignGenotypes || MathUtils.sum(newLikelihoods) > GATKVariantContextUtils.SUM_GL_THRESH_NOCALL)
                gb.alleles(GATKVariantContextUtils.noCallAlleles(ploidy));
            else
                assignGenotype(gb, newLikelihoods, allelesToUse, ploidy);
            newGTs.add(gb.make());
        }
    }

    return newGTs;

}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
private Map<String, Object> calculateIC(final VariantContext vc) {
    final GenotypesContext genotypes = (founderIds == null || founderIds.isEmpty()) ? vc.getGenotypes() : vc.getGenotypes(founderIds);
    if ( genotypes == null || genotypes.size() < MIN_SAMPLES || !vc.isVariant())
        return null;

    int idxAA = 0, idxAB = 1, idxBB = 2;

    if (!vc.isBiallelic()) {
        // for non-bliallelic case, do test with most common alt allele.
        // Get then corresponding indeces in GL vectors to retrieve GL of AA,AB and BB.
        int[] idxVector = vc.getGLIndecesOfAlternateAllele(vc.getAltAlleleWithHighestAlleleCount());
        idxAA = idxVector[0];
        idxAB = idxVector[1];
        idxBB = idxVector[2];
    }

    double refCount = 0.0;
    double hetCount = 0.0;
    double homCount = 0.0;
    int N = 0; // number of samples that have likelihoods
    for ( final Genotype g : genotypes ) {
        if ( g.isNoCall() || !g.hasLikelihoods() )
            continue;

        if (g.getPloidy() != 2) // only work for diploid samples
            continue;

        N++;
        final double[] normalizedLikelihoods = MathUtils.normalizeFromLog10( g.getLikelihoods().getAsVector() );
        refCount += normalizedLikelihoods[idxAA];
        hetCount += normalizedLikelihoods[idxAB];
        homCount += normalizedLikelihoods[idxBB];
    }

    if( N < MIN_SAMPLES ) {
        return null;
    }

    final double p = ( 2.0 * refCount + hetCount ) / ( 2.0 * (refCount + hetCount + homCount) ); // expected reference allele frequency
    final double q = 1.0 - p; // expected alternative allele frequency
    final double F = 1.0 - ( hetCount / ( 2.0 * p * q * (double)N ) ); // inbreeding coefficient

    Map<String, Object> map = new HashMap<String, Object>();
    map.put(getKeyNames().get(0), String.format("%.4f", F));
    return map;
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
@Override
@Requires("vc != null && allelesToUse != null")
public GenotypesContext subsetAlleles(VariantContext vc, int defaultPloidy, List<Allele> allelesToUse, boolean assignGenotypes) {
    // the genotypes with PLs
    final GenotypesContext oldGTs = vc.getGenotypes();

    // samples
    final List<String> sampleIndices = oldGTs.getSampleNamesOrderedByName();

    // the new genotypes to create
    final GenotypesContext newGTs = GenotypesContext.create();

    // we need to determine which of the alternate alleles (and hence the likelihoods) to use and carry forward
    final int numOriginalAltAlleles = vc.getAlternateAlleles().size();
    final int numNewAltAlleles = allelesToUse.size() - 1;


    // create the new genotypes
    for ( int k = 0; k < oldGTs.size(); k++ ) {
        final Genotype g = oldGTs.get(sampleIndices.get(k));
        final int declaredPloidy = g.getPloidy();
        final int ploidy = declaredPloidy <= 0 ? defaultPloidy : declaredPloidy;
        if ( !g.hasLikelihoods() ) {
            newGTs.add(GenotypeBuilder.create(g.getSampleName(),GaeaGvcfVariantContextUtils.noCallAlleles(ploidy)));
            continue;
        }

        // create the new likelihoods array from the alleles we are allowed to use
        final double[] originalLikelihoods = g.getLikelihoods().getAsVector();
        double[] newLikelihoods;

        // Optimization: if # of new alt alleles = 0 (pure ref call), keep original likelihoods so we skip normalization
        // and subsetting
        if ( numOriginalAltAlleles == numNewAltAlleles || numNewAltAlleles == 0) {
            newLikelihoods = originalLikelihoods;
        } else {
            newLikelihoods = GeneralPloidyGenotypeLikelihoods.subsetToAlleles(originalLikelihoods, ploidy, vc.getAlleles(), allelesToUse);

            // might need to re-normalize
            newLikelihoods = GvcfMathUtils.normalizeFromLog10(newLikelihoods, false, true);
        }

        // if there is no mass on the (new) likelihoods, then just no-call the sample
        if ( GvcfMathUtils.sum(newLikelihoods) > GaeaGvcfVariantContextUtils.SUM_GL_THRESH_NOCALL ) {
            newGTs.add(GenotypeBuilder.create(g.getSampleName(), GaeaGvcfVariantContextUtils.noCallAlleles(ploidy)));
        } else {
            final GenotypeBuilder gb = new GenotypeBuilder(g);
            final String sampleName = g.getSampleName();

            if ( numNewAltAlleles == 0 )
                gb.noPL();
            else
                gb.PL(newLikelihoods);

            // if we weren't asked to assign a genotype, then just no-call the sample
            if ( !assignGenotypes || GvcfMathUtils.sum(newLikelihoods) > GaeaGvcfVariantContextUtils.SUM_GL_THRESH_NOCALL )
                gb.alleles(GaeaGvcfVariantContextUtils.noCallAlleles(ploidy));
            else
                assignGenotype(gb, vc, sampleName, newLikelihoods, allelesToUse, ploidy);
            newGTs.add(gb.make());
        }
    }

    return GaeaGvcfVariantContextUtils.fixADFromSubsettedAlleles(newGTs, vc, allelesToUse);
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
public static Genotype removePLsAndAD(final Genotype g) {
	return (g.hasLikelihoods() || g.hasAD()) ? new GenotypeBuilder(g).noPL().noAD().make() : g;
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
/**
 * Get the genotype counts for A/A, A/B, and B/B where A is the reference and B is any alternate allele
 * @param vc
 * @param genotypes may be subset to just founders if a pedigree file is provided
 * @return may be null, otherwise length-3 double[] representing homRef, het, and homVar counts
 */
public double[] getGenotypeCountsForRefVsAllAlts(final VariantContext vc, final GenotypesContext genotypes) {
    if (genotypes == null || !vc.isVariant())
        return null;

    final boolean doMultiallelicMapping = !vc.isBiallelic();

    int idxAA = 0, idxAB = 1, idxBB = 2;

    double refCount = 0;
    double hetCount = 0;
    double homCount = 0;

    sampleCount = 0;
    for (final Genotype g : genotypes) {
        if (g.isCalled() && g.hasLikelihoods() && g.getPloidy() == 2)  // only work for diploid samples
            sampleCount++;
        else
            continue;

        //Need to round the likelihoods to deal with small numerical deviations due to normalizing
        final double[] normalizedLikelihoodsUnrounded = GvcfMathUtils.normalizeFromLog10(g.getLikelihoods().getAsVector());
        double[] normalizedLikelihoods = new double[normalizedLikelihoodsUnrounded.length];
        if (returnRounded) {
            for (int i = 0; i < normalizedLikelihoodsUnrounded.length; i++) {
                normalizedLikelihoods[i] = Math.round(normalizedLikelihoodsUnrounded[i]);
            }
        } else {
            normalizedLikelihoods = normalizedLikelihoodsUnrounded;
        }

        if (doMultiallelicMapping) {
            if (g.isHetNonRef()) {
                //all likelihoods go to homCount
                homCount++;
                continue;
            }

            if (!g.isHomRef()) {
                //get alternate allele for each sample
                final Allele a1 = g.getAllele(0);
                final Allele a2 = g.getAllele(1);
                final int[] idxVector = vc.getGLIndecesOfAlternateAllele(a2.isNonReference() ? a2 : a1);
                idxAA = idxVector[0];
                idxAB = idxVector[1];
                idxBB = idxVector[2];
            }
        }

        refCount += normalizedLikelihoods[idxAA];
        hetCount += normalizedLikelihoods[idxAB];
        homCount += normalizedLikelihoods[idxBB];
    }
    return new double[]{refCount, hetCount, homCount};
}