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

The following examples show how to use htsjdk.variant.variantcontext.VariantContext#hasAttribute() . 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: EvaluateInfoFieldConcordance.java    From gatk with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
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()));
        }
    }
}
 
Example 2
Source File: AS_StrandBiasTest.java    From gatk with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
/**
 * Parses the raw data stings of combined contingency matrix data and calls client methods calculateReducedData(myData)
 * implementation to generate double digest of provided allele information which is stored in '|' delineated lists.
 *
 * @param vc -- contains the final set of alleles, possibly subset by GenotypeGVCFs
 * @param originalVC -- used to get all the alleles for all gVCFs
 * @return
 */
@Override
public  Map<String, Object> finalizeRawData(final VariantContext vc, final VariantContext originalVC) {
    if (!vc.hasAttribute(getPrimaryRawKey())) {
        return new HashMap<>();
    }
    String rawContingencyTableData = vc.getAttributeAsString(getPrimaryRawKey(),null);
    if (rawContingencyTableData == null) {
        return new HashMap<>();
    }
    AlleleSpecificAnnotationData<List<Integer>> myData = new AlleleSpecificAnnotationData<>(originalVC.getAlleles(), rawContingencyTableData);
    parseRawDataString(myData);

    Map<Allele, Double> perAltRankSumResults = calculateReducedData(myData);

    String annotationString = makeReducedAnnotationString(vc, perAltRankSumResults);
    String rawAnnotationsString = StrandBiasUtils.makeRawAnnotationString(vc.getAlleles(), myData.getAttributeMap());
    Map<String, Object> returnMap = new HashMap<>();
    returnMap.put(getKeyNames().get(0), annotationString);
    returnMap.put(getPrimaryRawKey(), rawAnnotationsString);  //this is in case raw annotations are requested
    return returnMap;
}
 
Example 3
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 4
Source File: FilterVariantTranches.java    From gatk with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
@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());
}
 
Example 5
Source File: VariantSummary.java    From gatk with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
public void update2(VariantContext eval, VariantContext comp, ReferenceContext referenceContext, ReadsContext readsContext, FeatureContext featureContext) {
    if ( eval == null || (getWalker().ignoreAC0Sites() && eval.isMonomorphicInSamples()) )
        return;

    final Type type = getType(eval);
    if ( type == null )
        return;

    TypeSampleMap titvTable = null;

    // update DP, if possible
    if ( eval.hasAttribute(VCFConstants.DEPTH_KEY) )
        depthPerSample.inc(type, ALL);

    // update counts
    allVariantCounts.inc(type, ALL);

    // type specific calculations
    if ( type == Type.SNP && eval.isBiallelic() ) {
        titvTable = GATKVariantContextUtils.isTransition(eval) ? transitionsPerSample : transversionsPerSample;
        titvTable.inc(type, ALL);
    }

    // novelty calculation
    if ( comp != null || (type == Type.CNV && overlapsKnownCNV(eval, featureContext)))
        knownVariantCounts.inc(type, ALL);

    // per sample metrics
    for (final Genotype g : eval.getGenotypes()) {
        if ( ! g.isNoCall() && ! g.isHomRef() ) {
            countsPerSample.inc(type, g.getSampleName());

            // update transition / transversion ratio
            if ( titvTable != null ) titvTable.inc(type, g.getSampleName());

            if ( g.hasDP() )
                depthPerSample.inc(type, g.getSampleName());
        }
    }
}
 
Example 6
Source File: GetPileupSummaries.java    From gatk with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
private boolean alleleFrequencyInRange(final VariantContext vc) {
    if (!vc.hasAttribute(VCFConstants.ALLELE_FREQUENCY_KEY)) {
        if (!sawVariantsWithoutAlleleFrequency) {
            logger.warn(String.format("Variant context at %s:%d lacks allele frequency (AF) field.", vc.getContig(), vc.getStart()));
            sawVariantsWithoutAlleleFrequency = true;
        }
        return false;
    } else {
        sawVariantsWithAlleleFrequency = true;
        final double alleleFrequency = vc.getAttributeAsDouble(VCFConstants.ALLELE_FREQUENCY_KEY, -1.0);
        return minPopulationAlleleFrequency < alleleFrequency && alleleFrequency < maxPopulationAlleleFrequency;
    }
}
 
