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

/**
 * @param simple        simple variant derived from pair-iteration logic that is to be checked
 * @param attributes    source CPX variant attributes
 */
@VisibleForTesting
static boolean isConsistentWithCPX(final VariantContext simple,
                                   final RelevantAttributes attributes) {

    final String typeString = simple.getAttributeAsString(SVTYPE, "");

    if ( typeString.equals(SimpleSVType.SupportedType.DEL.name()) ) {
        final List<SimpleInterval> refSegments = attributes.referenceSegments;
        final List<String> altArrangement = attributes.altArrangements;

        final Tuple3<Set<SimpleInterval>, Set<Integer>, List<Integer>> missingAndPresentAndInvertedSegments =
                getMissingAndPresentAndInvertedSegments(refSegments, altArrangement);
        final Set<SimpleInterval> missingSegments = missingAndPresentAndInvertedSegments._1();

        return deletionConsistencyCheck(simple, missingSegments);
    } else if ( typeString.equals(SimpleSVType.SupportedType.INV.name()) ) {
        return false;
    } else
        return true;
}
 

/**
 * Convert short, i.e. duplicated range is < 50 bp, duplication call to insertion call.
 */
@VisibleForTesting
static VariantContext postProcessConvertShortDupToIns(final VariantContext simple) {
    final String type = simple.getAttributeAsString(SVTYPE, "");
    if ( type.equals(SimpleSVType.SupportedType.DUP.name()) ) {
        final SimpleInterval duplicatedRegion = new SimpleInterval(simple.getAttributeAsString(DUP_REPEAT_UNIT_REF_SPAN, ""));
        if (duplicatedRegion.size() > EVENT_SIZE_THRESHOLD) {
            return simple;
        } else {
            return new VariantContextBuilder(simple)
                    .alleles(Arrays.asList(simple.getReference(), altSymbAlleleIns))
                    .rmAttribute(SVTYPE)
                    .attribute(SVTYPE, SimpleSVType.SupportedType.INS.name())
                    .make();
        }
    } else
        return simple;
}
 

/**
 * Apply filters (that implements {@link StructuralVariantFilter}) given list of variants,
 * and write the variants to a single VCF file.
 * @param outputPrefixWithSampleName    prefix with sample name
 * @param variants                      variants to which filters are to be applied and written to file
 * @param svDiscoveryInputMetaData      metadata for use in filtering and file output
 */
public static void filterAndWriteMergedVCF(final String outputPrefixWithSampleName,
                                           final List<VariantContext> variants,
                                           final SvDiscoveryInputMetaData svDiscoveryInputMetaData) {
    final List<VariantContext> variantsWithFilterApplied = new ArrayList<>(variants.size());
    final List<StructuralVariantFilter> filters = Arrays.asList(
            new SVMappingQualityFilter(svDiscoveryInputMetaData.getDiscoverStageArgs().minMQ),
            new SVAlignmentLengthFilter(svDiscoveryInputMetaData.getDiscoverStageArgs().minAlignLength));
    for (final VariantContext variant : variants) {
        String svType = variant.getAttributeAsString(GATKSVVCFConstants.SVTYPE, "");
        if (svType.equals(GATKSVVCFConstants.SYMB_ALT_ALLELE_DEL) || svType.equals(GATKSVVCFConstants.SYMB_ALT_ALLELE_INS) || svType.equals(GATKSVVCFConstants.SYMB_ALT_ALLELE_DUP)) {
            if (Math.abs(variant.getAttributeAsInt(GATKSVVCFConstants.SVLEN, 0)) < StructuralVariationDiscoveryArgumentCollection.STRUCTURAL_VARIANT_SIZE_LOWER_BOUND )
                continue;
        }
        variantsWithFilterApplied.add(applyFilters(variant, filters));
    }

    final String out = outputPrefixWithSampleName + MERGED_VCF_FILE_NAME;
    SVVCFWriter.writeVCF(variantsWithFilterApplied, out,
            svDiscoveryInputMetaData.getReferenceData().getReferenceSequenceDictionaryBroadcast().getValue(),
            svDiscoveryInputMetaData.getDefaultToolVCFHeaderLines(),
            svDiscoveryInputMetaData.getToolLogger());
}
 

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

