Java Code Examples for htsjdk.variant.variantcontext.GenotypeBuilder#create()

private Genotype buildAndAnnotateTruthOverlappingGenotype(final String sample, final VariantContext truth, final List<VariantContext> calls,
                                                          final TargetCollection<Target> targets) {
    final Genotype truthGenotype = truth.getGenotype(sample);
    // if there is no truth genotype for that sample, we output the "empty" genotype.
    if (truthGenotype == null) {
        return GenotypeBuilder.create(sample, Collections.emptyList());
    }
    final int truthCopyNumber = GATKProtectedVariantContextUtils.getAttributeAsInt(truthGenotype,
            GS_COPY_NUMBER_FORMAT_KEY, truthNeutralCopyNumber);
    final CopyNumberTriStateAllele truthAllele = copyNumberToTrueAllele(truthCopyNumber);

    final List<Pair<VariantContext, Genotype>> allCalls = calls.stream()
            .map(vc -> new ImmutablePair<>(vc, vc.getGenotype(sample)))
            .filter(pair -> pair.getRight() != null)
            .filter(pair -> GATKProtectedVariantContextUtils.getAttributeAsString(pair.getRight(), XHMMSegmentGenotyper.DISCOVERY_KEY,
                    XHMMSegmentGenotyper.DISCOVERY_FALSE).equals(XHMMSegmentGenotyper.DISCOVERY_TRUE))
            .collect(Collectors.toList());

    final List<Pair<VariantContext, Genotype>> qualifiedCalls = composeQualifyingCallsList(targets, allCalls);

    return buildAndAnnotateTruthOverlappingGenotype(sample, targets, truthGenotype, truthCopyNumber,
                truthAllele, qualifiedCalls);
}
 

@Test
public void testGenotypeConcordanceDetermineStateDp() throws Exception {
    final List<Allele> allelesNormal = makeUniqueListOfAlleles(Aref, C);
    final Genotype gtNormal = GenotypeBuilder.create(TRUTH_SAMPLE_NAME, Arrays.asList(Aref, C));
    final VariantContext vcNormal = new VariantContextBuilder("test", snpLoc, snpLocStart, snpLocStop, allelesNormal).genotypes(gtNormal).make();

    final List<Allele> allelesLowDp = makeUniqueListOfAlleles(Aref, C);
    final Genotype gtLowDp = new GenotypeBuilder(TRUTH_SAMPLE_NAME, Arrays.asList(Aref, C)).DP(4).make();
    final VariantContext vcLowDp = new VariantContextBuilder("test", snpLoc, snpLocStart, snpLocStop, allelesLowDp).genotypes(gtLowDp).make();

    testGenotypeConcordanceDetermineState(vcLowDp, TruthState.LOW_DP, vcNormal, CallState.HET_REF_VAR1, 0, 20);
    testGenotypeConcordanceDetermineState(vcLowDp, TruthState.HET_REF_VAR1, vcLowDp, CallState.HET_REF_VAR1, 0, 2);

    testGenotypeConcordanceDetermineState(vcNormal, TruthState.HET_REF_VAR1, vcLowDp, CallState.LOW_DP, 0, 20);
    testGenotypeConcordanceDetermineState(vcNormal, TruthState.HET_REF_VAR1, vcLowDp, CallState.HET_REF_VAR1, 0, 2);

    testGenotypeConcordanceDetermineState(vcLowDp, TruthState.LOW_DP, vcLowDp, CallState.LOW_DP, 0, 20);
    testGenotypeConcordanceDetermineState(vcLowDp, TruthState.HET_REF_VAR1, vcLowDp, CallState.HET_REF_VAR1, 0, 2);
}
 