Example 7
Source File: Mutect2FilteringEngine.java    From gatk-protected with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
private static void applyInsufficientEvidenceFilter(final M2FiltersArgumentCollection MTFAC, final VariantContext vc, final Collection<String> filters) {
    if (vc.hasAttribute(GATKVCFConstants.TUMOR_LOD_KEY)) {
        final double[] tumorLods = getArrayAttribute(vc, GATKVCFConstants.TUMOR_LOD_KEY);

        if (MathUtils.arrayMax(tumorLods) < MTFAC.TUMOR_LOD_THRESHOLD) {
            filters.add(GATKVCFConstants.TUMOR_LOD_FILTER_NAME);
        }
    }
}
 
Example 8
Source File: GetPileupSummaries.java    From gatk-protected with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
private boolean alleleFrequencyInRange(final VariantContext vc) {
    if (!vc.hasAttribute(VCFConstants.ALLELE_FREQUENCY_KEY)) {
        return false;
    } else {
        final double alleleFrequency = vc.getAttributeAsDouble(VCFConstants.ALLELE_FREQUENCY_KEY, -1.0);
        return minPopulationAlleleFrequency < alleleFrequency && alleleFrequency < maxPopulationAlleleFrequency;
    }
}
 
Example 9
Source File: SvDiscoverFromLocalAssemblyContigAlignmentsSpark.java    From gatk with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
@Override
public boolean test(final VariantContext variantContext) {
    if ( !variantContext.hasAttribute(GATKSVVCFConstants.CONTIG_NAMES) )
        return true;

    final int alnLength = variantContext.getAttributeAsInt(attributeKey, 0);
    return alnLength >= threshold;
}
 
Example 10
Source File: CosmicAnnotationFactory.java    From hmftools with GNU General Public License v3.0 5 votes vote down vote up
@NotNull
public static List<CosmicAnnotation> fromContext(@NotNull final VariantContext context) {
    if (context.hasAttribute(COSMIC_IDENTIFIER)) {
        return context.getAttributeAsStringList(COSMIC_IDENTIFIER, "")
                .stream()
                .map(x -> x.trim().split(FIELD_SEPARATOR))
                .map(CosmicAnnotationFactory::fromParts)
                .collect(Collectors.toList());
    }
    return Collections.emptyList();
}
 
Example 11
Source File: StructuralVariantFactory.java    From hmftools with GNU General Public License v3.0 5 votes vote down vote up
public static Double unadjustedAllelicFrequency(@NotNull VariantContext context) {
    if (context.hasAttribute(TAF)) {
        return context.getAttributeAsDoubleList(TAF, 0.0).get(0);
    }

    if (context.hasAttribute(BPI_AF)) {
        return context.getAttributeAsDoubleList(BPI_AF, 0.0).get(0);
    }

    return null;
}
 
Example 12
Source File: SVContextUnitTest.java    From gatk with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
/**
 * Tests {@link SVContext#getStructuralVariantLength()}.
 * @param vc input variant context.
 */
@Test(dataProvider="validVariantContexts", dependsOnMethods = {"testCreate"}, groups = "sv")
public void testLength(final VariantContext vc, @SuppressWarnings("unused") final ReferenceMultiSparkSource reference) {
    final SVContext svc = SVContext.of(vc);
    final int length = svc.getStructuralVariantLength();
    if (!vc.hasAttribute(GATKSVVCFConstants.SVLEN)) {
        Assert.assertEquals(length, -1);
    } else {
        Assert.assertEquals(length, Math.abs(vc.getAttributeAsInt(GATKSVVCFConstants.SVLEN, 0)));
    }
}
 
Example 13
Source File: ApplyVQSR.java    From gatk with BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
/**
 * Calculate the allele-specific filter status of vc
 * @param vc
 * @param recals
 * @param builder   is modified by adding attributes
 * @return a String with the filter status for this site
 */
