public boolean sameGenotype(final Genotype other, boolean ignorePhase) { if (getPloidy() != other.getPloidy()) return false; // gotta have the same number of allele to be equal // By default, compare the elements in the lists of alleles, element-by-element Collection<Allele> thisAlleles = this.getAlleles(); Collection<Allele> otherAlleles = other.getAlleles(); if (ignorePhase) { // do not care about order, only identity of Alleles thisAlleles = new TreeSet<Allele>(thisAlleles); //implemented Allele.compareTo() otherAlleles = new TreeSet<Allele>(otherAlleles); } return thisAlleles.equals(otherAlleles); }
public boolean sameGenotype(final Genotype other, boolean ignorePhase) { if (getPloidy() != other.getPloidy()) return false; // gotta have the same number of allele to be equal // By default, compare the elements in the lists of alleles, element-by-element Collection<Allele> thisAlleles = this.getAlleles(); Collection<Allele> otherAlleles = other.getAlleles(); if (ignorePhase) { // do not care about order, only identity of Alleles thisAlleles = new TreeSet<Allele>(thisAlleles); //implemented Allele.compareTo() otherAlleles = new TreeSet<Allele>(otherAlleles); } return thisAlleles.equals(otherAlleles); }
public boolean sameGenotype(final Genotype other, boolean ignorePhase) { if (getPloidy() != other.getPloidy()) return false; // gotta have the same number of allele to be equal // By default, compare the elements in the lists of alleles, element-by-element Collection<Allele> thisAlleles = this.getAlleles(); Collection<Allele> otherAlleles = other.getAlleles(); if (ignorePhase) { // do not care about order, only identity of Alleles thisAlleles = new TreeSet<Allele>(thisAlleles); //implemented Allele.compareTo() otherAlleles = new TreeSet<Allele>(otherAlleles); } return thisAlleles.equals(otherAlleles); }
/** * Calculates the total ploidy of a variant context as the sum of all ploidies across genotypes. * @param vc the target variant context. * @param defaultPloidy the default ploidy to be assume when there is no ploidy information for a genotype. * @return never {@code null}. */ public static int totalPloidy(final VariantContext vc, final int defaultPloidy) { if (vc == null) throw new IllegalArgumentException("the vc provided cannot be null"); if (defaultPloidy < 0) throw new IllegalArgumentException("the default ploidy must 0 or greater"); int result = 0; for (final Genotype genotype : vc.getGenotypes()) { final int declaredPloidy = genotype.getPloidy(); result += declaredPloidy <= 0 ? defaultPloidy : declaredPloidy; } return result; }
/** * What is the max ploidy among all samples? Returns defaultPloidy if no genotypes are present * * @param defaultPloidy the default ploidy, if all samples are no-called * @return */ public int getMaxPloidy(final int defaultPloidy) { if ( defaultPloidy < 0 ) throw new IllegalArgumentException("defaultPloidy must be greater than or equal to 0"); if ( maxPloidy == -1 ) { maxPloidy = 0; // necessary in the case where there are no genotypes for ( final Genotype g : getGenotypes() ) { maxPloidy = Math.max(g.getPloidy(), maxPloidy); } // everything is no called so we return the default ploidy if ( maxPloidy == 0 ) maxPloidy = defaultPloidy; } return maxPloidy; }
/** * Utility that returns a list of allele strings corresponding to the alleles in this sample * @return */ protected List<String> getAlleleStrings() { final List<String> al = new ArrayList<String>(getPloidy()); for ( Allele a : getAlleles() ) al.add(a.getBaseString()); return al; }
/** * What is the max ploidy among all samples? Returns defaultPloidy if no genotypes are present * * @param defaultPloidy the default ploidy, if all samples are no-called * @return */ public int getMaxPloidy(final int defaultPloidy) { if ( defaultPloidy < 0 ) throw new IllegalArgumentException("defaultPloidy must be greater than or equal to 0"); if ( maxPloidy == -1 ) { maxPloidy = 0; // necessary in the case where there are no genotypes for ( final Genotype g : getGenotypes() ) { maxPloidy = Math.max(g.getPloidy(), maxPloidy); } // everything is no called so we return the default ploidy if ( maxPloidy == 0 ) maxPloidy = defaultPloidy; } return maxPloidy; }
/** * Utility that returns a list of allele strings corresponding to the alleles in this sample * @return */ protected List<String> getAlleleStrings() { final List<String> al = new ArrayList<String>(getPloidy()); for ( Allele a : getAlleles() ) al.add(a.getBaseString()); return al; }
/** * Utility that returns a list of allele strings corresponding to the alleles in this sample * @return */ protected List<String> getAlleleStrings() { final List<String> al = new ArrayList<String>(getPloidy()); for ( Allele a : getAlleles() ) al.add(a.getBaseString()); return al; }
public static boolean isUnfilteredCalledDiploidGenotype(Genotype gt) { return (! gt.isFiltered() && gt.isCalled() && gt.getPloidy() == 2); }
protected List<Haplotype> getAllHaplotypes() { int numSites = genotypes.length; int[] genotypeCards = new int[numSites]; for (int i = 0; i < numSites; i++) genotypeCards[i] = genotypes[i].getPloidy(); LinkedList<Haplotype> allHaps = new LinkedList<Haplotype>(); CardinalityCounter alleleCounter = new CardinalityCounter(genotypeCards); for (int[] alleleInds : alleleCounter) { byte[] hapBases = new byte[numSites]; for (int i = 0; i < numSites; i++) { Allele alleleI = genotypes[i].getAllele(alleleInds[i]); hapBases[i] = SNPallelePair.getSingleBase(alleleI); } allHaps.add(new Haplotype(hapBases)); } return allHaps; }
public AllelePair(Genotype gt) { if (gt.getPloidy() != 2) throw new ReviewedGATKException("AllelePair must have ploidy of 2! incoming gt was"+gt.toBriefString()); this.top = gt.getAllele(0); this.bottom = gt.getAllele(1); }
@Override public void addGenotype(final BCF2Encoder encoder, final VariantContext vc, final Genotype g) throws IOException { final int samplePloidy = g.getPloidy(); for ( int i = 0; i < nValuesPerGenotype; i++ ) { if ( i < samplePloidy ) { // we encode the actual allele final Allele a = g.getAllele(i); final int offset = getAlleleOffset(a); final int encoded = ((offset+1) << 1) | (g.isPhased() ? 0x01 : 0x00); encoder.encodeRawBytes(encoded, encodingType); } else { // we need to pad with missing as we have ploidy < max for this sample encoder.encodeRawBytes(encodingType.getMissingBytes(), encodingType); } } }
@Override public void addGenotype(final BCF2Encoder encoder, final VariantContext vc, final Genotype g) throws IOException { final int samplePloidy = g.getPloidy(); for ( int i = 0; i < nValuesPerGenotype; i++ ) { if ( i < samplePloidy ) { // we encode the actual allele final Allele a = g.getAllele(i); final int offset = getAlleleOffset(a); final int encoded = ((offset+1) << 1) | ((g.isPhased() && i!=0) ? 0x01 : 0x00); encoder.encodeRawBytes(encoded, encodingType); } else { // we need to pad with missing as we have ploidy < max for this sample encoder.encodeRawBytes(encodingType.getMissingBytes(), encodingType); } } }
@Override public void addGenotype(final BCF2Encoder encoder, final VariantContext vc, final Genotype g) throws IOException { final int samplePloidy = g.getPloidy(); for ( int i = 0; i < nValuesPerGenotype; i++ ) { if ( i < samplePloidy ) { // we encode the actual allele final Allele a = g.getAllele(i); final int offset = getAlleleOffset(a); final int encoded = ((offset+1) << 1) | ((g.isPhased() && i!=0) ? 0x01 : 0x00); encoder.encodeRawBytes(encoded, encodingType); } else { // we need to pad with missing as we have ploidy < max for this sample encoder.encodeRawBytes(encodingType.getMissingBytes(), encodingType); } } }
private boolean genotypeCanBeMergedInCurrentBlock(final Genotype g) { return currentBlock != null && currentBlock.withinBounds(capToMaxGQ(g.getGQ())) && currentBlock.getPloidy() == g.getPloidy() && (currentBlock.getMinPLs() == null || !g.hasPL() || (currentBlock.getMinPLs().length == g.getPL().length)); }
@Test(dataProvider = "SubsetAllelesData") public void testSubsetAlleles(final VariantContext inputVC, final List<Allele> allelesToUse, final List<Genotype> expectedGenotypes) { // initialize cache of allele anyploid indices for (final Genotype genotype : inputVC.getGenotypes()) { GenotypeLikelihoods.initializeAnyploidPLIndexToAlleleIndices(inputVC.getNAlleles() - 1, genotype.getPloidy()); } final GenotypesContext actual = GATKVariantContextUtils.subsetAlleles(inputVC, allelesToUse, GATKVariantContextUtils.GenotypeAssignmentMethod.USE_PLS_TO_ASSIGN); Assert.assertEquals(actual.size(), expectedGenotypes.size()); for ( final Genotype expected : expectedGenotypes ) { final Genotype actualGT = actual.get(expected.getSampleName()); Assert.assertNotNull(actualGT); assertGenotypesAreEqual(actualGT, expected); } }
private double[] effectiveAlleleCounts(final VariantContext vc, final double[] log10AlleleFrequencies) { final int numAlleles = vc.getNAlleles(); Utils.validateArg(numAlleles == log10AlleleFrequencies.length, "number of alleles inconsistent"); final double[] log10Result = new double[numAlleles]; Arrays.fill(log10Result, Double.NEGATIVE_INFINITY); for (final Genotype g : vc.getGenotypes()) { if (!g.hasLikelihoods()) { continue; } final GenotypeLikelihoodCalculator glCalc = GL_CALCS.getInstance(g.getPloidy(), numAlleles); final double[] log10GenotypePosteriors = log10NormalizedGenotypePosteriors(g, glCalc, log10AlleleFrequencies); new IndexRange(0, glCalc.genotypeCount()).forEach(genotypeIndex -> glCalc.genotypeAlleleCountsAt(genotypeIndex).forEachAlleleIndexAndCount((alleleIndex, count) -> log10Result[alleleIndex] = MathUtils.log10SumLog10(log10Result[alleleIndex], log10GenotypePosteriors[genotypeIndex] + MathUtils.Log10Cache.get(count)))); } return MathUtils.applyToArrayInPlace(log10Result, x -> Math.pow(10.0, x)); }
/** * Return a VCF-like string representation for the alleles of this genotype. * * If ignoreRefState is true, will not append the reference * marker on the alleles. * * @return a string representing the genotypes, or null if the type is unavailable. */ @Ensures("result != null || ! isAvailable()") public String getGenotypeString(boolean ignoreRefState) { if ( getPloidy() == 0 ) return "NA"; // Notes: // 1. Make sure to use the appropriate separator depending on whether the genotype is phased // 2. If ignoreRefState is true, then we want just the bases of the Alleles (ignoring the '*' indicating a ref Allele) // 3. So that everything is deterministic with regards to integration tests, we sort Alleles (when the genotype isn't phased, of course) return ParsingUtils.join(isPhased() ? PHASED_ALLELE_SEPARATOR : UNPHASED_ALLELE_SEPARATOR, ignoreRefState ? getAlleleStrings() : (isPhased() ? getAlleles() : ParsingUtils.sortList(getAlleles()))); }
@Test(dataProvider = "updatePLsSACsAndADData") public void testUpdatePLsAndADData(final VariantContext originalVC, final VariantContext selectedVC, final List<Genotype> expectedGenotypes) { // initialize cache of allele anyploid indices for (final Genotype genotype : originalVC.getGenotypes()) { GenotypeLikelihoods.initializeAnyploidPLIndexToAlleleIndices(originalVC.getNAlleles() - 1, genotype.getPloidy()); } final VariantContext selectedVCwithGTs = new VariantContextBuilder(selectedVC).genotypes(originalVC.getGenotypes()).make(); final GenotypesContext actual = GATKVariantContextUtils.updatePLsSACsAD(selectedVCwithGTs, originalVC); Assert.assertEquals(actual.size(), expectedGenotypes.size()); for ( final Genotype expected : expectedGenotypes ) { final Genotype actualGT = actual.get(expected.getSampleName()); Assert.assertNotNull(actualGT); assertGenotypesAreEqual(actualGT, expected); } }