@Test
public void testGenotypeConcordanceDetermineStateGq() throws Exception {
    final List<Allele> allelesNormal = makeUniqueListOfAlleles(Aref, C);
    final Genotype gtNormal = GenotypeBuilder.create(TRUTH_SAMPLE_NAME, Arrays.asList(Aref, C));
    final VariantContext vcNormal = new VariantContextBuilder("test", snpLoc, snpLocStart, snpLocStop, allelesNormal).genotypes(gtNormal).make();

    final List<Allele> allelesLowGq = makeUniqueListOfAlleles(Aref, C);
    final Genotype gtLowGq = new GenotypeBuilder(TRUTH_SAMPLE_NAME, Arrays.asList(Aref, C)).GQ(4).make();
    final VariantContext vcLowGq = new VariantContextBuilder("test", snpLoc, snpLocStart, snpLocStop, allelesLowGq).genotypes(gtLowGq).make();

    testGenotypeConcordanceDetermineState(vcLowGq, TruthState.LOW_GQ, vcNormal, CallState.HET_REF_VAR1, 20, 0);
    testGenotypeConcordanceDetermineState(vcLowGq, TruthState.HET_REF_VAR1, vcLowGq, CallState.HET_REF_VAR1, 2, 0);

    testGenotypeConcordanceDetermineState(vcNormal, TruthState.HET_REF_VAR1, vcLowGq, CallState.LOW_GQ, 20, 0);
    testGenotypeConcordanceDetermineState(vcNormal, TruthState.HET_REF_VAR1, vcLowGq, CallState.HET_REF_VAR1, 2, 0);

    testGenotypeConcordanceDetermineState(vcLowGq, TruthState.LOW_GQ, vcLowGq, CallState.LOW_GQ, 20, 0);
    testGenotypeConcordanceDetermineState(vcLowGq, TruthState.HET_REF_VAR1, vcLowGq, CallState.HET_REF_VAR1, 2, 0);
}
 

@DataProvider(name="basicTransitions")
public Object [] [] basicTransitions() {

    final List<Allele> ctoTAlleles = new ArrayList<>();
    ctoTAlleles.add(Allele.create("C", true));
    ctoTAlleles.add(Allele.create("T", false));

    final List<Allele> gtoTAlleles = new ArrayList<>();
    gtoTAlleles.add(Allele.create("G", true));
    gtoTAlleles.add(Allele.create("T", false));

    final List<Allele> gtoAAlleles = new ArrayList<>();
    gtoAAlleles.add(Allele.create("G", true));
    gtoAAlleles.add(Allele.create("A", false));

    final List<Allele> gtoCAlleles = new ArrayList<>();
    gtoCAlleles.add(Allele.create("G", true));
    gtoCAlleles.add(Allele.create("C", false));

    return new Object[] [] {
            // Genotype, transition, is in transition, is in transition or transition rc
            {GenotypeBuilder.create("DUMMYSAMPLE", ctoTAlleles), Transition.CtoT, true, true },
            {GenotypeBuilder.create("DUMMYSAMPLE", gtoTAlleles), Transition.CtoT, false, false },
            {GenotypeBuilder.create("DUMMYSAMPLE", gtoAAlleles), Transition.CtoT, false, true },
            {GenotypeBuilder.create("DUMMYSAMPLE", gtoCAlleles), Transition.CtoT, false, false },
            {GenotypeBuilder.create("DUMMYSAMPLE", ctoTAlleles), Transition.AtoT, false, false },
            {GenotypeBuilder.create("DUMMYSAMPLE", gtoTAlleles), Transition.AtoT, false, false },
            {GenotypeBuilder.create("DUMMYSAMPLE", gtoAAlleles), Transition.AtoT, false, false },
            {GenotypeBuilder.create("DUMMYSAMPLE", gtoCAlleles), Transition.AtoT, false, false }
    };
}
 

public Genotype getGenotype(final String sampleName,
                            final InfiniumGTCRecord infiniumGtcRecord,
                            final IlluminaManifestRecord record,
                            final Allele A,
                            final Allele B) {

    // The Sample Alleles
    final List<Allele> alleles;

    if (infiniumGtcRecord.genotype == InfiniumGTCFile.NO_CALL) alleles = NO_CALL_ALLELES;
    else if (infiniumGtcRecord.genotype == InfiniumGTCFile.AA_CALL) alleles = Arrays.asList(A, A);
    else if (infiniumGtcRecord.genotype == InfiniumGTCFile.AB_CALL) alleles = Arrays.asList(A, B);
    else if (infiniumGtcRecord.genotype == InfiniumGTCFile.BB_CALL) alleles = Arrays.asList(B, B);
    else {
        throw new PicardException("Unexpected genotype call [" + infiniumGtcRecord.genotype + "]" + " for SNP: " + record.getName());
    }

    final Map<String, Object> attributes = new HashMap<>();
    attributes.put(InfiniumVcfFields.IGC, formatFloatForVcf(infiniumGtcRecord.genotypeScore));
    attributes.put(InfiniumVcfFields.X, infiniumGtcRecord.rawXIntensity);
    attributes.put(InfiniumVcfFields.Y, infiniumGtcRecord.rawYIntensity);
    attributes.put(InfiniumVcfFields.NORMX, formatFloatForVcf(infiniumGtcRecord.normalizedXIntensity));
    attributes.put(InfiniumVcfFields.NORMY, formatFloatForVcf(infiniumGtcRecord.normalizedYIntensity));
    attributes.put(InfiniumVcfFields.R, formatFloatForVcf(infiniumGtcRecord.RIlmn));
    attributes.put(InfiniumVcfFields.THETA, formatFloatForVcf(infiniumGtcRecord.thetaIlmn));
    attributes.put(InfiniumVcfFields.BAF, formatFloatForVcf(infiniumGtcRecord.bAlleleFreq));
    attributes.put(InfiniumVcfFields.LRR, formatFloatForVcf(infiniumGtcRecord.logRRatio));

    return GenotypeBuilder.create(sampleName, alleles, attributes);
}
 

