HaplotypeCaller and Mutect2 in 4.1.6.0 failing with Smith-Waterman cigar reference length mismatch

[droazen]

As reported in this forum thread: https://gatk.broadinstitute.org/hc/en-us/community/posts/360060174372-Haplotype-Caller-4-1-6-0-java-lang-IllegalStateException-Smith-Waterman-alignment-failure- as well as in the closed ticket #6490 (comment), there appears to be a regression in GATK 4.1.6.0 causing HaplotypeCaller/Mutect2 to fail with an error like the following:

java.lang.IllegalStateException: Smith-Waterman alignment failure. Cigar = 275M with reference length 275 but expecting referenc
e length of 303 ref = GGACTGGTCATCAGCTACCCCGAGACCAGGACTGGTCATCAGCTACCCCGAGACCAGGACTGGTCATCAGCTACCCCGAGACCAGGACTGGTCATCAGCTACCCCG
AGACCAGGACTGGTCATCAGCTACCCCGAGAACAGGACTGGTCATCAGCTACCCCGAGACCAGGACTGGTCATCAGCTACCCCGAGACCAGGACTGGTCATCAGCTACCCCGAGACCAGGACTGGTCA
TCAGCTACTCCGAGACCAGCATGGAGAGGTTTGCTGATGGTTGGCTGACTGCTAGTGTGAGCACTTGTC path GGACTGGTCATCAGCTACCCCGAGACCAGGACTGGTCATCAGCTACCCCGAGA
CCAGGACTGGTCATCAGCTACCCCGAGACCAGGACTGGTCATCAGCTACCCCGAGAACAGGACTGGTCATCAGCTACCCCGAGAACAGGACTGGTCATCAGCTACCCCGAGACCAGGACTGGTCATCA
GCTACCCCGAGACCAGGACTGGTCATCAGCTACTCCGAGACCAGCATGGAGAGGTTTGCTGATGGTTGGCTGACTGCTAGTGTGAGCACTTGTC
at org.broadinstitute.hellbender.tools.walkers.haplotypecaller.readthreading.ReadThreadingAssembler.findBestPaths(ReadTh
readingAssembler.java:354)
at org.broadinstitute.hellbender.tools.walkers.haplotypecaller.readthreading.ReadThreadingAssembler.assembleKmerGraphsAn
dHaplotypeCall(ReadThreadingAssembler.java:196)
at org.broadinstitute.hellbender.tools.walkers.haplotypecaller.readthreading.ReadThreadingAssembler.runLocalAssembly(Rea
dThreadingAssembler.java:146)
at org.broadinstitute.hellbender.tools.walkers.haplotypecaller.AssemblyBasedCallerUtils.assembleReads(AssemblyBasedCalle
rUtils.java:269)
at org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCallerEngine.callRegion(HaplotypeCallerEngine.ja
va:541)
at org.broadinstitute.hellbender.tools.walkers.haplotypecaller.HaplotypeCaller.apply(HaplotypeCaller.java:212)
at org.broadinstitute.hellbender.engine.AssemblyRegionWalker.processReadShard(AssemblyRegionWalker.java:200)
at org.broadinstitute.hellbender.engine.AssemblyRegionWalker.traverse(AssemblyRegionWalker.java:173)
at org.broadinstitute.hellbender.engine.GATKTool.doWork(GATKTool.java:1048)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.runTool(CommandLineProgram.java:139)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMainPostParseArgs(CommandLineProgram.java:191)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:210)
at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:163)
at org.broadinstitute.hellbender.Main.mainEntry(Main.java:206)
at org.broadinstitute.hellbender.Main.main(Main.java:292)

[droazen]

@davidbenjamin I've spun this off into a separate ticket. I think this might be serious/widespread enough that we'll want a 4.1.6.1 release once we have a working fix. I recommend that we manually do full-scale test runs of both HaplotypeCaller and Mutect2 to confirm that the issue is resolved rather than relying on the integration tests, which are not large-scale enough to catch edge cases like this, unfortunately.

[davidbenjamin]

It seems like the patch in 4.1.6 didn't go far enough and that exception needs to be replaced with a continue in all cases. This seems to be occurring for haplotypes with long indels inside their assembly padding that don't have enough spanning sequence to resolve. Since the variation is inside the padding, it seems safe to ignore. Increasing padding resolves the issue, alhtough this is at the cost of runtime and should not be necessary.

For example, suppose we have a ref haplotype ABCDD, where A, B, and C represent sequences of, say, 100 bases and D is a sequence of 50 bases. Suppose further that A and DD are the padding. Then the cigar of an alt haplotype ABCD gets aligned as a 350 base match that doesn't span the full padded reference region, leading to the error.

I still need to figure out why this didn't happen in 4.1.4 (my guess is that elsewhere the code effectively skipped these haplotypes before the exception).

[davidbenjamin]

It seems like the patch in 4.1.6 didn't go far enough and that exception needs to be replaced with a continue in all cases.

That would work, but I see where I caused the regression upstream. I chopped leading and trailing deletions from haplotype cigars, same as for read cigars, but for haplotypes we want to keep these deletions because the start and end positions need to remain pegged to the reference start and end.

I have a fix + regression test branch, which is running on every M2 validation.

[LakshmanTejaG]

gatk 4.1.4.1 version working fine for
HaplotypeCaller Smith-Waterman cigar reference length mismatch

[droazen]

@davidbenjamin How's the patch coming? Did the M2 validation tests pass on your branch? We'll definitely try to expedite the code review, but I'll think we'll want some additional heavy-duty testing prior to release.

[davidbenjamin]

@droazen The PR is #6544. James has reviewed and requested a few more tests. It's working fine on validations including ~30 exomes and ~15 genomes.

[jamesemery]

We have asked the green team to run their pipeline tests on this branch to at least limit the risk of more full sample failures. It will probably be a few more days before we have those results. @gbggrant

[droazen]

Fixed by #6544 -- closing.