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

/**
 * Use the given metadata to create funcotations from variant context attributes (and alt alleles)
 * @param vc Never {@code null}
 * @param metadata Fields that should be present in the funcotations.  Can be a superset of the fields in the
 *                 funcotations.  Never {@code null}
 * @param datasourceName Name to appear in all funcotations.  Never {@code null}
 * @return Instances of {@link Funcotation} for each field in the metadata x alternate allele in the variant context.
 * If a field is not present in the variant context attributes, the field will ave value empty string ("") in all output
 * funcotations.  Fields will be the same names and values for each alternate allele in the funcotations.
 */
static List<Funcotation> createFuncotationsFromMetadata(final VariantContext vc, final FuncotationMetadata metadata, final String datasourceName) {

    Utils.nonNull(vc);
    Utils.nonNull(metadata);
    Utils.nonNull(datasourceName);

    final List<String> fields = metadata.retrieveAllHeaderInfo().stream().map(VCFInfoHeaderLine::getID).collect(Collectors.toList());
    final List<Funcotation> result = new ArrayList<>();
    for (final Allele allele: vc.getAlternateAlleles()) {

        // We must have fields for everything in the metadata.
        final List<String> funcotationFieldValues = new ArrayList<>();
        for (final String funcotationFieldName : fields) {
            funcotationFieldValues.add(vc.getAttributeAsString(funcotationFieldName, ""));
        }

        result.add(TableFuncotation.create(fields, funcotationFieldValues, allele, datasourceName, metadata));
    }

    return result;
}
 

/** Determine whether a given variant context represents a segment.
 *
 * Dev note: this is done by examining the length and the alt allele of the segment.
 *
 * @param vc Never {@code null}
 * @param minSizeForSegment Minimum size for a segment to be valid.
 * @return Boolean whether the given variant context could represent a copy number segment.
 */
public static boolean isSegmentVariantContext(final VariantContext vc, final int minSizeForSegment) {
    Utils.nonNull(vc);
    final List<String> ACCEPTABLE_ALT_ALLELES = Stream.concat(
            Stream.of(SimpleSVType.SupportedType.values())
                    .map(s -> SimpleSVType.createBracketedSymbAlleleString(s.toString())),
            Stream.of(AnnotatedIntervalToSegmentVariantContextConverter.COPY_NEUTRAL_ALLELE.getDisplayString(), Allele.UNSPECIFIED_ALTERNATE_ALLELE_STRING)
    ).collect(Collectors.toList());

    boolean acceptableAlternateAllele = false;
    for (final Allele a: vc.getAlternateAlleles()) {
        final String representation = a.getDisplayString();
        if (ACCEPTABLE_ALT_ALLELES.contains(representation)) {
            acceptableAlternateAllele = true;
        }
    }

    return acceptableAlternateAllele && (VariantContextUtils.getSize(vc) > minSizeForSegment);
}
 

@Override
protected List<Funcotation> createDefaultFuncotationsOnVariant( final VariantContext variant, final ReferenceContext referenceContext ) {

    final List<Funcotation> funcotationList = new ArrayList<>();

    // Add our tally for all alternate alleles in this variant:
    for ( final Allele altAllele : variant.getAlternateAlleles() ) {
        funcotationList.add(
                TableFuncotation.create(
                        new ArrayList<>(supportedFields),
                        new ArrayList<>(Collections.singletonList("")),
                        altAllele,
                        name, null
                )
        );
    }

    return funcotationList;
}
 

private List<Funcotation> createDefaultFuncotationsOnVariantHelper( final VariantContext variant, final ReferenceContext referenceContext, final Set<Allele> annotatedAltAlleles  ) {

        final List<Funcotation> funcotationList = new ArrayList<>();

        if ( supportedFieldNames.size() != 0 ) {

            final List<Allele> alternateAlleles = variant.getAlternateAlleles();

            for ( final Allele altAllele : alternateAlleles ) {
                if ( !annotatedAltAlleles.contains(altAllele) ) {
                    funcotationList.add(createDefaultFuncotation(altAllele));
                }
            }
        }

        return funcotationList;
    }
 

