Java Code Examples for htsjdk.variant.variantcontext.Allele#isSymbolic()

/**
 *
 * @param referenceAllele allele to treat as reference.  Create with {@link Allele#create(String, boolean)}, where
 *                  second parameter is {@code true}.  Never {@code null}.   If the reference is symbolic, exception will be thrown.
 * @param alternateAlleles List of alleles to treat as the alternates.  Easy to create with {@link Allele#create(String, boolean)}, where
 *                  second parameter is {@code false}.  Never {@code null}
 * @param minBaseQualityCutoff minimum base quality for the bases that match the allele in order to be counted.
 *                             Must be positive or zero.
 */
public AllelePileupCounter(final Allele referenceAllele, final List<Allele> alternateAlleles, int minBaseQualityCutoff) {

    this.referenceAllele = Utils.nonNull(referenceAllele);
    this.alternateAlleles = Utils.nonNull(alternateAlleles);

    // Additional checks
    if (referenceAllele.isSymbolic()) {
        throw new UserException.BadInput("A symbolic reference allele was specified.");
    }

    Utils.validateArg(!referenceAllele.isNonReference(), "Reference allele was non-reference: " + referenceAllele);
    Utils.validateArg(alternateAlleles.stream().allMatch(a -> a.isNonReference()),
            "One or more alternate alleles were reference: " + alternateAlleles.stream().map(a-> a.toString()).collect(Collectors.joining(", ")));

    this.minBaseQualityCutoff = ParamUtils.isPositiveOrZero(minBaseQualityCutoff, "Minimum base quality must be positive or zero.");;

    alternateAlleles.forEach(a -> countMap.put(a, new MutableInt(0)));
    countMap.put(referenceAllele, new MutableInt(0));
}
 

/**
 * @param genotype The genotype to test whether we can attempt to validate.  Never {@code null}
 * @param referenceAllele The reference allele (from parent VariantContext).  Never {@code null}
 * @return whether this class can even attempt a validation of the genotype in question.
 */
public static boolean isAbleToValidateGenotype(final Genotype genotype, final Allele referenceAllele) {
    Utils.nonNull(genotype);
    Utils.nonNull(referenceAllele);

    // In order to proceed, we have some assumptions:
    //  - The genotype has a ploidy of 2
    //  - The genotype is a simple indel or a xNP
    //  - The first allele of the genotype is reference
    final boolean isDiploid = genotype.getAlleles().size() == 2;
    final boolean doesGenotypeHaveReference = genotype.getAllele(0).equals(referenceAllele);
    final boolean isReferenceNotSymbolic = !referenceAllele.isSymbolic();
    final VariantContext.Type variantType = GATKVariantContextUtils.typeOfVariant(genotype.getAllele(0), genotype.getAllele(1));
    final boolean isValidatableVariantType = VALIDATABLE_TYPES.contains(variantType)
                    && !GATKVariantContextUtils.isComplexIndel(genotype.getAllele(0), genotype.getAllele(1));
    final boolean hasKnownCoverage = genotype.hasAD() && (genotype.getAD().length == 2);
    final boolean isValidateable = (isDiploid && doesGenotypeHaveReference && isValidatableVariantType && hasKnownCoverage && isReferenceNotSymbolic);
    if (!isValidateable) {
        logger.info("Cannot validate genotype: " + genotype + "  ploidy2: " + isDiploid +
                "  genotypeHasReferenceAllele: " + doesGenotypeHaveReference + "   validatableVariant: " + isValidatableVariantType +
                "  hasCompleteAD field: " + hasKnownCoverage + "  isNonSymbolicReference: " + isReferenceNotSymbolic);
    }

    return isValidateable;
}
 

/**
 * Checks the given variant alt allele for whether it's a special case.
 * If so, will create the appropriate funcotation for that case.
 * @param variant {@link VariantContext} for the given {@code altAllele}.
 * @param altAllele The alternate {@link Allele} to be checked for special cases.
 * @param reference {@link ReferenceContext} supporting the given {@code variant}.
 * @param transcript {@link GencodeGtfTranscriptFeature} containing the given {@code variant}.
 * @return A {@link GencodeFuncotation} for the given special case alt allele, or {@code null} if the allele is not a special case.
 */
