Search in sources :

Example 1 with GenotypesContext

use of htsjdk.variant.variantcontext.GenotypesContext in project gatk-protected by broadinstitute.

the class CalculateGenotypePosteriors method apply.

@Override
public void apply(final VariantContext variant, final ReadsContext readsContext, final ReferenceContext referenceContext, final FeatureContext featureContext) {
    final Collection<VariantContext> vcs = featureContext.getValues(getDrivingVariantsFeatureInput());
    final Collection<VariantContext> otherVCs = featureContext.getValues(supportVariants);
    final int missing = supportVariants.size() - otherVCs.size();
    for (final VariantContext vc : vcs) {
        final VariantContext vc_familyPriors;
        final VariantContext vc_bothPriors;
        //do family priors first (if applicable)
        final VariantContextBuilder builder = new VariantContextBuilder(vc);
        //only compute family priors for biallelelic sites
        if (!skipFamilyPriors && vc.isBiallelic()) {
            final GenotypesContext gc = famUtils.calculatePosteriorGLs(vc);
            builder.genotypes(gc);
        }
        VariantContextUtils.calculateChromosomeCounts(builder, false);
        vc_familyPriors = builder.make();
        if (!skipPopulationPriors) {
            vc_bothPriors = PosteriorProbabilitiesUtils.calculatePosteriorProbs(vc_familyPriors, otherVCs, missing * numRefIfMissing, globalPrior, !ignoreInputSamples, defaultToAC, useACoff);
        } else {
            final VariantContextBuilder builder2 = new VariantContextBuilder(vc_familyPriors);
            VariantContextUtils.calculateChromosomeCounts(builder, false);
            vc_bothPriors = builder2.make();
        }
        vcfWriter.add(vc_bothPriors);
    }
}
Also used : GenotypesContext(htsjdk.variant.variantcontext.GenotypesContext) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder) VariantContext(htsjdk.variant.variantcontext.VariantContext)

Example 2 with GenotypesContext

use of htsjdk.variant.variantcontext.GenotypesContext in project gatk by broadinstitute.

the class QualByDepth method annotate.

@Override
public Map<String, Object> annotate(final ReferenceContext ref, final VariantContext vc, final ReadLikelihoods<Allele> likelihoods) {
    Utils.nonNull(vc);
    if (!vc.hasLog10PError()) {
        return Collections.emptyMap();
    }
    final GenotypesContext genotypes = vc.getGenotypes();
    if (genotypes == null || genotypes.isEmpty()) {
        return Collections.emptyMap();
    }
    int depth = 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) MathUtils.sum(AD);
            if (totalADdepth != 0) {
                if (totalADdepth - AD[0] > 1) {
                    ADrestrictedDepth += totalADdepth;
                }
                depth += totalADdepth;
            }
        } else if (likelihoods != null) {
            depth += likelihoods.sampleReadCount(likelihoods.indexOfSample(genotype.getSampleName()));
        } else if (genotype.hasDP()) {
            depth += 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) {
        depth = ADrestrictedDepth;
    }
    if (depth == 0) {
        return Collections.emptyMap();
    }
    final double qual = -10.0 * vc.getLog10PError();
    double QD = qual / depth;
    // Hack: see note in the fixTooHighQD method below
    QD = fixTooHighQD(QD);
    return Collections.singletonMap(getKeyNames().get(0), String.format("%.2f", QD));
}
Also used : GenotypesContext(htsjdk.variant.variantcontext.GenotypesContext) Genotype(htsjdk.variant.variantcontext.Genotype)

Example 3 with GenotypesContext

use of htsjdk.variant.variantcontext.GenotypesContext in project gatk by broadinstitute.

the class AFCalculatorProvider method getInstance.

/**
     * Returns a AF calculator capable to handle a particular variant-context.
     * @param variantContext the target context build.
     * @param defaultPloidy the assumed ploidy in case that there is no a GT call present to determine it.
     * @return never {@code null}
     */
