Some reads are almost entirely soft-clipped

[marcelm]

Strobealign sometimes produces alignments where only a couple of bases are aligned and the rest is soft clipped. For example, running the tests/compare-baseline.sh script with option -s (single ends) currently produces the file baseline/bam/acc4cffe5ac2c4db266c58d00b7b6462c6b4189c.se.bam, where read SRR6055476.83000 is mapped with CIGAR 2M149S and mapping quality 60, which doesn’t make sense. This particular case seems to be caused by a false positive hit.

It would be better to mark such extremely short hits as unmapped. The question is which minimum number of aligned bases we require. It should definitely be at least $k$ but could be higher, perhaps $2k$?

In this dataset, there are 40 alignments shorter than $k=20$ and 664 shorter than $2k=40$.

[marcelm]

Alternatively, we could filter out false positive hits.

[ksahlin]

I think we removed the 'check' that the strobe and reference are identical (i.e., checking for false positive hits). Perhaps we should reintroduce it?

[marcelm]

I did a couple of experiments on single-end reads and here are some intermediate results.

We already have reverse_nam_if_needed, which will fix NAMs that are in the wrong orientation, and which returns false if the NAM is inconsistent (if the beginning matches in a different orientation than the end). I attempted to just discard NAMs that are inconsistent, but this decreases mapping accuracy (by up to 0.2 percentage points).

I then added the criterion that only inconsistent NAMs based upon a single match should be discarded. This also decreases accuracy, but not quite as much.

The mapping rate decreases by almost the same percentage, so most of the alignments we lose were actually correct.

So we kind of need these 'false positives' (which aren’t false at all).

I wanted to tackle this problem at an earlier stage than just removing too short alignments, but maybe that is the easiest solution here. I’ll test that as well.

I think we removed the 'check' that the strobe and reference are identical (i.e., checking for false positive hits). Perhaps we should reintroduce it?

Do you remember when we did this and where? If you mean the reverse_nam_if_needed thing, then that is still there.

[ksahlin]

Do you remember when we did this and where? If you mean the reverse_nam_if_needed thing, then that is still there.

Yes, I meant that function.

I just looked at reverse_nam_if_needed and I think it may have a flaw. A NAM can consist of several hits. One hit could be good and another one could be a 'false positive'. Currently the function only checks either if both are 'good'

if (ref_start_kmer == read_start_kmer && ref_end_kmer == read_end_kmer) {

or if both ends are "false positives".

It could, perhaps, be that the fully softclipped cases we are observing is when one end is correct and the other one is a false positive. So it means four possible combinations to test.

[ksahlin]

It could, perhaps, be that the fully softclipped cases we are observing is when one end is correct and the other one is a false positive. So it means four possible combinations to test.

Maybe this could also explain the discrepancy between the mapping accuracy and the score accuracy that we saw. It would be great to start using the score accuracy to the same extent as we rely on mapping accuracy. I am hopeful :)

[ksahlin]

The optimism didn't last long :) Looking at the code again I doubt that what I wrote above (ends in different orientations which can happen for the matches within a NAM) actually affects the checks in reverse_nam_if_needed.

My guess would then be that these almost entirely softclipped reads come from true hash collisions (thus, not matching the reference in either orientation), but your experiments show that they should be mapped to that location. This leaves me clueless on why this happens.

Is there a coordinate extraction bug somewhere?

Without understanding the underlying cause of the issue, is the remedy to align the read in both directions at the location?