use of htsjdk.variant.variantcontext.Allele in project gatk by broadinstitute.
the class EvaluateCopyNumberTriStateCalls method buildAndAnnotateTruthOverlappingGenotype.
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();
}
use of htsjdk.variant.variantcontext.Allele in project gatk by broadinstitute.
the class ValidateVariants method apply.
@Override
public void apply(final VariantContext vc, final ReadsContext readsContext, final ReferenceContext ref, final FeatureContext featureContext) {
if (DO_NOT_VALIDATE_FILTERED && vc.isFiltered()) {
return;
}
// get the true reference allele
final Allele reportedRefAllele = vc.getReference();
final int refLength = reportedRefAllele.length();
final Allele observedRefAllele = hasReference() ? Allele.create(Arrays.copyOf(ref.getBases(), refLength)) : null;
final Set<String> rsIDs = getRSIDs(featureContext);
for (final ValidationType t : validationTypes) {
try {
applyValidationType(vc, reportedRefAllele, observedRefAllele, rsIDs, t);
} catch (TribbleException e) {
if (WARN_ON_ERROR) {
logger.warn("***** " + e.getMessage() + " *****");
} else {
throw new UserException.FailsStrictValidation(drivingVariantFile, t, e.getMessage());
}
}
}
}
use of htsjdk.variant.variantcontext.Allele in project gatk by broadinstitute.
the class SelectVariants method addAnnotations.
/*
* Add annotations to the new VC
*
* @param builder the new VC to annotate
* @param originalVC the original VC
* @param selectedSampleNames the post-selection list of sample names
*/
private void addAnnotations(final VariantContextBuilder builder, final VariantContext originalVC, final Set<String> selectedSampleNames) {
if (fullyDecode) {
// TODO -- annotations are broken with fully decoded data
return;
}
if (keepOriginalChrCounts) {
final int[] indexOfOriginalAlleleForNewAllele;
final List<Allele> newAlleles = builder.getAlleles();
final int numOriginalAlleles = originalVC.getNAlleles();
// if the alleles already match up, we can just copy the previous list of counts
if (numOriginalAlleles == newAlleles.size()) {
indexOfOriginalAlleleForNewAllele = null;
} else // otherwise we need to parse them and select out the correct ones
{
indexOfOriginalAlleleForNewAllele = new int[newAlleles.size() - 1];
Arrays.fill(indexOfOriginalAlleleForNewAllele, -1);
// note that we don't care about the reference allele at position 0
for (int newI = 1; newI < newAlleles.size(); newI++) {
final Allele newAlt = newAlleles.get(newI);
for (int oldI = 0; oldI < numOriginalAlleles - 1; oldI++) {
if (newAlt.equals(originalVC.getAlternateAllele(oldI), false)) {
indexOfOriginalAlleleForNewAllele[newI - 1] = oldI;
break;
}
}
}
}
if (originalVC.hasAttribute(VCFConstants.ALLELE_COUNT_KEY)) {
builder.attribute(GATKVCFConstants.ORIGINAL_AC_KEY, getReorderedAttributes(originalVC.getAttribute(VCFConstants.ALLELE_COUNT_KEY), indexOfOriginalAlleleForNewAllele));
}
if (originalVC.hasAttribute(VCFConstants.ALLELE_FREQUENCY_KEY)) {
builder.attribute(GATKVCFConstants.ORIGINAL_AF_KEY, getReorderedAttributes(originalVC.getAttribute(VCFConstants.ALLELE_FREQUENCY_KEY), indexOfOriginalAlleleForNewAllele));
}
if (originalVC.hasAttribute(VCFConstants.ALLELE_NUMBER_KEY)) {
builder.attribute(GATKVCFConstants.ORIGINAL_AN_KEY, originalVC.getAttribute(VCFConstants.ALLELE_NUMBER_KEY));
}
}
VariantContextUtils.calculateChromosomeCounts(builder, false);
if (keepOriginalDepth && originalVC.hasAttribute(VCFConstants.DEPTH_KEY)) {
builder.attribute(GATKVCFConstants.ORIGINAL_DP_KEY, originalVC.getAttribute(VCFConstants.DEPTH_KEY));
}
boolean sawDP = false;
int depth = 0;
for (final String sample : selectedSampleNames) {
final Genotype g = originalVC.getGenotype(sample);
if (!g.isFiltered()) {
if (g.hasDP()) {
depth += g.getDP();
sawDP = true;
}
}
}
if (sawDP) {
builder.attribute(VCFConstants.DEPTH_KEY, depth);
}
}
use of htsjdk.variant.variantcontext.Allele in project gatk by broadinstitute.
the class XHMMSegmentGenotyperIntegrationTest method assertVariantsInfoFieldsAreConsistent.
private void assertVariantsInfoFieldsAreConsistent(final List<VariantContext> variants) {
for (final VariantContext variant : variants) {
final int expectedAN = variant.getGenotypes().size();
final int[] expectedAC = new int[CopyNumberTriState.values().length];
final List<Allele> alleles = variant.getAlleles();
for (final Genotype genotype : variant.getGenotypes()) {
Assert.assertEquals(genotype.getAlleles().size(), 1);
final int alleleIndex = alleles.indexOf(genotype.getAllele(0));
Assert.assertTrue(alleleIndex >= 0);
expectedAC[alleleIndex]++;
}
Assert.assertEquals(variant.getAlleles(), CopyNumberTriStateAllele.ALL_ALLELES);
Assert.assertTrue(variant.hasAttribute(XHMMSegmentGenotyper.NUMBER_OF_TARGETS_KEY));
Assert.assertTrue(variant.hasAttribute(VCFConstants.ALLELE_COUNT_KEY));
Assert.assertTrue(variant.hasAttribute(VCFConstants.END_KEY));
Assert.assertTrue(variant.hasAttribute(VCFConstants.ALLELE_FREQUENCY_KEY));
Assert.assertTrue(variant.hasAttribute(VCFConstants.ALLELE_NUMBER_KEY));
final int expectedANAnotherWay = IntStream.of(expectedAC).sum();
Assert.assertEquals(expectedANAnotherWay, expectedAN);
Assert.assertEquals(variant.getAttributeAsInt(VCFConstants.ALLELE_NUMBER_KEY, -1), expectedAN);
final double[] expectedAF = IntStream.of(expectedAC).mapToDouble(c -> c / (double) expectedAN).toArray();
final double[] observedAFWithoutRef = variant.getAttributeAsList(VCFConstants.ALLELE_FREQUENCY_KEY).stream().mapToDouble(o -> Double.parseDouble(String.valueOf(o))).toArray();
Assert.assertEquals(observedAFWithoutRef.length, expectedAF.length - 1);
for (int i = 0; i < observedAFWithoutRef.length; i++) {
Assert.assertEquals(observedAFWithoutRef[i], expectedAF[i + 1], 0.001);
}
final int[] observedACWithoutRef = variant.getAttributeAsList(VCFConstants.ALLELE_COUNT_KEY).stream().mapToInt(o -> Integer.parseInt(String.valueOf(o))).toArray();
Assert.assertEquals(observedACWithoutRef.length, expectedAC.length - 1);
for (int i = 0; i < observedACWithoutRef.length; i++) {
Assert.assertEquals(observedACWithoutRef[i], expectedAC[i + 1]);
}
Assert.assertEquals(variant.getAttributeAsInt(XHMMSegmentGenotyper.NUMBER_OF_TARGETS_KEY, -1), XHMMSegmentCallerBaseIntegrationTest.REALISTIC_TARGETS.targetCount(variant));
}
}
use of htsjdk.variant.variantcontext.Allele in project gatk by broadinstitute.
the class ReadThreadingAssemblerUnitTest method testAssembleRefAndInsertion.
@Test(dataProvider = "AssembleIntervalsWithVariantData")
public void testAssembleRefAndInsertion(final ReadThreadingAssembler assembler, final SimpleInterval loc, final int nReadsToUse, final int variantSite) {
final byte[] refBases = seq.getSubsequenceAt(loc.getContig(), loc.getStart(), loc.getEnd()).getBases();
for (int insertionLength = 1; insertionLength < 10; insertionLength++) {
final Allele refBase = Allele.create(refBases[variantSite], false);
final Allele altBase = Allele.create(new String(refBases).substring(variantSite, variantSite + insertionLength + 1), true);
final VariantContextBuilder vcb = new VariantContextBuilder("x", loc.getContig(), variantSite, variantSite + insertionLength, Arrays.asList(refBase, altBase));
testAssemblyWithVariant(assembler, refBases, loc, nReadsToUse, vcb.make());
}
}
Aggregations