Java Code Examples for htsjdk.variant.variantcontext.VariantContext#isIndel()

@Override
protected void secondPassApply(VariantContext variant, ReadsContext readsContext, ReferenceContext referenceContext, FeatureContext featureContext) {
    final VariantContextBuilder builder = new VariantContextBuilder(variant);

    if (removeOldFilters) {
        builder.unfiltered();
    }

    if (variant.hasAttribute(infoKey)) {
        final double score = Double.parseDouble((String) variant.getAttribute(infoKey));
        if (variant.isSNP() && snpCutoffs.size() != 0 && isTrancheFiltered(score, snpCutoffs)) {
            builder.filter(filterStringFromScore(SNPString, score, snpTranches, snpCutoffs));
            filteredSnps++;
        } else if (variant.isIndel() && indelCutoffs.size() != 0 && isTrancheFiltered(score, indelCutoffs)) {
            builder.filter(filterStringFromScore(INDELString, score, indelTranches, indelCutoffs));
            filteredIndels++;
        }
    }

    if (builder.getFilters() == null || builder.getFilters().size() == 0){
        builder.passFilters();
    }
    
    vcfWriter.add(builder.make());
}
 

public void add(final VariantContext vc) {
    if (lastIndel == null && vc.isIndel()) {
        buffer.add(vc);
        lastIndel = vc;
    } else if (nearby(lastIndel, vc)) {
        if (vc.isIndel()) {
            emitAllNonIndels(); // throw out {@code lastIndel} and {@code vc}
        } else {
            buffer.add(vc);
        }
    } else {
        emitAllVariants();
        buffer.add(vc);
    }

    lastIndel = vc.isIndel() ? vc : lastIndel;
}
 

@Override
public Map<String, Object> annotate(final ReferenceContext ref,
                                    final VariantContext vc,
                                    final AlleleLikelihoods<GATKRead, Allele> likelihoods) {
    Utils.nonNull(vc);
    if ( !vc.isIndel()) {
        return Collections.emptyMap();
    }

    final Pair<List<Integer>,byte[]> result = getNumTandemRepeatUnits(ref, vc);
    if (result == null) {
        return Collections.emptyMap();
    }

    final byte[] repeatUnit = result.getRight();
    final List<Integer> numUnits = result.getLeft();

    final Map<String, Object> map = new LinkedHashMap<>();
    map.put(GATKVCFConstants.STR_PRESENT_KEY, true);
    map.put(GATKVCFConstants.REPEAT_UNIT_KEY, new String(repeatUnit));
    map.put(GATKVCFConstants.REPEATS_PER_ALLELE_KEY, numUnits);
    return Collections.unmodifiableMap(map);
}
 

private void infoDifference(final VariantContext eval, final VariantContext truth) {
    if(eval.hasAttribute(this.evalInfoKey) && truth.hasAttribute(truthInfoKey)) {
        final double evalVal = Double.valueOf((String) eval.getAttribute(this.evalInfoKey));
        final double truthVal = Double.valueOf((String) truth.getAttribute(this.truthInfoKey));
        final double delta = evalVal - truthVal;
        final double deltaSquared = delta * delta;
        if (eval.isSNP()) {
            this.snpSumDelta += Math.sqrt(deltaSquared);
            this.snpSumDeltaSquared += deltaSquared;
        } else if (eval.isIndel()) {
            this.indelSumDelta += Math.sqrt(deltaSquared);
            this.indelSumDeltaSquared += deltaSquared;
        }
        if (warnBigDifferences && Math.abs(delta) > this.epsilon) {
            this.logger.warn(String.format("Difference (%f) greater than epsilon (%f) at %s:%d %s:", delta, this.epsilon, eval.getContig(), eval.getStart(), eval.getAlleles().toString()));
            this.logger.warn(String.format("\t\tTruth info: " + truth.getAttributes().toString()));
            this.logger.warn(String.format("\t\tEval info: " + eval.getAttributes().toString()));
        }
    }
}
 

