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

The following examples show how to use htsjdk.variant.variantcontext.VariantContextBuilder#make() . You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example 1
Source File: VariantEvalUtils.java    From gatk with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
public VariantContext ensureAnnotations(final VariantContext vc, final VariantContext vcsub) {
    final int originalAlleleCount = vc.getHetCount() + 2 * vc.getHomVarCount();
    final int newAlleleCount = vcsub.getHetCount() + 2 * vcsub.getHomVarCount();
    final boolean isSingleton = originalAlleleCount == newAlleleCount && newAlleleCount == 1;
    final boolean hasChrCountAnnotations = vcsub.hasAttribute(VCFConstants.ALLELE_COUNT_KEY) &&
            vcsub.hasAttribute(VCFConstants.ALLELE_FREQUENCY_KEY) &&
            vcsub.hasAttribute(VCFConstants.ALLELE_NUMBER_KEY);

    if ( ! isSingleton && hasChrCountAnnotations ) {
        // nothing to update
        return vcsub;
    } else {
        // have to do the work
        VariantContextBuilder builder = new VariantContextBuilder(vcsub);

        if ( isSingleton )
            builder.attribute(VariantEval.IS_SINGLETON_KEY, true);

        if ( ! hasChrCountAnnotations )
            VariantContextUtils.calculateChromosomeCounts(builder, true);

        return builder.make();
    }
}
 
Example 2
Source File: AnnotatedIntervalToSegmentVariantContextConverter.java    From gatk with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
/**
 * Convert a segment (as annotated interval) to a {@link VariantContext}.  The variant context will have a ref allele
 *  of the starting base of the segment.  The end position will be in the attribute {@link VCFConstants#END_KEY}.
 *  Any remaining annotations in the annotated interval will be converted to string attributes.  The segment call
 *  will be the alternate allele.  Currently, the alternate allele will be one of neutral, amplification, or deletion.
 *
 * The variant context will use symbolic alleles for the calls (copy neutral, amplification, or deletion).
 * See {@link CalledCopyRatioSegment.Call} for insertions and and deletions.
 *
 * @param segment Never {@code null}
 * @param referenceContext Never {@code null}
 * @return Never {@code null}
 */
public static VariantContext convert(final AnnotatedInterval segment, final ReferenceContext referenceContext) {
    Utils.nonNull(segment);
    Utils.nonNull(referenceContext);
    final SimpleInterval startPointInterval = new SimpleInterval(segment.getContig(), segment.getStart(), segment.getStart());
    final Allele refAlleleAtStart = Allele.create(referenceContext.getBases(startPointInterval), true);

    final CalledCopyRatioSegment.Call call = retrieveCall(segment);

    final VariantContextBuilder variantContextBuilder = new VariantContextBuilder().chr(segment.getContig())
            .start(segment.getStart())
            .stop(segment.getEnd())
            .alleles(Arrays.asList(
                    Allele.create(refAlleleAtStart, false),
                    convertActualSegCallToAllele(call)
            ))
            .attribute(VCFConstants.END_KEY, segment.getEnd());

    // Grab every field in the annotated interval and make it a string attribute
    segment.getAnnotations().keySet().stream().forEach(k -> variantContextBuilder.attribute(k, segment.getAnnotationValue(k)));

    return variantContextBuilder
            .make();
}
 
Example 3
Source File: SelectVariants.java    From gatk with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
private VariantContext buildVariantContextWithDroppedAnnotationsRemoved(final VariantContext vc) {
    if (infoAnnotationsToDrop.isEmpty() && genotypeAnnotationsToDrop.isEmpty()) {
        return vc;
    }
    final VariantContextBuilder rmAnnotationsBuilder = new VariantContextBuilder(vc);
    for (String infoField : infoAnnotationsToDrop) {
        rmAnnotationsBuilder.rmAttribute(infoField);
    }
    if (!genotypeAnnotationsToDrop.isEmpty()) {
        final ArrayList<Genotype> genotypesToWrite = new ArrayList<>();
        for (Genotype genotype : vc.getGenotypes()) {
            final GenotypeBuilder genotypeBuilder = new GenotypeBuilder(genotype).noAttributes();
            final Map<String, Object> attributes = new HashMap<>(genotype.getExtendedAttributes());
            for (String genotypeAnnotation : genotypeAnnotationsToDrop) {
                attributes.remove(genotypeAnnotation);
            }
            genotypeBuilder.attributes(attributes);
            genotypesToWrite.add(genotypeBuilder.make());
        }
        rmAnnotationsBuilder.genotypes(GenotypesContext.create(genotypesToWrite));
    }
    return rmAnnotationsBuilder.make();
}
 
Example 4
Source File: CosmicFuncotationFactoryUnitTest.java    From gatk with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
private VariantContext createVariantContext(final String contig,
                                            final int start,
                                            final int end,
                                            final String refString,
                                            final String altString) {

    final Allele refAllele = Allele.create(refString, true);
    final Allele altAllele = Allele.create(altString);

    final VariantContextBuilder variantContextBuilder = new VariantContextBuilder(
            FuncotatorReferenceTestUtils.retrieveHg19Chr3Ref(),
            contig,
            start,
            end,
            Arrays.asList(refAllele, altAllele)
    );

    return variantContextBuilder.make();
}
 
Example 5
Source File: AmberVCF.java    From hmftools with GNU General Public License v3.0 5 votes vote down vote up
@NotNull
private VariantContext create(@NotNull final TumorBAF tumorBaf, final List<Genotype> genotypes) {

    final List<Allele> alleles = alleles(tumorBaf);

    final VariantContextBuilder builder = new VariantContextBuilder().chr(tumorBaf.chromosome())
            .start(tumorBaf.position())
            .computeEndFromAlleles(alleles, (int) tumorBaf.position())
            .genotypes(genotypes)
            .alleles(alleles);

    final VariantContext context = builder.make();
    context.getCommonInfo().setLog10PError(tumorBaf.tumorAltQuality() / -10d);
    return context;
}
 
