Add repetitive hits that are close to non-repetitive hits with a linear sweep

[ksahlin]

Implemented a sweep between non-repetitive hits (anchors) and filtered hits to recover relevant filtered hits before chaining. Do not add all filtered reads, but only select largely non-overlapping ones with a minimum hit-count with a Weighted Interval Scheduling (WIS) algorithm.

When adding nearly all seeds with the sweep algorithm (commit 1eefd33), the accuracy further increases, but does so at a large runtime cost (1.2-2.5x compared to main for 50-250nt reads, SE and PE).

The WIS selection of seeds reduces runtime overhead to be, on average, about 10-20% slower than main (but sometimes also 5-10% faster). However, it improves accuracy on nearly all instances in both extension and mapping-only modes. In particular, the mapping-only accuracy is substantially improved; see the benchmark below.

General idea

After building the standard (non-repetitive) anchor set, filtered hit's reference positions are added back if they are on the same rough diagonal as at least one existing non-repetitive anchor, within a tolerance of k bp.

The diagonal of a seed hit is ref_start - query_start. Co-linear seeds share the same
diagonal, so this filter admits exactly the repetitive positions that are consistent with
"trusted" non-repetitive anchors.

Instead of adding back all filtered seeds (runtime), we add those with the lowest count that are also not overlapping too much (on the query), so this selects a subset of seeds to add back. We do this with WIS.

Implementation

Three changes are made:

1. RefRandstrobe field reorder (src/index.rs)

The order of the fields in RefRandstrobe controls how the index is sorted. The field order is changed from (hash_offset, position, ref_index) to (hash_offset, ref_index, position), so within a hash group entries are now sorted by (ref_id, ref_start). This enables an efficient pure two-pointer sweep, as described below. The .sti file format version is bumped from 7 to 8; existing cached index files must be regenerated.

(Note: the algorithm also works with the old order (ref_start, ref_id), but then one must maintain a hash table keyed on ref_id with separate sweep pointers per chromosome. This adds complexity and would perform poorly on genomes with thousands of contigs such as metagenomic assemblies. The field reorder eliminates this entirely.)

2. add_repetitive_via_sweep in src/chainer.rs

Called after hits_to_anchors() and before chaining. The algorithm:

Step 1 — Build the eligible diagonal set.
Collect all anchors from non-repetitive (non-filtered) hits, compute each anchor's diagonal d = ref_start − query_start, and sort the resulting anchor_diags: Vec<(ref_id, diagonal)>. Co-linear seeds share the same diagonal value, so this set encodes trusted mapping regions.

Step 2 — Select filtered hits via Weighted Interval Scheduling (WIS).
Rather than sweeping all filtered hits, we first select a non-overlapping tiling on the query. Hits are treated as intervals [query_start, query_end] with weight CAP − count + 1 (higher weight for less-repetitive seeds; all weights are positive so the DP always selects at least one hit). The standard WIS DP runs in O(n log n) and recovers the maximum-weight subset of hits whose query intervals do not overlap by more than 10 bp. Seeds with more than 10,000 genomic positions are excluded before selection. This tiling step is critical: without it, overlapping filtered hits produce nearly-identical anchors that inflate the anchor set and double chaining time without improving accuracy.

Step 3 — Two-pointer sweep.
For each selected filtered hit at query position q, iterate its positions in the index (now sorted by (ref_id, ref_start) within the hash group). Maintain a single pointer into anchor_diags. Because both sequences increase in the same (ref_id, value) order, the pointer never goes backwards — a pure O(N + A) two-pointer sweep per selected hit, where N = repetitive positions and A = non-repetitive anchors. Any repetitive position whose diagonal falls within tolerance = k bp of an anchor diagonal is added to the anchor set.

This sweep is closely related to the paired-chain sweeping approach of @NicolasBuchin used in strobealign's PE mapping mode.

Complexity per read: O(A log A) for sorting anchor_diags once; O(n log n) for WIS selection over the typically 2–5 filtered hits; O(N + A) amortised per selected hit for the sweep. An early-exit skips the entire function when no filtered hits exist (the common case for reads in unique regions).

3. Removal of local rescue

The existing local rescue mechanism (rescuing filtered seeds in stretches longer than rescue_distance) is removed. The sweep now covers what rescue previously handled for reads in repetitive regions. A minimal global rescue is retained only for reads with fewer than 5 non-repetitive hits.

4. Added --filter-cutoff command-line parameter

Allows directly overriding the filter cutoff at runtime without recompiling, which proved useful for testing the effect of including more or fewer repetitive seeds.

Breaking change

STI_FILE_FORMAT_VERSION is bumped from 7 to 8 due to the RefRandstrobe field reorder.
Existing .sti index files will be rejected with a clear version-mismatch error and must
be regenerated with strobealign --create-index.

Benchmarks

WIS solution percent point (pp) change in accuracy (left) and percent change in runtime (right) compared to baseline. Upper panels are extension, lower mapping-only. Red lines for PE blue for SE.

SIM1 (--snp-rate 0.001 --small-indel-rate 0.0001 --max-small-indel-size 50)

CHM13

Maize

SIM5 (--snp-rate 0.01 --small-indel-rate 0.002 --max-small-indel-size 100)

CHM13

Maize

[ksahlin]

I have generated both SIM1 and SIM5-like variations and added plots instead of tables above, as the number of settings became quite large (some experiment data points are missing from the plots above). The commit is ready to be tested on our proper benchmark. I am sure we could still tweak some settings/parameters if needed (e.g., a 65% runtime overhead for maize 75bp low-variation may not be acceptable), but the overall picture looks quite favorable.

[marcelm]

I’m runnning the benchmarks now.

I’ve fixed the merge conflicts. Ok if I push to this branch or are you still working on this at the moment?

[ksahlin]

I am not working on this branch anymore. Feel free to push away.

Here are my two notes that I had on this branch that could be done/tested:

1. Remove the unused "local rescue function"? (Since this commit mostly adds lines, would be nice to remove some)
2. Test to save only 3 initial anchors instead of 5 in global rescue. This will, in practice, not affect reads >150nt, but may decrease runtime noticeably for short reads 50-100, where we currently have the largest runtime increase for the method.