/**
 * Make sure that the order of the output file is identical to the order
 * of the input files by iterating through the output, making sure that,
 * if the context is an indel (snp), the next genomic position in the indel
 * (snp) queue is the same. Also make sure that the context is in the order
 * specified by the input files.
 */
private void validateSnpAndIndelResults(final File output, final Queue<String> indelContigPositions, final Queue<String> snpContigPositions) {
    final VCFFileReader outputReader = new VCFFileReader(output, false);
    final VariantContextComparator outputComparator = outputReader.getFileHeader().getVCFRecordComparator();
    VariantContext last = null;
    final CloseableIterator<VariantContext> iterator = outputReader.iterator();
    while (iterator.hasNext()) {
        final VariantContext outputContext = iterator.next();
        if (outputContext.isIndel()) Assert.assertEquals(getContigPosition(outputContext), indelContigPositions.poll());
        if (outputContext.isSNP()) Assert.assertEquals(getContigPosition(outputContext), snpContigPositions.poll());
        if (last != null) Assert.assertTrue(outputComparator.compare(last, outputContext) <= 0);
        last = outputContext;
    }
    iterator.close();

    // We should have polled everything off the indel (snp) queues
    Assert.assertEquals(indelContigPositions.size(), 0);
    Assert.assertEquals(snpContigPositions.size(), 0);
}
 

@Override
public void firstPassApply(final VariantContext variant, final ReadsContext readsContext, final ReferenceContext referenceContext, final FeatureContext featureContext) {
    if (!variant.hasAttribute(infoKey)){
        return;
    } else if (variant.isSNP()){
        scoredSnps++;
    } else if (variant.isIndel()){
        scoredIndels++;
    }

    for (FeatureInput<VariantContext> featureSource : resources) {
        for (VariantContext v : featureContext.getValues(featureSource)) {
            for (final Allele a : variant.getAlternateAlleles()) {
                if ((variant.getStart() == v.getStart()) && GATKVariantContextUtils.isAlleleInList(variant.getReference(), a, v.getReference(), v.getAlternateAlleles())) {
                    if (variant.isSNP()) {
                        resourceSNPScores.add(Double.parseDouble((String) variant.getAttribute(infoKey)));
                        return;
                    } else {
                        resourceIndelScores.add(Double.parseDouble((String)variant.getAttribute(infoKey)));
                        return;
                    }
                }
            }
        }
    }
}
 

private static void applySTRFilter(final VariantContext vc, final Collection<String> filters) {
    // STR contractions, such as ACTACTACT -> ACTACT, are overwhelmingly false positives so we hard filter by default
    if (vc.isIndel()) {
        final int[] rpa = vc.getAttributeAsList(GATKVCFConstants.REPEATS_PER_ALLELE_KEY).stream()
                .mapToInt(o -> Integer.parseInt(String.valueOf(o))).toArray();
        final String ru = vc.getAttributeAsString(GATKVCFConstants.REPEAT_UNIT_KEY, "");
        if (rpa != null && rpa.length > 1 && ru.length() > 1) {
            final int refCount = rpa[0];
            final int altCount = rpa[1];

            if (refCount - altCount == 1) {
                filters.add(GATKVCFConstants.STR_CONTRACTION_FILTER_NAME);
            }
        }
    }
}
 

public void add(final VariantContext vc) {
    if (lastIndel == null && vc.isIndel()) {
        buffer.add(vc);
        lastIndel = vc;
    } else if (nearby(lastIndel, vc)) {
        if (vc.isIndel()) {
            emitAllNonIndels(); // throw out {@code lastIndel} and {@code vc}
        } else {
            buffer.add(vc);
        }
    } else {
        emitAllVariants();
        buffer.add(vc);
    }

    lastIndel = vc.isIndel() ? vc : lastIndel;
}
 

