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

private void assertVariantSegmentsAreCovered(final List<VariantContext> variants,
                                             final List<HiddenStateSegmentRecord<CopyNumberTriState, Target>> variantSegments) {
    for (final HiddenStateSegmentRecord<CopyNumberTriState, Target> variantSegment : variantSegments) {
        final Optional<VariantContext> match = variants.stream()
                .filter(vc -> new SimpleInterval(vc).equals(variantSegment.getSegment().getInterval()))
                .findFirst();
        Assert.assertTrue(match.isPresent());
        final VariantContext matchedVariant = match.get();
        final Genotype genotype = matchedVariant.getGenotype(variantSegment.getSampleName());
        final String discovery = genotype.getAnyAttribute(XHMMSegmentGenotyper.DISCOVERY_KEY).toString();
        Assert.assertTrue(discovery.equals(XHMMSegmentGenotyper.DISCOVERY_TRUE));
        final CopyNumberTriState call = variantSegment.getSegment().getCall();

        final List<Allele> gt = genotype.getAlleles();
        Assert.assertEquals(gt.size(), 1);
        // The call may not be the same for case where the event-quality is relatively low.
        if (variantSegment.getSegment().getEventQuality() > 10) {
            Assert.assertEquals(CopyNumberTriStateAllele.valueOf(gt.get(0)).state, call, genotype.toString());
        }
        final String[] SQ = genotype.getAnyAttribute(XHMMSegmentGenotyper.SOME_QUALITY_KEY).toString().split(VCFConstants.INFO_FIELD_ARRAY_SEPARATOR);
        final double someQual = variantSegment.getSegment().getSomeQuality();
        Assert.assertEquals(Double.parseDouble(SQ[call == CopyNumberTriState.DELETION ? 0 : 1]), someQual, XHMMSegmentGenotyper.PHRED_SCORE_PRECISION, variantSegment.getSampleName() + " => " + genotype.toString());

        final String[] LQ = genotype.getAnyAttribute(XHMMSegmentGenotyper.START_QUALITY_KEY).toString().split(VCFConstants.INFO_FIELD_ARRAY_SEPARATOR);
        final double startQuality = variantSegment.getSegment().getStartQuality();
        Assert.assertEquals(Double.parseDouble(LQ[call == CopyNumberTriState.DELETION ? 0 : 1]), startQuality, XHMMSegmentGenotyper.PHRED_SCORE_PRECISION, variantSegment.getSampleName() + " => " + genotype.toString());

        final String[] RQ = genotype.getAnyAttribute(XHMMSegmentGenotyper.END_QUALITY_KEY).toString().split(VCFConstants.INFO_FIELD_ARRAY_SEPARATOR);
        final double endQuality = variantSegment.getSegment().getEndQuality();
        Assert.assertEquals(Double.parseDouble(RQ[call == CopyNumberTriState.DELETION ? 0 : 1]), endQuality, XHMMSegmentGenotyper.PHRED_SCORE_PRECISION, variantSegment.getSampleName() + " => " + genotype.toString());

        // Check the PL.
        final int[] PL = genotype.getPL();
        final int observedGQFromPL = Math.min(XHMMSegmentGenotyper.MAX_GQ, PL[CopyNumberTriStateAllele.REF.index()] - PL[CopyNumberTriStateAllele.valueOf(call).index()]);
        final double expectedCallPL = GATKProtectedMathUtils.roundPhred(QualityUtils.phredScaleErrorRate(QualityUtils.qualToProb(variantSegment.getSegment().getExactQuality())), HMMPostProcessor.PHRED_SCORE_PRECISION);
        final double expectedRefPL = GATKProtectedMathUtils.roundPhred(QualityUtils.phredScaleCorrectRate(QualityUtils.qualToProb(variantSegment.getSegment().getEventQuality())), HMMPostProcessor.PHRED_SCORE_PRECISION);
        final int expectedGQFromPL = Math.min(XHMMSegmentGenotyper.MAX_GQ, (int) Math.round(expectedRefPL - expectedCallPL));
        Assert.assertTrue(Math.abs(observedGQFromPL - expectedGQFromPL) <= 1, genotype.toString() + " " + variantSegment.getSegment().getEventQuality());
    }
}
 

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

