Java GenotypeBuilder.PL属性代码示例

/**
 * Convert a HomRefBlock into a VariantContext
 *
 * @param block the block to convert
 * @return a VariantContext representing the gVCF encoding for this block.
 * It will return {@code null} if input {@code block} is {@code null}, indicating that there
 * is no variant-context to be output into the VCF.
 */
private VariantContext blockToVCF(final HomRefBlock block) {
    if ( block == null ) return null;

    final VariantContextBuilder vcb = new VariantContextBuilder(block.getStartingVC());
    vcb.attributes(new HashMap<String, Object>(2)); // clear the attributes
    vcb.stop(block.getStop());
    vcb.attribute(VCFConstants.END_KEY, block.getStop());

    // create the single Genotype with GQ and DP annotations
    final GenotypeBuilder gb = new GenotypeBuilder(sampleName, GATKVariantContextUtils.homozygousAlleleList(block.getRef(),block.getPloidy()));
    gb.noAD().noPL().noAttributes(); // clear all attributes
    gb.GQ(block.getMedianGQ());
    gb.DP(block.getMedianDP());
    gb.attribute(MIN_DP_FORMAT_FIELD, block.getMinDP());
    gb.PL(block.getMinPLs());

    // This annotation is no longer standard
    //gb.attribute(MIN_GQ_FORMAT_FIELD, block.getMinGQ());

    return vcb.genotypes(gb.make()).make();
}
 

/**
 * From a given genotype, subset the PLs and SACs
 * @param g                                 genotype to subset
 * @param allelesToUse                      alleles to subset
 * @param vc                                variant context with alleles and genotypes
 * @param ploidy                            number of chromosomes
 * @param assignGenotypes                   true: assign hard genotypes, false: leave as no-call
 * @param newLikelihoods                    the PL values
 * @return genotype with the subsetted PLsL and SACs
 */
private Genotype subsetGenotypeAllelesWithLikelihoods(final Genotype g, final List<Allele> allelesToUse, final VariantContext vc, int ploidy,
                                                      final boolean assignGenotypes, final double[] newLikelihoods) {

    final GenotypeBuilder gb = new GenotypeBuilder(g);
    final String sampleName = g.getSampleName();

    // add likelihoods
    gb.PL(newLikelihoods);

    // get and add subsetted SACs
    final int[] newSACs = subsetSACAlleles(g, allelesToUse, vc);
    if (newSACs != null)
        gb.attribute(GaeaVCFConstants.STRAND_COUNT_BY_SAMPLE_KEY, newSACs);
    if (assignGenotypes)
        assignGenotype(gb, vc, sampleName, newLikelihoods, allelesToUse, ploidy);
    else
        gb.alleles(GaeaGvcfVariantContextUtils.noCallAlleles(ploidy));

    return gb.make();
}
 

private VariantContext makeCombinedAltAllelesVariantContext(final VariantContext vc) {
    final int nAltAlleles = vc.getNAlleles() - 1;

    if ( nAltAlleles == 1 )
        return vc;
    else {
        final VariantContextBuilder vcb = new VariantContextBuilder(vc);
        final Allele reference = vcb.getAlleles().get(0);
        vcb.alleles(Arrays.asList(reference, GATKVariantContextUtils.NON_REF_SYMBOLIC_ALLELE));
        final int genotypeCount = GenotypeLikelihoodCalculators.genotypeCount(2, vc.getNAlleles());
        final double[] hetLikelihoods = new double[vc.getNAlleles() - 1];
        final double[] homAltLikelihoods = new double[genotypeCount - hetLikelihoods.length - 1];
        final double[] newLikelihoods = new double[3];
        final List<Genotype> newGenotypes = new ArrayList<>(vc.getNSamples());
        for (final Genotype oldGenotype : vc.getGenotypes()) {
            final GenotypeBuilder gb = new GenotypeBuilder(oldGenotype);
            final List<Allele> oldAlleles = oldGenotype.getAlleles();
            if (oldAlleles != null) {
                final List<Allele> newAlleles = new ArrayList<>(oldAlleles.size());
                for (int i = 0; i < oldAlleles.size(); i++) {
                    final Allele oldAllele = oldAlleles.get(i);
                    if (oldAllele.isReference())
                        newAlleles.add(reference);
                    else if (oldAllele.isNoCall())
                        newAlleles.add(Allele.NO_CALL);
                    else
                        newAlleles.add(GATKVariantContextUtils.NON_REF_SYMBOLIC_ALLELE);
                }
                gb.alleles(newAlleles);
            }
            if (combineAltAlleleLikelihoods(oldGenotype, genotypeCount, newLikelihoods, hetLikelihoods, homAltLikelihoods))
                gb.PL(newLikelihoods);
            newGenotypes.add(gb.make());
        }
        return vcb.genotypes(newGenotypes).make();
    }
}
 

