Java Code Examples for htsjdk.variant.variantcontext.Genotype#getAllele()

private IlluminaGenotype getIlluminaGenotype(final Genotype genotype, final VariantContext context) {
    final IlluminaGenotype illuminaGenotype;
    if (genotype.isCalled()) {
        // Note that we remove the trailing '*' that appears in alleles that are reference.
        final String illuminaAlleleA = StringUtils.stripEnd(getStringAttribute(context, InfiniumVcfFields.ALLELE_A), "*");
        final String illuminaAlleleB = StringUtils.stripEnd(getStringAttribute(context, InfiniumVcfFields.ALLELE_B), "*");
        if (genotype.getAlleles().size() != 2) {
            throw new PicardException("Unexpected number of called alleles in variant context " + context + " found alleles: " + genotype.getAlleles());
        }
        final Allele calledAllele1 = genotype.getAllele(0);
        final Allele calledAllele2 = genotype.getAllele(1);

        if (calledAllele1.basesMatch(illuminaAlleleA)) {
            if (calledAllele2.basesMatch(illuminaAlleleA)) {
                illuminaGenotype = picard.arrays.illumina.IlluminaGenotype.AA;
            } else if (calledAllele2.basesMatch(illuminaAlleleB)) {
                illuminaGenotype = picard.arrays.illumina.IlluminaGenotype.AB;
            } else {
                throw new PicardException("Error matching called alleles to Illumina alleles.  Context: " + context);
            }
        } else if (calledAllele1.basesMatch(illuminaAlleleB)) {
            if (calledAllele2.basesMatch(illuminaAlleleA)) {
                illuminaGenotype = picard.arrays.illumina.IlluminaGenotype.AB;
            } else if (calledAllele2.basesMatch(illuminaAlleleB)) {
                illuminaGenotype = picard.arrays.illumina.IlluminaGenotype.BB;
            } else {
                throw new PicardException("Error matching called alleles to Illumina alleles.  Context: " + context);
            }
        } else {
            // We didn't match up the illumina alleles to called alleles
            throw new PicardException("Error matching called alleles to Illumina alleles.  Context: " + context);
        }
    } else {
        illuminaGenotype = picard.arrays.illumina.IlluminaGenotype.NN;
    }
    return illuminaGenotype;
}
 

/** Tests whether the alleles of the offspring are possible given the alleles of the parents. */
private boolean isMendelianViolation(final Genotype p1, final Genotype p2, final Genotype off) {
    final Allele offAllele1 = off.getAllele(0);
    final Allele offAllele2 = off.getAllele(1);

    if(p1.getAlleles().contains(offAllele1) && p2.getAlleles().contains(offAllele2)) return false;
    if(p2.getAlleles().contains(offAllele1) && p1.getAlleles().contains(offAllele2)) return false;

    return true;
}
 

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

/** Perform basic somatic pileup validation and return a result instance.
 *
 * NOTE: This only looks at the first alternate allele.  Multiallelics are not supported.
 *
 * @param genotype given genotype to validate.  Never {@code null}
 * @param referenceAllele Never {@code null}
 * @param validationNormalPileup Pileup for the coresponding location of the genotype in the validation normal.  Never {@code null}
 * @param validationTumorAltCount Alt read count for the validation tumor sample.  Must be positive or zero.
 * @param validationTumorTotalCount Total read count for the validation tumor sample.  Must be positive or zero.
 * @param minBaseQualityCutoff Minimum base quality threshold to count any read.  Must be positive or zero.
 * @param interval location represented in the genotype.  Never {@code null}
 * @param filters additional filters besides the ones on the genotype.  Typically,
 *                the parent variant context filters is provided here.  Never {@code null}
 * @return validation result.  {@code null} if the genotype cannot be validated
 */
public static BasicValidationResult calculateBasicValidationResult(final Genotype genotype, final Allele referenceAllele,
                                                                   final ReadPileup validationNormalPileup,
                                                                   int validationTumorAltCount, int validationTumorTotalCount,
                                                                   int minBaseQualityCutoff, final SimpleInterval interval, final String filters) {

    if (!isAbleToValidateGenotype(genotype, referenceAllele) || validationNormalPileup == null){
        return null;
    }

    Utils.nonNull(referenceAllele);
    Utils.nonNull(genotype);
    Utils.nonNull(interval);
    Utils.nonNull(filters);
    ParamUtils.isPositiveOrZero(validationTumorAltCount, "Validation alt count must be >= 0");
    ParamUtils.isPositiveOrZero(validationTumorTotalCount, "Validation total count must be >= 0");
    ParamUtils.isPositiveOrZero(minBaseQualityCutoff, "Minimum base quality cutoff must be >= 0");

    final double maxAltRatioSeenInNormalValidation = PowerCalculationUtils.calculateMaxAltRatio(validationNormalPileup,
            referenceAllele, minBaseQualityCutoff);
    final int discoveryTumorAltCount = genotype.getAD()[1];
    final int discoveryTumorTotalCount = genotype.getAD()[0] + discoveryTumorAltCount;

    if (Double.isNaN(maxAltRatioSeenInNormalValidation) || discoveryTumorTotalCount == 0) {
        return null;
    }

    final int minCountForSignal = PowerCalculationUtils.calculateMinCountForSignal(validationTumorTotalCount, maxAltRatioSeenInNormalValidation);
    final double power = PowerCalculationUtils.calculatePower(validationTumorTotalCount, discoveryTumorAltCount, discoveryTumorTotalCount, minCountForSignal);

    boolean isNotNoise = (validationTumorAltCount >= minCountForSignal);
    boolean isEnoughValidationCoverageToValidate = validationTumorAltCount >= 2;

    final String genotypeFilters = genotype.getFilters() == null ? "" : genotype.getFilters();
    final long numReadsSupportingAltInValidationNormal = PowerCalculationUtils.calculateNumReadsSupportingAllele(validationNormalPileup,
            genotype.getAllele(0), genotype.getAllele(1), minBaseQualityCutoff);
    return new BasicValidationResult(interval, minCountForSignal, isEnoughValidationCoverageToValidate,
            isNotNoise, power, validationTumorAltCount, validationTumorTotalCount-validationTumorAltCount,
            discoveryTumorAltCount, discoveryTumorTotalCount-discoveryTumorAltCount, referenceAllele,
            genotype.getAllele(1), filters + genotypeFilters, numReadsSupportingAltInValidationNormal);
}