private String makeReducedAnnotationString(final VariantContext vc, final Map<Allele,Double> perAltRankSumResults) {
    final StringBuilder annotationString = new StringBuilder();
    for (final Allele a : vc.getAlternateAlleles()) {
        if (annotationString.length() != 0) {
            annotationString.append(AnnotationUtils.ALLELE_SPECIFIC_REDUCED_DELIM);
        }
        if (!perAltRankSumResults.containsKey(a)) {
            logger.warn("VC allele not found in annotation alleles -- maybe there was trimming?  Allele " + a.getDisplayString() + " will be marked as missing.");
            annotationString.append(VCFConstants.MISSING_VALUE_v4); //add the missing value or the array size and indexes won't check out
        } else {
            final Double numericAlleleValue = perAltRankSumResults.get(a);
            final String perAlleleValue = numericAlleleValue != null ? String.format("%.3f", numericAlleleValue) : VCFConstants.MISSING_VALUE_v4;
            annotationString.append(perAlleleValue);
        }
    }
    return annotationString.toString();
}
 

/**
 Allocate and fill a 2x2 strand contingency table.  In the end, it'll look something like this:
 *             fw      rc
 *   allele1   #       #
 *   allele2   #       #
 * @return a 2x2 contingency table
 */
public static int[][] getContingencyTable( final AlleleLikelihoods<GATKRead, Allele> likelihoods,
                                           final VariantContext vc,
                                           final int minCount,
                                           final Collection<String> samples) {
    if( likelihoods == null || vc == null) {
        return null;
    }

    final Allele ref = vc.getReference();
    final List<Allele> allAlts = vc.getAlternateAlleles();

    final int[][] table = new int[ARRAY_DIM][ARRAY_DIM];
    for (final String sample : samples) {
        final int[] sampleTable = new int[ARRAY_SIZE];
        likelihoods.bestAllelesBreakingTies(sample).stream()
                .filter(ba -> ba.isInformative())
                .forEach(ba -> updateTable(sampleTable, ba.allele, ba.evidence, ref, allAlts));
        if (passesMinimumThreshold(sampleTable, minCount)) {
            copyToMainTable(sampleTable, table);
        }
    }

    return table;
}
 

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

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

private AlleleLikelihoods<GATKRead, Allele> makeLikelihoods(final VariantContext vc, final ReadsContext readsContext) {
    final List<GATKRead> reads = Utils.stream(readsContext).collect(Collectors.toList());
    final AlleleLikelihoods<GATKRead, Allele> result = new AlleleLikelihoods<>(variantSamples,
            new IndexedAlleleList<>(vc.getAlleles()), AssemblyBasedCallerUtils.splitReadsBySample(variantSamples, getHeaderForReads(), reads));

    final ReadPileup readPileup = new ReadPileup(vc, readsContext);
    final Map<String, ReadPileup> pileupsBySample = readPileup.splitBySample(getHeaderForReads(), "__UNKNOWN__");

    final Allele refAllele = vc.getReference();
    final List<Allele> altAlleles = vc.getAlternateAlleles();
    final int numAlleles = vc.getNAlleles();

    // manually fill each sample's likelihoods one read at a time, assigning probability 1 (0 in log space) to the matching allele, if present
    for (final Map.Entry<String, ReadPileup> samplePileups : pileupsBySample.entrySet()) {
        final LikelihoodMatrix<GATKRead, Allele> sampleMatrix = result.sampleMatrix(result.indexOfSample(samplePileups.getKey()));

        final MutableInt readIndex = new MutableInt(0);
        for (final PileupElement pe : samplePileups.getValue()) {
            // initialize all alleles as equally unlikely
            IntStream.range(0, numAlleles).forEach(a -> sampleMatrix.set(a, readIndex.intValue(), Double.NEGATIVE_INFINITY));

            // if present set the matching allele as likely
            final Allele pileupAllele = GATKVariantContextUtils.chooseAlleleForRead(pe, refAllele, altAlleles, minBaseQualityScore);
            if (pileupAllele != null) {
                sampleMatrix.set(sampleMatrix.indexOfAllele(pileupAllele), readIndex.intValue(), 0);
            }


            readIndex.increment();
        }
    }

    return result;
}
 

private static Object splitAltAlleles(final VariantContext vc) {
    final int numAltAlleles = vc.getAlternateAlleles().size();
    if ( numAltAlleles == 1 ) {
        return vc.getAlternateAllele(0);
    }

    return vc.getAlternateAlleles();
}
 

