Search in sources :

Example 26 with VariantContext

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

the class ApplyVQSR method doSiteSpecificFiltering.

/**
     * Calculate the filter status for a given VariantContext using the combined data from all alleles at a site
     * @param vc
     * @param recals
     * @param builder   is modified by adding attributes
     * @return a String with the filter status for this site
     */
private String doSiteSpecificFiltering(final VariantContext vc, final List<VariantContext> recals, final VariantContextBuilder builder) {
    VariantContext recalDatum = getMatchingRecalVC(vc, recals, null);
    if (recalDatum == null) {
        throw new UserException("Encountered input variant which isn't found in the input recal file. Please make sure VariantRecalibrator and ApplyRecalibration were run on the same set of input variants. First seen at: " + vc);
    }
    final String lodString = recalDatum.getAttributeAsString(GATKVCFConstants.VQS_LOD_KEY, null);
    if (lodString == null) {
        throw new UserException("Encountered a malformed record in the input recal file. There is no lod for the record at: " + vc);
    }
    final double lod;
    try {
        lod = Double.valueOf(lodString);
    } catch (NumberFormatException e) {
        throw new UserException("Encountered a malformed record in the input recal file. The lod is unreadable for the record at: " + vc);
    }
    builder.attribute(GATKVCFConstants.VQS_LOD_KEY, lod);
    builder.attribute(GATKVCFConstants.CULPRIT_KEY, recalDatum.getAttribute(GATKVCFConstants.CULPRIT_KEY));
    if (recalDatum != null) {
        if (recalDatum.hasAttribute(GATKVCFConstants.POSITIVE_LABEL_KEY))
            builder.attribute(GATKVCFConstants.POSITIVE_LABEL_KEY, true);
        if (recalDatum.hasAttribute(GATKVCFConstants.NEGATIVE_LABEL_KEY))
            builder.attribute(GATKVCFConstants.NEGATIVE_LABEL_KEY, true);
    }
    return generateFilterString(lod);
}
Also used : VariantContext(htsjdk.variant.variantcontext.VariantContext) UserException(org.broadinstitute.hellbender.exceptions.UserException)

Example 27 with VariantContext

use of htsjdk.variant.variantcontext.VariantContext in project gatk 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 28 with VariantContext

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

the class GVCFWriter method add.

/**
     * Add a VariantContext to this writer for emission
     *
     * Requires that the VC have exactly one genotype
     *
     * @param vc a non-null VariantContext
     */
@Override
public void add(VariantContext vc) {
    Utils.nonNull(vc);
    Utils.validateArg(vc.hasGenotypes(), "GVCF assumes that the VariantContext has genotypes");
    Utils.validateArg(vc.getGenotypes().size() == 1, () -> "GVCF assumes that the VariantContext has exactly one genotype but saw " + vc.getGenotypes().size());
    if (sampleName == null) {
        sampleName = vc.getGenotype(0).getSampleName();
    }
    if (currentBlock != null && !currentBlock.isContiguous(vc)) {
        // we've made a non-contiguous step (across interval, onto another chr), so finalize
        emitCurrentBlock();
    }
    final Genotype g = vc.getGenotype(0);
    if (g.isHomRef() && vc.hasAlternateAllele(GATKVCFConstants.NON_REF_SYMBOLIC_ALLELE) && vc.isBiallelic()) {
        // create bands
        final VariantContext maybeCompletedBand = addHomRefSite(vc, g);
        if (maybeCompletedBand != null) {
            underlyingWriter.add(maybeCompletedBand);
        }
    } else {
        // g is variant, so flush the bands and emit vc
        emitCurrentBlock();
        nextAvailableStart = vc.getEnd();
        contigOfNextAvailableStart = vc.getContig();
        underlyingWriter.add(vc);
    }
}
Also used : Genotype(htsjdk.variant.variantcontext.Genotype) VariantContext(htsjdk.variant.variantcontext.VariantContext)

Example 29 with VariantContext

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

the class GVCFWriter method addHomRefSite.

/**
     * Add hom-ref site from vc to this gVCF hom-ref state tracking, emitting any pending states if appropriate
     *
     * @param vc a non-null VariantContext
     * @param g  a non-null genotype from VariantContext
     * @return a VariantContext to be emitted, or null if non is appropriate
     */