private VariantContext makeCombinedAltAllelesVariantContext(final VariantContext vc) {
    final int nAltAlleles = vc.getNAlleles() - 1;

    if ( nAltAlleles == 1 )
        return vc;
    else {
        final VariantContextBuilder vcb = new VariantContextBuilder(vc);
        final Allele reference = vcb.getAlleles().get(0);
        vcb.alleles(Arrays.asList(reference, GaeaVCFConstants.NON_REF_SYMBOLIC_ALLELE));
        final int genotypeCount = GenotypeLikelihoodCalculators.genotypeCount(2, vc.getNAlleles());
        final double[] hetLikelihoods = new double[vc.getNAlleles() - 1];
        final double[] homAltLikelihoods = new double[genotypeCount - hetLikelihoods.length - 1];
        final double[] newLikelihoods = new double[3];
        final List<Genotype> newGenotypes = new ArrayList<>(vc.getNSamples());
        for (final Genotype oldGenotype : vc.getGenotypes()) {
            final GenotypeBuilder gb = new GenotypeBuilder(oldGenotype);
            final List<Allele> oldAlleles = oldGenotype.getAlleles();
            if (oldAlleles != null) {
                final List<Allele> newAlleles = new ArrayList<>(oldAlleles.size());
                for (int i = 0; i < oldAlleles.size(); i++) {
                    final Allele oldAllele = oldAlleles.get(i);
                    if (oldAllele.isReference())
                        newAlleles.add(reference);
                    else if (oldAllele.isNoCall())
                        newAlleles.add(Allele.NO_CALL);
                    else
                        newAlleles.add(GaeaVCFConstants.NON_REF_SYMBOLIC_ALLELE);
                }
                gb.alleles(newAlleles);
            }
            if (oldGenotype.isNonInformative())
                gb.PL(BIALLELIC_NON_INFORMATIVE_PLS);
            else if (combineAltAlleleLikelihoods(oldGenotype, genotypeCount, newLikelihoods, hetLikelihoods, homAltLikelihoods))
                gb.PL(newLikelihoods);
            newGenotypes.add(gb.make());
        }
        return vcb.genotypes(newGenotypes).make();
    }
}
 

private Genotype createHomRefGenotype(String sampleName) {
    final GenotypeBuilder gb = new GenotypeBuilder(sampleName, Collections.nCopies(getPloidy(), getRef()));
    gb.noAD().noPL().noAttributes(); // clear all attributes

    final int[] minPLs = getMinPLs();
    gb.PL(minPLs);
    gb.GQ(GATKVariantContextUtils.calculateGQFromPLs(minPLs));
    gb.DP(getMedianDP());
    gb.attribute(GATKVCFConstants.MIN_DP_FORMAT_KEY, getMinDP());

    return gb.make();
}
 

@Override
public Genotype entryToObject(TupleInput in)
	{
	GenotypeBuilder gb=new GenotypeBuilder(in.readString());
	if(in.readBoolean()) gb.DP(in.readInt());
	if(in.readBoolean()) gb.AD(arrayOfIntToEntry(in));
	if(in.readBoolean()) gb.GQ(in.readInt());
	if(in.readBoolean()) gb.PL(arrayOfIntToEntry(in));
	
	/* ALLELES ======================================== */
	int n=in.readInt();
	List<Allele> alleles=new ArrayList<Allele>(n);
	for(int i=0;i< n;++i)
		{
		alleles.add(this.alleleBinding.entryToObject(in));
		}
	gb.alleles(alleles);
	/* ATTRIBUTES ===================================== */
	n=in.readInt();
	for(int i=0;i< n;++i)
		{
		String key=in.readString();
		gb.attribute(key, super.readAttribute(in));
		}
	 
	return gb.make();
	}
 

