[BUG] IndexError: list index out of range in v3.7 on BAMs with heavy soft-clipping (regression from v3.6)

[AchiadChenzi]

Bug Report: Regression in v3.7 on soft-clipped reads

What I Ran:
I executed REDItools3 analyze on a BAM file containing paired-end RNA-seq reads where overlapping pairs were soft-clipped using bamutil clipOverlap (resulting in reads with large leading soft-clips, e.g., 80S20M).

python -m reditools analyze \
    --bam clipped.bam \
    --reference genome.fa \
    --strand 2 \
    --threads 16 \
    --window 1000000 \
    --min-read-depth 20 \
    --variants all \
    -o out.tsv

What I Got:
The command crashed during the run, terminating the joblib worker pool with the following error output:

[ERROR] (<class 'IndexError'>) list index out of range
[ERROR] (<class 'IndexError'>) list index out of range
[ERROR] Killing job

This error occurs consistently across samples with short insert sizes (which contain a high fraction of heavily soft-clipped reads). Running the same command on the original, unclipped STAR BAM files completes successfully.

Comparison with v3.6:
This is a regression. When running the exact same command on the exact same clipped BAM files using REDItools3 v3.6, the execution completes successfully without any errors.

[ahanden]

Can you run the 3.7 analysis again with the --debug flag and post the traceback here?

[ahanden]

I think I may have found the problem. I'm able to replicate the error if I have a sequence that has been entirely soft clipped (e.g. 100S). It's an edge case I hadn't considered when writing REDItools since such a read would never come out of an aligner directly. I will have to decide how to address this.

May I ask, why are you using clipOverlap for your project?

[AchiadChenzi]

The reason we use clipOverlap is that our upstream pipeline prepares BAM files to be compatible with deNovo-Detect, which explicitly requires overlapping paired-end reads to be soft-clipped to eliminate algorithm-dependent biases and false positives.
Since we are running both deNovo-Detect and REDItools3 on the same processed datasets, these completely soft-clipped reads are present in our inputs.

[ahanden]

I have found the offending code in REDItools. The change from 3.6 to 3.7 was to reduce memory consumption. I can add a conditional to check for fully clipped reads, but I'm reluctant to slow down the main loop for an edge case.

Have you tried running clipOverlap with the --unmapped option?