private String doAlleleSpecificFiltering(final VariantContext vc, final List<VariantContext> recals, final VariantContextBuilder builder) {
    double bestLod = VariantRecalibratorEngine.MIN_ACCEPTABLE_LOD_SCORE;
    final List<String> culpritStrings = new ArrayList<>();
    final List<String> lodStrings = new ArrayList<>();
    final List<String> AS_filterStrings = new ArrayList<>();

    String[] prevCulpritList = null;
    String[] prevLodList = null;
    String[] prevASfiltersList = null;

    //get VQSR annotations from previous run of ApplyRecalibration, if applicable
    if(foundINDELTranches || foundSNPTranches) {
        final String prevCulprits = vc.getAttributeAsString(GATKVCFConstants.AS_CULPRIT_KEY,"");
        prevCulpritList = prevCulprits.isEmpty()? new String[0] : prevCulprits.split(listPrintSeparator);
        final String prevLodString = vc.getAttributeAsString(GATKVCFConstants.AS_VQS_LOD_KEY,"");
        prevLodList = prevLodString.isEmpty()? new String[0] : prevLodString.split(listPrintSeparator);
        final String prevASfilters = vc.getAttributeAsString(GATKVCFConstants.AS_FILTER_STATUS_KEY,"");
        prevASfiltersList = prevASfilters.isEmpty()? new String[0] : prevASfilters.split(listPrintSeparator);
    }

    //for each allele in the current VariantContext
    for (int altIndex = 0; altIndex < vc.getNAlleles()-1; altIndex++) {
        final Allele allele = vc.getAlternateAllele(altIndex);

        //if the current allele is not part of this recalibration mode, add its annotations to the list and go to the next allele
        if (!VariantDataManager.checkVariationClass(vc, allele, MODE)) {
            updateAnnotationsWithoutRecalibrating(altIndex, prevCulpritList, prevLodList, prevASfiltersList, culpritStrings, lodStrings, AS_filterStrings);
            continue;
        }

        //if the current allele does need to have recalibration applied...

        //initialize allele-specific VQSR annotation data with values for spanning deletion
        String alleleLodString = emptyFloatValue;
        String alleleFilterString = emptyStringValue;
        String alleleCulpritString = emptyStringValue;

        //if it's not a spanning deletion, replace those allele strings with the real values
        if (!GATKVCFConstants.isSpanningDeletion(allele)) {
            VariantContext recalDatum = getMatchingRecalVC(vc, recals, allele);
            if (recalDatum == null) {
                throw new UserException("Encountered input allele which isn't found in the input recal file. Please make sure VariantRecalibrator and ApplyRecalibration were run on the same set of input variants with flag -AS. First seen at: " + vc);
            }

            //compare VQSLODs for all alleles in the current mode for filtering later
            final double lod = recalDatum.getAttributeAsDouble(GATKVCFConstants.VQS_LOD_KEY, VariantRecalibratorEngine.MIN_ACCEPTABLE_LOD_SCORE);
            if (lod > bestLod)
                bestLod = lod;

            alleleLodString = String.format("%.4f", lod);
            alleleFilterString = generateFilterString(lod);
            alleleCulpritString = recalDatum.getAttributeAsString(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);
            }
        }

        //append per-allele VQSR annotations
        lodStrings.add(alleleLodString);
        AS_filterStrings.add(alleleFilterString);
        culpritStrings.add(alleleCulpritString);
    }

    // Annotate the new record with its VQSLOD, AS_FilterStatus, and the worst performing annotation
    if(!AS_filterStrings.isEmpty() )
        builder.attribute(GATKVCFConstants.AS_FILTER_STATUS_KEY, AnnotationUtils.encodeStringList(AS_filterStrings));
    if(!lodStrings.isEmpty())
        builder.attribute(GATKVCFConstants.AS_VQS_LOD_KEY, AnnotationUtils.encodeStringList(lodStrings));
    if(!culpritStrings.isEmpty())
        builder.attribute(GATKVCFConstants.AS_CULPRIT_KEY, AnnotationUtils.encodeStringList(culpritStrings));

    return generateFilterStringFromAlleles(vc, bestLod);
}
 