/**
 * Create the new GenotypesContext with the subsetted PLs and ADs
 *
 * @param originalGs               the original GenotypesContext
 * @param vc                       the original VariantContext
 * @param allelesToUse             the actual alleles to use with the new Genotypes
 * @param likelihoodIndexesToUse   the indexes in the PL to use given the allelesToUse (@see #determineLikelihoodIndexesToUse())
 * @param assignGenotypes          assignment strategy for the (subsetted) PLs
 * @return a new non-null GenotypesContext
 */
private static GenotypesContext createGenotypesWithSubsettedLikelihoods(final GenotypesContext originalGs,
                                                                        final VariantContext vc,
                                                                        final List<Allele> allelesToUse,
                                                                        final List<Integer> likelihoodIndexesToUse,
                                                                        final GenotypeAssignmentMethod assignGenotypes) {
    // the new genotypes to create
    final GenotypesContext newGTs = GenotypesContext.create(originalGs.size());

    // make sure we are seeing the expected number of likelihoods per sample
    final int expectedNumLikelihoods = GenotypeLikelihoods.numLikelihoods(vc.getNAlleles(), 2);

    // the samples
    final List<String> sampleIndices = originalGs.getSampleNamesOrderedByName();

    // create the new genotypes
    for ( int k = 0; k < originalGs.size(); k++ ) {
        final Genotype g = originalGs.get(sampleIndices.get(k));
        final GenotypeBuilder gb = new GenotypeBuilder(g);

        // create the new likelihoods array from the alleles we are allowed to use
        double[] newLikelihoods;
        if ( !g.hasLikelihoods() ) {
            // we don't have any likelihoods, so we null out PLs and make G ./.
            newLikelihoods = null;
            gb.noPL();
        } else {
            final double[] originalLikelihoods = g.getLikelihoods().getAsVector();
            if ( likelihoodIndexesToUse == null ) {
                newLikelihoods = originalLikelihoods;
            } else if ( originalLikelihoods.length != expectedNumLikelihoods ) {
                newLikelihoods = null;
            } else {
                newLikelihoods = new double[likelihoodIndexesToUse.size()];
                int newIndex = 0;
                for ( final int oldIndex : likelihoodIndexesToUse )
                    newLikelihoods[newIndex++] = originalLikelihoods[oldIndex];

                // might need to re-normalize
                newLikelihoods = MathUtils.normalizeFromLog10(newLikelihoods, false, true);
            }

            if ( newLikelihoods == null || likelihoodsAreUninformative(newLikelihoods) )
                gb.noPL();
            else
                gb.PL(newLikelihoods);
        }

        updateGenotypeAfterSubsetting(g.getAlleles(), gb, assignGenotypes, newLikelihoods, allelesToUse);
        newGTs.add(gb.make());
    }

    return fixADFromSubsettedAlleles(newGTs, vc, allelesToUse);
}
 

/**
 * From a given variant context, extract a given subset of alleles, and update genotype context accordingly,
 * including updating the PL's, and assign genotypes accordingly
 *
 * @param vc              variant context with alleles and genotype likelihoods
 * @param defaultPloidy   ploidy to assume in case that {@code vc} does not contain that information
 *                        for a sample.
 * @param allelesToUse    alleles to subset
 * @param assignGenotypes true: assign hard genotypes, false: leave as no-call
 * @return GenotypesContext with new PLs
 */