@NotNull
public static VariantContext enrich(@NotNull final VariantContext context) {
    if (!context.isIndel() && !context.isSNP() ) {
        return context;
    }

    final VariantContextBuilder contextBuilder = new VariantContextBuilder(context).noGenotypes();

    final List<Allele> alleles = context.getAlleles();
    final Function<Allele, String> alleleKey = alleleKey(context);

    final List<Genotype> updatedGenotypes = Lists.newArrayList();
    for (Genotype genotype : context.getGenotypes()) {
        if (!genotype.hasAD() && hasRequiredAttributes(genotype, alleles, alleleKey)) {
            final GenotypeBuilder genotypeBuilder = new GenotypeBuilder(genotype).AD(readAD(genotype, alleles, alleleKey));
            updatedGenotypes.add(genotypeBuilder.make());
        } else {
            updatedGenotypes.add(genotype);
        }
    }

    return contextBuilder.genotypes(updatedGenotypes).make();
}
 

@Override
public List<Object> getRelevantStates(ReferenceContext referenceContext, ReadsContext readsContext, FeatureContext featureContext, VariantContext comp, String compName, VariantContext eval, String evalName, String sampleName, String FamilyName) {
    if (eval != null && eval.isIndel()) {
        for ( int l : eval.getIndelLengths() )
            if ( Math.abs(l) > 1 )
                return TWO_PLUS_BP; // someone is too long
        return ONE_BP; // all lengths are one
    } else
        return ALL;
}
 

@VisibleForTesting
static int qualityScore(@NotNull final VariantContext variant) {
    if (variant.isSNP()) {
        return getIntField(variant, SNP_QUAL_FIELD);
    } else if (variant.isIndel()) {
        return getIntField(variant, INDEL_QUAL_FIELD);
    } else {
        throw new IllegalStateException("record is not indel or snp: " + variant);
    }
}
 

/**
 * Returns a trimming result object from which the variant trimmed region and flanking non-variant sections
 * can be recovered latter.
 *
 * @param region the genome location range to trim.
 * @param variants list of variants contained in the trimming location. Variants therein
 *                                        not overlapping with {@code region} are simply ignored.
 * @param referenceContext
 * @return never {@code null}.
 */
public Result trim(final AssemblyRegion region, final SortedSet<VariantContext> variants, ReferenceContext referenceContext) {
    final List<VariantContext> variantsInRegion = variants.stream().filter(region::overlaps).collect(Collectors.toList());

    if ( variantsInRegion.isEmpty() ) {
        return noVariation(region);
    }

    int minStart = variantsInRegion.stream().mapToInt(VariantContext::getStart).min().getAsInt();
    int maxEnd = variantsInRegion.stream().mapToInt(VariantContext::getEnd).max().getAsInt();
    final SimpleInterval variantSpan = new SimpleInterval(region.getContig(), minStart, maxEnd).intersect(region);

    for (final VariantContext vc : variantsInRegion) {
        int padding = assemblyRegionArgs.snpPaddingForGenotyping;
        if (vc.isIndel()) {
            padding = assemblyRegionArgs.indelPaddingForGenotyping;
            final Pair<List<Integer>, byte[]> numRepeatsAndUnit = TandemRepeat.getNumTandemRepeatUnits(referenceContext, vc);
            if (numRepeatsAndUnit != null && numRepeatsAndUnit.getRight() != null) {
                final int repeatLength = numRepeatsAndUnit.getRight() == null ? 0 : numRepeatsAndUnit.getRight().length;
                final int mostRepeats = numRepeatsAndUnit.getLeft().stream().max(Integer::compareTo).orElse(0);
                final int longestSTR = mostRepeats * repeatLength;
                padding = assemblyRegionArgs.strPaddingForGenotyping + longestSTR;
            }
        }

        minStart = Math.min(minStart, Math.max(vc.getStart() - padding,1));
        maxEnd = Math.max(maxEnd, vc.getEnd() + padding);
    }

    final SimpleInterval paddedVariantSpan = new SimpleInterval(region.getContig(), minStart, maxEnd).intersect(region.getPaddedSpan());

    return new Result(region, variantSpan, paddedVariantSpan);
}
 