Example 6
Source File: VariantHotspotEnrichment.java    From hmftools with GNU General Public License v3.0 5 votes vote down vote up
@NotNull
private VariantContext enrich(@NotNull final VariantContext context) {
    VariantContextBuilder builder =
            new VariantContextBuilder(context).attribute(HOTSPOT_FLAG, false).attribute(NEAR_HOTSPOT_FLAG, false);

    if (hotspotEnrichment.isOnHotspot(context)) {
        return builder.attribute(HOTSPOT_FLAG, true).make();
    } else if (hotspotEnrichment.isNearHotspot(context)) {
        return builder.attribute(NEAR_HOTSPOT_FLAG, true).make();
    }

    return builder.make();

}
 
Example 7
Source File: SageHotspotAnnotation.java    From hmftools with GNU General Public License v3.0 5 votes vote down vote up
@NotNull
private VariantContext annotate(@NotNull final VariantContext context) {
    VariantContextBuilder builder =
            new VariantContextBuilder(context).attribute(HOTSPOT_FLAG, false).attribute(NEAR_HOTSPOT_FLAG, false);

    if (hotspotEnrichment.isOnHotspot(context)) {
        return builder.attribute(HOTSPOT_FLAG, true).make();
    } else if (hotspotEnrichment.isNearHotspot(context)) {
        return builder.attribute(NEAR_HOTSPOT_FLAG, true).make();
    }

    return builder.make();
}
 
Example 8
Source File: GATKToolUnitTest.java    From gatk with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
private void writeHeaderAndBadVariant(final VariantContextWriter writer) {
    final VariantContextBuilder vcBuilder = new VariantContextBuilder(
            "chr1","1", 1, 1, Arrays.asList(Allele.create("A", true)));
    vcBuilder.attribute("fake", new Object());
    final VariantContext vc = vcBuilder.make();
    final VCFHeader vcfHeader = new VCFHeader();
    writer.writeHeader(vcfHeader);
    writer.add(vc);
}
 
Example 9
Source File: AnnotatedVariantProducer.java    From gatk with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
public static VariantContext produceAnnotatedVcFromEvidenceTargetLink(final EvidenceTargetLink evidenceTargetLink,
                                                                      final SvType svType) {
    final PairedStrandedIntervals pairedStrandedIntervals = evidenceTargetLink.getPairedStrandedIntervals();
    final StrandedInterval strandedIntervalLeft = pairedStrandedIntervals.getLeft();
    final StrandedInterval strandedIntervalRight = pairedStrandedIntervals.getRight();
    final int start = strandedIntervalLeft.getInterval().midpoint();
    final int end = strandedIntervalRight.getInterval().midpoint();
    final VariantContextBuilder builder = svType
            .getBasicInformation()
            .attribute(GATKSVVCFConstants.CIPOS, produceCIInterval(start, strandedIntervalLeft.getInterval()))
            .attribute(GATKSVVCFConstants.CIEND, produceCIInterval(end, strandedIntervalRight.getInterval()))
            .attribute(GATKSVVCFConstants.READ_PAIR_SUPPORT, evidenceTargetLink.getReadPairs())
            .attribute(GATKSVVCFConstants.SPLIT_READ_SUPPORT, evidenceTargetLink.getSplitReads());
    return builder.make();
}
 
Example 10
Source File: SegmentedCpxVariantSimpleVariantExtractorUnitTest.java    From gatk with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
private VariantContext makeTestComplexVariant(final SimpleInterval affectedRefRegion, final int svLen,
                                              final String referenceBases, final String altSeqBases,
                                              final List<String> contigNames,
                                              final List<SimpleInterval> referenceSegments, final List<String> altArrangement) {
    final int maxAlignmentLength = random.nextInt(4000) + 1;
    int evidenceContigCount = contigNames.size();
    StringBuilder mqs = new StringBuilder();
    for (int i = 0; i <= evidenceContigCount; ++i) mqs.append(random.nextInt(61)).append(",");
    String mq = mqs.toString();

    final VariantContextBuilder builder = new VariantContextBuilder()
            .chr(affectedRefRegion.getContig()).start(affectedRefRegion.getStart()).stop(affectedRefRegion.getEnd())
            .alleles(Arrays.asList(Allele.create(referenceBases, true),
                    Allele.create(SimpleSVType.createBracketedSymbAlleleString(CPX_SV_SYB_ALT_ALLELE_STR))))
            .id(CPX_SV_SYB_ALT_ALLELE_STR + INTERVAL_VARIANT_ID_FIELD_SEPARATOR + affectedRefRegion.toString())
            .attribute(VCFConstants.END_KEY, affectedRefRegion.getEnd())
            .attribute(SVLEN, svLen)
            .attribute(SVTYPE, CPX_SV_SYB_ALT_ALLELE_STR)
            .attribute(CPX_EVENT_ALT_ARRANGEMENTS, String.join(VCFConstants.INFO_FIELD_ARRAY_SEPARATOR, altArrangement))
            .attribute(CONTIG_NAMES, String.join(VCFConstants.INFO_FIELD_ARRAY_SEPARATOR, contigNames))
            .attribute(SEQ_ALT_HAPLOTYPE, altSeqBases)
            .attribute(MAX_ALIGN_LENGTH, maxAlignmentLength)
            .attribute(MAPPING_QUALITIES, mq.substring(0, mq.length() - 1)); // drop last coma
    if (referenceSegments.isEmpty())
        return builder.make();
    else
        return builder
                .attribute(CPX_SV_REF_SEGMENTS, String.join(VCFConstants.INFO_FIELD_ARRAY_SEPARATOR,
                        referenceSegments.stream().map(SimpleInterval::toString).collect(Collectors.toList())))
                .make();
}
 
Example 11
Source File: FuncotatorEngine.java    From gatk with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
/**
 * Create a new {@link VariantContext} which will match the given Reference if there is a mismatch for input between the B37 reference and the HG19 reference.
 * @param variant A {@link VariantContext} object containing the variant to convert.
 * @return A {@link VariantContext} whose contig has been transformed to HG19 if requested by the user.  Otherwise, an identical variant.
 */
private VariantContext getCorrectVariantContextForReference(final VariantContext variant) {
    if ( mustConvertInputContigsToHg19 ) {
        final VariantContextBuilder vcb = new VariantContextBuilder(variant);
        vcb.chr(FuncotatorUtils.convertB37ContigToHg19Contig(variant.getContig()));
        return vcb.make();
    }
    else {
        return variant;
    }
}
 