public GenotypesContext subsetAlleles(final VariantContext vc, final int defaultPloidy,
                                      final List<Allele> allelesToUse,
                                      final boolean assignGenotypes) {
    // the genotypes with PLs
    final GenotypesContext oldGTs = vc.getGenotypes();

    // samples
    final List<String> sampleIndices = oldGTs.getSampleNamesOrderedByName();

    // the new genotypes to create
    final GenotypesContext newGTs = GenotypesContext.create();

    // we need to determine which of the alternate alleles (and hence the likelihoods) to use and carry forward
    final int numOriginalAltAlleles = vc.getAlternateAlleles().size();
    final int numNewAltAlleles = allelesToUse.size() - 1;


    // create the new genotypes
    for (int k = 0; k < oldGTs.size(); k++) {
        final Genotype g = oldGTs.get(sampleIndices.get(k));
        final int declaredPloidy = g.getPloidy();
        final int ploidy = declaredPloidy <= 0 ? defaultPloidy : declaredPloidy;
        if (!g.hasLikelihoods()) {
            newGTs.add(GenotypeBuilder.create(g.getSampleName(), GATKVariantContextUtils.noCallAlleles(ploidy)));
            continue;
        }

        // create the new likelihoods array from the alleles we are allowed to use
        final double[] originalLikelihoods = g.getLikelihoods().getAsVector();
        double[] newLikelihoods;

        // Optimization: if # of new alt alleles = 0 (pure ref call), keep original likelihoods so we skip normalization
        // and subsetting
        if (numOriginalAltAlleles == numNewAltAlleles || numNewAltAlleles == 0) {
            newLikelihoods = originalLikelihoods;
        } else {
            newLikelihoods = GeneralPloidyGenotypeLikelihoods.subsetToAlleles(originalLikelihoods, ploidy, vc.getAlleles(), allelesToUse);

            // might need to re-normalize
            newLikelihoods = MathUtils.normalizeFromLog10(newLikelihoods, false, true);
        }

        // if there is no mass on the (new) likelihoods, then just no-call the sample
        if (MathUtils.sum(newLikelihoods) > GATKVariantContextUtils.SUM_GL_THRESH_NOCALL) {
            newGTs.add(GenotypeBuilder.create(g.getSampleName(), GATKVariantContextUtils.noCallAlleles(ploidy)));
        } else {
            final GenotypeBuilder gb = new GenotypeBuilder(g);

            if (numNewAltAlleles == 0)
                gb.noPL();
            else
                gb.PL(newLikelihoods);

            // if we weren't asked to assign a genotype, then just no-call the sample
            if (!assignGenotypes || MathUtils.sum(newLikelihoods) > GATKVariantContextUtils.SUM_GL_THRESH_NOCALL)
                gb.alleles(GATKVariantContextUtils.noCallAlleles(ploidy));
            else
                assignGenotype(gb, newLikelihoods, allelesToUse, ploidy);
            newGTs.add(gb.make());
        }
    }

    return newGTs;

}
 

