Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix case where hom-ref "variant" with no data had wrong sized PLs and #7644

Merged
merged 2 commits into from
Jan 31, 2022
Merged

Fix case where hom-ref "variant" with no data had wrong sized PLs and #7644

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
5154d3f
Fix case where hom-ref "variant" with no data had wrong sized PLs and
ldgauthier Jan 25, 2022
6d5e9ab
Address review comments
ldgauthier Jan 28, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 6 additions & 1 deletion src/main/java/org/broadinstitute/hellbender/tools/walkers/variantutils/ReblockGVCF.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -507,9 +507,14 @@ protected GenotypeBuilder changeCallToHomRefVersusNonRef(final VariantContext lo
} else {
final List<Allele> bestAlleles = AlleleSubsettingUtils.calculateMostLikelyAlleles(lowQualVariant, genotype.getPloidy(), 1);
final Allele bestAlt = bestAlleles.stream().filter(a -> !a.isReference()).findFirst().orElse(Allele.NON_REF_ALLELE); //allow span dels
//we care about the best alt even though it's getting removed because NON_REF should get the best likelihoods
//it shouldn't matter that we're passing in different alt alleles since the GenotypesContext only knows
// the called alleles and this is a reference genotype that will stay hom-ref
final GenotypesContext context = AlleleSubsettingUtils.subsetAlleles(lowQualVariant.getGenotypes(),
genotype.getPloidy(), lowQualVariant.getAlleles(), Arrays.asList(inputRefAllele, bestAlt),
null, GenotypeAssignmentMethod.BEST_MATCH_TO_ORIGINAL, lowQualVariant.getAttributeAsInt(VCFConstants.DEPTH_KEY, 0), false); //BEST_MATCH to avoid no-calling low qual genotypes
null, GenotypeAssignmentMethod.BEST_MATCH_TO_ORIGINAL, //BEST_MATCH to avoid no-calling low qual genotypes
lowQualVariant.getAttributeAsInt(VCFConstants.DEPTH_KEY, 0),
true); //emitEmptyPLs = true to make sure we always subset
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure what a true here results into looking at the subsetAllele code. If it is what you need and that behavior from subsetAlleles is stable then I don't a problem with that.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The alternative is to drop the PLs, which I do not want. I've changed the parameter name to be a little more clear.

final Genotype subsetG = context.get(0);
gb = new GenotypeBuilder(subsetG).noAttributes(); //remove attributes because hom ref blocks shouldn't have posteriors
//subsetting may strip GQ and PLs for low qual genotypes
Expand Down
13 changes: 13 additions & 0 deletions ...st/java/org/broadinstitute/hellbender/tools/walkers/variantutils/ReblockGVCFUnitTest.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -110,6 +110,12 @@ public void testLowQualVariantToGQ0HomRef() {
Assert.assertTrue(!modified.filtersWereApplied());
Assert.assertEquals(modified.getLog10PError(), VariantContext.NO_LOG10_PERROR);

final Genotype longPls = VariantContextTestUtils.makeG("sample1", Allele.NO_CALL, Allele.NO_CALL, 0,0,0,0,0,0);
final VariantContext lotsOfZeroPls = makeNoDepthVC("lowQualVar", Arrays.asList(LONG_REF, DELETION, Allele.NON_REF_ALLELE), LONG_REF.length(), longPls);
final VariantContext properlySubset = reblocker.lowQualVariantToGQ0HomRef(lotsOfZeroPls).make();
Assert.assertEquals(properlySubset.getGenotype(0).getPL().length, 3);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you check the content of PL also is correnct? is it all 0s or '.'?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Assert.assertEquals(...getPL(), new int [] { 0, 0, 0}) may just work.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It does work, thanks. This is done.



//No-calls were throwing NPEs. Now they're not.
final Genotype g2 = VariantContextTestUtils.makeG("sample1", Allele.NO_CALL,Allele.NO_CALL);
final VariantContext noData = makeDeletionVC("noData", Arrays.asList(LONG_REF, DELETION, Allele.NON_REF_ALLELE), LONG_REF.length(), g2);
Expand Down Expand Up @@ -425,6 +431,13 @@ private VariantContext makeDeletionVC(final String source, final List<Allele> al
.genotypes(Arrays.asList(genotypes)).unfiltered().log10PError(-3.0).attribute(VCFConstants.DEPTH_KEY, EXAMPLE_DP).make();
}

private VariantContext makeNoDepthVC(final String source, final List<Allele> alleles, final int refLength, final Genotype... genotypes) {
final int start = DEFAULT_START;
final int stop = start+refLength-1;
return new VariantContextBuilder(source, "1", start, stop, alleles)
.genotypes(Arrays.asList(genotypes)).unfiltered().log10PError(-3.0).make();
}

private Genotype addAD(final Genotype g, final int... ads) {
return new GenotypeBuilder(g).AD(ads).make();
}
Expand Down