@NotNull
private static Function<Allele, String> alleleKey(@NotNull final VariantContext variant) {
    if (variant.isSNP()) {
        return StrelkaAllelicDepth::snpAlleleKey;
    } else if (variant.isIndel()) {
        return StrelkaAllelicDepth::indelAlleleKey;
    } else {
        throw new IllegalStateException("record is neither indel nor snp: " + variant);
    }
}
 

/** 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);
}
 

@NotNull
private Optional<RepeatContext> getRepeatContext(@NotNull final VariantContext variant, int relativePosition,
        @NotNull final String sequence) {
    if (variant.isIndel()) {
        return RepeatContextFactory.repeats(relativePosition + 1, sequence);
    } else if (variant.isSNP() || variant.isMNP()) {
        int altLength = variant.getAlternateAllele(0).getBaseString().length();

        Optional<RepeatContext> priorRepeat = RepeatContextFactory.repeats(relativePosition - 1, sequence);
        Optional<RepeatContext> postRepeat = RepeatContextFactory.repeats(relativePosition + altLength, sequence);
        return max(priorRepeat, postRepeat);
    } else {
        return Optional.empty();
    }
}
 

protected static boolean checkVariationClass( final VariantContext evalVC, final VariantRecalibratorArgumentCollection.Mode mode ) {
    switch( mode ) {
        case SNP:
            return evalVC.isSNP() || evalVC.isMNP();
        case INDEL:
            return evalVC.isStructuralIndel() || evalVC.isIndel() || evalVC.isMixed() || evalVC.isSymbolic();
        case BOTH:
            return true;
        default:
            throw new IllegalStateException( "Encountered unknown recal mode: " + mode );
    }
}
 

private static boolean isPassIndel(@NotNull final VariantContext context) {
    int altLength = alt(context).length();
    int refLength = context.getReference().getBaseString().length();

    return context.isIndel() && refLength < MAX_REF_ALT_LENGTH && altLength < MAX_REF_ALT_LENGTH && PASS.test(context);
}
 

/** 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<VariantType>> bitSetsToVariantTypes,
                            final IntervalList intervals,
                            final Optional<Log> log) {

    final Optional<ProgressLogger> progress = log.map(l -> new ProgressLogger(l, (int) 1e5, "Read", "variants"));
    final VCFFileReader variantReader = new VCFFileReader(dbSnpFile, intervals != null);
    final Iterator<VariantContext> variantIterator;
    if (intervals != null) {
        variantIterator = new ByIntervalListVariantContextIterator(variantReader, intervals);
    }
    else {
        variantIterator = variantReader.iterator();
    }

    while (variantIterator.hasNext()) {
        final VariantContext kv = variantIterator.next();

        for (final Map.Entry<DbSnpBitSetUtil, Set<VariantType>> tuple : bitSetsToVariantTypes.entrySet()) {
            final DbSnpBitSetUtil bitset            = tuple.getKey();
            final Set<VariantType> 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(VariantType.SNP)) ||
                    (kv.isIndel() && variantsToMatch.contains(VariantType.insertion)) ||
                    (kv.isIndel() && variantsToMatch.contains(VariantType.deletion))) {

                for (int i = kv.getStart(); i <= kv.getEnd(); i++)  bits.set(i, true);
            }
        }
        progress.map(p -> p.record(kv.getContig(), kv.getStart()));
    }

    CloserUtil.close(variantReader);
}
 

public static GATKVariant sparkVariantAdapter(VariantContext vc) {
    return new MinimalVariant(new SimpleInterval(vc.getContig(),vc.getStart(),vc.getEnd()), vc.isSNP(), vc.isIndel());
}