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

The following examples show how to use htsjdk.variant.variantcontext.Genotype#isHomRef() . You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example 1
Source File: MendelianViolation.java    From gatk with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up 
/**
 * Evaluate the genotypes of mom, dad, and child to detect Mendelian violations
 *
 * @param gMom
 * @param gDad
 * @param gChild
 * @return true if the three genotypes represent a Mendelian violation; false otherwise
 */
public static boolean isViolation(final Genotype gMom, final Genotype gDad, final Genotype gChild) {
    if (gChild.isNoCall()){ //cannot possibly be a violation is child is no call
        return false;
    }
    if(gMom.isHomRef() && gDad.isHomRef() && gChild.isHomRef()) {
        return false;
    }

    //1 parent is no "call
    if(!gMom.isCalled()){
        return (gDad.isHomRef() && gChild.isHomVar()) || (gDad.isHomVar() && gChild.isHomRef());
    }
    else if(!gDad.isCalled()){
        return (gMom.isHomRef() && gChild.isHomVar()) || (gMom.isHomVar() && gChild.isHomRef());
    }
    //Both parents have genotype information
    final Allele childRef = gChild.getAlleles().get(0);
    return !(gMom.getAlleles().contains(childRef) && gDad.getAlleles().contains(gChild.getAlleles().get(1)) ||
        gMom.getAlleles().contains(gChild.getAlleles().get(1)) && gDad.getAlleles().contains(childRef));
}
Example 2
Source File: SampleList.java    From gatk with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up 
@Override
public Map<String, Object> annotate(final ReferenceContext ref,
                                    final VariantContext vc,
                                    final AlleleLikelihoods<GATKRead, Allele> likelihoods) {
    Utils.nonNull(vc);
    if ( vc.isMonomorphicInSamples() || !vc.hasGenotypes() ) {
        return Collections.emptyMap();
    }

    final StringBuilder samples = new StringBuilder();
    for ( final Genotype genotype : vc.getGenotypesOrderedByName() ) {
        if ( genotype.isCalled() && !genotype.isHomRef() ){
            if ( samples.length() > 0 ) {
                samples.append(",");
            }
            samples.append(genotype.getSampleName());
        }
    }

    if ( samples.length() == 0 ) {
        return Collections.emptyMap();
    }

    return Collections.singletonMap(getKeyNames().get(0), samples.toString());
}
Example 3
Source File: ReadOrientationFilter.java    From gatk with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up 
@VisibleForTesting
double artifactProbability(final ReferenceContext referenceContext, final VariantContext vc, final Genotype g) {
    // As of June 2018, genotype is hom ref iff we have the normal sample, but this may change in the future
    // TODO: handle MNVs
    if (g.isHomRef() || (!vc.isSNP() && !vc.isMNP()) ){
        return 0;
    } else if (!artifactPriorCollections.containsKey(g.getSampleName())) {
        return 0;
    }

    final double[] tumorLods = VariantContextGetters.getAttributeAsDoubleArray(vc, GATKVCFConstants.TUMOR_LOG_10_ODDS_KEY, () -> null, -1);
    final int indexOfMaxTumorLod = MathUtils.maxElementIndex(tumorLods);
    final Allele altAllele = vc.getAlternateAllele(indexOfMaxTumorLod);
    final byte[] altBases = altAllele.getBases();

    // for MNVs, treat each base as an independent substitution and take the maximum of all error probabilities
    return IntStream.range(0, altBases.length).mapToDouble(n -> {
        final Nucleotide altBase = Nucleotide.valueOf(new String(new byte[] {altBases[n]}));

        return artifactProbability(referenceContext, vc.getStart() + n, g, indexOfMaxTumorLod, altBase);
    }).max().orElse(0.0);

}
Example 4
Source File: MendelianViolationDetector.java    From picard with MIT License 5 votes vote down vote up 
/** Tests to ensure that a locus is bi-allelic within the given set of samples. */
private boolean isVariant(final Genotype... gts) {
    for (final Genotype gt : gts) {
        if (gt.isCalled() && !gt.isHomRef()) return true;
    }

    return false;
}
Example 5
Source File: GVCFBlockCombiner.java    From gatk with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up 
@Override
public void submit(VariantContext vc) {
    Utils.nonNull(vc);
    Utils.validateArg(vc.hasGenotypes(), "GVCF assumes that the VariantContext has genotypes");
    Utils.validateArg(vc.getGenotypes().size() == 1, () -> "GVCF assumes that the VariantContext has exactly one genotype but saw " + vc.getGenotypes().size());

    if (sampleName == null) {
        sampleName = vc.getGenotype(0).getSampleName();
    }

    if (currentBlock != null && !currentBlock.isContiguous(vc)) {
        // we've made a non-contiguous step (across interval, onto another chr), so finalize
        emitCurrentBlock();
    }

    final Genotype g = vc.getGenotype(0);
    if (g.isHomRef() && vc.hasAlternateAllele(Allele.NON_REF_ALLELE) && vc.isBiallelic()) {
        // create bands
        final VariantContext maybeCompletedBand = addHomRefSite(vc, g);
        if (maybeCompletedBand != null) {
            toOutput.add(maybeCompletedBand);
        }
    } else {
        // g is variant, so flush the bands and emit vc
        emitCurrentBlock();
        nextAvailableStart = vc.getEnd();
        contigOfNextAvailableStart = vc.getContig();
        toOutput.add(vc);
    }
}
Example 6
Source File: VariantSummary.java    From gatk with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up 
public void update2(VariantContext eval, VariantContext comp, ReferenceContext referenceContext, ReadsContext readsContext, FeatureContext featureContext) {
    if ( eval == null || (getWalker().ignoreAC0Sites() && eval.isMonomorphicInSamples()) )
        return;

    final Type type = getType(eval);
    if ( type == null )
        return;

    TypeSampleMap titvTable = null;

    // update DP, if possible
    if ( eval.hasAttribute(VCFConstants.DEPTH_KEY) )
        depthPerSample.inc(type, ALL);

    // update counts
    allVariantCounts.inc(type, ALL);

    // type specific calculations
    if ( type == Type.SNP && eval.isBiallelic() ) {
        titvTable = GATKVariantContextUtils.isTransition(eval) ? transitionsPerSample : transversionsPerSample;
        titvTable.inc(type, ALL);
    }

    // novelty calculation
    if ( comp != null || (type == Type.CNV && overlapsKnownCNV(eval, featureContext)))
        knownVariantCounts.inc(type, ALL);

    // per sample metrics
    for (final Genotype g : eval.getGenotypes()) {
        if ( ! g.isNoCall() && ! g.isHomRef() ) {
            countsPerSample.inc(type, g.getSampleName());

            // update transition / transversion ratio
            if ( titvTable != null ) titvTable.inc(type, g.getSampleName());

            if ( g.hasDP() )
                depthPerSample.inc(type, g.getSampleName());
        }
    }
}
Example 7
Source File: MendelianViolation.java    From gatk with BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up 
protected void updateViolations(final String familyId, final String motherId, final String fatherId, final String childId, final VariantContext vc){
    final Genotype gMom = vc.getGenotype(motherId);
    final Genotype gDad = vc.getGenotype(fatherId);
    final Genotype gChild = vc.getGenotype(childId);

    if (gMom == null || gDad == null || gChild == null){
        if(abortOnSampleNotFound) {
            throw new IllegalArgumentException(String.format("Variant %s:%d: Missing genotypes for family %s: mom=%s dad=%s family=%s", vc.getContig(), vc.getStart(), familyId, motherId, fatherId, childId));
        }
        else {
            return;
        }
    }

    //Count No calls
    if(allCalledOnly && (!gMom.isCalled() || !gDad.isCalled() || !gChild.isCalled())){
        noCall++;
    }
    else if (!gMom.isCalled() && !gDad.isCalled() || !gChild.isCalled()){
        noCall++;
    }
    //Count lowQual. Note that if min quality is set to 0, even values with no quality associated are returned
    else if (minGenotypeQuality > 0 && (
        gMom.getGQ()   < minGenotypeQuality ||
        gDad.getGQ()   < minGenotypeQuality ||
        gChild.getGQ() < minGenotypeQuality )) {
        //no call
        lowQual++;
    }
    else {
        //Count all families per loci called
        familyCalled++;

        if (!(gMom.isHomRef() && gDad.isHomRef() && gChild.isHomRef())) {
            varFamilyCalled++;
        }

        if(isViolation(gMom, gDad, gChild)){
            violationFamilies.add(familyId);
            violations_total++;
        }
        final int count = inheritance.get(gMom.getType()).get(gDad.getType()).get(gChild.getType());
        inheritance.get(gMom.getType()).get(gDad.getType()).put(gChild.getType(),count+1);
    }
}
Example 8
Source File: SelectVariants.java    From gatk with BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up 
private boolean sampleHasVariant(final Genotype g) {
    return (g !=null && !g.isHomRef() && (g.isCalled() || (g.isFiltered() && !XLfiltered)));
}
Example 9
Source File: RMSMappingQuality.java    From gatk with BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up 
private static boolean hasReferenceDepth(Genotype gt) {
    return gt.isHomRef() || (gt.isNoCall() && gt.hasPL() && gt.getPL()[0] == 0 && gt.getPL()[1] != 0);
}
Example 10
Source File: RMSMappingQuality.java    From gatk with BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up 
/**
 *
 * @return the number of reads at the given site, trying first {@Link GATKVCFConstants.RAW_MAPPING_QUALITY_WITH_DEPTH_KEY},
 * falling back to calculating the value as InfoField {@link VCFConstants#DEPTH_KEY} minus the
 * format field {@link GATKVCFConstants#MIN_DP_FORMAT_KEY} or DP of each of the HomRef genotypes at that site.
 * If neither of those is possible, will fall back to calculating the reads from the likelihoods data if provided.
 * @throws UserException.BadInput if the {@link VCFConstants#DEPTH_KEY} is missing or if the calculated depth is <= 0
 */