public VariantContext getLikelihoods(final RefMetaDataTracker tracker,
                                         final ReferenceContext ref,
                                         final Map<String, AlignmentContext> contexts,
                                         final AlignmentContextUtils.ReadOrientation contextType,
                                         final List<Allele> allAllelesToUse,
                                         final boolean useBAQedPileup,
                                         final GenomeLocParser locParser,
                                         final Map<String, PerReadAlleleLikelihoodMap> perReadAlleleLikelihoodMap) {

        GenomeLoc loc = ref.getLocus();
//        if (!ref.getLocus().equals(lastSiteVisited)) {
        if (contextType == AlignmentContextUtils.ReadOrientation.COMPLETE) {
            // starting a new site: clear allele list
            haplotypeMap.clear();
            perReadAlleleLikelihoodMap.clear(); // clean mapping sample-> per read, per allele likelihoods
            alleleList = getInitialAlleleList(tracker, ref, contexts, contextType, UAC, ignoreSNPAllelesWhenGenotypingIndels);
            if (alleleList.isEmpty())
                return null;
        }

        getHaplotypeMapFromAlleles(alleleList, ref, loc, haplotypeMap); // will update haplotypeMap adding elements
        if (haplotypeMap == null || haplotypeMap.isEmpty())
            return null;

        // start making the VariantContext
        // For all non-snp VC types, VC end location is just startLocation + length of ref allele including padding base.
        final int endLoc = loc.getStart() + alleleList.get(0).length() - 1;
        final int eventLength = getEventLength(alleleList);

        final VariantContextBuilder builder = new VariantContextBuilder("UG_call", loc.getContig(), loc.getStart(), endLoc, alleleList);

        // create the genotypes; no-call everyone for now
        GenotypesContext genotypes = GenotypesContext.create();
        final int ploidy = UAC.genotypeArgs.samplePloidy;
        final List<Allele> noCall = GATKVariantContextUtils.noCallAlleles(ploidy);

        // For each sample, get genotype likelihoods based on pileup
        // compute prior likelihoods on haplotypes, and initialize haplotype likelihood matrix with them.

        for (Map.Entry<String, AlignmentContext> sample : contexts.entrySet()) {
            AlignmentContext context = AlignmentContextUtils.stratify(sample.getValue(), contextType);

            if (!perReadAlleleLikelihoodMap.containsKey(sample.getKey())){
                // no likelihoods have been computed for this sample at this site
                perReadAlleleLikelihoodMap.put(sample.getKey(), new PerReadAlleleLikelihoodMap());
            }
            final ReadBackedPileup pileup = context.getBasePileup();
            if (pileup != null) {
                final GenotypeBuilder b = new GenotypeBuilder(sample.getKey());
                final double[] genotypeLikelihoods = pairModel.computeDiploidReadHaplotypeLikelihoods(pileup, haplotypeMap, ref, eventLength, perReadAlleleLikelihoodMap.get(sample.getKey()), UAC.getSampleContamination().get(sample.getKey()));
                b.PL(genotypeLikelihoods);
                b.alleles(noCall);
                b.DP(getFilteredDepth(pileup));
                genotypes.add(b.make());

                if (DEBUG) {
                    System.out.format("Sample:%s Alleles:%s GL:", sample.getKey(), alleleList.toString());
                    for (int k = 0; k < genotypeLikelihoods.length; k++)
                        System.out.format("%1.4f ", genotypeLikelihoods[k]);
                    System.out.println();
                }
            }
        }

        return builder.genotypes(genotypes).make();
    }
 

@Override
@Requires("vc != null && allelesToUse != null")
public GenotypesContext subsetAlleles(VariantContext vc, int defaultPloidy, List<Allele> allelesToUse, boolean assignGenotypes) {
    // the genotypes with PLs
    final GenotypesContext oldGTs = vc.getGenotypes();

    // samples
    final List<String> sampleIndices = oldGTs.getSampleNamesOrderedByName();

    // the new genotypes to create
    final GenotypesContext newGTs = GenotypesContext.create();

    // we need to determine which of the alternate alleles (and hence the likelihoods) to use and carry forward
    final int numOriginalAltAlleles = vc.getAlternateAlleles().size();
    final int numNewAltAlleles = allelesToUse.size() - 1;


    // create the new genotypes
    for ( int k = 0; k < oldGTs.size(); k++ ) {
        final Genotype g = oldGTs.get(sampleIndices.get(k));
        final int declaredPloidy = g.getPloidy();
        final int ploidy = declaredPloidy <= 0 ? defaultPloidy : declaredPloidy;
        if ( !g.hasLikelihoods() ) {
            newGTs.add(GenotypeBuilder.create(g.getSampleName(),GaeaGvcfVariantContextUtils.noCallAlleles(ploidy)));
            continue;
        }

        // create the new likelihoods array from the alleles we are allowed to use
        final double[] originalLikelihoods = g.getLikelihoods().getAsVector();
        double[] newLikelihoods;

        // Optimization: if # of new alt alleles = 0 (pure ref call), keep original likelihoods so we skip normalization
        // and subsetting
        if ( numOriginalAltAlleles == numNewAltAlleles || numNewAltAlleles == 0) {
            newLikelihoods = originalLikelihoods;
        } else {
            newLikelihoods = GeneralPloidyGenotypeLikelihoods.subsetToAlleles(originalLikelihoods, ploidy, vc.getAlleles(), allelesToUse);

            // might need to re-normalize
            newLikelihoods = GvcfMathUtils.normalizeFromLog10(newLikelihoods, false, true);
        }

        // if there is no mass on the (new) likelihoods, then just no-call the sample
        if ( GvcfMathUtils.sum(newLikelihoods) > GaeaGvcfVariantContextUtils.SUM_GL_THRESH_NOCALL ) {
            newGTs.add(GenotypeBuilder.create(g.getSampleName(), GaeaGvcfVariantContextUtils.noCallAlleles(ploidy)));
        } else {
            final GenotypeBuilder gb = new GenotypeBuilder(g);
            final String sampleName = g.getSampleName();

            if ( numNewAltAlleles == 0 )
                gb.noPL();
            else
                gb.PL(newLikelihoods);

            // if we weren't asked to assign a genotype, then just no-call the sample
            if ( !assignGenotypes || GvcfMathUtils.sum(newLikelihoods) > GaeaGvcfVariantContextUtils.SUM_GL_THRESH_NOCALL )
                gb.alleles(GaeaGvcfVariantContextUtils.noCallAlleles(ploidy));
            else
                assignGenotype(gb, vc, sampleName, newLikelihoods, allelesToUse, ploidy);
            newGTs.add(gb.make());
        }
    }

    return GaeaGvcfVariantContextUtils.fixADFromSubsettedAlleles(newGTs, vc, allelesToUse);
}
 

