public boolean isHet() { return vcfGenotype.isHet(); }
/** * Decides whether or not to choose the reference haplotype, depending on the given genotype * * @param g the genotype of the given sample * @param pReferenceGivenHet probability of choosing reference for hets * * @return true if one should use the reference haplotype, false otherwise */ private boolean chooseRefHaplotype(final Genotype g, final double pReferenceGivenHet) { final double refP; if ( g.isHomRef() ) refP = 1; else if ( g.isHet() ) refP = pReferenceGivenHet; else refP = 0.0; return ran.nextDouble() < refP; }
private void updatePairMetricsCounters(Genotype parent, Genotype child, int mvCount, HashMap<Byte,Integer> counters){ //Increment metrics counters if(parent.isCalled() && child.isCalled()){ counters.put(NUM_PAIR_GENOTYPES_CALLED,counters.get(NUM_PAIR_GENOTYPES_CALLED)+1); if(parent.isPhased()) counters.put(NUM_PAIR_GENOTYPES_PHASED,counters.get(NUM_PAIR_GENOTYPES_PHASED)+1); else{ counters.put(NUM_PAIR_VIOLATIONS,counters.get(NUM_PAIR_VIOLATIONS)+mvCount); if(parent.isHet() && child.isHet()) counters.put(NUM_PAIR_HET_HET,counters.get(NUM_PAIR_HET_HET)+1); } }else{ counters.put(NUM_PAIR_GENOTYPES_NOCALL,counters.get(NUM_PAIR_GENOTYPES_NOCALL)+1); } }
private void updateTrioMetricsCounters(Genotype mother, Genotype father, Genotype child, int mvCount, HashMap<Byte,Integer> counters){ //Increment metrics counters if(mother.isCalled() && father.isCalled() && child.isCalled()){ counters.put(NUM_TRIO_GENOTYPES_CALLED,counters.get(NUM_TRIO_GENOTYPES_CALLED)+1); if(mother.isPhased()) counters.put(NUM_TRIO_GENOTYPES_PHASED,counters.get(NUM_TRIO_GENOTYPES_PHASED)+1); else{ if(mvCount > 0){ if(mvCount >1) counters.put(NUM_TRIO_DOUBLE_VIOLATIONS,counters.get(NUM_TRIO_DOUBLE_VIOLATIONS)+1); else counters.put(NUM_TRIO_VIOLATIONS,counters.get(NUM_TRIO_VIOLATIONS)+1); } else if(mother.isHet() && father.isHet() && child.isHet()) counters.put(NUM_TRIO_HET_HET_HET,counters.get(NUM_TRIO_HET_HET_HET)+1); } }else{ counters.put(NUM_TRIO_GENOTYPES_NOCALL,counters.get(NUM_TRIO_GENOTYPES_NOCALL)+1); } }
/** * @return true if variantContext is to be kept, otherwise false * Assumes that this.sample is a sample in the variantContext, if not null, * otherwise looks for the first genotype (and assumes it exists). * @param variantContext the record to examine for heterozygosity */ @Override public boolean test(final VariantContext variantContext) { final Genotype gt = (sample == null) ? variantContext.getGenotype(0) : variantContext.getGenotype(sample); if (gt == null) { throw new IllegalArgumentException((sample == null) ? "Cannot find any genotypes in VariantContext: " + variantContext : "Cannot find sample requested: " + sample); } //XOR operator to reverse behaviour if keepHets is true. return gt.isHet() ^ !keepHets; } }
/** * @return true if variantContext is to be kept, otherwise false * Assumes that this.sample is a sample in the variantContext, if not null, * otherwise looks for the first genotype (and assumes it exists). * @param variantContext the record to examine for heterozygosity */ @Override public boolean test(final VariantContext variantContext) { final Genotype gt = (sample == null) ? variantContext.getGenotype(0) : variantContext.getGenotype(sample); if (gt == null) { throw new IllegalArgumentException((sample == null) ? "Cannot find any genotypes in VariantContext: " + variantContext : "Cannot find sample requested: " + sample); } //XOR operator to reverse behaviour if keepHets is true. return gt.isHet() ^ !keepHets; } }
/** * @return Sample name if there is only one sample that contains alternate allele(s), else null if either multiple samples that * are not homref, or no samples that are not homref. */ protected static String getSingletonSample(final VariantContext vc) { // peek can only change effectively final variables...workaround final String[] sampleName = new String[1]; if (vc.getGenotypes() .stream() // look at het or homVar genotypes .filter(genotype -> genotype.isHet() || genotype.isHomVar()) // two such genotypes will be enough .limit(2) //get any of the sample names .peek(genotype -> sampleName[0] = genotype.getSampleName()) //map to the number of variant chromosomes .mapToInt(genotype -> genotype.isHet() ? 1 : 2) //add them up .reduce(Integer::sum) // compare to 1 with 0 as default .orElse(0) == 1) { return sampleName[0]; } else { return null; } }
/** * @return Sample name if there is only one sample that contains alternate allele(s), else null if either multiple samples that * are not homref, or no samples that are not homref. */ protected static String getSingletonSample(final VariantContext vc) { // peek can only change effectively final variables...workaround final String[] sampleName = new String[1]; if (vc.getGenotypes() .stream() // look at het or homVar genotypes .filter(genotype -> genotype.isHet() || genotype.isHomVar()) // two such genotypes will be enough .limit(2) //get any of the sample names .peek(genotype -> sampleName[0] = genotype.getSampleName()) //map to the number of variant chromosomes .mapToInt(genotype -> genotype.isHet() ? 1 : 2) //add them up .reduce(Integer::sum) // compare to 1 with 0 as default .orElse(0) == 1) { return sampleName[0]; } else { return null; } }
private byte getStandardEncoding(Genotype g, int offset) { byte b; if ( ! checkGQIsGood(g) ) { b = NO_CALL; } else if ( g.isHomRef() ) { b = HOM_REF; } else if ( g.isHomVar() ) { b = HOM_VAR; } else if ( g.isHet() ) { b = HET; } else { b = NO_CALL; } return (byte) (b << (2*offset)); }
private byte getFlippedEncoding(Genotype g, int offset) { byte b; if ( ! checkGQIsGood(g) ) { b = NO_CALL; } else if ( g.isHomRef() ) { b = HOM_VAR; } else if ( g.isHomVar() ) { b = HOM_REF; } else if ( g.isHet() ) { b = HET; } else { b = NO_CALL; } return (byte) (b << (2*offset)); }
private void checkCoverage(Genotype gt, Builder<String> builder, GenotypeFilterImpl impl) { if (gt.isHet()) { if (impl.getCoverage(gt) < options.getMinGtCovHet()) builder.add(ThresholdFilterHeaderExtender.FILTER_GT_MIN_COV_HET); } else if (!gt.isHomRef()) { if (impl.getCoverage(gt) < options.getMinGtCovHomAlt()) builder.add(ThresholdFilterHeaderExtender.FILTER_GT_MIN_COV_HOM_ALT); } if (impl.getCoverage(gt) > options.getMaxCov()) builder.add(ThresholdFilterHeaderExtender.FILTER_GT_MAX_COV); }
@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 checkAaf(Genotype gt, Builder<String> builder, GenotypeFilterImpl impl) { if (gt.isHet()) { if (impl.getAlternativeAlleleFraction(gt) < options.getMinGtAafHet()) builder.add(ThresholdFilterHeaderExtender.FILTER_GT_MIN_AAF_HET); if (impl.getAlternativeAlleleFraction(gt) > options.getMaxGtAafHet()) builder.add(ThresholdFilterHeaderExtender.FILTER_GT_MAX_AAF_HET); } else if (gt.isHomRef()) { if (impl.getAlternativeAlleleFraction(gt) > options.getMaxGtAafHomRef()) builder.add(ThresholdFilterHeaderExtender.FILTER_GT_MAX_AAF_HOM_REF); } else { if (impl.getAlternativeAlleleFraction(gt) < options.getMinGtAafHomAlt()) builder.add(ThresholdFilterHeaderExtender.FILTER_GT_MIN_AAF_HOM_ALT); } }
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; } } }
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; } } }
public void annotate(final RefMetaDataTracker tracker, final AnnotatorCompatible walker, final ReferenceContext ref, final AlignmentContext stratifiedContext, final VariantContext vc, final Genotype g, final GenotypeBuilder gb, final PerReadAlleleLikelihoodMap alleleLikelihoodMap){ // We need a heterozygous genotype and either a context or non-empty alleleLikelihoodMap if ( g == null || !g.isCalled() || !g.isHet() || ( stratifiedContext == null && (alleleLikelihoodMap == null || alleleLikelihoodMap.isEmpty())) ) return; // If we have a <NON_REF> allele the SNP is biallelic if there are 3 alleles and both the reference and first alt allele are length 1. final boolean biallelicSNP = vc.hasAllele(GATKVCFConstants.NON_REF_SYMBOLIC_ALLELE) ? vc.getAlleles().size() == 3 && vc.getReference().length() == 1 && vc.getAlternateAllele(0).length() == 1 : vc.isSNP() && vc.isBiallelic(); if ( !biallelicSNP ) return; final Double ratio = (alleleLikelihoodMap != null && !alleleLikelihoodMap.isEmpty()) ? annotateWithLikelihoods(alleleLikelihoodMap, vc) : annotateWithPileup(stratifiedContext, vc); if (ratio == null) return; gb.attribute(getKeyNames().get(0), Double.valueOf(String.format("%.2f", ratio))); }
/** * Returns the IUPAC encoding for the given genotype or the reference base if not possible * * @param genotype the genotype to encode * @param ref the reference base * @return non-null, non-empty String of bases */ private String getIUPACbase(final Genotype genotype, final String ref) { if ( genotype == null ) throw new IllegalStateException("The genotype is null for sample " + iupacSample); // If have a spanning deletion, if both alleles are spanning deletions, use the empty allele. Otherwise, use the allele is not a // spanning deletion. if ( genotype.getAlleles().contains(Allele.SPAN_DEL) ) { if ( genotype.isHomVar() ) { return EMPTY_BASE; } else { return genotype.getAllele(0).equals(Allele.SPAN_DEL) ? genotype.getAllele(1).getBaseString() : genotype.getAllele(0).getBaseString(); } } if ( !genotype.isHet() ) return genotype.isHom() ? genotype.getAllele(0).getBaseString() : ref; final byte allele1 = genotype.getAllele(0).getBases()[0]; final byte allele2 = genotype.getAllele(1).getBases()[0]; return new String(new byte[] {BaseUtils.basesToIUPAC(allele1, allele2)}); } }
@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); }
public void addFormatGenotype(Genotype gt) { if(gt == null){ this.setText(' '); } else if(gt.isHomRef()){ this.setText('.'); this.setBgColor(Config.get(ConfigKey.feature_background_no_strand)); } else if(gt.isHomVar()){ this.setText('O'); this.setBgColor(Config.get(ConfigKey.feature_background_negative_strand)); } else if(gt.isHet()){ this.setText('E'); this.setBgColor(Config.get(ConfigKey.feature_background_positive_strand)); } else { this.setText('?'); } }