Example 12
Source File: CalculateGenotypePosteriors.java    From gatk with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
@Override
public void apply(final VariantContext variant,
                  final ReadsContext readsContext,
                  final ReferenceContext referenceContext,
                  final FeatureContext featureContext) {

    final Collection<VariantContext> otherVCs = featureContext.getValues(supportVariants);

    //If no resource contains a matching variant, then add numRefIfMissing as a pseudocount to the priors
    List<VariantContext> resourcesWithMatchingStarts = otherVCs.stream()
            .filter(vc -> variant.getStart() == vc.getStart()).collect(Collectors.toList());

    //do family priors first (if applicable)
    final VariantContextBuilder builder = new VariantContextBuilder(variant);
    //only compute family priors for biallelelic sites
    if (!skipFamilyPriors && variant.isBiallelic()){
        final GenotypesContext gc = famUtils.calculatePosteriorGLs(variant);
        builder.genotypes(gc);
    }
    VariantContextUtils.calculateChromosomeCounts(builder, false);
    final VariantContext vc_familyPriors = builder.make();

    final VariantContext vc_bothPriors;
    if (!skipPopulationPriors) {
        vc_bothPriors = PosteriorProbabilitiesUtils.calculatePosteriorProbs(vc_familyPriors, resourcesWithMatchingStarts,
                resourcesWithMatchingStarts.isEmpty() ? numRefIfMissing : 0, options);
    } else {
        vc_bothPriors = vc_familyPriors;
    }
    vcfWriter.add(vc_bothPriors);
}
 
Example 13
Source File: AlleleFilterUtils.java    From gatk with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
/**
 * Adds the new allele filters to the existing allele filters in the vc. Computes whether there are
 * new site filters and updates the filter in the vc. If there are no site filters, sets filters to pass
 * Sets the filters for each allele and calculates the intersection of the allele filters to set on the variant.
 * PASS if the intersection is empty.
 * @param vc The variant context to add the filters to, both at the allele and site level
 * @param newAlleleFilters filters to be applied to each allele, the intersection of these filters are applied at the site level
 * @param invalidatePreviousFilters whether existing filters should be removed
 * @return The updated variant context
 */
public static VariantContext addAlleleAndSiteFilters(VariantContext vc, List<Set<String>> newAlleleFilters, boolean invalidatePreviousFilters) {
    if (newAlleleFilters.isEmpty()) {
        return vc;
    }
    List<List<String>> currentAlleleFilters = decodeASFilters(vc);
    if (!currentAlleleFilters.isEmpty() && newAlleleFilters.size() != currentAlleleFilters.size()) {
        // log an error
        return vc;
    }

    if (currentAlleleFilters.isEmpty() || invalidatePreviousFilters) {
        currentAlleleFilters = new ArrayList<>(Collections.nCopies(newAlleleFilters.size(), Collections.singletonList(GATKVCFConstants.SITE_LEVEL_FILTERS)));
    }
    ListIterator<List<String>> currentAlleleFiltersIt = currentAlleleFilters.listIterator();
    List<List<String>> updatedAlleleFilters = newAlleleFilters.stream().map(newfilters -> addAlleleFilters(currentAlleleFiltersIt.next(), newfilters.stream().collect(Collectors.toList()))).collect(Collectors.toList());
    String encodedFilters = encodeASFilters(updatedAlleleFilters);
    VariantContextBuilder vcb = new VariantContextBuilder(vc).attribute(GATKVCFConstants.AS_FILTER_STATUS_KEY, encodedFilters);

    if (invalidatePreviousFilters) {
        vcb.unfiltered();
    }
    Set<String> siteFiltersToAdd = newAlleleFilters.stream().skip(1)
            .collect(()->new HashSet<>(newAlleleFilters.get(0)), Set::retainAll, Set::retainAll);
    siteFiltersToAdd.stream().forEach(filter -> vcb.filter(filter));
    if ((vcb.getFilters() == null || vcb.getFilters().isEmpty()) && !invalidatePreviousFilters) {
        vcb.passFilters();
    }
    return vcb.make();
}
 
Example 14
Source File: FuncotateSegments.java    From gatk with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
private VariantContext transformAttributesToStandardNames(final VariantContext vc) {

        final Map<String,Object> transformedAttributes = vc.getAttributes().entrySet().stream()
                .collect(Collectors.toMap(e-> aliasToKeyMapping.getOrDefault(e.getKey(), e.getKey()), e -> e.getValue()));
        final VariantContextBuilder vcb = new VariantContextBuilder(vc).attributes(transformedAttributes);
        return vcb.make();
    }
 