private void validateFuncotationsOnVcf(final Iterable<VariantContext> vcfIterable, final String[] funcotationFieldNames) {
    for (final VariantContext vc : vcfIterable ) {
        final String funcotation = vc.getAttributeAsString(VcfOutputRenderer.FUNCOTATOR_VCF_FIELD_NAME, "");

        Assert.assertNotEquals(funcotation, "");

        final String rawFuncotations = funcotation.substring(1,funcotation.length()-1);

        Assert.assertEquals(StringUtils.countMatches(rawFuncotations, VcfOutputRenderer.FIELD_DELIMITER), funcotationFieldNames.length - 1);

        // This is here to make sure we can create the FuncotationMap object without exploding.
        // It serves as a secondary check.
        final FuncotationMap funkyMap = FuncotationMap.createAsAllTableFuncotationsFromVcf(FuncotationMap.NO_TRANSCRIPT_AVAILABLE_KEY, funcotationFieldNames,
                funcotation, vc.getAlternateAllele(0), "VCF");
    }
}
 

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

/**
 * Parses the raw data histograms generated for each allele and outputs the median value for each histogram
 * as a representative double value.
 *
 * @param vc -- contains the final set of alleles, possibly subset by GenotypeGVCFs
 * @param originalVC -- used to get all the alleles for all gVCFs
 * @return the finalized key and value as well as the raw key and value
 */
public  Map<String, Object> finalizeRawData(final VariantContext vc, final VariantContext originalVC) {
    if (!vc.hasAttribute(getPrimaryRawKey())) {
        return new HashMap<>();
    }
    final String rawRankSumData = vc.getAttributeAsString(getPrimaryRawKey(),null);
    if (rawRankSumData == null) {
        return new HashMap<>();
    }
    final Map<String,Object> annotations = new HashMap<>();
    final AlleleSpecificAnnotationData<Histogram> myData = new AlleleSpecificAnnotationData<>(originalVC.getAlleles(), rawRankSumData);
    parseRawDataString(myData);

    final Map<Allele, Double> perAltRankSumResults = calculateReducedData(myData.getAttributeMap(), myData.getRefAllele());
    //shortcut for no ref values
    if (perAltRankSumResults.isEmpty()) {
        return annotations;
    }
    final String annotationString = makeReducedAnnotationString(vc, perAltRankSumResults);
    annotations.put(getKeyNames().get(0), annotationString);
    annotations.put(getPrimaryRawKey(), makeCombinedAnnotationString(vc.getAlleles(), myData.getAttributeMap()));
    return annotations;
}
 

/**
 * Takes combined raw annotation data sums, and calculates per allele the average root mean squared from the raw data
 * using expected Allele Depth counts data.
 *
 * Will output delineated doubles in the format: sqrt(TotalAllele1RMS/Allele1Depth)|sqrt(TotalAllele1RMS/Allele1Depth)|...
 */
@Override
@SuppressWarnings({"unchecked", "rawtypes"})//FIXME generics here blow up
public Map<String, Object> finalizeRawData(final VariantContext vc, final VariantContext originalVC) {
    if (!vc.hasAttribute(getPrimaryRawKey())) {
        return new HashMap<>();
    }
    final String rawMQdata = vc.getAttributeAsString(getPrimaryRawKey(),null);
    if (rawMQdata == null) {
        return new HashMap<>();
    }

    final Map<String,Object> annotations = new HashMap<>();
    final ReducibleAnnotationData myData = new AlleleSpecificAnnotationData<Double>(originalVC.getAlleles(), rawMQdata);
    parseRawDataString(myData);

    final String annotationString = makeFinalizedAnnotationString(vc, myData.getAttributeMap());
    annotations.put(getKeyNames().get(0), annotationString);
    annotations.put(getPrimaryRawKey(), makeRawAnnotationString(vc.getAlleles(), myData.getAttributeMap()));
    return annotations;
}
 

private static boolean isReverseCandidate(@NotNull final VariantContext context) {
    final boolean altMatch =
            context.getAlternateAlleles().stream().anyMatch(x -> x.getBaseString().equals("C") || x.getBaseString().equals("A"));

    final String triContext = context.getAttributeAsString(SomaticRefContextEnrichment.TRINUCLEOTIDE_FLAG, Strings.EMPTY);
    final boolean triMatch = triContext.endsWith("GA");

    return isNotFiltered(context) && triMatch && altMatch;
}
 