@Test(dataProvider = "genotypeConcordanceDetermineStateDataProvider")
public void testGenotypeConcordanceDetermineState(final Allele truthAllele1, final Allele truthAllele2, final TruthState expectedTruthState,
                                                  final Allele callAllele1, final Allele callAllele2, final CallState expectedCallState) throws Exception {
    final List<Allele> truthAlleles = makeUniqueListOfAlleles(truthAllele1, truthAllele2);
    final Genotype truthGt = GenotypeBuilder.create(TRUTH_SAMPLE_NAME, Arrays.asList(truthAllele1, truthAllele2));

    final VariantContext truthVariantContext = new VariantContextBuilder("test", snpLoc, snpLocStart, snpLocStop, truthAlleles).genotypes(truthGt).make();

    final List<Allele> callAlleles = makeUniqueListOfAlleles(callAllele1, callAllele2);
    final Genotype callGt = GenotypeBuilder.create(CALL_SAMPLE_NAME, Arrays.asList(callAllele1, callAllele2));
    final VariantContext callVariantContext = new VariantContextBuilder("test", snpLoc, snpLocStart, snpLocStop, callAlleles).genotypes(callGt).make();

    testGenotypeConcordanceDetermineState(truthVariantContext, expectedTruthState, callVariantContext, expectedCallState, 0, 0);
}
 

@Test
public void testGenotypeConcordanceDetermineStateNull() throws Exception {
    final List<Allele> alleles = makeUniqueListOfAlleles(Aref, C);
    final Genotype gt1 = GenotypeBuilder.create(TRUTH_SAMPLE_NAME, Arrays.asList(Aref, C));
    final VariantContext vc1 = new VariantContextBuilder("test", snpLoc, snpLocStart, snpLocStop, alleles).genotypes(gt1).make();

    testGenotypeConcordanceDetermineState(null, TruthState.MISSING, null, CallState.MISSING, 0, 0);
    testGenotypeConcordanceDetermineState(vc1, TruthState.HET_REF_VAR1, null, CallState.MISSING, 0, 0);
    testGenotypeConcordanceDetermineState(null, TruthState.MISSING, vc1, CallState.HET_REF_VAR1, 0, 0);
}
 

@Test
public void testGenotypeConcordanceDetermineStateFilter() throws Exception {
    final Set<String> filters = new HashSet<String>(Arrays.asList("BAD!"));

    // Filtering on the variant context
    final List<Allele> alleles1 = makeUniqueListOfAlleles(Aref, C);
    final Genotype gt1 = GenotypeBuilder.create(TRUTH_SAMPLE_NAME, Arrays.asList(Aref, C));
    final VariantContext vcFiltered = new VariantContextBuilder("test", snpLoc, snpLocStart, snpLocStop, alleles1).genotypes(gt1).filters(filters).make();

    final List<Allele> alleles2 = makeUniqueListOfAlleles(Aref, T);
    final Genotype gt2 = GenotypeBuilder.create(TRUTH_SAMPLE_NAME, Arrays.asList(Aref, T));
    final VariantContext vcNotFiltered = new VariantContextBuilder("test", snpLoc, snpLocStart, snpLocStop, alleles2).genotypes(gt2).make();

    testGenotypeConcordanceDetermineState(vcFiltered, TruthState.VC_FILTERED, vcNotFiltered, CallState.HET_REF_VAR1, 0, 0);
    testGenotypeConcordanceDetermineState(vcNotFiltered, TruthState.HET_REF_VAR1, vcFiltered, CallState.VC_FILTERED, 0, 0);
    testGenotypeConcordanceDetermineState(vcFiltered, TruthState.VC_FILTERED, vcFiltered, CallState.VC_FILTERED, 0, 0);

    // Filtering on the genotype
    final List<String> gtFilters = new ArrayList<String>(Arrays.asList("WICKED"));
    final List<Allele> alleles3 = makeUniqueListOfAlleles(Aref, C);
    final Genotype gt3 = new GenotypeBuilder(TRUTH_SAMPLE_NAME, Arrays.asList(Aref, C)).filters(gtFilters).make();
    final VariantContext vcGtFiltered = new VariantContextBuilder("test", snpLoc, snpLocStart, snpLocStop, alleles3).genotypes(gt3).make();

    testGenotypeConcordanceDetermineState(vcGtFiltered, TruthState.GT_FILTERED, vcNotFiltered, CallState.HET_REF_VAR1, 0, 0);
    testGenotypeConcordanceDetermineState(vcNotFiltered, TruthState.HET_REF_VAR1, vcGtFiltered, CallState.GT_FILTERED, 0, 0);
    testGenotypeConcordanceDetermineState(vcGtFiltered, TruthState.GT_FILTERED, vcGtFiltered, CallState.GT_FILTERED, 0, 0);
}
 