private String makeFinalizedAnnotationString(final VariantContext vc, final Map<Allele, Double> perAlleleValues) {
    final Map<Allele, Integer> variantADs = getADcounts(vc);
    String annotationString = "";
    for (final Allele current : vc.getAlternateAlleles()) {
        if (!annotationString.isEmpty()) {
            annotationString += ",";
        }
        if (perAlleleValues.containsKey(current)) {
            annotationString += String.format(printFormat, Math.sqrt((double) perAlleleValues.get(current) / variantADs.get(current)));
        } else {
            allele_logger.warn("ERROR: VC allele is not found in annotation alleles -- maybe there was trimming?");
        }
    }
    return annotationString;
}
 

protected String makeReducedAnnotationString(VariantContext vc, Map<Allele,Double> perAltsStrandCounts) {
    String annotationString = "";
    for (Allele a : vc.getAlternateAlleles()) {
        if (!annotationString.isEmpty()) {
            annotationString += REDUCED_DELIM;
        }
        if (!perAltsStrandCounts.containsKey(a)) {
            logger.warn("ERROR: VC allele not found in annotation alleles -- maybe there was trimming?");
        } else {
            annotationString += String.format("%.3f", perAltsStrandCounts.get(a));
        }
    }
    return annotationString;
}
 

private void addVariantDatum(final VariantContext vc, final boolean isInput, final FeatureContext context ) {
    if( vc != null && ( IGNORE_ALL_FILTERS || vc.isNotFiltered() || ignoreInputFilterSet.containsAll(vc.getFilters()) ) ) {
        if( VariantDataManager.checkVariationClass( vc, VRAC.MODE ) && !VRAC.useASannotations) {
            addDatum(reduceSum, isInput, context, vc, null, null);
        }
        else if( VRAC.useASannotations ) {
            for (final Allele allele : vc.getAlternateAlleles()) {
                if (!GATKVCFConstants.isSpanningDeletion(allele) && VariantDataManager.checkVariationClass(vc, allele, VRAC.MODE)) {
                    //note that this may not be the minimal representation for the ref and alt allele
                    addDatum(reduceSum, isInput, context, vc, vc.getReference(), allele);
                }
            }
        }
    }
}
 

/**
 * Helper method for testMaxAlternateAlleles
 *
 * @param vc VariantContext to check
 * @return number of alt alleles in vc, excluding NON_REF (if present)
 */
private int getNumAltAllelesExcludingNonRef( final VariantContext vc ) {
    final List<Allele> altAlleles = vc.getAlternateAlleles();
    int numAltAllelesExcludingNonRef = 0;

    for ( final Allele altAllele : altAlleles ) {
        if ( ! altAllele.equals(Allele.NON_REF_ALLELE) ) {
            ++numAltAllelesExcludingNonRef;
        }
    }

    return numAltAllelesExcludingNonRef;
}
 

public List<Object> getRelevantStates(ReferenceContext referenceContext, ReadsContext readsContext, FeatureContext featureContext, VariantContext comp, String compName, VariantContext eval, String evalName, String sampleName, String familyName) {
    if (eval != null) {
        int AC = 0; // by default, the site is considered monomorphic

        try {
            if ( eval.isBiallelic() ) {
                if ( eval.hasAttribute(GATKVCFConstants.MLE_ALLELE_COUNT_KEY) ) {
                    // the MLEAC is allowed to be larger than the AN (e.g. in the case of all PLs being 0, the GT is ./. but the exact model may arbitrarily choose an AC>1)
                    AC = Math.min(eval.getAttributeAsInt(GATKVCFConstants.MLE_ALLELE_COUNT_KEY, 0), nchrom);
                } else if ( eval.hasAttribute(VCFConstants.ALLELE_COUNT_KEY) ) {
                    AC = eval.getAttributeAsInt(VCFConstants.ALLELE_COUNT_KEY, 0);
                }
            }
        } catch ( ClassCastException e ) {
            // protect ourselves from bad inputs
            // TODO -- fully decode VC
        }

        if ( AC == 0 && eval.isVariant() ) {
            // fall back to the direct calculation
            for (Allele allele : eval.getAlternateAlleles())
                AC = Math.max(AC, eval.getCalledChrCount(allele));
        }

        // make sure that the AC isn't invalid
        if ( AC > nchrom )
            throw new UserException(String.format("The AC value (%d) at position %s:%d " +
                    "is larger than the number of chromosomes over all samples (%d)", AC,
                    eval.getContig(), eval.getStart(), nchrom));

        return Collections.singletonList((Object) AC);
    } else {
        return Collections.emptyList();
    }
}
 

/**
 * Attempts to treat the given features as {@link GencodeGtfFeature} objects in order to
 * create funcotations for the given variant and reference.
 *
 * This is the entry point into the Factory
 */