private GencodeFuncotation checkForAndCreateSpecialAlleleFuncotation(final VariantContext variant,
                                                                     final Allele altAllele,
                                                                     final ReferenceContext reference,
                                                                     final GencodeGtfTranscriptFeature transcript) {
    // If the alt allele is a spanning deletion or a symbolic allele, create an unknown funcotation:
    if (altAllele.isSymbolic() || altAllele.equals(Allele.SPAN_DEL) ) {
        return createFuncotationForSymbolicAltAllele(variant, altAllele, transcript, reference);
    }

    // If the reference allele or alternate allele contains any masked bases:
    if ( variant.getReference().getBaseString().contains(FuncotatorConstants.MASKED_ANY_BASE_STRING) || altAllele.getBaseString().contains(FuncotatorConstants.MASKED_ANY_BASE_STRING) ) {
        return createFuncotationForMaskedBases(variant, altAllele, transcript, reference);
    }

    return null;
}
 

/**
 * Determines the variant type based on the given reference allele and alternate allele.
 * @param refAllele The reference {@link Allele} for this variant.
 * @param altAllele The alternate {@link Allele} for this variant.
 * @return A {@link GencodeFuncotation.VariantType} representing the variation type between the given reference and alternate {@link Allele}.
 * Spanning deletions and no calls will get a type of {@link org.broadinstitute.hellbender.tools.funcotator.dataSources.gencode.GencodeFuncotation.VariantType#NA}
 */
@VisibleForTesting
static GencodeFuncotation.VariantType getVariantType( final Allele refAllele, final Allele altAllele ) {

    if ( altAllele.length() > refAllele.length() ) {
        return GencodeFuncotation.VariantType.INS;
    }
    else if (altAllele.equals(Allele.SPAN_DEL) || altAllele.equals(Allele.NO_CALL) || altAllele.isSymbolic() ) {
        return GencodeFuncotation.VariantType.NA;
    }
    else if (altAllele.length() < refAllele.length()) {
        return GencodeFuncotation.VariantType.DEL;
    }
    else {
        // We know they are the same length, now we just need to check one of them:
        switch (refAllele.length()) {
            case 1:  return GencodeFuncotation.VariantType.SNP;
            case 2:  return GencodeFuncotation.VariantType.DNP;
            case 3:  return GencodeFuncotation.VariantType.TNP;
            default: return GencodeFuncotation.VariantType.ONP;
        }
    }
}
 

private static Allele reverseComplement(final Allele oldAllele) {
    if (oldAllele.isSymbolic() || oldAllele.isNoCall() || oldAllele.equals(Allele.SPAN_DEL)) {
        return oldAllele;
    } else {
        return Allele.create(SequenceUtil.reverseComplement(oldAllele.getBaseString()), oldAllele.isReference());
    }
}
 

protected static boolean checkVariationClass( final VariantContext evalVC, final Allele allele, final VariantRecalibratorArgumentCollection.Mode mode ) {
    switch( mode ) {
        case SNP:
            //note that spanning deletions are considered SNPs by this logic
            return evalVC.getReference().length() == allele.length();
        case INDEL:
            return (evalVC.getReference().length() != allele.length()) || allele.isSymbolic();
        case BOTH:
            return true;
        default:
            throw new IllegalStateException( "Encountered unknown recal mode: " + mode );
    }
}
 

/**
 * Creates a {@link GencodeFuncotation}s based on the given {@link Allele} with type
 * {@link GencodeFuncotation.VariantClassification#IGR}.
 * Reports reference bases as if they are on the {@link Strand#POSITIVE} strand.
 * @param variant The {@link VariantContext} associated with this annotation.
 * @param altAllele The alternate {@link Allele} to use for this {@link GencodeFuncotation}.
 * @param reference The {@link ReferenceContext} in which the given {@link Allele}s appear.
 * @return An IGR funcotation for the given allele.
 */