Example 15
Source File: TestUtilsCpxVariantInference.java    From gatk with BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
private static PreprocessedAndAnalysisReadyContigWithExpectedResults buildForContig_2428_0(final Set<String> canonicalChromosomes,
                                                                                           final SAMSequenceDictionary refSeqDict) {
    final AssemblyContigWithFineTunedAlignments analysisReadyContig =
            makeContigAnalysisReadyForCpxBranch("asm002428:tig00000\t16\tchr2\t4450935\t60\t1216M3I256M1030S\t*\t0\t0\tTTCCTCTGTGGAGGACACAGCATTCAAGGCGCCATCTTGGAATCAGAATCACTAACCACATCAAAGCCTAATACTGGACTTCCCAACCTCCTAAACTGTAAGCAAATAATTTTGTTCATAAGTTACCCAGCCTGTGCTATTCTTCTATACCAGCACATATGGACTAAGACACTACGTATTCTAATAAGTGTGTGAAAACCCAAGATAAAATAGGCGTTATTTCCCATTTTACAGATGAGAAAACTGAAGCTTGAGGAATTTGTGAATTTGTTTAATCCAAATCTGTTAGAAAATAGTAGACCTCACCTCTTCAATAGTTGAGCATGGATAATTATGATACTAATTTTAGACTCATATTTCTTCATTTGTACATGATTTAAATCTCAGAAAAAATTTAGAATTTCTCTTTACTGTCTGTTGGTATTATGAAATGTGATAGCAGTATGTCTAGGCACAGATTTTTTTTTTTTTTTTCCTTAGCAGGGGTAGTTATTTTTTTCTTTCTTAGCAAGGGTGGCTTTTTCAGTCTGAAAGCTGAGACATACACTTTTCTCTGAGAAGAAATATTTTCTATTATTTGAATGTTTCACCTCTCCTTTTTCCTATTCTCTCTTTCTGAGACATACATTGATTCTCCAAGTTATCATCTTTCTCTCTTATCTCTATGTTTGTCTTTGGTTCTCTTTTTTCTCTGTATTCTGGGAAATTTTCATAGCTTTATTCAAGTTCATTTTTCAAGCTCTAAGAACTCTCTTTTTGTCTGATTTTACCATTTTTATTATACTAGTTGTTACAGTTTTGATAGATAAACTATCTTCCAAAATTTCTAAACGTATGGGTTTGGGAAATTGCTTTTAAAATCCTTAGGGATATAGCCAAATTACAGCATGCCAAGAAATGACAGAAATGTATTTTTTAATTAAAAAAAGGAGTATTAGTAAAGCTTCAACTTTAGAACAGGGTCTATGTCAGAGAAACACCACAGAGCAATGAAAACTGTGTTCTCCAAATTATTTAGAACGTAGGAGGATTTCCAAAATATTTATTTATTTTGAAAGTAGAATGATTAACATCTGCTGTCTATAAGCAACACTTTTGAAAAAAGGAATTTACTCCCTCCAATTTGCTCCATAGCCTACAGTTTCTTCTTTTATCTACCTCTTCCTCCTCCTCCTCCACCACCTTCTTTTCTTGTTGTTGTTGTTGTTCTTCTTCTTCCTCCTCCTCCTCTTCCTCTTCTTCTTCTTCTCCTTCTTCTTCCCCTCCCTTCCTTCCTTCCTTCTTTCCTTCCTTCCTTCCTCTTGTTCTTCTTCTTCCTCTTATTCTCCTTCTCCTTCTTCTTCTTCTCCTGCTTCTTCCTCTTCTCCTTCTTCTTCTTCCTCCTCCTCCTTCTTCTCTTCTTCTGCTTTTCTTCTCTTCTTCTTCTCTTTTCCTTCCTTCCTTCTTTCCTTCCTCTTCTTTCTCTTCTTCTTCCTCTTATTCTCCTTCTTCCTCCTCTTCTCCTTCTTCCTTCTCCTCCTCTTCCTCCTCGTTGTTCTTGTCATTGTTGTTTTTGTTGTTCTTCCTCCTCTTCTCTCTCCTCCTTCTCCTCCTCCTCCTCCTTCTTCTTCTCCTTGTCCTCCTCCTCTGTCTCTGTCTTCTTCTTTTTCTTCTTGTTCCTCTTCTTCTTCCTCTTATTCTCCTTCTTCTTCCTCCCCTTCTCCCTCTTCCGTCTTCTCCTCCTCCTCCCTCTTCCTTCTTCTCCTCCTCCTCCTTCTTTTTCTTTCTTCTTCTTTCTTCTTGTTCTTGTTCTTCTTCTTCTTCCTCCTCCTCATCCTTCTTCTCTTCTTCTGCTTTTCTTCTGTTCTTCTTCTTCTTTCCCTCCCTCCCTCCCTCCCTCCCTTCCTTCCTTCCTTCCTTCCTCTTCTTTCTCCTCTTCTTCTTCTTGTTCTTGTTCTTCTTCTTCCTCTTATTCTCCTTCTTCCTCCTCTTCTCCTTCTTCTTCCTTTTGCTCCTCTTCCTCCTTCTTGTTGTTGTTGTTCTTCTTCTTCTTCCTCTCATTCTTCTTCTTCCTCCTCTTCTCCCTCCTCCTTCTCCTCCTCCTCCTTCTTCTTCTCCTTCTCCTCCTCCTTTTTCTTCTCCTTCTCCTCCTCTTCCATCTCTGTCTTTGCCTTCTTCTTCTTGTTCCTCATCTTCTTCTTCCTCTTATTCTCCTTCTTCTTCCTTCCCTTCTCCCTCTTCCTTTTTCTCCTCCTCCTCCTTCTTTCTTCTTTCTTCCTTCTTCTTCTTTCTCTTCTTCTTCTTCTCCTTCTTCTTCTTCCTCCTCCTCCTCCCCTTCCCCCTTCCCCTCCTCCTCCTCCTTCTCCTCCTTCTCCTCCTCCTTCTTCTTCTTCCTCTTCTTCTTCTTCTTCTCCTTCTTCTTCTTCTCCTTCTTCTTCTTCTTCCTCCTCCTCCTCCCCTTCTCCCTTCCCCTCCTCCTCCTCCTCCTTCTCCTCCTCCTTCTTCTTCTTCTTCTTCTTCTTCTTC\t*\tSA:Z:chr2,4452399,-,2250H75M3D52M15I113M,60,25,157;chr9,34840411,-,2065H67M373H,33,2,57;chr10,75714830,-,1851H51M603H,0,0,51;chr2,4452298,-,1793H60M652H,60,3,45;chr6,15853372,-,1907H43M555H,60,0,43;\tMD:Z:259C382G365G125G337\tRG:Z:GATKSVContigAlignments\tNM:i:7\tAS:i:1433\tXS:i:0", canonicalChromosomes, refSeqDict);

    final CpxVariantInducingAssemblyContig.BasicInfo manuallyCalculatedBasicInfo =
            new CpxVariantInducingAssemblyContig.BasicInfo("chr2", false,
                    new SimpleInterval("chr2", 4450935, 4450935), new SimpleInterval("chr2", 4452641, 4452641));

    //////////
    final List<CpxVariantInducingAssemblyContig.Jump> manuallyCalculatedJumps =
            Arrays.asList(
                    new CpxVariantInducingAssemblyContig.Jump(new SimpleInterval("chr2", 4452399, 4452399), new SimpleInterval("chr9", 34840477, 34840477), StrandSwitch.NO_SWITCH, 118),
                    new CpxVariantInducingAssemblyContig.Jump(new SimpleInterval("chr9", 34840411, 34840411), new SimpleInterval("chr6", 15853414, 15853414), StrandSwitch.NO_SWITCH, 115),
                    new CpxVariantInducingAssemblyContig.Jump(new SimpleInterval("chr6", 15853372, 15853372), new SimpleInterval("chr2", 4452357, 4452357), StrandSwitch.NO_SWITCH, 54),
                    new CpxVariantInducingAssemblyContig.Jump(new SimpleInterval("chr2", 4452298, 4452298), new SimpleInterval("chr2", 4452406, 4452406), StrandSwitch.NO_SWITCH, 318)
            );
    final List<SimpleInterval> manuallyCalculatedEventPrimaryChromosomeSegmentingLocations =
            Arrays.asList(
                    new SimpleInterval("chr2", 4452298, 4452298),
                    new SimpleInterval("chr2", 4452357, 4452357),
                    new SimpleInterval("chr2", 4452399, 4452399),
                    new SimpleInterval("chr2", 4452406, 4452406)
            );
    final CpxVariantInducingAssemblyContig manuallyCalculatedCpxVariantInducingAssemblyContig =
            new CpxVariantInducingAssemblyContig(analysisReadyContig,
                    manuallyCalculatedBasicInfo,
                    manuallyCalculatedJumps,
                    manuallyCalculatedEventPrimaryChromosomeSegmentingLocations);

    //////////
    final SimpleInterval manuallyCalculatedAffectedRefRegion =
            new SimpleInterval("chr2", 4452298, 4452406);
    final List<SimpleInterval> manuallyCalculatedSegments =
            Arrays.asList(
                    new SimpleInterval("chr2", 4452298, 4452357),
                    new SimpleInterval("chr2", 4452357, 4452399),
                    new SimpleInterval("chr2", 4452399, 4452406));
    final List<String> manuallyCalculatedAltArrangementDescription =
            Arrays.asList("1", "2", "3", "UINS-318", "1", "UINS-54", "chr6:15853372-15853414", "UINS-115", "chr9:34840411-34840477", "UINS-118", "3");
    final byte[] manuallyCalculatedAltSeq = Arrays.copyOfRange(analysisReadyContig.getContigSequence(), 247, 1139);
    SequenceUtil.reverseComplement(manuallyCalculatedAltSeq);
    final CpxVariantCanonicalRepresentation manuallyCalculatedCpxVariantCanonicalRepresentation =
            new CpxVariantCanonicalRepresentation(
                    manuallyCalculatedAffectedRefRegion,
                    manuallyCalculatedSegments,
                    manuallyCalculatedAltArrangementDescription,
                    manuallyCalculatedAltSeq);

    //////////
    final byte[] dummyRefSequence = "T".getBytes();
    final VariantContextBuilder variantContextBuilder = new VariantContextBuilder()
            .chr("chr2").start(4452298).stop(4452406)
            .alleles(Arrays.asList(Allele.create(dummyRefSequence, true), Allele.create(SimpleSVType.createBracketedSymbAlleleString(CPX_SV_SYB_ALT_ALLELE_STR), false)))
            .id("CPX_chr2:4452298-4452406")
            .attribute(VCFConstants.END_KEY, 4452406)
            .attribute(SVTYPE, GATKSVVCFConstants.CPX_SV_SYB_ALT_ALLELE_STR)
            .attribute(SVLEN, 783)
            .attribute(SEQ_ALT_HAPLOTYPE, new String(manuallyCalculatedAltSeq));
    variantContextBuilder.attribute(CPX_EVENT_ALT_ARRANGEMENTS,
            String.join(VCFConstants.INFO_FIELD_ARRAY_SEPARATOR, manuallyCalculatedAltArrangementDescription));
    variantContextBuilder.attribute(CPX_SV_REF_SEGMENTS,
            String.join(VCFConstants.INFO_FIELD_ARRAY_SEPARATOR, manuallyCalculatedSegments.stream().map(SimpleInterval::toString).collect(Collectors.toList())));
    final VariantContext manuallyCalculatedVariantContext = variantContextBuilder.make();

    //////////
    final Set<SimpleInterval> manuallyCalculatedTwoBaseBoundaries = new HashSet<>();
    manuallyCalculatedTwoBaseBoundaries.add(new SimpleInterval("chr2", 4452399, 4452400));
    manuallyCalculatedTwoBaseBoundaries.add(new SimpleInterval("chr2", 4452640, 4452641));
    manuallyCalculatedTwoBaseBoundaries.add(new SimpleInterval("chr2", 4452298, 4452299));
    manuallyCalculatedTwoBaseBoundaries.add(new SimpleInterval("chr2", 4452356, 4452357));
    manuallyCalculatedTwoBaseBoundaries.add(new SimpleInterval("chr2", 4450935, 4450936));
    manuallyCalculatedTwoBaseBoundaries.add(new SimpleInterval("chr2", 4452405, 4452406));

    //////////
    return new PreprocessedAndAnalysisReadyContigWithExpectedResults(
            manuallyCalculatedCpxVariantInducingAssemblyContig,
            manuallyCalculatedCpxVariantCanonicalRepresentation,
            manuallyCalculatedVariantContext,
            dummyRefSequence,
            manuallyCalculatedTwoBaseBoundaries);
}
 