@Override
protected List<Funcotation> createFuncotationsOnVariant(final VariantContext variant, final ReferenceContext referenceContext, final List<Feature> geneFeatureList) {
    final List<Funcotation> outputFuncotations = new ArrayList<>();

    final List<Feature> nonNullFeatures = geneFeatureList.stream().filter(g -> g != null).collect(Collectors.toList());

    // If we have features we need to annotate, go through them and create annotations:
    if ( nonNullFeatures.size() > 0 ) {
        for ( final Allele altAllele : variant.getAlternateAlleles() ) {

            // At this point we know the feature list is composed of GTF Gene Features
            final List<GencodeGtfGeneFeature> gencodeGtfGeneFeatures = convertFeaturesToGencodeGtfGeneFeatures(nonNullFeatures);

            // By this point we know the feature type is correct, so we cast it:
            final List<GencodeFuncotation> gencodeFuncotationList = createGencodeFuncotationsByAllTranscripts(variant, referenceContext, altAllele, gencodeGtfGeneFeatures);

            // Add the returned funcotations here:
            if ( ! gencodeFuncotationList.isEmpty() ) {
                outputFuncotations.addAll(gencodeFuncotationList);
            } else {
                // If the funcotation List is empty, it means that this variant was IGR for all transcripts,
                // or that there were no basic transcripts:
                outputFuncotations.add(createIgrFuncotation(variant, altAllele, referenceContext));
            }
        }
    } else {
        // This is an IGR. We have to have an annotation on all variants, so we must annotate it.
        outputFuncotations.addAll(createIgrFuncotations(variant, referenceContext));
    }

    return outputFuncotations;
}
 

private void calculateSNPPairwiseNovelty(VariantContext eval, VariantContext comp) {
    if ( comp == null )
        return;

    int knownAlleles = 0;
    for ( Allele alt : eval.getAlternateAlleles() ) {
        if ( comp.getAlternateAlleles().contains(alt) )
            knownAlleles++;
    }

    if ( knownAlleles == eval.getAlternateAlleles().size() )
        knownSNPsComplete++;
    else if ( knownAlleles > 0 )
        knownSNPsPartial++;
}
 

@Override
public void update1(final VariantContext eval, final ReferenceContext referenceContext, final ReadsContext readsContext, final FeatureContext featureContext) {
    if ( eval.isIndel() && ! eval.isComplexIndel() ) {
        if ( ! ( getWalker().ignoreAC0Sites() && eval.isMonomorphicInSamples() )) {
            // only if we are actually polymorphic in the subsetted samples should we count the allele
            for ( Allele alt : eval.getAlternateAlleles() ) {
                final int alleleSize = alt.length() - eval.getReference().length();
                if ( alleleSize == 0 ) throw new GATKException("Allele size not expected to be zero for indel: alt = " + alt + " ref = " + eval.getReference());
                updateLengthHistogram(eval.getReference(), alt);
            }
        }
    }
}
 

private List<Funcotation> createDefaultFuncotationsOnVariantHelper( final VariantContext variant, final ReferenceContext referenceContext, final Set<Allele> annotatedAltAlleles  ) {

        final List<Funcotation> funcotationList = new ArrayList<>();

        final List<Allele> alternateAlleles = variant.getAlternateAlleles();

        for ( final Allele altAllele : alternateAlleles ) {
            if ( !annotatedAltAlleles.contains(altAllele) ) {
                funcotationList.add(TableFuncotation.create(annotationColumnNames, emptyAnnotationList, altAllele, name, null));
            }
        }

        return funcotationList;
    }
 

