General Supplementary Alignment for large varations inside reads

CHANGES.md

@@ -17,6 +17,11 @@
 * #686: Avoid unbounded memory usage when data cannot be written to the
   output as fast as reads are mapped (that is, when the output file/stdout
   blocks).
+* #682: Report a supplementary alignment for a single-end read that crosses a
+  large insertion, deletion or tandem duplication, instead of soft-clipping the
+  part of the read that does not fit the primary alignment (options `--sm`,
+  `--ss`, `--sd`, `--si` and `--so`, see the README). Chains are now kept per read
+  interval (`--vs`, `--vn`), which also changes mapping for paired-end reads.
 
 ## v0.20.0 (2026-11-06)
 

README.md

@@ -281,6 +281,40 @@ controls the chaining look-back window heuristic, `--gd` sets the diagonal gap c
 `--gl` sets the gap length cost, `--vp` determines the best-chain score threshold, 
 and `--sg` controls the maximum skip distance on the reference.
 
+Chains are grouped by the read interval they cover, and only the best scoring
+chains of each interval are kept: `--vp` and a cap of `--vn` chains (1000 by
+default) apply per interval. A chain is also kept if it scores at least `--vs`
+(40 by default), even when it scores below `--vp` of the best chain found, so a
+chain that covers a part of the read the best chain misses is kept as well.
+
+### Supplementary alignments
+
+Strobealign searches for a single supplementary alignment when the primary
+alignment does not explain the whole read, which happens when a read crosses a
+large variation. Once the primary is extended, candidate chains are checked
+against two variation patterns: for an insertion, the reference carries on
+within `--si` bases (200 by default) of the primary while the read jumps the
+inserted bases, and for a deletion or a tandem duplication, the read carries on
+while the reference jumps up to `--sd` bases (100,000 by default). A candidate
+also has to be on the same contig and the same strand as the primary, so
+translocations and inversions are not detected. Candidate chains may overlap
+the primary on the read by a fraction `--so` (0.5 by default) of what the
+primary covers.
+
+Candidates are extended like the primary, but only over the read interval the
+primary does not cover, so that the end bonus `-L` pulls the extension towards
+the primary. Strobealign reports the best extension as a supplementary record
+if it scores at least `--sm` (80 by default) and meets the primary within
+`--ss` bases (50 by default), either on the read or on the reference. The
+supplementary gets its own mapping quality, computed from the second best
+candidate.
+
+Duplication detection is run on the supplementary alignment the same way it is
+run on the primary.
+
+The supplementary alignment search runs for single-end reads only. Use
+`--no-supplementary` to turn it off.
+
 ### Duplication detection
 
 Strobealign attempts to reinterpret insertions inside an alignment as tandem

src/chainer.rs

