From 8ee86e7c147399277a57f74bf751c66b3bda157f Mon Sep 17 00:00:00 2001
From: jamesemery <emeryj@broadinstitute.org>
Date: Tue, 12 Mar 2024 09:55:44 -0400
Subject: [PATCH] Fix to long deletions that overhang into the assembly window
 causing exceptions in HaplotypeCaller (#8731)

---
 .../HaplotypeCallerEngine.java                |   1 +
 .../HaplotypeCallerIntegrationTest.java       |   3 +-
 ...ode_givenAlleles_ExtremeLengthDeletion.vcf |  29 ++++++++++++++++++
 ...givenAlleles_ExtremeLengthDeletion.vcf.idx | Bin 0 -> 359 bytes
 4 files changed, 32 insertions(+), 1 deletion(-)
 create mode 100644 src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/testGenotypeGivenAllelesMode_givenAlleles_ExtremeLengthDeletion.vcf
 create mode 100644 src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/testGenotypeGivenAllelesMode_givenAlleles_ExtremeLengthDeletion.vcf.idx

diff --git a/src/main/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerEngine.java b/src/main/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerEngine.java
index 9e09eeba8a9..36918f1fb91 100644
--- a/src/main/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerEngine.java
+++ b/src/main/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerEngine.java
@@ -758,6 +758,7 @@ public List<VariantContext> callRegion(final AssemblyRegion region, final Featur
         final List<Event> givenAlleles = features.getValues(hcArgs.alleles).stream()
                 .filter(vc -> hcArgs.forceCallFiltered || vc.isNotFiltered())
                 .flatMap(vc -> GATKVariantContextUtils.splitVariantContextToEvents(vc, false, GenotypeAssignmentMethod.BEST_MATCH_TO_ORIGINAL, false).stream())
+                .filter(event -> event.getStart() >= region.getSpan().getStart()) // filter out events that do not start within the region. This approach works because events that begin upstream of the calling window cannot be called by this region calling code in the frist place. 
                 .collect(Collectors.toList());
 
         if( givenAlleles.isEmpty() && region.size() == 0 ) {
diff --git a/src/test/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java b/src/test/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java
index cddb87c5dac..1e059b444b5 100644
--- a/src/test/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java
+++ b/src/test/java/org/broadinstitute/hellbender/tools/walkers/haplotypecaller/HaplotypeCallerIntegrationTest.java
@@ -671,7 +671,8 @@ public void testFloorGVCFBlocks(final String inputFileName, final String referen
     public Object[][] getForceCallingInputs() {
         return new Object[][] {
                 {NA12878_20_21_WGS_bam, new File(TEST_FILES_DIR, "testGenotypeGivenAllelesMode_givenAlleles.vcf"), "20:10000000-10010000"},
-                {NA12878_20_21_WGS_bam, new File(toolsTestDir, "mutect/gga_mode.vcf"), "20:9998500-10010000"}
+                {NA12878_20_21_WGS_bam, new File(toolsTestDir, "mutect/gga_mode.vcf"), "20:9998500-10010000"},
+                {NA12878_20_21_WGS_bam, new File(TEST_FILES_DIR, "testGenotypeGivenAllelesMode_givenAlleles_ExtremeLengthDeletion.vcf"), "20:9998500-10010000"} // This is designed to test https://github.com/broadinstitute/gatk/issues/8675, which stemmed from an edge case in the force calling logic where a deletion allele that is longer than the assembly window padding spans into the assembly window. This tests that we do not see an exception in this case.
         };
     }
 
diff --git a/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/testGenotypeGivenAllelesMode_givenAlleles_ExtremeLengthDeletion.vcf b/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/testGenotypeGivenAllelesMode_givenAlleles_ExtremeLengthDeletion.vcf
new file mode 100644
index 00000000000..e36b2052864
--- /dev/null
+++ b/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/testGenotypeGivenAllelesMode_givenAlleles_ExtremeLengthDeletion.vcf
@@ -0,0 +1,29 @@
+##fileformat=VCFv4.2
+##FILTER=<ID=LowQual,Description="Low quality">
+##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
+##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
+##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">
+##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes, for each ALT allele, in the same order as listed">
+##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency, for each ALT allele, in the same order as listed">
+##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
+##INFO=<ID=BaseQRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt Vs. Ref base qualities">
+##INFO=<ID=ClippingRankSum,Number=1,Type=Float,Description="Z-score From Wilcoxon rank sum test of Alt vs. Ref number of hard clipped bases">
+##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth; some reads may have been filtered">
+##INFO=<ID=DS,Number=0,Type=Flag,Description="Were any of the samples downsampled?">
+##INFO=<ID=ExcessHet,Number=1,Type=Float,Description="Phred-scaled p-value for exact test of excess heterozygosity">
+##INFO=<ID=FS,Number=1,Type=Float,Description="Phred-scaled p-value using Fisher's exact test to detect strand bias">
+##INFO=<ID=InbreedingCoeff,Number=1,Type=Float,Description="Inbreeding coefficient as estimated from the genotype likelihoods per-sample when compared against the Hardy-Weinberg expectation">
+##INFO=<ID=MLEAC,Number=A,Type=Integer,Description="Maximum likelihood expectation (MLE) for the allele counts (not necessarily the same as the AC), for each ALT allele, in the same order as listed">
+##INFO=<ID=MLEAF,Number=A,Type=Float,Description="Maximum likelihood expectation (MLE) for the allele frequency (not necessarily the same as the AF), for each ALT allele, in the same order as listed">
+##INFO=<ID=MQ,Number=1,Type=Float,Description="RMS Mapping Quality">
+##INFO=<ID=MQRankSum,Number=1,Type=Float,Description="Z-score From Wilcoxon rank sum test of Alt vs. Ref read mapping qualities">
+##INFO=<ID=QD,Number=1,Type=Float,Description="Variant Confidence/Quality by Depth">
+##INFO=<ID=RAW_MQ,Number=1,Type=Float,Description="Raw data for RMS Mapping Quality">
+##INFO=<ID=ReadPosRankSum,Number=1,Type=Float,Description="Z-score from Wilcoxon rank sum test of Alt vs. Ref read position bias">
+##INFO=<ID=SOR,Number=1,Type=Float,Description="Symmetric Odds Ratio of 2x2 contingency table to detect strand bias">
+##contig=<ID=20,length=63025520>
+##contig=<ID=21,length=48129895>
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	SAMPLE1
+20	10000694	.	GAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAAGAAAAAA	A	.	.	.	GT	0|1
diff --git a/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/testGenotypeGivenAllelesMode_givenAlleles_ExtremeLengthDeletion.vcf.idx b/src/test/resources/org/broadinstitute/hellbender/tools/haplotypecaller/testGenotypeGivenAllelesMode_givenAlleles_ExtremeLengthDeletion.vcf.idx
new file mode 100644
index 0000000000000000000000000000000000000000..f469a4b8452bb36aea9d4360945adfaccb3d7ab2
GIT binary patch
literal 359
zcmZ8cO-sW-5M4iR;?<wgi6V&JLbJv~8=)q}+ma^JB)aaF-D$B0kDdfi9{ijB180kY
z;K00r_vYhGK1m)Agb+uZt4`w$02Cgb2Yf-ddx93VUZSbdK~vH*cvk_TCvfQP##QLS
zx(1fcmUU-5b+jQIPLZ|RLtC!3rQI5<QfuB7R-|b5RN8IP#7t}SGh3r-285EY<k-s;
zO+#(MJ=9q<+Ey#^b;_nWKEBW3GdvOEnDbV~c`=u2Dzh}dQ!^Q7;_~`x6#Z8o@<b(a
uAtIJSm+B`M@q7|zW4Tv)a(rK8xi}Z$)J4SbCJq^Y{!WMTe!?)Qy~SUqKxC=_

literal 0
HcmV?d00001