Example 16
Source File: Mutect2FilteringEngine.java    From gatk with BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
/**
 * Create a filtered variant and record statistics for the final pass of {@link FilterMutectCalls}
 */
public VariantContext applyFiltersAndAccumulateOutputStats(final VariantContext vc, final ReferenceContext referenceContext) {
    final VariantContextBuilder vcb = new VariantContextBuilder(vc).filters(new HashSet<>());

    final ErrorProbabilities errorProbabilities = new ErrorProbabilities(filters, vc, this, referenceContext);
    filteringOutputStats.recordCall(errorProbabilities, getThreshold() - EPSILON);

    // error probability must exceed threshold, and just in case threshold is bad, probabilities close to 1 must be filtered
    // and probabilities close to 0 must not be filtered
    double errorThreshold = Math.min(1 - EPSILON, Math.max(EPSILON, getThreshold()));

    Map<String, Double> siteFiltersWithErrorProb = new LinkedHashMap<>();

    // apply allele specific filters
    List<List<String>> alleleStatusByFilter =
            errorProbabilities.getProbabilitiesForAlleleFilters().entrySet().stream()
                    .filter(entry -> !entry.getValue().isEmpty())
                    .map(entry -> addFilterStrings(entry.getValue(), errorThreshold, entry.getKey().filterName())).collect(Collectors.toList());

    // for each allele, merge all allele specific filters
    List<List<String>> filtersByAllele = ErrorProbabilities.transpose(alleleStatusByFilter);
    List<List<String>> distinctFiltersByAllele = filtersByAllele.stream().map(this::getDistinctFiltersForAllele).collect(Collectors.toList());
    ListIterator<String> mergedFilterStringByAllele = distinctFiltersByAllele.stream().map(AnnotationUtils::encodeStringList).collect(Collectors.toList()).listIterator();

    List<String> orderedASFilterStrings = vc.getAlternateAlleles().stream().map(allele -> allele.isSymbolic() ?
            GATKVCFConstants.SITE_LEVEL_FILTERS : mergedFilterStringByAllele.next()).collect(Collectors.toList());
    String finalAttrString = AnnotationUtils.encodeAnyASListWithRawDelim(orderedASFilterStrings);

    vcb.putAttributes(Collections.singletonMap(GATKVCFConstants.AS_FILTER_STATUS_KEY, finalAttrString));


    // compute site-only filters
    // from allele specific filters
     alleleStatusByFilter.stream().forEachOrdered(alleleStatusForFilter -> {
        if (!alleleStatusForFilter.isEmpty() && alleleStatusForFilter.stream().distinct().count() == 1 && !alleleStatusForFilter.contains(GATKVCFConstants.SITE_LEVEL_FILTERS)) {
            siteFiltersWithErrorProb.put(alleleStatusForFilter.get(0), 1.0);
        }
    });


    // from variant filters
    errorProbabilities.getProbabilitiesForVariantFilters().entrySet().stream()
            .forEach(entry -> {
                entry.getKey().phredScaledPosteriorAnnotationName().ifPresent(annotation -> {
                    if (entry.getKey().requiredInfoAnnotations().stream().allMatch(vc::hasAttribute)) {
                        vcb.attribute(annotation, QualityUtils.errorProbToQual(entry.getValue()));
                    }
                });
                if (entry.getValue() > errorThreshold) {
                    siteFiltersWithErrorProb.put(entry.getKey().filterName(), entry.getValue());
                }

            });

    // if all alleles have been filtered out, but for different reasons, fail the site.
    // if the site is only ref and symbolic, no filters will be applied so don't fail
    if (siteFiltersWithErrorProb.isEmpty() && !distinctFiltersByAllele.stream().allMatch(List::isEmpty)) {
        List<List<String>> filtersNonSymbolicAlleles =  GATKVariantContextUtils.removeDataForSymbolicAltAlleles(vc, distinctFiltersByAllele);
        // if any allele passed, don't fail the site
        if (!filtersNonSymbolicAlleles.stream().anyMatch(filterList -> filterList.contains(GATKVCFConstants.SITE_LEVEL_FILTERS))) {
            // we know the allele level filters exceeded their threshold - so set this prob to 1
            siteFiltersWithErrorProb.put(GATKVCFConstants.FAIL, 1.0);
        }
    }

    // this code limits the number of filters specified for any variant to the highest probability filters
    // this will not change the status of whether a variant is actually filtered or not
    final double maxErrorProb = siteFiltersWithErrorProb.values().stream().mapToDouble(p->p).max().orElse(1);
    siteFiltersWithErrorProb.entrySet().stream().forEach(entry -> {
        if (entry.getValue() >= Math.min(maxErrorProb, MIN_REPORTABLE_ERROR_PROBABILITY)) {
            vcb.filter(entry.getKey());
        }
    });

    return vcb.make();
}
 