Example 14
Source File: SelectVariants.java    From gatk with BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
private void addAnnotations(final VariantContextBuilder builder, final VariantContext originalVC, final Set<String> selectedSampleNames) {
    if (fullyDecode) {
        return; // TODO -- annotations are broken with fully decoded data
    }

    if (keepOriginalChrCounts) {
        final int[] indexOfOriginalAlleleForNewAllele;
        final List<Allele> newAlleles = builder.getAlleles();
        final int numOriginalAlleles = originalVC.getNAlleles();

        // if the alleles already match up, we can just copy the previous list of counts
        if (numOriginalAlleles == newAlleles.size()) {
            indexOfOriginalAlleleForNewAllele = null;
        }
        // otherwise we need to parse them and select out the correct ones
        else {
            indexOfOriginalAlleleForNewAllele = new int[newAlleles.size() - 1];
            Arrays.fill(indexOfOriginalAlleleForNewAllele, -1);

            // note that we don't care about the reference allele at position 0
            for (int newI = 1; newI < newAlleles.size(); newI++) {
                final Allele newAlt = newAlleles.get(newI);
                for (int oldI = 0; oldI < numOriginalAlleles - 1; oldI++) {
                    if (newAlt.equals(originalVC.getAlternateAllele(oldI), false)) {
                        indexOfOriginalAlleleForNewAllele[newI - 1] = oldI;
                        break;
                    }
                }
            }
        }

        if (originalVC.hasAttribute(VCFConstants.ALLELE_COUNT_KEY)) {
            builder.attribute(GATKVCFConstants.ORIGINAL_AC_KEY,
                    getReorderedAttributes(originalVC.getAttribute(VCFConstants.ALLELE_COUNT_KEY), indexOfOriginalAlleleForNewAllele));
        }
        if (originalVC.hasAttribute(VCFConstants.ALLELE_FREQUENCY_KEY)) {
            builder.attribute(GATKVCFConstants.ORIGINAL_AF_KEY,
                    getReorderedAttributes(originalVC.getAttribute(VCFConstants.ALLELE_FREQUENCY_KEY), indexOfOriginalAlleleForNewAllele));
        }
        if (originalVC.hasAttribute(VCFConstants.ALLELE_NUMBER_KEY)) {
            builder.attribute(GATKVCFConstants.ORIGINAL_AN_KEY, originalVC.getAttribute(VCFConstants.ALLELE_NUMBER_KEY));
        }
    }

    VariantContextUtils.calculateChromosomeCounts(builder, false);

    if (keepOriginalDepth && originalVC.hasAttribute(VCFConstants.DEPTH_KEY)) {
        builder.attribute(GATKVCFConstants.ORIGINAL_DP_KEY, originalVC.getAttribute(VCFConstants.DEPTH_KEY));
    }

    boolean sawDP = false;
    int depth = 0;
    for (final String sample : selectedSampleNames ) {
        final Genotype g = originalVC.getGenotype(sample);
        if (!g.isFiltered()) {
            if (g.hasDP()) {
                depth += g.getDP();
                sawDP = true;
            }
        }
    }

    if (sawDP) {
        builder.attribute(VCFConstants.DEPTH_KEY, depth);
    }
}
 
Example 15
Source File: Mutect2FilteringEngine.java    From gatk-protected with BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
private static void applyPanelOfNormalsFilter(final M2FiltersArgumentCollection MTFAC, final VariantContext vc, final Collection<String> filters) {
    final boolean siteInPoN = vc.hasAttribute(SomaticGenotypingEngine.IN_PON_VCF_ATTRIBUTE);
    if (siteInPoN) {
        filters.add(GATKVCFConstants.PON_FILTER_NAME);
    }
}
 
Example 16
Source File: HaplotypeMap.java    From picard with MIT License 4 votes vote down vote up
/**
 * Constructs a HaplotypeMap from the provided file.
 */