private GencodeFuncotation createIgrFuncotation(final VariantContext variant,
                                                final Allele altAllele,
                                                final ReferenceContext reference){

    final GencodeFuncotationBuilder funcotationBuilder = new GencodeFuncotationBuilder();

    // Get GC Content:
    funcotationBuilder.setGcContent( calculateGcContent( variant.getReference(), altAllele, reference, gcContentWindowSizeBases ) );

    final String alleleString = altAllele.getBaseString().isEmpty() ? altAllele.toString() : altAllele.getBaseString();

    funcotationBuilder.setVariantClassification( GencodeFuncotation.VariantClassification.IGR )
                      .setRefAllele( variant.getReference() )
                      .setTumorSeqAllele2( alleleString )
                      .setStart(variant.getStart())
                      .setEnd(variant.getEnd())
                      .setVariantType(getVariantType(variant.getReference(), altAllele))
                      .setChromosome(variant.getContig())
                      .setAnnotationTranscript(FuncotationMap.NO_TRANSCRIPT_AVAILABLE_KEY);

    if ( (!altAllele.isSymbolic()) && (!altAllele.equals(Allele.SPAN_DEL)) ) {
        funcotationBuilder.setGenomeChange(getGenomeChangeString(variant, altAllele));
    }

    funcotationBuilder.setNcbiBuild(ncbiBuildVersion);

    // Set our reference context in the the FuncotatonBuilder:
    funcotationBuilder.setReferenceContext(FuncotatorUtils.createReferenceSnippet(variant.getReference(), altAllele, reference, Strand.POSITIVE, referenceWindow).getBaseString());

    // Set our version:
    funcotationBuilder.setVersion(version);

    // Set our data source name:
    funcotationBuilder.setDataSourceName(getName());

    return funcotationBuilder.build();
}
 

/**
 * Creates a {@link GencodeFuncotation}s based on the given {@link Allele} with type
 * {@link GencodeFuncotation.VariantClassification#FIVE_PRIME_FLANK} or
 * {@link GencodeFuncotation.VariantClassification#THREE_PRIME_FLANK}
 *
 * @param variant The {@link VariantContext} associated with this annotation.
 * @param altAllele The alternate {@link Allele} to use for this {@link GencodeFuncotation}.
 * @param transcript {@link GencodeGtfTranscriptFeature} associated with this flank funcotation.
 * @param reference The {@link ReferenceContext} in which the given {@link Allele}s appear.
 * @param flankType {@link GencodeFuncotation.VariantClassification#FIVE_PRIME_FLANK} or
 *                  {@link GencodeFuncotation.VariantClassification#THREE_PRIME_FLANK}
 * @return A flank funcotation for the given allele
 */
private GencodeFuncotation createFlankFuncotation(final VariantContext variant, final Allele altAllele, final GencodeGtfTranscriptFeature transcript, final ReferenceContext reference, final GencodeFuncotation.VariantClassification flankType) {

    // Create a symbolic allele funcotation for the flank type, if appropriate:
    if (altAllele.isSymbolic() || altAllele.equals(Allele.SPAN_DEL) ) {
        return createFuncotationForSymbolicAltAllele(variant, altAllele, transcript, reference, flankType);
    }

    final GencodeFuncotationBuilder funcotationBuilder = new GencodeFuncotationBuilder();

    // NOTE: The reference context is ALWAYS from the + strand
    final StrandCorrectedReferenceBases referenceBases = FuncotatorUtils.createReferenceSnippet(variant.getReference(), altAllele, reference, Strand.POSITIVE, referenceWindow);
    final String referenceBasesString = referenceBases.getBaseString(Strand.POSITIVE);

    final double gcContent = calculateGcContent(variant.getReference(), altAllele, reference, gcContentWindowSizeBases);

    funcotationBuilder.setHugoSymbol(transcript.getGeneName())
            .setChromosome(variant.getContig())
            .setStart(variant.getStart())
            .setEnd(variant.getEnd())
            .setStrand(transcript.getGenomicStrand())
            .setVariantClassification(flankType)
            .setVariantType(getVariantType(variant.getReference(), altAllele))
            .setRefAllele(variant.getReference())
            .setTumorSeqAllele2(altAllele.getBaseString())
            .setGenomeChange(getGenomeChangeString(variant, altAllele))
            .setAnnotationTranscript(transcript.getTranscriptId())
            .setReferenceContext(referenceBasesString)
            .setGcContent(gcContent)
            .setNcbiBuild(ncbiBuildVersion)
            .setGeneTranscriptType(transcript.getTranscriptType());

    // Set our version:
    funcotationBuilder.setVersion(version);

    // Set our data source name:
    funcotationBuilder.setDataSourceName(getName());

    return funcotationBuilder.build();
}
 

