Java Code Examples for htsjdk.variant.variantcontext.VariantContextBuilder#filters()

@Override
public Object onTraversalSuccess() {
    final String tumorSample = getHeaderForVariants().getMetaDataLine(Mutect2Engine.TUMOR_SAMPLE_KEY_IN_VCF_HEADER).getValue();
    final Mutect2FilteringEngine filteringEngine = new Mutect2FilteringEngine(MTFAC, tumorSample);
    // TODO: implement sophisticated filtering
    for (final VariantContext vc : unfilteredCalls) {
        final VariantContextBuilder vcb = new VariantContextBuilder(vc);
        vcb.filters(filteringEngine.calculateFilters(MTFAC, vc));
        vcfWriter.add(vcb.make());
    }
    return "SUCCESS";
}
 

@Override
public void apply(final VariantContext vc, final ReadsContext readsContext, final ReferenceContext ref, final FeatureContext featureContext) {

    final List<VariantContext> recals =  featureContext.getValues(recal, vc.getStart());
    final boolean evaluateThisVariant = useASannotations || VariantDataManager.checkVariationClass( vc, MODE );

    //vc.isNotFiltered is true for PASS; vc.filtersHaveBeenApplied covers PASS and filters
    final boolean variantIsNotFiltered = IGNORE_ALL_FILTERS || vc.isNotFiltered() ||
            (!ignoreInputFilterSet.isEmpty() && ignoreInputFilterSet.containsAll(vc.getFilters()));

    if( evaluateThisVariant && variantIsNotFiltered) {
        String filterString;
        final VariantContextBuilder builder = new VariantContextBuilder(vc);
        if (!useASannotations) {
            filterString = doSiteSpecificFiltering(vc, recals, builder);
        }
        else {  //allele-specific mode
            filterString = doAlleleSpecificFiltering(vc, recals, builder);
        }

        //for both non-AS and AS modes:
        if( filterString.equals(VCFConstants.PASSES_FILTERS_v4) ) {
            builder.passFilters();
        } else if(filterString.equals(VCFConstants.UNFILTERED)) {
            builder.unfiltered();
        } else {
            builder.filters(filterString);
        }

        final VariantContext outputVC = builder.make();
        if( !EXCLUDE_FILTERED || outputVC.isNotFiltered() ) {
            vcfWriter.add( outputVC );
        }
    } else { // valid VC but not compatible with this mode, so just emit the variant untouched
        vcfWriter.add( vc );
    }
}
 

/**
 * Merges multiple VariantContexts all for the same locus into a single hybrid.
 *
 * @param variantContexts           list of VCs
 * @return new VariantContext       representing the merge of variantContexts
 */
public static VariantContext merge(final List<VariantContext> variantContexts) {
    if ( variantContexts == null || variantContexts.isEmpty() )
        return null;

    // establish the baseline info from the first VC
    final VariantContext first = variantContexts.get(0);
    final String name = first.getSource();

    final Set<String> filters = new HashSet<>();

    int depth = 0;
    double log10PError = CommonInfo.NO_LOG10_PERROR;
    boolean anyVCHadFiltersApplied = false;
    GenotypesContext genotypes = GenotypesContext.create();

    // Go through all the VCs, verify that the loci and ID and other attributes agree.
    final Map<String, Object> firstAttributes = first.getAttributes();
    for (final VariantContext vc : variantContexts ) {
        if ((vc.getStart() != first.getStart()) || (!vc.getContig().equals(first.getContig()))) {
            throw new PicardException("Mismatch in loci among input VCFs");
        }
        if (!vc.getID().equals(first.getID())) {
            throw new PicardException("Mismatch in ID field among input VCFs");
        }
        if (!vc.getReference().equals(first.getReference())) {
            throw new PicardException("Mismatch in REF allele among input VCFs");
        }
        checkThatAllelesMatch(vc, first);

        genotypes.addAll(vc.getGenotypes());

        // We always take the QUAL of the first VC with a non-MISSING qual for the combined value
        if ( log10PError == CommonInfo.NO_LOG10_PERROR )
            log10PError =  vc.getLog10PError();

        filters.addAll(vc.getFilters());
        anyVCHadFiltersApplied |= vc.filtersWereApplied();

        // add attributes
        // special case DP (add it up)
        if (vc.hasAttribute(VCFConstants.DEPTH_KEY))
            depth += vc.getAttributeAsInt(VCFConstants.DEPTH_KEY, 0);

        // Go through all attributes - if any differ (other than AC, AF, or AN) that's an error.  (recalc AC, AF, and AN later)
        for (final Map.Entry<String, Object> p : vc.getAttributes().entrySet()) {
            final String key = p.getKey();
            if ((!key.equals("AC")) && (!key.equals("AF")) && (!key.equals("AN")) && (!key.equals("devX_AB")) && (!key.equals("devY_AB"))) {
                final Object value = p.getValue();
                final Object extantValue = firstAttributes.get(key);
                if (extantValue == null) {
                    // attribute in one VCF but not another.  Die!
                    throw new PicardException("Attribute '" + key + "' not found in all VCFs");
                }
                else if (!extantValue.equals(value)) {
                    // Attribute disagrees in value between one VCF Die! (if not AC, AF, nor AN, skipped above)
                    throw new PicardException("Values for attribute '" + key + "' disagrees among input VCFs");
                }
            }
        }
    }

    if (depth > 0)
        firstAttributes.put(VCFConstants.DEPTH_KEY, String.valueOf(depth));

    final VariantContextBuilder builder = new VariantContextBuilder().source(name).id(first.getID());
    builder.loc(first.getContig(), first.getStart(), first.getEnd());
    builder.alleles(first.getAlleles());
    builder.genotypes(genotypes);

    builder.attributes(new TreeMap<>(firstAttributes));
    // AC, AF, and AN are recalculated here
    VariantContextUtils.calculateChromosomeCounts(builder, false);

    builder.log10PError(log10PError);
    if ( anyVCHadFiltersApplied ) {
        builder.filters(filters.isEmpty() ? filters : new TreeSet<>(filters));
    }

    return builder.make();
}
 

