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

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

/** 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);
}
 

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

@Override
public String filter(final VariantContext ctx, final Genotype gt) {
    if (gt.getGQ() < minGq) return "LowGQ";
    else return null;
}
 

/**
 * 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
}
 

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