private void fromVcf(final File file) {
    try ( final VCFFileReader reader = new VCFFileReader(file, false)) {

        final SAMSequenceDictionary dict = reader.getFileHeader().getSequenceDictionary();

        if (dict == null || dict.getSequences().isEmpty()) {
            throw new IllegalStateException("Haplotype map VCF file must contain header: " + file.getAbsolutePath());
        }

        initialize(new SAMFileHeader(dict));

        final Map<String, HaplotypeBlock> anchorToHaplotype = new HashMap<>();

        for (final VariantContext vc : reader) {

            if (vc.getNSamples() > 1) {
                throw new IllegalStateException("Haplotype map VCF file must contain at most one sample: " + file.getAbsolutePath());
            }

            final Genotype gc = vc.getGenotype(0); // may be null
            final boolean hasGc = gc != null;

            if (vc.getAlternateAlleles().size() != 1) {
                throw new IllegalStateException("Haplotype map VCF file must contain exactly one alternate allele per site: " + vc.toString());
            }

            if (!vc.isSNP()) {
                throw new IllegalStateException("Haplotype map VCF file must contain only SNPs: " + vc.toString());
            }

            if (!vc.hasAttribute(VCFConstants.ALLELE_FREQUENCY_KEY)) {
                throw new IllegalStateException("Haplotype map VCF Variants must have an '"+ VCFConstants.ALLELE_FREQUENCY_KEY + "' INFO field: " + vc.toString());
            }


            if (hasGc && gc.isPhased() && !gc.hasExtendedAttribute(VCFConstants.PHASE_SET_KEY)) {
                throw new IllegalStateException("Haplotype map VCF Variants' genotypes that are phased must have a PhaseSet (" + VCFConstants.PHASE_SET_KEY+")" + vc.toString());
            }

            if (hasGc && gc.isPhased() && !gc.isHet()) {
                throw new IllegalStateException("Haplotype map VCF Variants' genotypes that are phased must be HET" + vc.toString());
            }

            // Then parse them out
            final String chrom = vc.getContig();
            final int pos = vc.getStart();
            final String name = vc.getID();

            final byte ref = vc.getReference().getBases()[0];
            final byte var = vc.getAlternateAllele(0).getBases()[0];

            final double temp_maf = vc.getAttributeAsDouble(VCFConstants.ALLELE_FREQUENCY_KEY, 0D);
            final boolean swapped = hasGc && !gc.getAllele(0).equals(vc.getReference());

            final byte major, minor;
            final double maf;

            if (swapped) {
                major = var;
                minor = ref;
                maf = 1 - temp_maf;
            } else {
                major = ref;
                minor = var;
                maf = temp_maf;
            }

            final String anchor = anchorFromVc(vc);

            // If it's the anchor snp, start the haplotype
            if (!anchorToHaplotype.containsKey(anchor)) {
                final HaplotypeBlock newBlock = new HaplotypeBlock(maf);
                anchorToHaplotype.put(anchor, newBlock);
            }
            final HaplotypeBlock block = anchorToHaplotype.get(anchor);
            block.addSnp(new Snp(name, chrom, pos, major, minor, maf, null));
        }

        // And add them all now that they are all ready.
        fromHaplotypes(anchorToHaplotype.values());
    }
}
 
Example 17
Source File: Degeneracy.java    From gatk with BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up
public List<Object> getRelevantStates(ReferenceContext referenceContext, ReadsContext readsContext, FeatureContext featureContext, VariantContext comp, String compName, VariantContext eval, String evalName, String sampleName, String FamilyName) {
    ArrayList<Object> relevantStates = new ArrayList<Object>();

    relevantStates.add("all");

    if (eval != null && eval.isVariant()) {
        String type = null;
        String aa = null;
        Integer frame = null;

        if (eval.hasAttribute("refseq.functionalClass")) {
            aa = eval.getAttributeAsString("refseq.variantAA", null);
            frame = eval.getAttributeAsInt("refseq.frame", 0);
        } else if (eval.hasAttribute("refseq.functionalClass_1")) {
            int annotationId = 1;
            String key;

            do {
                key = String.format("refseq.functionalClass_%d", annotationId);

                String newtype = eval.getAttributeAsString(key, null);

                if ( newtype != null &&
                        ( type == null ||
                                ( type.equals("silent") && !newtype.equals("silent") ) ||
                                ( type.equals("missense") && newtype.equals("nonsense") ) )
                        ) {
                    type = newtype;

                    String aakey = String.format("refseq.variantAA_%d", annotationId);
                    aa = eval.getAttributeAsString(aakey, null);

                    if (aa != null) {
                        String framekey = String.format("refseq.frame_%d", annotationId);

                        if (eval.hasAttribute(framekey)) {
                            frame = eval.getAttributeAsInt(framekey, 0);
                        }
                    }
                }

                annotationId++;
            } while (eval.hasAttribute(key));
        }

        if (aa != null && degeneracies.containsKey(aa) && frame != null) {
            relevantStates.add(degeneracies.get(aa).get(frame));
        }
    }

    return relevantStates;
}
 
Example 18
Source File: LiftoverUtils.java    From picard with MIT License 4 votes vote down vote up
/**
 * 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();
}