protected static GenotypesContext fixGenotypes(final GenotypesContext originals, final List<Allele> originalAlleles, final List<Allele> newAlleles) {
    // optimization: if nothing needs to be fixed then don't bother
    if (originalAlleles.equals(newAlleles)) {
        return originals;
    }

    if (originalAlleles.size() != newAlleles.size()) {
        throw new IllegalStateException("Error in allele lists: the original and new allele lists are not the same length: " + originalAlleles.toString() + " / " + newAlleles.toString());
    }

    final Map<Allele, Allele> alleleMap = new HashMap<>();
    for (int idx = 0; idx < originalAlleles.size(); idx++) {
        alleleMap.put(originalAlleles.get(idx), newAlleles.get(idx));
    }

    final GenotypesContext fixedGenotypes = GenotypesContext.create(originals.size());
    for (final Genotype genotype : originals) {
        final List<Allele> fixedAlleles = new ArrayList<>();
        for (final Allele allele : genotype.getAlleles()) {
            if (allele.isNoCall()) {
                fixedAlleles.add(allele);
            } else {
                Allele newAllele = alleleMap.get(allele);
                if (newAllele == null) {
                    throw new IllegalStateException("Allele not found: " + allele.toString() + ", " + originalAlleles + "/ " + newAlleles);
                }
                fixedAlleles.add(newAllele);
            }
        }
        fixedGenotypes.add(new GenotypeBuilder(genotype).alleles(fixedAlleles).make());
    }
    return fixedGenotypes;
}
 

private boolean haveSameGenotypes(final Genotype g1, final Genotype g2) {
    if (g1 == null || g2 == null) {
        return false;
    }

    if ((g1.isCalled() && g2.isFiltered()) ||
            (g2.isCalled() && g1.isFiltered()) ||
            (g1.isFiltered() && g2.isFiltered() && XLfiltered)) {
        return false;
    }

    final List<Allele> a1s = g1.getAlleles();
    final List<Allele> a2s = g2.getAlleles();
    return (a1s.containsAll(a2s) && a2s.containsAll(a1s));
}
 

private static boolean hasSpanningDeletionAllele(final Genotype gt) {
    for(final Allele a : gt.getAlleles()) {
        boolean hasSpanningDeletion = GATKVCFConstants.isSpanningDeletion(a);
        if(hasSpanningDeletion) {
            return true;
        }
    }
    return false;
}
 

/**
 * method to swap the reference and alt alleles of a bi-allelic, SNP
 *
 * @param vc                   the {@link VariantContext} (bi-allelic SNP) that needs to have it's REF and ALT alleles swapped.
 * @param annotationsToReverse INFO field annotations (of double value) that will be reversed (x->1-x)
 * @param annotationsToDrop    INFO field annotations that will be dropped from the result since they are invalid when REF and ALT are swapped
 * @return a new {@link VariantContext} with alleles swapped, INFO fields modified and in the genotypes, GT, AD and PL corrected appropriately
 */