@Override
public double calculateErrorProbability(final VariantContext vc, final Mutect2FilteringEngine filteringEngine, ReferenceContext referenceContext) {

    final int[] rpa = vc.getAttributeAsList(GATKVCFConstants.REPEATS_PER_ALLELE_KEY).stream()
            .mapToInt(o -> Integer.parseInt(String.valueOf(o))).toArray();
    if (rpa.length < 2) {
        return 0;
    }
    final String ru = vc.getAttributeAsString(GATKVCFConstants.REPEAT_UNIT_KEY, "");

    final int referenceSTRBaseCount = ru.length() * rpa[0];
    final int numPCRSlips = rpa[0] - rpa[1];
    if (referenceSTRBaseCount >= minSlippageLength && Math.abs(numPCRSlips) == 1) {
        // calculate the p-value that out of n reads we would have at least k slippage reads
        // if this p-value is small we keep the variant (reject the PCR slippage hypothesis)
        final int[] ADs = filteringEngine.sumADsOverSamples(vc, true, false);
        if (ADs == null || ADs.length < 2) {
            return 0;
        }
        final int depth = (int) MathUtils.sum(ADs);

        final int altCount = (int) MathUtils.sum(ADs) - ADs[0];
        final double logSomaticLikelihood = filteringEngine.getSomaticClusteringModel().logLikelihoodGivenSomatic(depth, altCount);

        double likelihoodGivenSlippageArtifact;
        try {
            likelihoodGivenSlippageArtifact = Beta.regularizedBeta(slippageRate, ADs[1] + 1, ADs[0] + 1);
        } catch (final MaxCountExceededException e) {
            //if the special function can't be computed, use a binomial with fixed probability
            likelihoodGivenSlippageArtifact = new BinomialDistribution(null, depth, slippageRate).probability(ADs[1]);
        }

        final double logOdds = logSomaticLikelihood - Math.log(likelihoodGivenSlippageArtifact);
        return filteringEngine.posteriorProbabilityOfError(vc, logOdds, 0);
    } else {
        return 0;
    }
}
 

private AnnotatedInterval(final VariantContext vc) {
    sourceVC = vc;
    interval = new SimpleInterval( vc.getContig(), vc.getStart(), vc.getEnd());
    id = vc.getID();
    type = vc.getAttributeAsString(SVTYPE, "");
    svlen = vc.getAttributeAsInt(SVLEN, 0);
    alleles = vc.getAlleles();
}
 

/**
 * Calculate the filter status for a given VariantContext using the combined data from all alleles at a site
 * @param vc
 * @param recals
 * @param builder   is modified by adding attributes
 * @return a String with the filter status for this site
 */
private String doSiteSpecificFiltering(final VariantContext vc, final List<VariantContext> recals, final VariantContextBuilder builder) {
    VariantContext recalDatum = getMatchingRecalVC(vc, recals, null);
    if( recalDatum == null ) {
        throw new UserException("Encountered input variant which isn't found in the input recal file. Please make sure VariantRecalibrator and ApplyRecalibration were run on the same set of input variants. First seen at: " + vc );
    }

    final String lodString = recalDatum.getAttributeAsString(GATKVCFConstants.VQS_LOD_KEY, null);
    if( lodString == null ) {
        throw new UserException("Encountered a malformed record in the input recal file. There is no lod for the record at: " + vc );
    }
    final double lod;
    try {
        lod = Double.valueOf(lodString);
    } catch (NumberFormatException e) {
        throw new UserException("Encountered a malformed record in the input recal file. The lod is unreadable for the record at: " + vc );
    }

    builder.attribute(GATKVCFConstants.VQS_LOD_KEY, lod);
    builder.attribute(GATKVCFConstants.CULPRIT_KEY, recalDatum.getAttribute(GATKVCFConstants.CULPRIT_KEY));
    if(recalDatum != null) {
        if (recalDatum.hasAttribute(GATKVCFConstants.POSITIVE_LABEL_KEY))
            builder.attribute(GATKVCFConstants.POSITIVE_LABEL_KEY, true);
        if (recalDatum.hasAttribute(GATKVCFConstants.NEGATIVE_LABEL_KEY))
            builder.attribute(GATKVCFConstants.NEGATIVE_LABEL_KEY, true);
    }

    return generateFilterString(lod);
}
 