@Override
protected List<Funcotation> createFuncotationsOnVariant(final VariantContext variant, final ReferenceContext referenceContext, final List<Feature> featureList) {
    final List<Funcotation> outputFuncotations = new ArrayList<>();

    final List<Allele> alternateAlleles = variant.getAlternateAlleles();

    // Create a set to put our annotated Alternate alleles in.
    // We'll use this to determine if the alt allele has been annotated.
    final Set<Allele> annotatedAltAlleles = new HashSet<>(alternateAlleles.size());

    if ( !featureList.isEmpty() ) {
        for ( final Feature feature : featureList ) {

            // Get the kind of feature we want here:
            if ( feature != null ) {

                // By this point we know the feature type is correct, so we cast it:
                final XsvTableFeature tableFeature = (XsvTableFeature) feature;

                // Make sure we are annotating the correct XSVTableFeature:
                // NOTE: This is probably redundant to the name check in DataSourceFuncotationFactory
                if ( tableFeature.getDataSourceName().equals(name) ) {

                    // Now we create one funcotation for each Alternate allele:
                    for ( final Allele altAllele : alternateAlleles ) {
                        outputFuncotations.add(TableFuncotation.create(tableFeature, altAllele, name, null));
                        annotatedAltAlleles.add(altAllele);
                    }
                }

                // TODO: Must break the loop now to prevent multiple entries messing up the number of fields in the funcotation (issue #4930 - https://github.com/broadinstitute/gatk/issues/4930)
                break;
            }
        }
    }

    // If we didn't add funcotations for an allele, we should add in blank funcotations to that allele for each field that can be produced
    // by this LocatableXsvFuncotationFactory:
    if ( annotatedAltAlleles.size() != alternateAlleles.size() ) {
        outputFuncotations.addAll(createDefaultFuncotationsOnVariantHelper(variant, referenceContext, annotatedAltAlleles));
    }

    return outputFuncotations;
}
 

@Override
public void write(final VariantContext variant, final FuncotationMap txToFuncotationMap) {

    // Create a new variant context builder:
    final VariantContextBuilder variantContextOutputBuilder = new VariantContextBuilder(variant);

    final StringBuilder funcotatorAnnotationStringBuilder = new StringBuilder();

    // Get the old VCF Annotation field and append the new information to it:
    final Object existingAnnotation = variant.getAttribute(FUNCOTATOR_VCF_FIELD_NAME, null);
    final List<String> existingAlleleAnnotations;
    if ( existingAnnotation != null) {
        existingAlleleAnnotations = Utils.split(existingAnnotation.toString(), ',');
    }
    else {
        existingAlleleAnnotations = Collections.emptyList();
    }

    // Go through each allele and add it to the writer separately:
    final List<Allele> alternateAlleles = variant.getAlternateAlleles();
    for ( int alleleIndex = 0; alleleIndex < alternateAlleles.size() ; ++alleleIndex ) {

        final Allele altAllele = alternateAlleles.get(alleleIndex);

        if ( alleleIndex < existingAlleleAnnotations.size() ) {
            funcotatorAnnotationStringBuilder.append( existingAlleleAnnotations.get(alleleIndex) );
            funcotatorAnnotationStringBuilder.append(FIELD_DELIMITER);
        }

        for (final String txId : txToFuncotationMap.getTranscriptList()) {
            funcotatorAnnotationStringBuilder.append(START_TRANSCRIPT_DELIMITER);
            final List<Funcotation> funcotations = txToFuncotationMap.get(txId);
            final Funcotation manualAnnotationFuncotation = createManualAnnotationFuncotation(altAllele);

            funcotatorAnnotationStringBuilder.append(
                    Stream.concat(funcotations.stream(), Stream.of(manualAnnotationFuncotation))
                            .filter(f -> f.getAltAllele().equals(altAllele))
                            .filter(f -> f.getFieldNames().size() > 0)
                            .filter(f -> !f.getDataSourceName().equals(FuncotatorConstants.DATASOURCE_NAME_FOR_INPUT_VCFS))
                            .map(VcfOutputRenderer::adjustIndelAlleleInformation)
                            .map(f -> FuncotatorUtils.renderSanitizedFuncotationForVcf(f, finalFuncotationFieldNames))
                            .collect(Collectors.joining(FIELD_DELIMITER))
            );

            funcotatorAnnotationStringBuilder.append(END_TRANSCRIPT_DELIMITER + ALL_TRANSCRIPT_DELIMITER);
        }
        // We have a trailing "#" - we need to remove it:
        funcotatorAnnotationStringBuilder.deleteCharAt(funcotatorAnnotationStringBuilder.length()-1);
        funcotatorAnnotationStringBuilder.append(VCFConstants.INFO_FIELD_ARRAY_SEPARATOR);
    }

    // We have a trailing "," - we need to remove it:
    funcotatorAnnotationStringBuilder.deleteCharAt(funcotatorAnnotationStringBuilder.length()-1);

    // Add our new annotation:
    variantContextOutputBuilder.attribute(FUNCOTATOR_VCF_FIELD_NAME, funcotatorAnnotationStringBuilder.toString());

    // Add the genotypes from the variant:
    variantContextOutputBuilder.genotypes( variant.getGenotypes() );

    // Render and add our VCF line:
    vcfWriter.add( variantContextOutputBuilder.make() );
}