From 5d95d60cd9e29e2cd3ba7bbabd6f996e377b72a0 Mon Sep 17 00:00:00 2001 From: cwhelan Date: Tue, 23 Oct 2018 11:15:16 -0400 Subject: [PATCH] Fix error in genotype given alleles mode when input alleles have genotypes (#5341) * fix issue of alleles not being replaced with star alleles when the given allele variants have genotypes * add unit test to replaceWithSpanDelVC --- .../HaplotypeCallerGenotypingEngine.java | 9 ++-- ...plotypeCallerGenotypingEngineUnitTest.java | 39 ++++++++++++++++++ ...ted.testGenotypeGivenAllelesMode.gatk4.vcf | 1 + ...tGenotypeGivenAllelesMode_givenAlleles.vcf | 17 ++++---- ...otypeGivenAllelesMode_givenAlleles.vcf.idx | Bin 448 -> 433 bytes 5 files changed, 55 insertions(+), 11 deletions(-) diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerGenotypingEngine.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerGenotypingEngine.java index 0335af96c4b..e990f484f68 100644 --- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerGenotypingEngine.java +++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerGenotypingEngine.java @@ -216,18 +216,21 @@ public CalledHaplotypes assignGenotypeLikelihoods(final List haplotyp activeRegionWindow,tracker,activeAllelesToGenotype,emitReferenceConfidence,maxMnpDistance,header,false); } - private List replaceSpanDels(final List eventsAtThisLoc, final Allele refAllele, final int loc) { + @VisibleForTesting + static List replaceSpanDels(final List eventsAtThisLoc, final Allele refAllele, final int loc) { return eventsAtThisLoc.stream().map(vc -> replaceWithSpanDelVC(vc, refAllele, loc)).collect(Collectors.toList()); } - private VariantContext replaceWithSpanDelVC(final VariantContext variantContext, final Allele refAllele, final int loc) { + @VisibleForTesting + static VariantContext replaceWithSpanDelVC(final VariantContext variantContext, final Allele refAllele, final int loc) { if (variantContext.getStart() == loc) { return variantContext; } else { VariantContextBuilder builder = new VariantContextBuilder(variantContext) .start(loc) .stop(loc) - .alleles(Arrays.asList(refAllele, Allele.SPAN_DEL)); + .alleles(Arrays.asList(refAllele, Allele.SPAN_DEL)) + .genotypes(GenotypesContext.NO_GENOTYPES); return builder.make(); } diff --git a/src/test/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerGenotypingEngineUnitTest.java b/src/test/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerGenotypingEngineUnitTest.java index 3a04f3a5999..a74474081a8 100644 --- a/src/test/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerGenotypingEngineUnitTest.java +++ b/src/test/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerGenotypingEngineUnitTest.java @@ -6,6 +6,7 @@ import org.apache.commons.lang3.tuple.Pair; import org.broadinstitute.gatk.nativebindings.smithwaterman.SWOverhangStrategy; import org.broadinstitute.gatk.nativebindings.smithwaterman.SWParameters; +import org.broadinstitute.hellbender.testutils.VariantContextTestUtils; import org.broadinstitute.hellbender.utils.QualityUtils; import org.broadinstitute.hellbender.utils.SimpleInterval; import org.broadinstitute.hellbender.utils.Utils; @@ -559,4 +560,42 @@ public void testRemoveExcessiveAltAlleleFromVC(){ Assert.assertEquals(reducedVC.getNAlleles(), 3); Assert.assertTrue(reducedVC.getAlleles().containsAll(Arrays.asList(Allele.create("A", true), Allele.create("T", false), Allele.create("C", false)))); } + + @Test + public void testReplaceWithSpanDelVC() { + final VariantContext snp = new VariantContextBuilder("source", "1", 1000000, 1000000, + Arrays.asList(Allele.create("A", true), + Allele.create("G", false))).make(); + final VariantContext snpReplacement = + HaplotypeCallerGenotypingEngine.replaceWithSpanDelVC(snp, Allele.create("A", true), 1000000); + + Assert.assertEquals(snpReplacement, snp); + + final VariantContext spanDel = new VariantContextBuilder("source", "1", 999995, 1000005, + Arrays.asList(Allele.create("AAAAAAAAAAA", true), + Allele.create("A", false))).make(); + + final VariantContext spanDelReplacement = + HaplotypeCallerGenotypingEngine.replaceWithSpanDelVC(spanDel, Allele.create("A", true), 1000000); + + final VariantContext expectedSpanDelReplacement = new VariantContextBuilder("source", "1", 1000000, 1000000, + Arrays.asList(Allele.create("A", true), + Allele.SPAN_DEL)).make(); + + VariantContextTestUtils.assertVariantContextsAreEqual(spanDelReplacement, + expectedSpanDelReplacement, Collections.emptyList()); + + final VariantContext spanDelWithGt = new VariantContextBuilder("source", "1", 999995, 1000005, + Arrays.asList(Allele.create("AAAAAAAAAAA", true), + Allele.create("A", false))) + .genotypes(new GenotypeBuilder("S1", + Arrays.asList(spanDel.getAlleles().get(0), spanDel.getAlleles().get(1))).make()).make(); + + final VariantContext spanDelWithGTReplacement = + HaplotypeCallerGenotypingEngine.replaceWithSpanDelVC(spanDelWithGt, Allele.create("A", true), 1000000); + + VariantContextTestUtils.assertVariantContextsAreEqual(spanDelWithGTReplacement, + expectedSpanDelReplacement, Collections.emptyList()); + + } } diff --git a/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/expected.testGenotypeGivenAllelesMode.gatk4.vcf b/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/expected.testGenotypeGivenAllelesMode.gatk4.vcf index 842236c0226..fdf87b8f6d7 100644 --- a/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/expected.testGenotypeGivenAllelesMode.gatk4.vcf +++ b/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/expected.testGenotypeGivenAllelesMode.gatk4.vcf @@ -31,3 +31,4 @@ 20 10004769 . TAAAACTATGC T 980.73 . AC=1;AF=0.500;AN=2;BaseQRankSum=1.603;DP=79;ExcessHet=3.0103;FS=3.758;MLEAC=1;MLEAF=0.500;MQ=52.37;MQRankSum=-6.625;QD=15.32;ReadPosRankSum=1.879;SOR=1.306 GT:AD:DP:GQ:PL 0/1:37,27:64:99:1018,0,1471 20 10004771 . A *,T 998.77 . AC=1,0;AF=0.500,0.00;AN=2;BaseQRankSum=-1.270;DP=74;ExcessHet=3.0103;FS=2.397;MLEAC=1,0;MLEAF=0.500,0.00;MQ=51.81;MQRankSum=-6.325;QD=16.93;ReadPosRankSum=2.185;SOR=1.115 GT:AD:DP:GQ:PL 0/1:32,27,0:59:99:1027,0,1345,1162,1305,2852 20 10006819 . AAAAC T 0 LowQual AC=0;AF=0.00;AN=2;DP=75;ExcessHet=3.0103;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=61.14;SOR=0.572 GT:AD:DP:GQ:PL 0/0:75,0:75:99:0,244,2147483647 +20 10006823 . C *,G 0 LowQual AC=0,0;AF=0.00,0.00;AN=2;DP=73;ExcessHet=3.0103;FS=0.000;MLEAC=0,0;MLEAF=0.00,0.00;MQ=61.18;SOR=0.069 GT:AD:DP:GQ:PL 0/0:17,0,0:17:51:0,229,2147483647,51,819,590 diff --git a/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/testGenotypeGivenAllelesMode_givenAlleles.vcf b/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/testGenotypeGivenAllelesMode_givenAlleles.vcf index c9c28c26d07..48ed82e665e 100644 --- a/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/testGenotypeGivenAllelesMode_givenAlleles.vcf +++ b/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/testGenotypeGivenAllelesMode_givenAlleles.vcf @@ -25,11 +25,12 @@ ##INFO= ##contig= ##contig= -#CHROM POS ID REF ALT QUAL FILTER INFO -20 10000694 . G A . . . -20 10001436 . A AAGGCT . . . -20 10001661 . T C . . . -20 10004094 . A T . . . -20 10004769 . TAAAACTATGC T . . . -20 10004771 . A T . . . -20 10006819 . AAAAC T . . . +#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT SAMPLE1 +20 10000694 . G A . . . GT 0|1 +20 10001436 . A AAGGCT . . . GT 0|1 +20 10001661 . T C . . . GT 0|1 +20 10004094 . A T . . . GT 0|1 +20 10004769 . TAAAACTATGC T . . . GT 0|1 +20 10004771 . A T . . . GT 0|1 +20 10006819 . AAAAC T . . . GT 0|1 +20 10006823 . C G . . . GT 0|1 diff --git a/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/testGenotypeGivenAllelesMode_givenAlleles.vcf.idx b/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/testGenotypeGivenAllelesMode_givenAlleles.vcf.idx index 36d7758150f5b72a804d6f2b33ff0a22e785e721..f0c2511184ce4fdbecc5ccf767cafbfce6ae71c2 100644 GIT binary patch delta 151 zcmX@WypegrYTgMv3}C=vvUN_{#2pIEW_kt_UurN~OcrL;5Hx3Sb9D>}4RVcl_7C+7 XVX&C&!zdBM0o2jT4bev}8s-iFylxf+ delta 166 zcmdnUe1LhvYF