public AFCalculator getInstance(final VariantContext variantContext, final int defaultPloidy, final int maximumAltAlleles) {
    Utils.nonNull(variantContext, "variant context cannot be null");
    final int sampleCount = variantContext.getNSamples();
    if (sampleCount == 0) {
        return getInstance(defaultPloidy, maximumAltAlleles);
    }
    final GenotypesContext genotypes = variantContext.getGenotypes();
    final Genotype firstGenotype = genotypes.get(0);
    int ploidy = firstGenotype.getPloidy();
    if (ploidy <= 0) {
        ploidy = defaultPloidy;
    }
    for (int i = 1; i < sampleCount; i++) {
        final Genotype genotype = genotypes.get(i);
        final int declaredPloidy = genotype.getPloidy();
        final int actualPloidy = declaredPloidy <= 0 ? defaultPloidy : declaredPloidy;
        if (actualPloidy != ploidy) {
            ploidy = AFCalculatorImplementation.UNBOUND_PLOIDY;
            break;
        }
    }
    return getInstance(ploidy, Math.min(variantContext.getNAlleles() - 1, maximumAltAlleles));
}
Also used : GenotypesContext(htsjdk.variant.variantcontext.GenotypesContext) Genotype(htsjdk.variant.variantcontext.Genotype)

Example 4 with GenotypesContext

use of htsjdk.variant.variantcontext.GenotypesContext in project gatk by broadinstitute.

the class InbreedingCoeff method annotate.

@Override
public Map<String, Object> annotate(final ReferenceContext ref, final VariantContext vc, final ReadLikelihoods<Allele> likelihoods) {
    Utils.nonNull(vc);
    final GenotypesContext genotypes = (founderIds == null || founderIds.isEmpty()) ? vc.getGenotypes() : vc.getGenotypes(founderIds);
    if (genotypes == null || genotypes.size() < MIN_SAMPLES || !vc.isVariant()) {
        return Collections.emptyMap();
    }
    final Pair<Integer, Double> sampleCountCoeff = calculateIC(vc, genotypes);
    final int sampleCount = sampleCountCoeff.getLeft();
    final double F = sampleCountCoeff.getRight();
    if (sampleCount < MIN_SAMPLES) {
        logger.warn("Annotation will not be calculated, must provide at least " + MIN_SAMPLES + " samples");
        return Collections.emptyMap();
    }
    return Collections.singletonMap(getKeyNames().get(0), String.format("%.4f", F));
}
Also used : GenotypesContext(htsjdk.variant.variantcontext.GenotypesContext)

Example 5 with GenotypesContext

use of htsjdk.variant.variantcontext.GenotypesContext in project gatk by broadinstitute.

the class MappingQualityZeroUnitTest method makeVC.

private VariantContext makeVC() {
    final GenotypesContext testGC = GenotypesContext.create(2);
    final Allele refAllele = Allele.create("A", true);
    final Allele altAllele = Allele.create("T");
    return (new VariantContextBuilder()).alleles(Arrays.asList(refAllele, altAllele)).chr("1").start(15L).stop(15L).genotypes(testGC).make();
}
Also used : Allele(htsjdk.variant.variantcontext.Allele) GenotypesContext(htsjdk.variant.variantcontext.GenotypesContext) VariantContextBuilder(htsjdk.variant.variantcontext.VariantContextBuilder)

Aggregations

GenotypesContext (htsjdk.variant.variantcontext.GenotypesContext)14 Allele (htsjdk.variant.variantcontext.Allele)7 VariantContextBuilder (htsjdk.variant.variantcontext.VariantContextBuilder)7 Genotype (htsjdk.variant.variantcontext.Genotype)6 VariantContext (htsjdk.variant.variantcontext.VariantContext)5 VcfIterator (com.github.lindenb.jvarkit.util.vcf.VcfIterator)2 GenotypeBuilder (htsjdk.variant.variantcontext.GenotypeBuilder)2 VariantContextWriter (htsjdk.variant.variantcontext.writer.VariantContextWriter)2 VCFHeader (htsjdk.variant.vcf.VCFHeader)2 VCFHeaderLine (htsjdk.variant.vcf.VCFHeaderLine)2 File (java.io.File)2 Collectors (java.util.stream.Collectors)2 Parameter (com.beust.jcommander.Parameter)1 Launcher (com.github.lindenb.jvarkit.util.jcommander.Launcher)1 Program (com.github.lindenb.jvarkit.util.jcommander.Program)1 Logger (com.github.lindenb.jvarkit.util.log.Logger)1 GenomicSequence (com.github.lindenb.jvarkit.util.picard.GenomicSequence)1 SAMSequenceDictionaryProgress (com.github.lindenb.jvarkit.util.picard.SAMSequenceDictionaryProgress)1 SAMSequenceDictionary (htsjdk.samtools.SAMSequenceDictionary)1 LiftOver (htsjdk.samtools.liftover.LiftOver)1