index out of bound error related to hard clipping STILL an issue

[rebber]

Hi,

I have an issue that is similar to #283 and #254 but not fixed by the fixes for these (PRs #293 and #284 ). I commented on #283 last year, but maybe no-one saw it since that issue was already closed. So here it comes again.

The error is:

There was Exception while processing record on KH20201014:17728542 on region 17:41264305-41264924. The processing will be continued from the next record.
java.lang.StringIndexOutOfBoundsException: String index out of range: 71
        at java.lang.String.charAt(String.java:658)
        at com.astrazeneca.vardict.modules.CigarParser.parseCigar(CigarParser.java:523)
        at com.astrazeneca.vardict.modules.CigarParser.process(CigarParser.java:88)
        at java.util.concurrent.CompletableFuture.uniApply(CompletableFuture.java:602)
        at java.util.concurrent.CompletableFuture.uniApplyStage(CompletableFuture.java:614)
        at java.util.concurrent.CompletableFuture.thenApply(CompletableFuture.java:1983)
        at com.astrazeneca.vardict.modes.AbstractMode.pipeline(AbstractMode.java:48)
        at com.astrazeneca.vardict.modes.SomaticMode.processBothBamsInPipeline(SomaticMode.java:111)
        at com.astrazeneca.vardict.modes.SomaticMode.notParallel(SomaticMode.java:49)
        at com.astrazeneca.vardict.VarDictLauncher.start(VarDictLauncher.java:65)
        at com.astrazeneca.vardict.Main.main(Main.java:15)

and repeated for a number of reads with the same issue.

The program just skips the read with error and then continues to the next read, so we have just ignored these errors. However, we are now sequencing to higher depths, and due to that often get >10 reads with this error in the same region. Then the Configuration.MAX_EXCEPTION_COUNT is hit and the whole program is stopped with error like:

VarDictJava fails (there were 11 continued exceptions during the run).
Critical exception occurs on region: 17:41264305-41264924, program will be stopped.
java.util.concurrent.CompletionException: java.lang.RuntimeException: java.lang.StringIndexOutOfBoundsException: String index out of range: 79

Since the program won't complete to the end of the bed file for these cases this has become an increasing problem for us.

The command I run is:

vardict-java -G human_g1k_v37_decoy.fasta -f 0.01 -N tumor -b "indel_HC_error.tumor.bam|indel_HC_error.normal.bam" -mfreq 0.01 -nmfreq 0.01  -c 1 -S 2 -E 3 -g 4 -Q 10 problem_region.bed 

I have VarDictJava v 1.8.2 (installed with conda, thus the direct "vardict-java" call). The error happens only when local realignment is used. With option "-k 0" it doesn't happen, but we really want to use that functionality. I attach a zip with minimal bams and bed which trigger the error (although <10 times, so the program is not stopped by it). minimal_bam_github.zip

I have paired reads (mapped to GRCh37) where the overlapping parts in the pairs have been trimmed with fgbio ClipBam, leaving hard clipped bases in the CIGAR. The issue happens at a very specific pattern in the CIGAR: xM3D8M3IyH (with x as num mapped and y as num hard clipped). I have tried to follow in the source code what happens to these reads, and has found the following flow for parsing the CIGAR:

• In src/main/java/com/astrazeneca/vardict/modules/CigarParser.java, CIGAR is run through function parseCigar():
  • with local realignment turned on (option -k = 1, default) the CIGAR is modified in src/main/java/com/astrazeneca/vardict/modules/CigarModifier.java with function cigarModifier.modifyCigar()
    • CIGAR of type xM3D8M3IyH matches "NOTDIG_DIG_I_DIG_M_DIG_DI_DIGI" (regex "(\d+)D(\d+)M(\d+)([DI])(\d+I)?") and is thus modified with function combineToCloseToCorrect() (line 209), which turns the 3D8M3I pattern into 11D11I, so modified CIGAR is xM11D11IyH
      • pull request #254 (fix for issue 253) makes sure to avoid CIGAR modification with similar function combineToCloseToOne() when the CIGAR starts with hard clipped bases, but no such check for this case
  • back in CigarParser.java, still within parseCigar, there is a function cleanupCigar() supposed to remove e.g. hard-clipped bases at the ends, but this operates on the original CIGAR (rec.getCigar(), rec.setCigar()), while the modified cigar in variable "cigar" is what is used going forward, so we continue with xM11D11IyH
  • in loop processCigar,
    • the first matching element is parsed, and then (line 490) it checks if next element is a deletion, with function isCloserThenVextAndGoodBase(),
      • this returns true, since we have deletion element after matching element (xM11D11IyH)
        • pull request #293 (fix for issue 283) makes sure to return a false from this function, if the second next element is hard clipping, but that doesn't happen for this case since there is also an insertion after the deletion, and then the hard clipping comes
    • since isCloserThenVextAndGoodBase() in line 490 gives true, it's checked if element after the deletion is an insertion with function isTwoInsertionsAhead() (line 515)
      • this returns true, since we have insertion element after deletion element (xM11D11IyH)
    • since isTwoInsertionsAhead() in line 515 returns true, the base qualities are checked up until one position after the insertion ends (line 523)
      • this fails since we have no more bases (just hard clippings), and we get the "String index out of range" error
      • it is not expected that this insertion is at end of the actual read since earlier checks should have set any insertions in the ends to soft clips, but this case wasn't discovered since there is hard clipping at the end, thus we get the error

I see multiple possible solutions:

• Avoid CIGAR modification by combineToCloseToCorrect() by checking that there isn't hard clipping following the deletion - matching - deletion/insertion - (insertion) pattern, e.g. by extending the "NOTDIG_DIG_I_DIG_M_DIG_DI_DIGI" regex to "(\d+)D(\d+)M(\d+)([DI])(\d+I)?(\d+H)?$" and make a similar exception as for combineToCloseToOne(), but for hard clipping at the end
• Change the exception in isCloserThenVextAndGoodBase() so that false is returned whenever there is no mapping element between indel elements and any possible final hard clipping (and not only when the first deletion is immediately followed by hard clipping)
• Fix so that cleanupCigar() runs on the modified CIGAR

Is it possible that you could add a fix for this? I am not very experienced in writing java (just reading it), and also doesn't know how my suggested solutions would go with other parts of the logic, therefore I don't want to make a pull request with any of my solutions. But hopefully you could do this in a good way.

Many thanks! Rebecka

[rebber]

@asmlgkj I'm not sure it's the exact same error you are having. But anyway, we solved it temporarily by making a local clone of vardict and in there increasing the value of Configuration.MAX_EXCEPTION_COUNT to a value high enough so that the number of exceptions per site wouldn't reach it, with our typical read depth. So maybe you could try something like that.

[ghost]

As of March 2022, there are still problems with this part of the code (CigarParser.java), from what I can see mostly related to not checking whether the index created is within the bounds of the string, eg.

        // $q1 quality of last base before deletion
        char qualityOfLastSegmentBeforeDel = queryQuality.charAt(readPositionIncludingSoftClipped - 1);

there is no check that readPositionIncludingSoftClipped is > 0 (and indeed sometimes is not). There are at least two or three more places like that.