Example 17
Source File: CpxVariantCanonicalRepresentationUnitTest.java    From gatk with BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
@Test
public void testSpecialCtor() {
    final AlignedContig alignedContig = TestUtilsForAssemblyBasedSVDiscovery.fromPrimarySAMRecordString("asm000308:tig00000\t0\tchr1\t14491391\t60\t1276M530S\t*\t0\t0\tTTGCTGCACCCATCAATCCGTCGTCTACATTAGGTATTTCTCCTAATGCTATCCCTCCCCTAGTCCCCTACCCGCCGACAGGTCCCGGTGTGTGATATTCCCCTCCCTGTGTCCATGTTACTCTTTTGATATTACCAGGGACACCTGGATTTCTACTGATTTTAATGAGAATACCTTCTGTATTCACCATTAAATATGATGGTAGTTGCTGGTTTTAGCTCGATATTATTTGTCATGCTGATTAAGTGGACTTGCATTCCTAGCTTTCAAAGAGGTTTTCTTTCCTTTTTAATAAGGAGTGGGTGTTGCATGTTATCTAATACATTGTCAGTGTGTGGCTATTTATAAGTTCTGTGTGATTATACCATTTATCTATATAAATGTCCCCTTTGTTCTATTTAATAATGCTTTTCCCCCTTCTGAATTCCACTTTGAATTTGAATTCTACTTTGTCTGAAATAGATCCTGCCACCCCTGCTTTTAAAAAGAAAAAAATCTTTTTGCTTGTATTATTTAACTTTTTTGCCTATCCCTCCCTTTTTAATCTTTTCATACCATTGCTTTTCAGTGTCTCGAGCAGTAAGACATTTAACAATTATCAGCCCCATGCTTACTTTGTGCCAGACACTGGATTAAACAAAAATGGAAAAAGAGGATAGAATGTGCTGGAAGGGGTACATTCAAACCCAGTCTGAACTGGCCACTGCTGTGAGCAGGTTTGGGGACAGCAGTAGATCCTAGAAGGGCCTGACCAGCTGGGGAAACTGGCCAGGCTGTCCAGAGGTGACAAGAGGATTGTCACCCAGACTTGCCCAAGAAGAGTGAATCTGAGTCTTGGAGAGAACAGGAGTTTGGGTTCTTCTGGGCCCAGATGGCCTCAGGGCTCCCTGGAATTTGGGGACCCCACAGTTGGTCGCCACCATGAATTGAGGAGCCTTGCTTCTCTCCACACTGTCTTTTCCCTGCCTCCTCGTGGCTTCTGCTTCACTCATTCACTCATTCTGTCAGTGAATGATTCTTCAGCACCTGCCCTGCATAGGATGCCATTGTAGGTGCTGGGAAATCAACGGGAAGAAGATGGAAAACGAGACTTCCCTTATGAAGCTTCTGTTCTACAGAGGTGGGCAGACATGGCCAGAAAAAGCACAAGGCCATTTCCAATGGTGGAAGGGCCAGGACTGCTGCCCTTTCTGATAGCTTCTCTTTACACTTAGGAGAAAATTCAGGGCCCCATAATCCCTAGGCCCTACATAACCACACATGCACACACCACAAACCACACACACACACCACACACACCACACACCACACGCTACACACACCATGCACATACCACAAACCACACACACACCACACACACACCACACATTACACACACCACACAGACACCACACACCACACAAACATCACACACCACAGACACATCACACACCACACACACACCACACATACACAGCACACAGCTCACGCATACACAGCACACACATCACACATACACATACCACATACACACCACACACACACCACAAACCACATATACACAGCACACACATCACACAAACACATACCACATACACACCACACACCGTACATACACAGCATACACATTACACATACACACACGACACACCACACACAGACCACACACCACACAGATACAGCACAGAGACACTACACATACCACATACACACTACACACCACACACACACCACACATACACAGCGCACATACACACACCACACACACCACATACAAACCACATACCACATACACACCACACATATACCACACAGACACCATACATA\t*\tSA:Z:chr1,14492666,+,1684S71M51S,60,0,71;chr4,8687087,-,422S46M1338S,60,1,41;\tMD:Z:144T1131\tRG:Z:GATKSVContigAlignments\tNM:i:1\tAS:i:1271\tXS:i:70", true);
    final AssemblyContigWithFineTunedAlignments preprocessedTig = new AssemblyContigWithFineTunedAlignments(alignedContig);
    final CpxVariantInducingAssemblyContig analysisReadyContig = new CpxVariantInducingAssemblyContig(preprocessedTig, TestUtilsForAssemblyBasedSVDiscovery.bareBoneHg38SAMSeqDict);
    final SimpleInterval manuallyCalculatedAffectedRefRegion = new SimpleInterval("chr1", 14492666, 14492666);
    final byte[] manuallyCalculatedAltSeq = Arrays.copyOfRange(alignedContig.getContigSequence(), 1275, 1685);
    final List<String> manuallyCalculatedAltArrangements = Arrays.asList("1", "UINS-62", "-chr4:8687087-8687132", "UINS-300", "1");
    final List<SimpleInterval> manuallyCalculatedSegments = Collections.singletonList(manuallyCalculatedAffectedRefRegion);

    final CpxVariantCanonicalRepresentation cpxVariantCanonicalRepresentation = new CpxVariantCanonicalRepresentation(analysisReadyContig);

    Assert.assertEquals(cpxVariantCanonicalRepresentation.getAffectedRefRegion(),
                        manuallyCalculatedAffectedRefRegion);
    Assert.assertEquals(cpxVariantCanonicalRepresentation.getReferenceSegments(),
                        manuallyCalculatedSegments);
    Assert.assertEquals(cpxVariantCanonicalRepresentation.getEventDescriptions(),
                        manuallyCalculatedAltArrangements);
    Assert.assertTrue(Arrays.equals(cpxVariantCanonicalRepresentation.getAltSeq(),
                                    manuallyCalculatedAltSeq));

    final byte[] dummyRefSequence = "TCGA".getBytes();
    final VariantContextBuilder variantContextBuilder = new VariantContextBuilder()
            .chr("chr1").start(14492666).stop(14492666)
            .alleles(Arrays.asList(Allele.create(dummyRefSequence, true), Allele.create(SimpleSVType.createBracketedSymbAlleleString(CPX_SV_SYB_ALT_ALLELE_STR), false)))
            .id("CPX_chr1:14492666-14492666")
            .attribute(VCFConstants.END_KEY, 14492666)
            .attribute(SVTYPE, GATKSVVCFConstants.CPX_SV_SYB_ALT_ALLELE_STR)
            .attribute(SVLEN, 409)
            .attribute(SEQ_ALT_HAPLOTYPE, new String(manuallyCalculatedAltSeq));
    variantContextBuilder.attribute(CPX_EVENT_ALT_ARRANGEMENTS,
            String.join(VCFConstants.INFO_FIELD_ARRAY_SEPARATOR, manuallyCalculatedAltArrangements));
    variantContextBuilder.attribute(CPX_SV_REF_SEGMENTS,
            String.join(VCFConstants.INFO_FIELD_ARRAY_SEPARATOR, manuallyCalculatedSegments.stream().map(SimpleInterval::toString).collect(Collectors.toList())));
    final VariantContext manuallyCalculatedVariantContext = variantContextBuilder.make();

    VariantContextTestUtils.assertVariantContextsAreEqual(
            cpxVariantCanonicalRepresentation.toVariantContext(dummyRefSequence).make(),
            manuallyCalculatedVariantContext,
            Collections.emptyList(), Collections.emptyList());
}
 