@VisibleForTesting
static long getNumOfReads(final VariantContext vc,
                         final AlleleLikelihoods<GATKRead, Allele> likelihoods) {
    if(vc.hasAttribute(GATKVCFConstants.RAW_MAPPING_QUALITY_WITH_DEPTH_KEY)) {
        List<Long> mqTuple = VariantContextGetters.getAttributeAsLongList(vc, GATKVCFConstants.RAW_MAPPING_QUALITY_WITH_DEPTH_KEY,0L);
        if (mqTuple.get(TOTAL_DEPTH_INDEX) > 0) {
            return mqTuple.get(TOTAL_DEPTH_INDEX);
        }
    }

    long numOfReads = 0;
    if (vc.hasAttribute(VCFConstants.DEPTH_KEY)) {
        numOfReads = VariantContextGetters.getAttributeAsLong(vc, VCFConstants.DEPTH_KEY, -1L);
        if(vc.hasGenotypes()) {
            for(final Genotype gt : vc.getGenotypes()) {
                if(gt.isHomRef()) {
                    //site-level DP contribution will come from MIN_DP for gVCF-called reference variants or DP for BP resolution
                    if (gt.hasExtendedAttribute(GATKVCFConstants.MIN_DP_FORMAT_KEY)) {
                        numOfReads -= Long.parseLong(gt.getExtendedAttribute(GATKVCFConstants.MIN_DP_FORMAT_KEY).toString());
                    } else if (gt.hasDP()) {
                        numOfReads -= gt.getDP();
                    }
                }
            }
        }

    // If there is no depth key, try to compute from the likelihoods
    } else if (likelihoods != null && likelihoods.numberOfAlleles() != 0) {
        for (int i = 0; i < likelihoods.numberOfSamples(); i++) {
            for (GATKRead read : likelihoods.sampleEvidence(i)) {
                if (read.getMappingQuality() != QualityUtils.MAPPING_QUALITY_UNAVAILABLE) {
                    numOfReads++;
                }
            }
        }
    }
    if (numOfReads <= 0) {
        numOfReads = -1;  //return -1 to result in a NaN
    }
    return numOfReads;
}