/**
 * Merge into the context a new genotype represented by the given
 * VariantContext for the provided list of target alleles. This method
 * assumes that none of the alleles in the VC overlaps with any of the
 * alleles in the set.
 */
private static void mergeRefConfidenceGenotypes(final GenotypesContext mergedGenotypes, final VariantContext vc,
		final List<Allele> remappedAlleles, final List<Allele> targetAlleles, final boolean samplesAreUniquified,
		final boolean shouldComputePLs) {
	final int maximumPloidy = vc.getMaxPloidy(GaeaGvcfVariantContextUtils.DEFAULT_PLOIDY);
	// the map is different depending on the ploidy, so in order to keep
	// this method flexible (mixed ploidies)
	// we need to get a map done (lazily inside the loop) for each ploidy,
	// up to the maximum possible.
	final int[][] genotypeIndexMapsByPloidy = new int[maximumPloidy + 1][];
	final int maximumAlleleCount = Math.max(remappedAlleles.size(), targetAlleles.size());

	for (final Genotype g : vc.getGenotypes()) {
		final String name;
		if (samplesAreUniquified)
			name = g.getSampleName() + "." + vc.getSource();
		else
			name = g.getSampleName();
		final int ploidy = g.getPloidy();
		final GenotypeBuilder genotypeBuilder = new GenotypeBuilder(g)
				.alleles(GaeaGvcfVariantContextUtils.noCallAlleles(g.getPloidy())).noPL();
		genotypeBuilder.name(name);

		final boolean doPLs = shouldComputePLs && g.hasPL();
		final boolean hasAD = g.hasAD();
		final boolean hasSAC = g.hasExtendedAttribute(GaeaVCFConstants.STRAND_COUNT_BY_SAMPLE_KEY);
		if (doPLs || hasSAC || hasAD) {
			final int[] perSampleIndexesOfRelevantAlleles = getIndexesOfRelevantAlleles(remappedAlleles,
					targetAlleles, vc.getStart(), g);
			if (doPLs) {
				// lazy initialization of the genotype index map by ploidy.

				final int[] genotypeIndexMapByPloidy = genotypeIndexMapsByPloidy[ploidy] == null
						? GenotypeLikelihoodCalculators.getInstance(ploidy, maximumAlleleCount).genotypeIndexMap(
								perSampleIndexesOfRelevantAlleles)
						: genotypeIndexMapsByPloidy[ploidy];
				final int[] PLs = generatePL(g, genotypeIndexMapByPloidy);
				genotypeBuilder.PL(PLs);
			}
			if (hasAD) {
				genotypeBuilder.AD(generateAD(g.getAD(), perSampleIndexesOfRelevantAlleles));
			}
			if (hasSAC) {
				final List<Integer> sacIndexesToUse = adaptToSACIndexes(perSampleIndexesOfRelevantAlleles);
				final int[] SACs = GaeaGvcfVariantContextUtils.makeNewSACs(g, sacIndexesToUse);
				genotypeBuilder.attribute(GaeaVCFConstants.STRAND_COUNT_BY_SAMPLE_KEY, SACs);
			}
		}
		mergedGenotypes.add(genotypeBuilder.make());
	}
}