Java Genotype.getAD方法代码示例

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包/类
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包/类
/**
 * 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包/类
/**
 * Provide a centralized method of getting the number of reads per allele, 
 * depending on the input given.  Will use the following (in order of preference):
 * - genotype.getAD()
 * - reads from an AlignmentContext
 * - reads from a PerReadAlleleLikelihoodMap (Not yet implemented) 
 *
 *
 * @param genotype The genotype of interest
 * @param stratifiedContexts A mapping 
 * @param vc
 * @return
 */
private int[] getCounts(final Genotype genotype,
                        final Map<String, AlignmentContext> stratifiedContexts,
                        final VariantContext vc){

    // Can't do anything without a genotype here
    if(genotype == null)
        return null;

    int[] retVal = genotype.getAD();
    AlignmentContext context;

    if ( retVal == null && stratifiedContexts != null &&
            (context = stratifiedContexts.get(genotype.getSampleName())) != null){
        // If we get to this point, the getAD() function returned no information
        // about AlleleDepth by Sample - perhaps it wasn't annotated?
        // In that case, let's try to build it up using the algorithm that
        // was here in v 3.1-1 and earlier
        // Also, b/c of the assignment check in the if statement above,
        // we know we have a valid AlignmentContext for this sample!

        final ReadBackedPileup pileup = context.getBasePileup();
        final String bases = new String(pileup.getBases());
        List<Allele> alleles = vc.getAlleles();
        final int n_allele = alleles.size();
        retVal = new int[n_allele];

        // Calculate the depth for each allele, under the assumption that
        // the allele is a single base
        int i=0;
        for(Allele a : alleles){
            retVal[i] = MathUtils.countOccurrences(a.toString().charAt(0), bases);
            i++;
        }

    }

    return retVal;

}
 

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包/类
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包/类
@NotNull
public static AlleleFrequencyData determineAlleleFrequencies(@NotNull final Genotype genotype) {
    Preconditions.checkArgument(genotype.hasAD());

    int[] adFields = genotype.getAD();
    int totalReadCount = 0;
    final int alleleReadCount = adFields[1];
    for (final int afField : adFields) {
        totalReadCount += afField;
    }

    return new AlleleFrequencyData(alleleReadCount, totalReadCount);
}
 

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;
}