Java GenotypeBuilder.make方法代码示例

import htsjdk.variant.variantcontext.GenotypeBuilder; //导入方法依赖的package包/类
/**
 * Fix the AD for the given Genotype
 *
 * @param genotype              the original Genotype
 * @param alleleIndexesToUse    a bitset describing whether or not to keep a given index
 * @param nAllelesToUse         how many alleles we are keeping
 * @return a non-null Genotype
 */
private static Genotype fixAD(final Genotype genotype, final boolean[] alleleIndexesToUse, final int nAllelesToUse) {
    // if it ain't broke don't fix it
    if ( !genotype.hasAD() )
        return genotype;

    final GenotypeBuilder builder = new GenotypeBuilder(genotype);

    final int[] oldAD = genotype.getAD();
    if ( oldAD.length != alleleIndexesToUse.length ) {
        builder.noAD();
    } else {
        final int[] newAD = new int[nAllelesToUse];
        int currentIndex = 0;
        for ( int i = 0; i < oldAD.length; i++ ) {
            if ( alleleIndexesToUse[i] )
                newAD[currentIndex++] = oldAD[i];
        }
        builder.AD(newAD);
    }
    return builder.make();
}
 

import htsjdk.variant.variantcontext.GenotypeBuilder; //导入方法依赖的package包/类
/**
 * 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();
}
 

import htsjdk.variant.variantcontext.GenotypeBuilder; //导入方法依赖的package包/类
private Genotype composeNonTruthOverlappingGenotype(final VariantContext enclosingContext, final Genotype genotype) {
    final GenotypeBuilder builder = new GenotypeBuilder(genotype.getSampleName());
    if (genotype.isCalled()) {
        GATKProtectedVariantContextUtils.setGenotypeQualityFromPLs(builder, genotype);
        final int[] PL = genotype.getPL();
        final int callAlleleIndex = GATKProtectedMathUtils.minIndex(PL);
        final double quality = callQuality(genotype);
        builder.alleles(Collections.singletonList(enclosingContext.getAlleles().get(callAlleleIndex)));
        builder.attribute(VariantEvaluationContext.CALL_QUALITY_KEY, quality);
        final boolean discovered = XHMMSegmentGenotyper.DISCOVERY_TRUE.equals(
                GATKProtectedVariantContextUtils.getAttributeAsString(genotype, XHMMSegmentGenotyper.DISCOVERY_KEY,
                        XHMMSegmentGenotyper.DISCOVERY_FALSE));
        if (callAlleleIndex != 0 && discovered) {
            builder.attribute(VariantEvaluationContext.EVALUATION_CLASS_KEY, EvaluationClass.UNKNOWN_POSITIVE.acronym);
        }
        if (quality < filterArguments.minimumCalledSegmentQuality) {
            builder.filter(EvaluationFilter.LowQuality.acronym);
        } else {
            builder.filter(EvaluationFilter.PASS);
        }
    } else { /* assume it is REF */
        /* TODO this is a hack to make Andrey's CODEX vcf work; and in general, VCFs that only include discovered
         * variants and NO_CALL (".") on other samples. The idea is to force the evaluation tool to take it call
         * as REF on all other samples. Otherwise, the effective allele frequency of the variant will be erroneously
         * high and will be filtered. */
        builder.alleles(Collections.singletonList(CopyNumberTriStateAllele.REF));
        builder.attribute(VariantEvaluationContext.CALL_QUALITY_KEY, 100000);
        builder.filter(EvaluationFilter.PASS);
    }
    return builder.make();
}
 

import htsjdk.variant.variantcontext.GenotypeBuilder; //导入方法依赖的package包/类
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();
}
 

import htsjdk.variant.variantcontext.GenotypeBuilder; //导入方法依赖的package包/类
/** Augment genotype with the given filters and return modified GenotypeBuilder */
public Genotype gtWithAppliedFilters(Genotype gt) {
	GenotypeBuilder gtBuilder = new GenotypeBuilder(gt);

	ArrayList<String> filters = new ArrayList<>();
	if (gt.isFiltered())
		filters.add(gt.getFilters());
	filters.addAll(getFiltersFor(gt));
	gtBuilder.filters(filters);

	return gtBuilder.make();
}
 

import htsjdk.variant.variantcontext.GenotypeBuilder; //导入方法依赖的package包/类
public Genotype makeGenotype() {
GenotypeBuilder gb= new GenotypeBuilder(this.g);
String fmt = getReview(); if(fmt==null) fmt="";
fmt= fmt.trim().replaceAll("[ \n\t\\:_]+","_");
  	gb.attribute(IgvReview.this.reviewFormat.getID(), fmt);
  	return gb.make();
}
 

import htsjdk.variant.variantcontext.GenotypeBuilder; //导入方法依赖的package包/类
@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();
	}
 