@Override
public void apply(VariantContext discoveryVariantContext, ReadsContext readsContext, ReferenceContext referenceContext, FeatureContext featureContext) {

    final Genotype genotype = discoveryVariantContext.getGenotype(discoverySampleInVcf);

    // Skip any symbolic reference alleles
    final Allele referenceAllele = discoveryVariantContext.getReference();
    if (referenceAllele.isSymbolic()) {
        logger.warn("Skipping variant with symbolic reference allele: " + discoveryVariantContext);
        return;
    }

    // If we cannot validate this genotype, we should simply skip it.
    final boolean isAbleToValidate = BasicSomaticShortMutationValidator.isAbleToValidateGenotype(genotype, referenceAllele);
    if (!isAbleToValidate) {
        if (annotatedVcf != null) {
            vcfWriter.add(new VariantContextBuilder(discoveryVariantContext).attribute(JUDGMENT_INFO_FIELD_KEY, Judgment.SKIPPED).make());
        }
        return;
    }

    // We assume that there is only one alternate allele that we are interested in.
    // Multiallelics are not supported.
    final Allele altAllele = genotype.getAllele(1);

    final SAMFileHeader samFileHeader = getHeaderForReads();

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

    // This could happen when read realignment moves reads such that a particular locus has zero reads
    if (pileupsBySample.isEmpty()){
        return;
    }

    final ReadPileup validationNormalPileup = pileupsBySample.get(validationControlName);
    // TODO: handle this case more carefully
    if (validationNormalPileup == null){
        return;
    }

    final ReadPileup validationTumorPileup = pileupsBySample.get(validationCaseName);

    final Map<Allele, MutableInt> validationTumorAllelicCounts = new AllelePileupCounter(referenceAllele, discoveryVariantContext.getAlternateAlleles(), minBqCutoff, validationTumorPileup)
            .getCountMap();

    final int validationTumorAltCount = validationTumorAllelicCounts.getOrDefault(altAllele, new MutableInt(0)).intValue();
    final int validationTumorRefCount = validationTumorAllelicCounts.get(referenceAllele).intValue();

    final String filterString = discoveryVariantContext.getFilters().stream().sorted().collect(Collectors.joining(";"));
    final BasicValidationResult basicValidationResult = BasicSomaticShortMutationValidator.calculateBasicValidationResult(
            genotype, referenceAllele, validationNormalPileup, validationTumorAltCount,
            validationTumorRefCount + validationTumorAltCount, minBqCutoff,
            new SimpleInterval(discoveryVariantContext.getContig(), discoveryVariantContext.getStart(), discoveryVariantContext.getEnd()),
            filterString);

    if (basicValidationResult != null) {
        results.add(basicValidationResult);
    }

    final boolean normalArtifact = basicValidationResult.getNumAltSupportingReadsInNormal() > maxValidationNormalCount;
    final boolean validated = !normalArtifact && basicValidationResult != null && basicValidationResult.isOutOfNoiseFloor();
    final boolean powered = normalArtifact || (basicValidationResult != null && basicValidationResult.getPower() > minPower);

    if (discoveryVariantContext.isSNP()) {
        if (validated) {
            snpTruePositiveCount.increment();
        } else if (powered) {
            snpFalsePositiveCount.increment();
        }
    } else {
        if (validated) {
            indelTruePositiveCount.increment();
        } else if (powered) {
            indelFalsePositiveCount.increment();
        }
    }

    if (annotatedVcf != null) {
        vcfWriter.add(new VariantContextBuilder(discoveryVariantContext)
                .attribute(JUDGMENT_INFO_FIELD_KEY, validated ? Judgment.VALIDATED : Judgment.UNVALIDATED)
                .attribute(POWER_INFO_FIELD_KEY, basicValidationResult == null ? 0 : basicValidationResult.getPower())
                .attribute(VALIDATION_AD_INFO_FIELD_KEY, new int[] {validationTumorRefCount, validationTumorAltCount}).make());
    }
}