protected VariantContext addHomRefSite(final VariantContext vc, final Genotype g) {
    if (nextAvailableStart != -1) {
        // don't create blocks while the hom-ref site falls before nextAvailableStart (for deletions)
        if (vc.getStart() <= nextAvailableStart && vc.getContig().equals(contigOfNextAvailableStart)) {
            return null;
        }
        // otherwise, reset to non-relevant
        nextAvailableStart = -1;
        contigOfNextAvailableStart = null;
    }
    final VariantContext result;
    if (genotypeCanBeMergedInCurrentBlock(g)) {
        currentBlock.add(vc.getStart(), g);
        result = null;
    } else {
        result = currentBlock != null ? currentBlock.toVariantContext(sampleName) : null;
        currentBlock = createNewBlock(vc, g);
    }
    return result;
}
Also used : VariantContext(htsjdk.variant.variantcontext.VariantContext)

Example 30 with VariantContext

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

the class DbSnpBitSetUtil method loadVcf.

/** Private helper method to read through the VCF and create one or more bit sets. */
private static void loadVcf(final File dbSnpFile, final SAMSequenceDictionary sequenceDictionary, final Map<DbSnpBitSetUtil, Set<DbSnpVariantType>> bitSetsToVariantTypes) {
    final VCFFileReader variantReader = new VCFFileReader(dbSnpFile);
    final CloseableIterator<VariantContext> variantIterator = variantReader.iterator();
    while (variantIterator.hasNext()) {
        final VariantContext kv = variantIterator.next();
        for (final Map.Entry<DbSnpBitSetUtil, Set<DbSnpVariantType>> tuple : bitSetsToVariantTypes.entrySet()) {
            final DbSnpBitSetUtil bitset = tuple.getKey();
            final Set<DbSnpVariantType> variantsToMatch = tuple.getValue();
            BitSet bits = bitset.sequenceToBitSet.get(kv.getContig());
            if (bits == null) {
                final int nBits;
                if (sequenceDictionary == null)
                    nBits = kv.getEnd() + 1;
                else
                    nBits = sequenceDictionary.getSequence(kv.getContig()).getSequenceLength() + 1;
                bits = new BitSet(nBits);
                bitset.sequenceToBitSet.put(kv.getContig(), bits);
            }
            if (variantsToMatch.isEmpty() || (kv.isSNP() && variantsToMatch.contains(DbSnpVariantType.SNP)) || (kv.isIndel() && variantsToMatch.contains(DbSnpVariantType.insertion)) || (kv.isIndel() && variantsToMatch.contains(DbSnpVariantType.deletion))) {
                for (int i = kv.getStart(); i <= kv.getEnd(); i++) bits.set(i, true);
            }
        }
    }
    CloserUtil.close(variantIterator);
    CloserUtil.close(variantReader);
}
Also used : VCFFileReader(htsjdk.variant.vcf.VCFFileReader) VariantContext(htsjdk.variant.variantcontext.VariantContext)

Aggregations

VariantContext (htsjdk.variant.variantcontext.VariantContext)237 Test (org.testng.annotations.Test)119 File (java.io.File)98 BaseTest (org.broadinstitute.hellbender.utils.test.BaseTest)68 Allele (htsjdk.variant.variantcontext.Allele)55 Genotype (htsjdk.variant.variantcontext.Genotype)49 VariantContextBuilder (htsjdk.variant.variantcontext.VariantContextBuilder)43 Collectors (java.util.stream.Collectors)43 SimpleInterval (org.broadinstitute.hellbender.utils.SimpleInterval)43 CommandLineProgramTest (org.broadinstitute.hellbender.CommandLineProgramTest)41 FeatureDataSource (org.broadinstitute.hellbender.engine.FeatureDataSource)36 StandardArgumentDefinitions (org.broadinstitute.hellbender.cmdline.StandardArgumentDefinitions)35 java.util (java.util)33 IOException (java.io.IOException)25 Assert (org.testng.Assert)25 SAMSequenceDictionary (htsjdk.samtools.SAMSequenceDictionary)24 VCFFileReader (htsjdk.variant.vcf.VCFFileReader)22 UserException (org.broadinstitute.hellbender.exceptions.UserException)22 Target (org.broadinstitute.hellbender.tools.exome.Target)22 DataProvider (org.testng.annotations.DataProvider)22