Almost perfect secondary matches not found

[tprodanov]

Good afternoon,

I encountered a problem (present in the latest commit) where almost perfect secondary alignments don't get reported no matter the -S and -M values.

Specifically, using the following flies:
haplotypes.fa.gz
read.fq.gz

strobealign --no-progress -M 10000 -N 10000 -S 0.1 -f 0.0001 -k 15 --eqx -r 150 -t 8 \
    db/loci/PARVB/haplotypes.fa.gz read.fq.gz

produces 116 alignments but I know for sure that there should be more. For example, the given read aligns to haplotypes HG02451.1 or NA19240.1 with just a single mismatch, so I hoped they would be reported. I tried different values of -M, -N, -S and -f but could not make it work. Could you please look into it and perhaps provide a flag to disable whatever heuristics that stops the read alignment when a perfect alignment is already found?

Best wishes,
Timofey

[ksahlin]

Hi Timofey,

Thanks for the detailed report and for providing the test files.

Is there a clean definition of "enough alignments" in a haplotype-aware context like yours, where many sequences are nearly identical by design? The current heuristics were tuned for reference genome mapping, not multi-haplotype panels.

What is the edit distance / alignment score of the 116 alignments that are reported -- are they all perfect?

On the question of a flag to disable early stopping: this is something we can consider, but it requires some thought on what would then be the stopping criterion - edit distance, alignment score dropoff to the best match, or similar? We want to make sure such a flag is clearly scoped, e.g. "report all alignments within -S of the best" rather than unbounded.

[tprodanov]

Thank you for looking into this!

I understand that my use-case is rather unusual, so no worries if this is hard to solve. I am writing some code right now that would transfer alignments from one haplotype to another by combining CIGARs, but my worry is that I may still miss some alignments so it would be more reliable if Strobealign could generate them as well.

In general I am fine with having bad alignments, I hoped setting -N and -M to high values and -S to something small would generate these suboptimal alignments, but that did not happen. In this specific case all reported alignments have perfect scores and I guess that all the dropped ones have at least one mismatch, maybe more. In terms of the arguments I am fine with all your suggestions. I guess ideally it would be difference in edit distance, but I am also fine with score dropoff or anything that is calculated before actual alignments are constructed.

[funnell]

I'm also interested in this (for metagenomics). Would setting -S to 0.0 help?