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

The following examples show how to use htsjdk.variant.variantcontext.Genotype#getGQ() . 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: CallingMetricAccumulator.java    From picard with MIT License 6 votes vote down vote up 
private void updateDetailMetric(final VariantCallingDetailMetrics metric,
                                final Genotype genotype,
                                final VariantContext vc,
                                final boolean hasSingletonSample) {
    updateSummaryMetric(metric, genotype, vc, hasSingletonSample);

    if (genotype != null && !vc.isFiltered()) {
        if (genotype.getGQ() == 0) {
            ++metric.TOTAL_GQ0_VARIANTS;
        }
        if (genotype.isHet()) {
            ++metric.numHets;
        } else if (genotype.isHomVar()) {
            ++metric.numHomVar;
        }
    }
}
Example 2
Source File: CreateSnpIntervalFromVcf.java    From Drop-seq with MIT License 5 votes vote down vote up 
/** Are all the samples listed heterozygous for this variant, and pass GQ threshold?*/
private boolean sitePassesFilters(final VariantContext ctx, final Set<String> samples, final int GQThreshold, final boolean hetSNPsOnly) {

	boolean flag = true;
	for (String sample : samples) {
		Genotype g = ctx.getGenotype(sample);
		if (g.getGQ()<GQThreshold || (!g.isHet() && hetSNPsOnly))
		 return (false);
	}
	return (flag);
}
Example 3
Source File: HaplotypeCallerIntegrationTest.java    From gatk with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up 
@Test
public void testGVCFModeGenotypePosteriors() throws Exception {
    Utils.resetRandomGenerator();

    final String inputFileName = NA12878_20_21_WGS_bam;
    final String referenceFileName =b37_reference_20_21;

    final File output = createTempFile("testGVCFModeIsConsistentWithPastResults", ".g.vcf");

    final String[] args = {
            "-I", inputFileName,
            "-R", referenceFileName,
            "-L", "20:10000000-10100000",
            "-O", output.getAbsolutePath(),
            "--" + AssemblyBasedCallerArgumentCollection.EMIT_REF_CONFIDENCE_LONG_NAME, ReferenceConfidenceMode.GVCF.toString(),
            "--" + GenotypeCalculationArgumentCollection.SUPPORTING_CALLSET_LONG_NAME,
                largeFileTestDir + "1000G.phase3.broad.withGenotypes.chr20.10100000.vcf",
            "--" + GenotypeCalculationArgumentCollection.NUM_REF_SAMPLES_LONG_NAME, "2500",
            "-pairHMM", "AVX_LOGLESS_CACHING",
            "--" + AssemblyBasedCallerArgumentCollection.ALLELE_EXTENSION_LONG_NAME, "2",
            "--" + StandardArgumentDefinitions.ADD_OUTPUT_VCF_COMMANDLINE, "false"
    };

    runCommandLine(args);

    Pair<VCFHeader, List<VariantContext>> results = VariantContextTestUtils.readEntireVCFIntoMemory(output.getAbsolutePath());

    for (final VariantContext vc : results.getRight()) {
        final Genotype g = vc.getGenotype(0);
        if (g.hasDP() && g.getDP() > 0 && g.hasGQ() && g.getGQ() > 0) {
            Assert.assertTrue(g.hasExtendedAttribute(GATKVCFConstants.PHRED_SCALED_POSTERIORS_KEY));
        }
        if (isGVCFReferenceBlock(vc) ) {
            Assert.assertTrue(!vc.hasAttribute(GATKVCFConstants.GENOTYPE_PRIOR_KEY));
        }
        else if (!vc.getAlternateAllele(0).equals(Allele.NON_REF_ALLELE)){      //there are some variants that don't have non-symbolic alts
            Assert.assertTrue(vc.hasAttribute(GATKVCFConstants.GENOTYPE_PRIOR_KEY));
        }
    }
}
Example 4
Source File: GenotypeQualityFilter.java    From picard with MIT License 4 votes vote down vote up 
@Override
public String filter(final VariantContext ctx, final Genotype gt) {
    if (gt.getGQ() < minGq) return "LowGQ";
    else return null;
}
Example 5
Source File: HomRefBlock.java    From gatk with BSD 3-Clause "New" or "Revised" License 4 votes vote down vote up 
/**
 * Add a homRef block to the current block
 *
 * @param pos current genomic position
 * @param newEnd new calculated block end position
 * @param genotype A non-null Genotype with GQ and DP attributes
 */
@Override
public void add(final int pos, final int newEnd, final Genotype genotype) {
    Utils.nonNull(genotype, "genotype cannot be null");
    if ( ! genotype.hasPL() ) { throw new IllegalArgumentException("genotype must have PL field");}
    if ( pos != end + 1 ) { throw new IllegalArgumentException("adding genotype at pos " + pos + " isn't contiguous with previous end " + end); }
    if ( genotype.getPloidy() != ploidy) { throw new IllegalArgumentException("cannot add a genotype with a different ploidy: " + genotype.getPloidy() + " != " + ploidy); }
    // Make sure the GQ is within the bounds of this band. Treat GQs > 99 as 99.
    if ( !withinBounds(Math.min(genotype.getGQ(), VCFConstants.MAX_GENOTYPE_QUAL))) {
        throw new IllegalArgumentException("cannot add a genotype with GQ=" + genotype.getGQ() + " because it's not within bounds ["
                + this.getGQLowerBound() + ',' + this.getGQUpperBound() + ')');
    }

    if( minPLs == null ) {
        minPLs = genotype.getPL();
    }
    else { // otherwise take the min with the provided genotype's PLs
        final int[] pls = genotype.getPL();
        if (pls.length != minPLs.length) {
            throw new GATKException("trying to merge different PL array sizes: " + pls.length + " != " + minPLs.length);
        }
        for (int i = 0; i < pls.length; i++) {
            minPLs[i] = Math.min(minPLs[i], pls[i]);
        }
    }

    if( genotype.hasExtendedAttribute(GATKVCFConstants.PHRED_SCALED_POSTERIORS_KEY)) {
        if (minPPs == null ) {
            minPPs = PosteriorProbabilitiesUtils.parsePosteriorsIntoPhredSpace(genotype);
        }
        else { // otherwise take the min with the provided genotype's PLs
            final int[] pps = PosteriorProbabilitiesUtils.parsePosteriorsIntoPhredSpace(genotype);
            if (pps.length != minPPs.length) {
                throw new GATKException("trying to merge different PP array sizes: " + pps.length + " != " + minPPs.length);
            }
            for (int i = 0; i < pps.length; i++) {
                minPPs[i] = Math.min(minPPs[i], pps[i]);
            }
        }
    }

    end = newEnd;
    DPs.add(Math.max(genotype.getDP(), 0)); // DP must be >= 0
}
Example 6
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);
    }
}