Java Genotype.hasAD方法代码示例

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
/**
 * Fix the AD for the given Genotype
 *
 * @param genotype              the original Genotype
 * @param alleleIndexesToUse    a bitset describing whether or not to keep a given index
 * @param nAllelesToUse         how many alleles we are keeping
 * @return a non-null Genotype
 */
private static Genotype fixAD(final Genotype genotype, final boolean[] alleleIndexesToUse, final int nAllelesToUse) {
    // if it ain't broke don't fix it
    if ( !genotype.hasAD() )
        return genotype;

    final GenotypeBuilder builder = new GenotypeBuilder(genotype);

    final int[] oldAD = genotype.getAD();
    if ( oldAD.length != alleleIndexesToUse.length ) {
        builder.noAD();
    } else {
        final int[] newAD = new int[nAllelesToUse];
        int currentIndex = 0;
        for ( int i = 0; i < oldAD.length; i++ ) {
            if ( alleleIndexesToUse[i] )
                newAD[currentIndex++] = oldAD[i];
        }
        builder.AD(newAD);
    }
    return builder.make();
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
/**
 * Fix the AD for the given Genotype
 *
 * @param genotype
 *            the original Genotype
 * @param alleleIndexesToUse
 *            a bitset describing whether or not to keep a given index
 * @param nAllelesToUse
 *            how many alleles we are keeping
 * @return a non-null Genotype
 */
private static Genotype fixAD(final Genotype genotype, final BitSet alleleIndexesToUse) {
	// if it ain't broke don't fix it
	if (!genotype.hasAD())
		return genotype;

	final GenotypeBuilder builder = new GenotypeBuilder(genotype);

	final int[] oldAD = genotype.getAD();
	final int[] newAD = new int[alleleIndexesToUse.cardinality()];

       int currentIndex = 0;
       for ( int i = alleleIndexesToUse.nextSetBit(0); i >= 0; i = alleleIndexesToUse.nextSetBit(i+1) ) {
           newAD[currentIndex++] = oldAD[i];
       }
	
       return builder.AD(newAD).make();
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
/**
 * Merge into the context a new genotype represented by the given VariantContext for the provided list of target alleles.
 * This method assumes that none of the alleles in the VC overlaps with any of the alleles in the set.
 *
 * @param mergedGenotypes   the genotypes context to add to
 * @param VC                the Variant Context for the sample
 * @param remappedAlleles   the list of remapped alleles for the sample
 * @param targetAlleles     the list of target alleles
 */
private static void mergeRefConfidenceGenotypes(final GenotypesContext mergedGenotypes,
                                                final VariantContext VC,
                                                final List<Allele> remappedAlleles,
                                                final List<Allele> targetAlleles) {
    final int maximumPloidy = VC.getMaxPloidy(GATKVariantContextUtils.DEFAULT_PLOIDY);
    // the map is different depending on the ploidy, so in order to keep this method flexible (mixed ploidies)
    // we need to get a map done (lazily inside the loop) for each ploidy, up to the maximum possible.
    final int[][] genotypeIndexMapsByPloidy = new int[maximumPloidy + 1][];
    final int maximumAlleleCount = Math.max(remappedAlleles.size(),targetAlleles.size());
    final int[] indexesOfRelevantAlleles = getIndexesOfRelevantAlleles(remappedAlleles, targetAlleles, VC.getStart());

    for ( final Genotype g : VC.getGenotypes() ) {
        final String name = g.getSampleName();
        if ( mergedGenotypes.containsSample(name) )
            continue;
        final int ploidy = g.getPloidy();
        final GenotypeBuilder genotypeBuilder = new GenotypeBuilder(g).alleles(GATKVariantContextUtils.noCallAlleles(g.getPloidy()));
        if (g.hasPL()) {
            // lazy initialization of the genotype index map by ploidy.
            final int[] genotypeIndexMapByPloidy = genotypeIndexMapsByPloidy[ploidy] == null
                    ? GenotypeLikelihoodCalculators.getInstance(ploidy, maximumAlleleCount).genotypeIndexMap(indexesOfRelevantAlleles)
                    : genotypeIndexMapsByPloidy[ploidy];
            final int[] PLs = generatePL(g, genotypeIndexMapByPloidy);
            final int[] AD = g.hasAD() ? generateAD(g.getAD(), indexesOfRelevantAlleles) : null;
            genotypeBuilder.PL(PLs).AD(AD).noGQ();
        }
        mergedGenotypes.add(genotypeBuilder.make());
    }
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
protected int getDepth(final GenotypesContext genotypes) {
    int standardDepth = 0;
    int ADrestrictedDepth = 0;

    for ( final Genotype genotype : genotypes ) {

        // we care only about variant calls with likelihoods
        if ( !genotype.isHet() && !genotype.isHomVar() )
            continue;

        // if we have the AD values for this sample, let's make sure that the variant depth is greater than 1!
        if ( genotype.hasAD() ) {
            final int[] AD = genotype.getAD();
            final int totalADdepth = (int)GvcfMathUtils.sum(AD);
            if ( totalADdepth - AD[0] > 1 )
                ADrestrictedDepth += totalADdepth;
            standardDepth += totalADdepth;
            continue;
        }
        
        if ( genotype.hasDP() )
        	standardDepth += genotype.getDP();
    }

    // if the AD-restricted depth is a usable value (i.e. not zero), then we should use that one going forward
    if ( ADrestrictedDepth > 0 )
        standardDepth = ADrestrictedDepth;

    return standardDepth;
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
@VisibleForTesting
@NotNull
Optional<SomaticVariant> createVariant(@NotNull final String sample, @NotNull final VariantContext context) {
    if (filter.test(context)) {
        final Genotype genotype = context.getGenotype(sample);
        if (genotype.hasAD() && genotype.getAD().length > 1) {
            final AlleleFrequencyData frequencyData = VariantFactoryFunctions.determineAlleleFrequencies(genotype);
            SomaticVariantImpl.Builder builder = new SomaticVariantImpl.Builder().chromosome(context.getContig())
                    .annotations(Collections.emptyList())
                    .position(context.getStart())
                    .ref(context.getReference().getBaseString())
                    .alt(alt(context))
                    .alleleReadCount(frequencyData.alleleReadCount())
                    .totalReadCount(frequencyData.totalReadCount())
                    .totalReadCount(frequencyData.totalReadCount())
                    .mappability(context.getAttributeAsDouble(MAPPABILITY_TAG, 0));

            attachCallers(builder, context);
            attachAnnotations(builder, context);
            attachFilter(builder, context);
            attachID(builder, context);
            attachType(builder, context);
            return Optional.of(builder.build());
        }
    }
    return Optional.empty();
}
 

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包/类
/**
 * Cleans up genotype-level annotations that need to be updated.
 * 1. move MIN_DP to DP if present
 * 2. propagate DP to AD if not present
 * 3. remove SB if present
 * 4. change the PGT value from "0|1" to "1|1" for homozygous variant genotypes
 *
 * @param VC            the VariantContext with the Genotypes to fix
 * @param createRefGTs  if true we will also create proper hom ref genotypes since we assume the site is monomorphic
 * @return a new set of Genotypes
 */
private List<Genotype> cleanupGenotypeAnnotations(final VariantContext VC, final boolean createRefGTs) {
    final GenotypesContext oldGTs = VC.getGenotypes();
    final List<Genotype> recoveredGs = new ArrayList<>(oldGTs.size());
    for ( final Genotype oldGT : oldGTs ) {
        final Map<String, Object> attrs = new HashMap<>(oldGT.getExtendedAttributes());

        final GenotypeBuilder builder = new GenotypeBuilder(oldGT);
        int depth = oldGT.hasDP() ? oldGT.getDP() : 0;

        // move the MIN_DP to DP
        if ( oldGT.hasExtendedAttribute("MIN_DP") ) {
            depth = Integer.parseInt((String)oldGT.getAnyAttribute("MIN_DP"));
            builder.DP(depth);
            attrs.remove("MIN_DP");
        }

        // remove SB
        attrs.remove("SB");

        // update PGT for hom vars
        if ( oldGT.isHomVar() && oldGT.hasExtendedAttribute(HaplotypeCaller.HAPLOTYPE_CALLER_PHASING_GT_KEY) ) {
            attrs.put(HaplotypeCaller.HAPLOTYPE_CALLER_PHASING_GT_KEY, "1|1");
        }

        // create AD if it's not there
        if ( !oldGT.hasAD() && VC.isVariant() ) {
            final int[] AD = new int[VC.getNAlleles()];
            AD[0] = depth;
            builder.AD(AD);
        }

        if ( createRefGTs ) {
            final int ploidy = oldGT.getPloidy();
            final List<Allele> refAlleles = Collections.nCopies(ploidy,VC.getReference());

            //keep 0 depth samples as no-call
            if (depth > 0) {
                builder.alleles(refAlleles);
            }

            // also, the PLs are technically no longer usable
            builder.noPL();
        }

        recoveredGs.add(builder.noAttributes().attributes(attrs).make());
    }
    return recoveredGs;
}
 

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
public Map<String, Object> annotate(final RefMetaDataTracker tracker,
                                    final AnnotatorCompatible walker,
                                    final ReferenceContext ref,
                                    final Map<String, AlignmentContext> stratifiedContexts,
                                    final VariantContext vc,
                                    final Map<String, PerReadAlleleLikelihoodMap> perReadAlleleLikelihoodMap ) {
    if ( !vc.hasLog10PError() )
        return null;

    final GenotypesContext genotypes = vc.getGenotypes();
    if ( genotypes == null || genotypes.size() == 0 )
        return null;

    int standardDepth = 0;
    int ADrestrictedDepth = 0;

    for ( final Genotype genotype : genotypes ) {

        // we care only about variant calls with likelihoods
        if ( !genotype.isHet() && !genotype.isHomVar() )
            continue;

        // if we have the AD values for this sample, let's make sure that the variant depth is greater than 1!
        // TODO -- If we like how this is working and want to apply it to a situation other than the single sample HC pipeline,
        // TODO --  then we will need to modify the annotateContext() - and related - routines in the VariantAnnotatorEngine
        // TODO --  so that genotype-level annotations are run first (to generate AD on the samples) and then the site-level
        // TODO --  annotations must come afterwards (so that QD can use the AD).
        if ( genotype.hasAD() ) {
            final int[] AD = genotype.getAD();
            final int totalADdepth = (int) MathUtils.sum(AD);
            if ( totalADdepth - AD[0] > 1 )
                ADrestrictedDepth += totalADdepth;
            standardDepth += totalADdepth;
            continue;
        }

        if (stratifiedContexts!= null && !stratifiedContexts.isEmpty()) {
            final AlignmentContext context = stratifiedContexts.get(genotype.getSampleName());
            if ( context == null )
                continue;
            standardDepth += context.getBasePileup().depthOfCoverage();

        } else if (perReadAlleleLikelihoodMap != null) {
            final PerReadAlleleLikelihoodMap perReadAlleleLikelihoods = perReadAlleleLikelihoodMap.get(genotype.getSampleName());
            if (perReadAlleleLikelihoods == null || perReadAlleleLikelihoods.isEmpty())
                continue;

            standardDepth += perReadAlleleLikelihoods.getNumberOfStoredElements();
        } else if ( genotype.hasDP() ) {
            standardDepth += genotype.getDP();
        }
    }

    // if the AD-restricted depth is a usable value (i.e. not zero), then we should use that one going forward
    if ( ADrestrictedDepth > 0 )
        standardDepth = ADrestrictedDepth;

    if ( standardDepth == 0 )
        return null;

    final double altAlleleLength = GATKVariantContextUtils.getMeanAltAlleleLength(vc);
    // Hack: when refContext == null then we know we are coming from the HaplotypeCaller and do not want to do a
    //  full length-based normalization (because the indel length problem is present only in the UnifiedGenotyper)
    double QD = -10.0 * vc.getLog10PError() / ((double)standardDepth * indelNormalizationFactor(altAlleleLength, ref != null));

    // Hack: see note in the fixTooHighQD method below
    QD = fixTooHighQD(QD);

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

import htsjdk.variant.variantcontext.Genotype; //导入方法依赖的package包/类
private List<Genotype> cleanupGenotypeAnnotations(final VariantContext VC, final boolean createRefGTs) {
	final GenotypesContext oldGTs = VC.getGenotypes();
	final List<Genotype> recoveredGs = new ArrayList<>(oldGTs.size());
	for (final Genotype oldGT : oldGTs) {
		final Map<String, Object> attrs = new HashMap<>(oldGT.getExtendedAttributes());

		final GenotypeBuilder builder = new GenotypeBuilder(oldGT);
		int depth = oldGT.hasDP() ? oldGT.getDP() : 0;

		// move the MIN_DP to DP
		if (oldGT.hasExtendedAttribute(GaeaVCFConstants.MIN_DP_FORMAT_KEY)) {
			depth = Integer.parseInt((String) oldGT.getAnyAttribute(GaeaVCFConstants.MIN_DP_FORMAT_KEY));
			builder.DP(depth);
			attrs.remove(GaeaVCFConstants.MIN_DP_FORMAT_KEY);
		}

		// move the GQ to RGQ
		if (createRefGTs && oldGT.hasGQ()) {
			builder.noGQ();
			attrs.put(GaeaVCFConstants.REFERENCE_GENOTYPE_QUALITY, oldGT.getGQ());
		}

		// remove SB
		attrs.remove(GaeaVCFConstants.STRAND_BIAS_BY_SAMPLE_KEY);

		// update PGT for hom vars
		if (oldGT.isHomVar() && oldGT.hasExtendedAttribute(GaeaVCFConstants.HAPLOTYPE_CALLER_PHASING_GT_KEY)) {
			attrs.put(GaeaVCFConstants.HAPLOTYPE_CALLER_PHASING_GT_KEY, "1|1");
		}

		// create AD if it's not there
		if (!oldGT.hasAD() && VC.isVariant()) {
			final int[] AD = new int[VC.getNAlleles()];
			AD[0] = depth;
			builder.AD(AD);
		}

		if (createRefGTs) {
			final int ploidy = oldGT.getPloidy();
			final List<Allele> refAlleles = Collections.nCopies(ploidy, VC.getReference());

			// keep 0 depth samples and 0 GQ samples as no-call
			if (depth > 0 && oldGT.hasGQ() && oldGT.getGQ() > 0) {
				builder.alleles(refAlleles);
			}

			// also, the PLs are technically no longer usable
			builder.noPL();
		}

		recoveredGs.add(builder.noAttributes().attributes(attrs).make());
	}
	return recoveredGs;
}
 

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包/类
/**
 * Merge into the context a new genotype represented by the given
 * VariantContext for the provided list of target alleles. This method
 * assumes that none of the alleles in the VC overlaps with any of the
 * alleles in the set.
 */
private static void mergeRefConfidenceGenotypes(final GenotypesContext mergedGenotypes, final VariantContext vc,
		final List<Allele> remappedAlleles, final List<Allele> targetAlleles, final boolean samplesAreUniquified,
		final boolean shouldComputePLs) {
	final int maximumPloidy = vc.getMaxPloidy(GaeaGvcfVariantContextUtils.DEFAULT_PLOIDY);
	// the map is different depending on the ploidy, so in order to keep
	// this method flexible (mixed ploidies)
	// we need to get a map done (lazily inside the loop) for each ploidy,
	// up to the maximum possible.
	final int[][] genotypeIndexMapsByPloidy = new int[maximumPloidy + 1][];
	final int maximumAlleleCount = Math.max(remappedAlleles.size(), targetAlleles.size());

	for (final Genotype g : vc.getGenotypes()) {
		final String name;
		if (samplesAreUniquified)
			name = g.getSampleName() + "." + vc.getSource();
		else
			name = g.getSampleName();
		final int ploidy = g.getPloidy();
		final GenotypeBuilder genotypeBuilder = new GenotypeBuilder(g)
				.alleles(GaeaGvcfVariantContextUtils.noCallAlleles(g.getPloidy())).noPL();
		genotypeBuilder.name(name);

		final boolean doPLs = shouldComputePLs && g.hasPL();
		final boolean hasAD = g.hasAD();
		final boolean hasSAC = g.hasExtendedAttribute(GaeaVCFConstants.STRAND_COUNT_BY_SAMPLE_KEY);
		if (doPLs || hasSAC || hasAD) {
			final int[] perSampleIndexesOfRelevantAlleles = getIndexesOfRelevantAlleles(remappedAlleles,
					targetAlleles, vc.getStart(), g);
			if (doPLs) {
				// lazy initialization of the genotype index map by ploidy.

				final int[] genotypeIndexMapByPloidy = genotypeIndexMapsByPloidy[ploidy] == null
						? GenotypeLikelihoodCalculators.getInstance(ploidy, maximumAlleleCount).genotypeIndexMap(
								perSampleIndexesOfRelevantAlleles)
						: genotypeIndexMapsByPloidy[ploidy];
				final int[] PLs = generatePL(g, genotypeIndexMapByPloidy);
				genotypeBuilder.PL(PLs);
			}
			if (hasAD) {
				genotypeBuilder.AD(generateAD(g.getAD(), perSampleIndexesOfRelevantAlleles));
			}
			if (hasSAC) {
				final List<Integer> sacIndexesToUse = adaptToSACIndexes(perSampleIndexesOfRelevantAlleles);
				final int[] SACs = GaeaGvcfVariantContextUtils.makeNewSACs(g, sacIndexesToUse);
				genotypeBuilder.attribute(GaeaVCFConstants.STRAND_COUNT_BY_SAMPLE_KEY, SACs);
			}
		}
		mergedGenotypes.add(genotypeBuilder.make());
	}
}