Example 18
Source File: MergePedIntoVcf.java    From picard with MIT License 4 votes vote down vote up
/**
 * Writes out the VariantContext objects in the order presented to the supplied output file
 * in VCF format.
 */
private void writeVcf(final CloseableIterator<VariantContext> variants,
                      final File output,
                      final SAMSequenceDictionary dict,
                      final VCFHeader vcfHeader, ZCallPedFile zCallPedFile) {

    try(VariantContextWriter writer = new VariantContextWriterBuilder()
            .setOutputFile(output)
            .setReferenceDictionary(dict)
            .setOptions(VariantContextWriterBuilder.DEFAULT_OPTIONS)
            .build()) {
        writer.writeHeader(vcfHeader);
        while (variants.hasNext()) {
            final VariantContext context = variants.next();

            final VariantContextBuilder builder = new VariantContextBuilder(context);
            if (zCallThresholds.containsKey(context.getID())) {
                final String[] zThresh = zCallThresholds.get(context.getID());
                builder.attribute(InfiniumVcfFields.ZTHRESH_X, zThresh[0]);
                builder.attribute(InfiniumVcfFields.ZTHRESH_Y, zThresh[1]);
            }
            final Genotype originalGenotype = context.getGenotype(0);
            final Map<String, Object> newAttributes = originalGenotype.getExtendedAttributes();
            final VCFEncoder vcfEncoder = new VCFEncoder(vcfHeader, false, false);
            final Map<Allele, String> alleleMap = vcfEncoder.buildAlleleStrings(context);

            final String zCallAlleles = zCallPedFile.getAlleles(context.getID());
            if (zCallAlleles == null) {
                throw new PicardException("No zCall alleles found for snp " + context.getID());
            }
            final List<Allele> zCallPedFileAlleles = buildNewAllelesFromZCall(zCallAlleles, context.getAttributes());
            newAttributes.put(InfiniumVcfFields.GTA, alleleMap.get(originalGenotype.getAllele(0)) + UNPHASED + alleleMap.get(originalGenotype.getAllele(1)));
            newAttributes.put(InfiniumVcfFields.GTZ, alleleMap.get(zCallPedFileAlleles.get(0)) + UNPHASED + alleleMap.get(zCallPedFileAlleles.get(1)));

            final Genotype newGenotype = GenotypeBuilder.create(originalGenotype.getSampleName(), zCallPedFileAlleles,
                    newAttributes);
            builder.genotypes(newGenotype);
            logger.record("0", 0);
            // AC, AF, and AN are recalculated here
            VariantContextUtils.calculateChromosomeCounts(builder, false);
            final VariantContext newContext = builder.make();
            writer.add(newContext);
        }
    }
}
 