/**
 * This will take an input VariantContext and lift to the provided interval.
 * If this interval is in the opposite orientation, all alleles and genotypes will be reverse complemented
 * and indels will be left-aligned.  Currently this is only able to invert biallelic indels, and null will be
 * returned for any unsupported VC.
 *
 * @param source                The VariantContext to lift
 * @param target                The target interval
 * @param refSeq                The reference sequence, which should match the target interval
 * @param writeOriginalPosition If true, INFO field annotations will be added to store the original position and contig
 * @param writeOriginalAlleles  If true, an INFO field annotation will be added to store the original alleles in order.  This can be useful in the case of more complex liftovers, like reverse complements, left aligned indels or swapped REF/ALT
 * @return The lifted VariantContext.  This will be null if the input VariantContext could not be lifted.
 */
public static VariantContext liftVariant(final VariantContext source, final Interval target, final ReferenceSequence refSeq, final boolean writeOriginalPosition, final boolean writeOriginalAlleles) {
    if (target == null) {
        return null;
    }

    final VariantContextBuilder builder;
    if (target.isNegativeStrand()) {
        builder = reverseComplementVariantContext(source, target, refSeq);
    } else {
        builder = liftSimpleVariantContext(source, target);
    }

    if (builder == null) {
        return null;
    }

    builder.filters(source.getFilters());
    builder.log10PError(source.getLog10PError());

    // If any of the source alleles is symbolic, do not populate the END tag (protecting things like <NON_REF> from
    // getting screwed up in reverse-complemented variants
    if (source.hasAttribute(VCFConstants.END_KEY) && builder.getAlleles().stream().noneMatch(Allele::isSymbolic)) {
        builder.attribute(VCFConstants.END_KEY, (int) builder.getStop());
    } else {
        builder.rmAttribute(VCFConstants.END_KEY);
    }

    // make sure that the variant isn't mistakenly set as "SwappedAlleles"
    builder.rmAttribute(SWAPPED_ALLELES);
    if (target.isNegativeStrand()) {
        builder.attribute(REV_COMPED_ALLELES, true);
    } else {
        // make sure that the variant isn't mistakenly set as "ReverseComplementedAlleles" (from a previous liftover, say)
        builder.rmAttribute(REV_COMPED_ALLELES);
    }

    builder.id(source.getID());

    if (writeOriginalPosition) {
        builder.attribute(LiftoverVcf.ORIGINAL_CONTIG, source.getContig());
        builder.attribute(LiftoverVcf.ORIGINAL_START, source.getStart());
    }

    if (writeOriginalAlleles && !source.getAlleles().equals(builder.getAlleles())) {
        builder.attribute(LiftoverVcf.ORIGINAL_ALLELES, allelesToStringList(source.getAlleles()));
    }

    return builder.make();
}