@@ -54,6 +54,8 @@ pub struct ChainingParameters {
     pub diag_diff_penalty: f32,
     pub gap_length_penalty: f32,
     pub valid_score_threshold: f32,
+    pub min_chain_score: f32,
+    pub max_chains_per_interval: usize,
     pub max_ref_gap: Option<usize>,
     pub matches_weight: f32,
     pub max_diagonal_ratio: f32,
@@ -66,6 +68,8 @@ impl Default for ChainingParameters {
             diag_diff_penalty: 0.1,
             gap_length_penalty: 0.05,
             valid_score_threshold: 0.7,
+            min_chain_score: 40.0,
+            max_chains_per_interval: 1000,
             max_ref_gap: None,
             matches_weight: 0.01,
             max_diagonal_ratio: 10.0,
@@ -452,6 +456,70 @@ fn hits_to_anchors(hits: &Vec<Hit>, index: &StrobemerIndex) -> Vec<Anchor> {
     anchors
 }
 
+/// A part of the query covered by chains, with the best chain score and the
+/// number of chains admitted for it
+struct QueryInterval {
+    start: usize,
+    end: usize,
+    best_score: f32,
+    n_chains: usize,
+}
+
+/// Decides whether a chain covering `start..end` on the query is kept, and
+/// records it in `intervals`.
+///
+/// The chain belongs to the first interval that overlaps it by at least half
+/// of the shorter of the two spans. If there is none, a new interval is added
+/// for the chain and the chain is kept.
+///
+/// Within its interval, a chain scoring higher than the interval's best is
+/// kept, and the interval takes over the chain's span and score. Otherwise the
+/// chain is rejected if it scores below `valid_score_threshold` times the
+/// interval's best, or if the interval already holds `max_chains_per_interval`
+/// chains. Each kept chain increments the interval's chain count.
+///
+/// Returns whether the chain is kept.
+fn admit_chain(
+    intervals: &mut Vec<QueryInterval>,
+    start: usize,
+    end: usize,
+    score: f32,
+    parameters: &ChainingParameters,
+) -> bool {
+    for interval in intervals.iter_mut() {
+        let overlap = interval
+            .end
+            .min(end)
+            .saturating_sub(interval.start.max(start));
+        if overlap * 2 < (end - start).min(interval.end - interval.start) {
+            continue;
+        }
+
+        if score > interval.best_score {
+            interval.start = start;
+            interval.end = end;
+            interval.best_score = score;
+        } else if score < interval.best_score * parameters.valid_score_threshold
+            || interval.n_chains >= parameters.max_chains_per_interval
+        {
+            return false;
+        }
+
+        interval.n_chains += 1;
+
+        return true;
+    }
+
+    intervals.push(QueryInterval {
+        start,
+        end,
+        best_score: score,
+        n_chains: 1,
+    });
+
+    true
+}
+
 impl ChainingResult {
     fn extract_chains(
         &self,
@@ -461,7 +529,8 @@ impl ChainingResult {
         chains: &mut Vec<Chain>,
     ) {
         let n = self.anchors.len();
-        let valid_score = self.best_score * self.parameters.valid_score_threshold;
+        let valid_score = (self.best_score * self.parameters.valid_score_threshold)
+            .min(self.parameters.min_chain_score);
 
         let mut candidates = vec![];
         for i in 0..n {
@@ -473,12 +542,12 @@ impl ChainingResult {
         candidates.sort_by(|a, b| b.1.total_cmp(&a.1));
 
         let mut used = vec![false; n];
-        // Each candidate yields at most one chain
-        chains.reserve(candidates.len());
-        for (i, score) in candidates {
+        let mut intervals = Vec::new();
+        for (i, mut score) in candidates {
             if used[i] {
                 continue;
             }
+            used[i] = true;
 
             let mut j = i;
             let mut overlaps = false;
@@ -488,11 +557,12 @@ impl ChainingResult {
             let mut ref_coverage = self.anchors[i].ref_start();
 
             while self.predecessors[j] != usize::MAX {
-                j = self.predecessors[j];
-                if used[j] {
+                let predecessor = self.predecessors[j];
+                if used[predecessor] {
                     overlaps = true;
                     break;
                 }
+                j = predecessor;
                 n_chain_anchors += 1;
                 used[j] = true;
 
@@ -503,6 +573,23 @@ impl ChainingResult {
             }
 
             if overlaps {
+                score = self.dp[i] - self.dp[j] + k as f32;
+                if score < valid_score {
+                    continue;
+                }
+            }
+
+            let first = &self.anchors[j];
+            let last = &self.anchors[i];
+            let query_start = first.query_start();
+            let query_end = last.query_start() + k;
+            if !admit_chain(
+                &mut intervals,
+                query_start,
+                query_end,
+                score,
+                &self.parameters,
+            ) {
                 continue;
             }
 
@@ -515,13 +602,10 @@ impl ChainingResult {
                 current = self.predecessors[current];
             }
 
-            let first = &self.anchors[j];
-            let last = &self.anchors[i];
-
             chains.push(Chain {
                 id: chains.len(),
-                query_start: first.query_start(),
-                query_end: last.query_start() + k,
+                query_start,
+                query_end,
                 ref_contig_start: contig_starts.ref_contig_start(first.ref_start()),
                 ref_start: first.ref_start(),
                 ref_end: last.ref_start() + k,
@@ -600,6 +684,49 @@ pub mod test {
         assert!(chaining_result.best_score > 42.0);
     }
 
+    #[test]
+    fn recover_lower_scoring_sub_chain_for_supplementary() {
+        let chainer = Chainer::new(20, ChainingParameters::default());
+        let starts = ContigStarts::new(vec![0], 100_000);
+        let extract = |anchors| {
+            let mut chains = vec![];
+            chainer
+                .collinear_chaining(anchors, &starts, 2000)
+                .extract_chains(20, false, &starts, &mut chains);
+            chains
+        };
+
+        #[rustfmt::skip]
+        let chains = extract(anchors![
+            (0, 0), (20, 20), (40, 40), (60, 60), (80, 80), (100, 100),
+            (1090, 200), (1110, 220), (1130, 240),
+        ]);
+        let alone = extract(anchors![(1090, 200), (1110, 220), (1130, 240)]);
+
+        assert_eq!(chains.len(), 2);
+        assert_eq!(alone.len(), 1);
+        let (cut, whole) = (&chains[1], &alone[0]);
+        assert_eq!(cut.anchors, whole.anchors);
+        assert_eq!(
+            (
+                cut.query_start,
+                cut.query_end,
+                cut.ref_start,
+                cut.ref_end,
+                cut.matching_bases
+            ),
+            (
+                whole.query_start,
+                whole.query_end,
+                whole.ref_start,
+                whole.ref_end,
+                whole.matching_bases
+            ),
+        );
+
+        assert!((cut.score - whole.score).abs() < 1e-4);
+    }
+
     #[test]
     fn linear_score_adjacent_anchors() {
         let chainer = Chainer::new(20, ChainingParameters::default());

src/details.rs

@@ -71,7 +71,20 @@ pub struct AlignmentDetails {
     /// No. of insertions re-interpreted as tandem duplications
     pub duplications: usize,
 
+    /// No. of chains beside a primary that fit a variant pattern and were kept as candidates
+    pub supplementary_candidates: usize,
+
+    /// No. of candidates whose extension produced an alignment
+    pub supplementary_extended: usize,
+
+    /// No. of supplementary alignments reported
+    pub supplementary: usize,
+
     pub time_extend: f64,
+
+    pub time_extend_primary: f64,
+
+    pub time_extend_supplementary: f64,
 }
 
 impl ops::AddAssign<AlignmentDetails> for AlignmentDetails {
@@ -84,7 +97,12 @@ impl ops::AddAssign<AlignmentDetails> for AlignmentDetails {
         self.best_alignments += rhs.best_alignments;
         self.duplication_candidates += rhs.duplication_candidates;
         self.duplications += rhs.duplications;
+        self.supplementary_candidates += rhs.supplementary_candidates;
+        self.supplementary_extended += rhs.supplementary_extended;
+        self.supplementary += rhs.supplementary;
         self.time_extend += rhs.time_extend;
+        self.time_extend_primary += rhs.time_extend_primary;
+        self.time_extend_supplementary += rhs.time_extend_supplementary;
     }
 }