@Test
public void testFastGenotypeIsSerializable() {
    Genotype genotype = GenotypeBuilder.create("sample1", Collections.nCopies(2, Allele.create("C", false)));
    SparkTestUtils.roundTripInKryo(genotype, genotype.getClass(), SparkContextFactory.getTestSparkContext().getConf());
}
 

/**
 * Writes out the VariantContext objects in the order presented to the supplied output file
 * in VCF format.
 */
private void writeVcf(final CloseableIterator<VariantContext> variants,
                      final File output,
                      final SAMSequenceDictionary dict,
                      final VCFHeader vcfHeader, ZCallPedFile zCallPedFile) {

    try(VariantContextWriter writer = new VariantContextWriterBuilder()
            .setOutputFile(output)
            .setReferenceDictionary(dict)
            .setOptions(VariantContextWriterBuilder.DEFAULT_OPTIONS)
            .build()) {
        writer.writeHeader(vcfHeader);
        while (variants.hasNext()) {
            final VariantContext context = variants.next();

            final VariantContextBuilder builder = new VariantContextBuilder(context);
            if (zCallThresholds.containsKey(context.getID())) {
                final String[] zThresh = zCallThresholds.get(context.getID());
                builder.attribute(InfiniumVcfFields.ZTHRESH_X, zThresh[0]);
                builder.attribute(InfiniumVcfFields.ZTHRESH_Y, zThresh[1]);
            }
            final Genotype originalGenotype = context.getGenotype(0);
            final Map<String, Object> newAttributes = originalGenotype.getExtendedAttributes();
            final VCFEncoder vcfEncoder = new VCFEncoder(vcfHeader, false, false);
            final Map<Allele, String> alleleMap = vcfEncoder.buildAlleleStrings(context);

            final String zCallAlleles = zCallPedFile.getAlleles(context.getID());
            if (zCallAlleles == null) {
                throw new PicardException("No zCall alleles found for snp " + context.getID());
            }
            final List<Allele> zCallPedFileAlleles = buildNewAllelesFromZCall(zCallAlleles, context.getAttributes());
            newAttributes.put(InfiniumVcfFields.GTA, alleleMap.get(originalGenotype.getAllele(0)) + UNPHASED + alleleMap.get(originalGenotype.getAllele(1)));
            newAttributes.put(InfiniumVcfFields.GTZ, alleleMap.get(zCallPedFileAlleles.get(0)) + UNPHASED + alleleMap.get(zCallPedFileAlleles.get(1)));

            final Genotype newGenotype = GenotypeBuilder.create(originalGenotype.getSampleName(), zCallPedFileAlleles,
                    newAttributes);
            builder.genotypes(newGenotype);
            logger.record("0", 0);
            // AC, AF, and AN are recalculated here
            VariantContextUtils.calculateChromosomeCounts(builder, false);
            final VariantContext newContext = builder.make();
            writer.add(newContext);
        }
    }
}
 

@Test
public void testFastGenotypeIsSerializable() {
    Genotype genotype = GenotypeBuilder.create("sample1", Collections.nCopies(2, Allele.create("C", false)));
    SparkTestUtils.roundTripInKryo(genotype, genotype.getClass(), SparkContextFactory.getTestSparkContext().getConf());
}