public boolean isNoCall() { return vcfGenotype.isNoCall(); }
/** * Find the number of no-call genotypes * * @param vc the variant context * @return number of no-call genotypes */ private int numNoCallGenotypes(final VariantContext vc) { if (vc == null) return 0; int numFiltered = 0; for ( final Genotype g : vc.getGenotypes() ) { if ( g.isNoCall() ) numFiltered++; } return numFiltered; }
static boolean allGenotypesAreUnfilteredAndCalled(VariantContext vc) { for (final Genotype gt : vc.getGenotypes()) { if (gt.isNoCall() || gt.isFiltered()) return false; } return true; }
private double log10PLFromSamples(final VariantContext vc, final String sample, boolean calcRefP) { Genotype g = vc.getGenotype(sample); double log10pSample = -1000; if ( ! g.isNoCall() ) { final double[] gLikelihoods = MathUtils.normalizeFromLog10(g.getLikelihoods().getAsVector()); log10pSample = Math.log10(calcRefP ? gLikelihoods[0] : 1 - gLikelihoods[0]); log10pSample = Double.isInfinite(log10pSample) ? -10000 : log10pSample; } return log10pSample; }
private String formatJsScriptWithGenotype(String js, Genotype gt) { // Comments are copied from htsjdk if(js.contains("{HOM}")){ js= js.replace("{HOM}", Boolean.toString(gt.isHom())); } if(js.contains("{HET}")){ js= js.replace("{HET}", Boolean.toString(gt.isHet())); } if(js.contains("{HOM_REF}")){ js= js.replace("{HOM_REF}", Boolean.toString(gt.isHomRef())); } if(js.contains("{HOM_VAR}")){ js= js.replace("{HOM_VAR}", Boolean.toString(gt.isHomVar())); // true if all observed alleles are alt; if any alleles are no-calls, return false. } if(js.contains("{HET_NON_REF}")){ js= js.replace("{HET_NON_REF}", Boolean.toString(gt.isHetNonRef())); // true if we're het (observed alleles differ) and neither allele is reference; if the ploidy is less than 2 or if any alleles are no-calls, this method will return false. } if(js.contains("{CALLED}")){ js= js.replace("{CALLED}", Boolean.toString(gt.isCalled())); // true if this genotype is comprised of any alleles that are not no-calls (even if some are). } if(js.contains("{NO_CALL}")){ js= js.replace("{NO_CALL}", Boolean.toString(gt.isNoCall())); // true if this genotype is not actually a genotype but a "no call" (e.g. './.' in VCF); if any alleles are not no-calls (even if some are), this method will return false. } if(js.contains("{MIXED}")){ js= js.replace("{MIXED}", Boolean.toString(gt.isMixed())); // true if this genotype is comprised of both calls and no-calls. } return js; }
@Override public String filter(final VariantContext ctx) { if (ctx.getHetCount() == 0) return null; final Map<List<Allele>, Counts> countsMap = new HashMap<List<Allele>, Counts>(); for (final Genotype gt : ctx.getGenotypesOrderedByName()) { if (gt.isNoCall() || !gt.isHet() || !gt.hasAD()) continue; final List<Allele> alleles = gt.getAlleles(); Counts counts = countsMap.get(alleles); if (counts == null) { counts = new Counts(); countsMap.put(alleles, counts); } counts.allele1 += gt.getAD()[ctx.getAlleleIndex(alleles.get(0))]; counts.allele2 += gt.getAD()[ctx.getAlleleIndex(alleles.get(1))]; } for (final Counts counts : countsMap.values()) { final int total = counts.allele1 + counts.allele2; if (total > 0 && Math.min(counts.allele1, counts.allele2) / (double) total < this.hetAlleleBalance) return AB_FILTER; } return null; } }
@Override public String filter(final VariantContext ctx) { if (ctx.getHetCount() == 0) return null; final Map<List<Allele>, Counts> countsMap = new HashMap<List<Allele>, Counts>(); for (final Genotype gt : ctx.getGenotypesOrderedByName()) { if (gt.isNoCall() || !gt.isHet() || !gt.hasAD()) continue; final List<Allele> alleles = gt.getAlleles(); Counts counts = countsMap.get(alleles); if (counts == null) { counts = new Counts(); countsMap.put(alleles, counts); } counts.allele1 += gt.getAD()[ctx.getAlleleIndex(alleles.get(0))]; counts.allele2 += gt.getAD()[ctx.getAlleleIndex(alleles.get(1))]; } for (final Counts counts : countsMap.values()) { final int total = counts.allele1 + counts.allele2; if (total > 0 && Math.min(counts.allele1, counts.allele2) / (double) total < this.hetAlleleBalance) return AB_FILTER; } return null; } }
private void putContig(VariantContext vc, String sampleName) { final Statistics stats = perSampleStats.get(sampleName); if (sampleName == null) { stats.putContig(vc.getContig()); } else { final Genotype gt = vc.getGenotype(sampleName); if (!gt.isHomRef() && !gt.isNoCall()) stats.putContig(vc.getContig()); } }
if (!genotype.isNoCall()) {
/** * Get de novo allele in <code>sampleName</code> or <code>null</code> if there is none * * @param vc {@link VariantContext} to query * @param sampleName Name of the sample * @return De novo allele */ private Allele getDeNovoAllele(VariantContext vc, String sampleName) { final Person person = this.pedigree.getNameToMember().get(sampleName).getPerson(); if (person.getFather() == null || person.getMother() == null) return null; // cannot make any judgement final Genotype gtPerson = vc.getGenotype(sampleName); final Genotype gtFather = vc.getGenotype(person.getFather().getName()); final Genotype gtMother = vc.getGenotype(person.getMother().getName()); if (gtPerson == null || gtPerson.isNoCall() || gtFather == null || gtFather.isNoCall() || gtMother == null || gtMother.isNoCall()) return null; // cannot make any judgement if (!gtPerson.isHet()) return null; // impossible or too unlikely // Count non-reference alleles not yet seen in parents. Should be exactly one. final Set<Allele> personAlleles = new HashSet<>(gtPerson.getAlleles()); personAlleles.remove(vc.getReference()); personAlleles.removeAll(gtFather.getAlleles()); personAlleles.removeAll(gtMother.getAlleles()); if (personAlleles.size() == 1) { return personAlleles.iterator().next(); } else { return null; } }
if( g.isNoCall() || ! g.hasAnyAttribute(GATKVCFConstants.STRAND_BIAS_BY_SAMPLE_KEY) ) continue;
int homCount = 0; for ( final Genotype g : genotypes ) { if ( g.isNoCall() ) continue;
@Requires({"eval != null","truth != null"}) public void update(final VariantContext eval, final VariantContext truth) { overallSiteConcordance.update(eval,truth); final Set<Allele> truthAlleles = new HashSet<>(truth.getAlleles()); for ( final String sample : perSampleGenotypeConcordance.keySet() ) { final Genotype evalGenotype = eval.getGenotype(sample); final Genotype truthGenotype = truth.getGenotype(sample); // ensure genotypes are either no-call ("."), missing (empty alleles), or diploid if ( ( ! evalGenotype.isNoCall() && evalGenotype.getPloidy() != 2 && evalGenotype.getPloidy() > 0) || ( ! truthGenotype.isNoCall() && truthGenotype.getPloidy() != 2 && truthGenotype.getPloidy() > 0) ) { throw new UserException(String.format("Concordance Metrics is currently only implemented for DIPLOID genotypes, found eval ploidy: %d, comp ploidy: %d",evalGenotype.getPloidy(),truthGenotype.getPloidy())); } final boolean allelesMatch = doAllelesMatch(evalGenotype, truthGenotype, truth.getReference(), truthAlleles); perSampleGenotypeConcordance.get(sample).update(allelesMatch, evalGenotype, truthGenotype); final boolean doPrint = overallGenotypeConcordance.update(allelesMatch, evalGenotype, truthGenotype); if(sitesFile != null && doPrint) sitesFile.println(eval.getChr() + ":" + eval.getStart() + "\t" + sample + "\t" + truthGenotype.getType() + "\t" + evalGenotype.getType()); } }
@Test public void testOneStartsBeforeTwoAndEndsAfterwards() throws Exception { final String cmd = baseTestString(" -L 1:69485-69509"); final WalkerTestSpec spec = new WalkerTestSpec(cmd, 1, Arrays.asList("")); spec.disableShadowBCF(); final File gVCF = executeTest("testOneStartsBeforeTwoAndEndsAfterwards", spec).first.get(0); final List<VariantContext> allVCs = GATKVCFUtils.readVCF(gVCF).getSecond(); Assert.assertEquals(allVCs.size(), 2, "Observed: " + allVCs); final VariantContext first = allVCs.get(0); Assert.assertEquals(first.getStart(), 69491); Assert.assertEquals(first.getEnd(), 69497); Assert.assertEquals(first.getGenotypes().size(), 2); Assert.assertTrue(first.getGenotype("NA1").isNoCall()); Assert.assertTrue(first.getGenotype("NA2").isNoCall()); final VariantContext second = allVCs.get(1); Assert.assertEquals(second.getStart(), 69498); Assert.assertEquals(second.getEnd(), 69506); Assert.assertEquals(second.getGenotypes().size(), 2); Assert.assertTrue(second.getGenotype("NA1").isNoCall()); Assert.assertTrue(second.getGenotype("NA2").isNoCall()); }
@Test public void testCreatingPartiallyCalledGenotype() { List<Allele> alleles = Arrays.asList(Aref, C); Genotype g = GenotypeBuilder.create("foo", Arrays.asList(C, Allele.NO_CALL)); VariantContext vc = new VariantContextBuilder("test", snpLoc, snpLocStart, snpLocStop, alleles).genotypes(g).make(); Assert.assertTrue(vc.isSNP()); Assert.assertEquals(vc.getNAlleles(), 2); Assert.assertTrue(vc.hasGenotypes()); Assert.assertFalse(vc.isMonomorphicInSamples()); Assert.assertTrue(vc.isPolymorphicInSamples()); Assert.assertEquals(vc.getGenotype("foo"), g); Assert.assertEquals(vc.getCalledChrCount(), 1); // we only have 1 called chromosomes, we exclude the NO_CALL one isn't called Assert.assertEquals(vc.getCalledChrCount(Aref), 0); Assert.assertEquals(vc.getCalledChrCount(C), 1); Assert.assertFalse(vc.getGenotype("foo").isHet()); Assert.assertFalse(vc.getGenotype("foo").isHom()); Assert.assertFalse(vc.getGenotype("foo").isNoCall()); Assert.assertFalse(vc.getGenotype("foo").isHom()); Assert.assertTrue(vc.getGenotype("foo").isMixed()); Assert.assertEquals(vc.getGenotype("foo").getType(), GenotypeType.MIXED); }
@Test(enabled = !DEBUG, dataProvider = "SplitBiallelics", dependsOnMethods = "testSplitBiallelicsNoGenotypes") public void testSplitBiallelicsGenotypes(final VariantContext vc, final List<VariantContext> expectedBiallelics) { final List<Genotype> genotypes = new ArrayList<Genotype>(); int sampleI = 0; for ( final List<Allele> alleles : Utils.makePermutations(vc.getAlleles(), 2, true) ) { genotypes.add(GenotypeBuilder.create("sample" + sampleI++, alleles)); } genotypes.add(GenotypeBuilder.createMissing("missing", 2)); final VariantContext vcWithGenotypes = new VariantContextBuilder(vc).genotypes(genotypes).make(); final List<VariantContext> biallelics = GATKVariantContextUtils.splitVariantContextToBiallelics(vcWithGenotypes); for ( int i = 0; i < biallelics.size(); i++ ) { final VariantContext actual = biallelics.get(i); Assert.assertEquals(actual.getNSamples(), vcWithGenotypes.getNSamples()); // not dropping any samples for ( final Genotype inputGenotype : genotypes ) { final Genotype actualGenotype = actual.getGenotype(inputGenotype.getSampleName()); Assert.assertNotNull(actualGenotype); if ( ! vc.isVariant() || vc.isBiallelic() ) Assert.assertEquals(actualGenotype, vcWithGenotypes.getGenotype(inputGenotype.getSampleName())); else Assert.assertTrue(actualGenotype.isNoCall()); } } }
private static GenotypeConcordanceStateCodes getStateCode(final VariantContext ctx, final String sample, final int minGq, final int minDp) { // Site level checks // Added potential to include missing sites as hom ref. if (ctx == null) return MISSING_CODE; else if (ctx.isMixed()) return IS_MIXED_CODE; else if (ctx.isFiltered()) return VC_FILTERED_CODE; else { // Genotype level checks final Genotype genotype = ctx.getGenotype(sample); if (genotype.isNoCall()) return NO_CALL_CODE; else if (genotype.isFiltered()) return GT_FILTERED_CODE; else if ((genotype.getGQ() != -1) && (genotype.getGQ() < minGq)) return LOW_GQ_CODE; else if ((genotype.getDP() != -1) && (genotype.getDP() < minDp)) return LOW_DP_CODE; // Note. Genotype.isMixed means that it is called on one chromosome and NOT on the other else if ((genotype.isMixed())) return NO_CALL_CODE; } return null; }
private static GenotypeConcordanceStateCodes getStateCode(final VariantContext ctx, final String sample, final int minGq, final int minDp) { // Site level checks // Added potential to include missing sites as hom ref. if (ctx == null) return MISSING_CODE; else if (ctx.isMixed()) return IS_MIXED_CODE; else if (ctx.isFiltered()) return VC_FILTERED_CODE; else { // Genotype level checks final Genotype genotype = ctx.getGenotype(sample); if (genotype.isNoCall()) return NO_CALL_CODE; else if (genotype.isFiltered()) return GT_FILTERED_CODE; else if ((genotype.getGQ() != -1) && (genotype.getGQ() < minGq)) return LOW_GQ_CODE; else if ((genotype.getDP() != -1) && (genotype.getDP() < minDp)) return LOW_DP_CODE; // Note. Genotype.isMixed means that it is called on one chromosome and NOT on the other else if ((genotype.isMixed())) return NO_CALL_CODE; } return null; }
} else { Assert.assertTrue(g.isNoCall());
if ( ! g.isNoCall() && ! g.isHomRef() ) { countsPerSample.inc(type, g.getSampleName());