public static VariantContext swapRefAlt(final VariantContext vc, final Collection<String> annotationsToReverse, final Collection<String> annotationsToDrop) {

    if (!vc.isBiallelic() || !vc.isSNP()) {
        throw new IllegalArgumentException("swapRefAlt can only process biallelic, SNPS, found " + vc.toString());
    }

    final VariantContextBuilder swappedBuilder = new VariantContextBuilder(vc);

    swappedBuilder.attribute(SWAPPED_ALLELES, true);

    // Use getBaseString() (rather than the Allele itself) in order to create new Alleles with swapped
    // reference and non-variant attributes
    swappedBuilder.alleles(Arrays.asList(vc.getAlleles().get(1).getBaseString(), vc.getAlleles().get(0).getBaseString()));

    final Map<Allele, Allele> alleleMap = new HashMap<>();

    // A mapping from the old allele to the new allele, to be used when fixing the genotypes
    alleleMap.put(vc.getAlleles().get(0), swappedBuilder.getAlleles().get(1));
    alleleMap.put(vc.getAlleles().get(1), swappedBuilder.getAlleles().get(0));

    final GenotypesContext swappedGenotypes = GenotypesContext.create(vc.getGenotypes().size());
    for (final Genotype genotype : vc.getGenotypes()) {
        final List<Allele> swappedAlleles = new ArrayList<>();
        for (final Allele allele : genotype.getAlleles()) {
            if (allele.isNoCall()) {
                swappedAlleles.add(allele);
            } else {
                swappedAlleles.add(alleleMap.get(allele));
            }
        }
        // Flip AD
        final GenotypeBuilder builder = new GenotypeBuilder(genotype).alleles(swappedAlleles);
        if (genotype.hasAD() && genotype.getAD().length == 2) {
            final int[] ad = ArrayUtils.clone(genotype.getAD());
            ArrayUtils.reverse(ad);
            builder.AD(ad);
        } else {
            builder.noAD();
        }

        //Flip PL
        if (genotype.hasPL() && genotype.getPL().length == 3) {
            final int[] pl = ArrayUtils.clone(genotype.getPL());
            ArrayUtils.reverse(pl);
            builder.PL(pl);
        } else {
            builder.noPL();
        }
        swappedGenotypes.add(builder.make());
    }
    swappedBuilder.genotypes(swappedGenotypes);

    for (final String key : vc.getAttributes().keySet()) {
        if (annotationsToDrop.contains(key)) {
            swappedBuilder.rmAttribute(key);
        } else if (annotationsToReverse.contains(key) && !vc.getAttributeAsString(key, "").equals(VCFConstants.MISSING_VALUE_v4)) {
            final double attributeToReverse = vc.getAttributeAsDouble(key, -1);

            if (attributeToReverse < 0 || attributeToReverse > 1) {
                log.warn("Trying to reverse attribute " + key +
                        " but found value that isn't between 0 and 1: (" + attributeToReverse + ") in variant " + vc + ". Results might be wrong.");
            }
            swappedBuilder.attribute(key, 1 - attributeToReverse);
        }
    }

    return swappedBuilder.make();
}
 

public void update1(VariantContext vc, final ReferenceContext referenceContext, final ReadsContext readsContext, final FeatureContext featureContext) {
    if (vc == null || !vc.isSNP() || (getWalker().ignoreAC0Sites() && vc.isMonomorphicInSamples())) {
        return;
    }

    //this maps allele to a count
    ConcurrentMap<String, Integer> alleleCounts = new ConcurrentHashMap<String, Integer>();

    int numHetsHere = 0;
    int numGenosHere = 0;
    int numIndsHere = 0;

    for (final Genotype genotype : vc.getGenotypes()) {
        numIndsHere++;
        if (!genotype.isNoCall()) {
            //increment stats for heterozygosity
            if (genotype.isHet()) {
                numHetsHere++;
            }

            numGenosHere++;
            //increment stats for pairwise mismatches

            for (Allele allele : genotype.getAlleles()) {
                if (allele.isCalled()) {
                    String alleleString = allele.toString();
                    alleleCounts.putIfAbsent(alleleString, 0);
                    alleleCounts.put(alleleString, alleleCounts.get(alleleString) + 1);
                }
            }
        }
    }
    if (numGenosHere > 0) {
        //only if have one called genotype at least
        this.numSites++;

        this.totalHet += numHetsHere / (double)numGenosHere;

        //compute based on num sites
        float harmonicFactor = 0;
        for (int i = 1; i <= numIndsHere; i++) {
            harmonicFactor += 1.0 / i;
        }
        this.thetaRegionNumSites += 1.0 / harmonicFactor;

        //now compute pairwise mismatches
        float numPairwise = 0;
        int numDiffs = 0;
        for (String allele1 : alleleCounts.keySet()) {
            int allele1Count = alleleCounts.get(allele1);

            for (String allele2 : alleleCounts.keySet()) {
                if (allele1.compareTo(allele2) < 0) {
                    continue;
                }
                if (allele1 .compareTo(allele2) == 0) {
                    numPairwise += allele1Count * (allele1Count - 1) * .5;

                }
                else {
                    int allele2Count = alleleCounts.get(allele2);
                    numPairwise += allele1Count * allele2Count;
                    numDiffs += allele1Count * allele2Count;
                }
            }
        }

        if (numPairwise > 0) {
            this.totalAvgDiffs += numDiffs / numPairwise;
        }
    }
}