Example 19
Source File: CombineGenotypingArrayVcfs.java    From picard with MIT License 4 votes vote down vote up
/**
 * 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();
}
 
Example 20
Source File: GencodeFuncotationFactoryUnitTest.java    From gatk with BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
@Test (dataProvider = "provideMuc16SnpDataForGetVariantClassification")
void testGetVariantClassificationForCodingRegions(final int chromosomeNumber,
                                  final int start,
                                  final int end,
                                  final GencodeFuncotation.VariantType variantType,
                                  final String ref,
                                  final String alt,
                                  final GencodeFuncotation.VariantClassification expectedVariantClassification) {

    // This test can only deal with variants in coding regions.
    // So we ignore any tests that are expected outside of coding regions.
    // i.e. expectedVariantClassification is one of:
    //     { INTRON, FIVE_PRIME_UTR, THREE_PRIME_UTR, IGR, FIVE_PRIME_FLANK, DE_NOVO_START_IN_FRAME, DE_NOVO_START_OUT_FRAME, RNA, LINCRNA }
    // We test these cases in another unit test.
    if ((expectedVariantClassification == GencodeFuncotation.VariantClassification.INTRON) ||
        (expectedVariantClassification == GencodeFuncotation.VariantClassification.FIVE_PRIME_UTR) ||
        (expectedVariantClassification == GencodeFuncotation.VariantClassification.THREE_PRIME_UTR) ||
        (expectedVariantClassification == GencodeFuncotation.VariantClassification.IGR) ||
        (expectedVariantClassification == GencodeFuncotation.VariantClassification.FIVE_PRIME_FLANK) ||
        (expectedVariantClassification == GencodeFuncotation.VariantClassification.DE_NOVO_START_IN_FRAME) ||
        (expectedVariantClassification == GencodeFuncotation.VariantClassification.DE_NOVO_START_OUT_FRAME) ||
        (expectedVariantClassification == GencodeFuncotation.VariantClassification.RNA) ||
        (expectedVariantClassification == GencodeFuncotation.VariantClassification.LINCRNA) )
    {
        return;
    }

    final String contig = "chr" + Integer.toString(chromosomeNumber);
    final SimpleInterval variantInterval = new SimpleInterval( contig, start, end );

    final Allele refAllele = Allele.create(ref, true);
    final Allele altAllele = Allele.create(alt);

    final VariantContextBuilder variantContextBuilder = new VariantContextBuilder(
            FuncotatorReferenceTestUtils.retrieveHg19Chr19Ref(),
            contig,
            start,
            end,
            Arrays.asList(refAllele, altAllele)
    );
    final VariantContext variantContext = variantContextBuilder.make();

    // Get our gene feature iterator:
    final CloseableTribbleIterator<GencodeGtfFeature> gtfFeatureIterator;
    try {
        gtfFeatureIterator = gencodeHg19FeatureReader.query(contig, start, end);
    }
    catch (final IOException ex) {
        throw new GATKException("Could not finish the test!", ex);
    }

    // Get the gene.
    // We know the first gene is the right one - the gene in question is the MUC16 gene:
    final GencodeGtfGeneFeature             gene = (GencodeGtfGeneFeature) gtfFeatureIterator.next();
    final GencodeGtfTranscriptFeature transcript = getMuc16Transcript(gene);
    final GencodeGtfExonFeature             exon = getExonForVariant( gene, variantInterval );

    final ReferenceContext referenceContext = new ReferenceContext( refDataSourceHg19Ch19, variantInterval );

    final List<? extends Locatable> exonPositionList = GencodeFuncotationFactory.getSortedCdsAndStartStopPositions(transcript);

    final ReferenceDataSource muc16TranscriptDataSource = ReferenceDataSource.of(new File(FuncotatorTestConstants.GENCODE_DATA_SOURCE_FASTA_PATH_HG19).toPath());
    final Map<String, GencodeFuncotationFactory.MappedTranscriptIdInfo> muc16TranscriptIdMap = GencodeFuncotationFactory. createTranscriptIdMap(muc16TranscriptDataSource);

    final SequenceComparison seqComp =
            GencodeFuncotationFactory.createSequenceComparison(
                    variantContext,
                    altAllele,
                    referenceContext,
                    transcript,
                    exonPositionList,
                    muc16TranscriptIdMap,
                    muc16TranscriptDataSource,
                    true);

    final GencodeFuncotation.VariantClassification varClass = GencodeFuncotationFactory.createVariantClassification(
            variantContext,
            altAllele,
            variantType,
            exon,
            transcript.getExons().size(),
            seqComp
    );

    Assert.assertEquals( varClass, expectedVariantClassification );
}