@Nullable
public static String mateId(@NotNull VariantContext context) {
    String mate = context.getAttributeAsString(MATE_ID, null);
    if (mate == null) {
        return context.getAttributeAsString(PAR_ID, null);
    }

    return mate;
}
 

private static boolean isForwardCandidate(@NotNull final VariantContext context) {
    final boolean altMatch =
            context.getAlternateAlleles().stream().anyMatch(x -> x.getBaseString().equals("T") || x.getBaseString().equals("G"));

    final String triContext = context.getAttributeAsString(SomaticRefContextEnrichment.TRINUCLEOTIDE_FLAG, Strings.EMPTY);
    final boolean triMatch = triContext.startsWith("TC");

    return isNotFiltered(context) && triMatch && altMatch;
}
 

@Override
@SuppressWarnings({"unchecked", "rawtypes"})//FIXME generics here blow up
public Map<String, Object> finalizeRawData(final VariantContext vc, final VariantContext originalVC) {
    String rawMQdata;
    if (vc.hasAttribute(getPrimaryRawKey())) {
        rawMQdata = vc.getAttributeAsString(getPrimaryRawKey(), null);
    }
    else if (vc.hasAttribute(getDeprecatedRawKeyName())) {
        if (!allowOlderRawKeyValues) {
            throw new UserException.BadInput("Presence of '-"+getDeprecatedRawKeyName()+"' annotation is detected. This GATK version expects key "
                    + getPrimaryRawKey() + " with a tuple of sum of squared MQ values and total reads over variant "
                    + "genotypes as the value. This could indicate that the provided input was produced with an older version of GATK. " +
                    "Use the argument '--"+RMS_MAPPING_QUALITY_OLD_BEHAVIOR_OVERRIDE_ARGUMENT+"' to override and attempt the deprecated MQ calculation. There " +
                    "may be differences in how newer GATK versions calculate DP and MQ that may result in worse MQ results. Use at your own risk.");
        }

        rawMQdata = vc.getAttributeAsString(getDeprecatedRawKeyName(), null);
        //the original version of ReblockGVCF produces a different MQ format -- try to handle that gracefully here just in case those files go through GenotypeGVCFs
        if (vc.hasAttribute(GATKVCFConstants.MAPPING_QUALITY_DEPTH_DEPRECATED)) {
            logger.warn("Presence of " + GATKVCFConstants.MAPPING_QUALITY_DEPTH_DEPRECATED + " key indicates that this tool may be running on an older output of ReblockGVCF " +
                    "that may not have compatible annotations with this GATK version. Attempting to reformat MQ data.");
            final String rawMQdepth = vc.getAttributeAsString(GATKVCFConstants.MAPPING_QUALITY_DEPTH_DEPRECATED,null);
            if (rawMQdepth == null) {
                throw new UserException.BadInput("MQ annotation data is not properly formatted. This version expects a " +
                        "long tuple of sum of squared MQ values and total reads over variant genotypes.");
            }
            rawMQdata = Math.round(Double.parseDouble(rawMQdata)) + "," + rawMQdepth;  //deprecated format was double so it needs to be converted to long
        }
        else {
            logger.warn("MQ annotation data is not properly formatted. This GATK version expects key "
                    + getPrimaryRawKey() + " with a tuple of sum of squared MQ values and total reads over variant "
                    + "genotypes as the value. Attempting to use deprecated MQ calculation.");
            final long numOfReads = getNumOfReads(vc, null);
            rawMQdata = Math.round(Double.parseDouble(rawMQdata)) + "," + numOfReads;   //deprecated format was double so it needs to be converted to long
        }
    }
    else {
        return new HashMap<>();
    }
    if (rawMQdata == null) {
        return new HashMap<>();
    }

    ReducibleAnnotationData<List<Long>> myData = new ReducibleAnnotationData(rawMQdata);
    parseRawDataString(myData);

    final String annotationString = makeFinalizedAnnotationString(myData.getAttribute(Allele.NO_CALL).get(TOTAL_DEPTH_INDEX), myData.getAttribute(Allele.NO_CALL).get(SUM_OF_SQUARES_INDEX));
    return Collections.singletonMap(getKeyNames().get(0), annotationString);
}
 

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

private static boolean imprecise(@NotNull VariantContext context) {
    final String impreciseStr = context.getAttributeAsString(IMPRECISE, "");
    boolean isPrecise = impreciseStr.isEmpty() || !impreciseStr.equals("true");
    return !isPrecise;
}