import htsjdk.variant.variantcontext.GenotypeBuilder; //导入方法依赖的package包/类
private Genotype buildAndAnnotateTruthOverlappingGenotype(final String sample,
                                                          final TargetCollection<Target> targets,
                                                          final Genotype truthGenotype,
                                                          final int truthCopyNumber,
                                                          final CopyNumberTriStateAllele truthAllele,
                                                          final List<Pair<VariantContext, Genotype>> calls) {
    final Set<CopyNumberTriStateAllele> calledAlleles = calls.stream()
            .map(pair -> CopyNumberTriStateAllele.valueOf(pair.getRight().getAllele(0)))
            .collect(Collectors.toSet());

    final Allele calledAllele = calledAlleles.size() == 1 ? calledAlleles.iterator().next() : Allele.NO_CALL;
    final GenotypeBuilder builder = new GenotypeBuilder(sample);

    // Set the call allele.
    builder.alleles(Collections.singletonList(calledAllele));

    // Set the truth allele.
    builder.attribute(VariantEvaluationContext.TRUTH_GENOTYPE_KEY, CopyNumberTriStateAllele.ALL_ALLELES.indexOf(truthAllele));

    // Annotate the genotype with the number of calls.
    builder.attribute(VariantEvaluationContext.CALLED_SEGMENTS_COUNT_KEY, calls.size());

    // When there is more than one qualified type of event we indicate how many.
    builder.attribute(VariantEvaluationContext.CALLED_ALLELE_COUNTS_KEY,
            CopyNumberTriStateAllele.ALL_ALLELES.stream()
                    .mapToInt(allele -> (int) calls.stream()
                            .filter(pair -> pair.getRight().getAllele(0).equals(allele, true))
                            .count())
                    .toArray());

    // Calculate the length in targets of the call as the sum across all calls.
    builder.attribute(VariantEvaluationContext.CALLED_TARGET_COUNT_KEY,
            calls.stream().mapToInt(pair -> getTargetCount(targets, pair.getLeft(), pair.getRight()))
                    .sum());

    // Calculate call quality-- if there is more than one overlapping call we take the maximum qual one.
    builder.attribute(VariantEvaluationContext.CALL_QUALITY_KEY,
            calls.stream().mapToDouble(pair -> GATKProtectedVariantContextUtils.calculateGenotypeQualityFromPLs(pair.getRight())).max().orElse(0.0));

    // Calculate the truth copy fraction.
    builder.attribute(VariantEvaluationContext.TRUTH_COPY_FRACTION_KEY,
            truthGenotype.getExtendedAttribute(GS_COPY_NUMBER_FRACTION_KEY));

    // Calculate the truth call quality.
    final double truthQuality = calculateTruthQuality(truthGenotype, truthCopyNumber);
    builder.attribute(VariantEvaluationContext.TRUTH_QUALITY_KEY, truthQuality);

    // Set genotype filters:
    final boolean truthPassQualityMinimum = truthQuality >= filterArguments.minimumTruthSegmentQuality;
    builder.filter(truthPassQualityMinimum ? EvaluationFilter.PASS
                : EvaluationFilter.LowQuality.acronym);

    // Calculate the evaluation class (TP, FN, etc.). Only if there is actually either a truth or a call that is not ref.
    if (calledAlleles.contains(CopyNumberTriStateAllele.DEL) || calledAlleles.contains(CopyNumberTriStateAllele.DUP) || truthAllele != CopyNumberTriStateAllele.REF) {
        final EvaluationClass evaluationClass;
        if (calledAlleles.isEmpty() || (calledAlleles.size() == 1 && calledAlleles.contains(CopyNumberTriStateAllele.REF))) {
            evaluationClass = EvaluationClass.FALSE_NEGATIVE;
        } else if (calledAlleles.size() == 1) {
            evaluationClass = calledAlleles.contains(truthAllele) ? EvaluationClass.TRUE_POSITIVE :
                    truthAllele == CopyNumberTriStateAllele.REF ? EvaluationClass.FALSE_POSITIVE :
                          /* else */ EvaluationClass.DISCORDANT_POSITIVE;
        } else {
            evaluationClass = truthAllele == CopyNumberTriStateAllele.REF ? EvaluationClass.FALSE_POSITIVE
                    : EvaluationClass.MIXED_POSITIVE;
        }
        builder.attribute(VariantEvaluationContext.EVALUATION_CLASS_KEY, evaluationClass.acronym);
    }
    return builder.make();
}
 

import htsjdk.variant.variantcontext.GenotypeBuilder; //导入方法依赖的package包/类
@Override
protected int doVcfToVcf(String inputName, VcfIterator in, VariantContextWriter out) 
	{
	final List<Allele> empty_g=new ArrayList<Allele>(2);
	empty_g.add(Allele.NO_CALL);
	empty_g.add(Allele.NO_CALL);
	VCFHeader h=in.getHeader();
	final List<String> vcf_samples=h.getSampleNamesInOrder();
	h.addMetaDataLine(new VCFFormatHeaderLine("GR", 1, VCFHeaderLineType.Integer, "(1) = Genotype was reset by "+getProgramName()));
	out.writeHeader(h);
	while(in.hasNext())
		{
		VariantContext ctx=in.next();
		ContigPosRef cap=new ContigPosRef(ctx);
		Set<String> samplesToReset=this.pos2sample.get(cap);
		if(samplesToReset!=null)
			{
			VariantContextBuilder vcb=new VariantContextBuilder(ctx);
			List<Genotype> genotypes=new ArrayList<Genotype>();
			for(String sample:vcf_samples)
				{
				Genotype g=ctx.getGenotype(sample);
				if(g==null) continue;
				GenotypeBuilder gb=new GenotypeBuilder(g);
				
				if(samplesToReset.contains(sample))
					{
					gb.alleles(empty_g);
					gb.attribute("GR",1);
					}
				else
					{
					gb.attribute("GR",0);
					}
				g=gb.make();
				genotypes.add(g);
				
				}
			vcb.genotypes(genotypes);
			ctx=VCFUtils.recalculateAttributes(vcb.make());
			}
		out.add(ctx);
		}
	return 0;
	}