General Supplementary Alignment for large varations inside reads by NicolasBuchin · Pull Request #682 · ksahlin/strobealign

NicolasBuchin · 2026-11-04T12:28:34Z

We currently only use supplementary alignments to report duplications in duplication detection. This PR introduces a new logic to search for general supplementary alignments that complete the primary alignment and explain the read by accounting for large variations.

By "general" I mean that these supplementary alignments are found independently of the primary: each one comes from its own chain and covers a part of the read the primary alignment does not explain, across a large insertion, deletion or duplication. This is different from the supplementary records produced by duplication detection, which only split an existing alignment into one record per copy and never align any uncovered part of the read.

Once the primary alignment is extended, we look for chains that fit a large variation pattern. Candidates must be on the same contig and the same strand (no translocations or inversions for now), and must fit one of two variation patterns: an insertion is close on the reference, within --si (200), and open on the read, a deletion or a duplication is close on the read and loose on the reference, up to --sd (100 kb). A candidate may also overlap the primary on the read by up to --so (0.5) of the primary's read span.

(This figure shows the patterns we look for in chain location compared to primary: left insertion, right duplication and deletion)

Candidates are extended in the same way as the primary alignment, and only over the interval not covered by the primary, with an end bonus to reach the primary on the read. The best extension is reported as the supplementary alignment if it has a score of at least --sm (default: 80) and meets the primary within --ss (50) bases, either on the read or on the reference. It gets its own MAPQ, computed from the second-best supplementary candidate. Duplication detection then runs over both primary and supplementary segments, and each copy it splits off keeps the MAPQ of the segment it came from.

Many of the candidates we need were not passed on after chaining, so extract_chains was updated. Chains are now grouped by the read interval they cover, --vp is applied to the best chain in each interval and --vn (1000) caps how many we keep per interval. A chain reaches that test if it scores at least --vs (40) or --vp of the read's best chain, whichever is lower, which means chains scoring below --vp are still kept if they cover a different read interval. This changes mapping for paired-end reads as well: no supplementary records there, but the chains we extend are not the same ones, which is why na:i and al:i move in the paired-end phix tests.

We also keep chains that share parents with a better chain. We cut these chains to retain only the anchors not used by the better chain. This helps recover small variations that could otherwise fit into a chain but reduce its score too much.

Two things I have not touched: the mapq for the primary alignment, sometimes the primary mapq is calculated with the second best score coming from the same chain as the supplementary alignment, we should change this in the future.

…iant Large SVs can't be fully spanned by chains, so we now report more than one alignment if it helps explain the read, either through a large insertion, deletion or duplication. The second alignment is always a supplementary of the primary. The duplication detection is done on each of the alignments: primary and supplementary. Is-new-baseline: yes

--vp used to be a fraction of the best chain on the whole read, and is now a fraction of the best chain over a read interval. Two chains fall in the same interval when they overlap by at least half of the shorter of the two. A chain also has to clear an absolute bar, the lower of --vs (default 40) and --vp of the best chain on the read, so that taking the fraction off the whole read does not admit every chain. Each interval has at most --vn (default 20) chains. Is-new-baseline: yes

Instead of extending a supplementary candidate on the full read, we now align only over the interval not covered by the primary, which saves runtime for supplementary extension. The extension stops at the primary, and has an end bonus to reach the primary's position on the read. Is-new-baseline: yes

A chain that ran into an anchor an earlier chain had already taken used to be dropped whole. It is now cut and scored on its own unique anchors and kept if that still clears the floor. Is-new-baseline: yes

Is-new-baseline: yes

marcelm

Here is a first round of review comments. There may be more, also because I haven’t understood everything, yet.

What I wonder about most is the introduction of so many new parameters. This type of supplementary alignments introduced in this PR is for very specific events, and I haven’t come to a conclusion for myself whether this is the right approach. I understand you want to be conservative now, and that is good, so I’m not arguing that anything needs to be changed for this PR. But it would be nice if a future, more general version required fewer parameters.

Reading through the code, I also noticed some refactoring that should definitely happen, but I’d be ok with doing this separately.

marcelm · 2026-11-05T13:59:25Z

        };
        sam_records.push(sam_output.make_mapped_record(
-            &segments[i],
+            &segments[i].alignment,


Suggestion: Let make_mapped_record accept a Segment (instead of Alignment plus mapq). Not sure how well this works at all call sites, though.

It only fits the primary/supplementary loop. The secondary alignments and the paired-end path (through make_record) work with plain Alignments, so they would have to convert each one into the struct just to call it. Maybe we can do some refactoring later, if we want paired-end to support supplementary for example.

marcelm · 2026-11-05T14:00:05Z

    let (mut ssw, mut pw) = if !mapping_parameters.use_ssw {
        (
-            extend_seed(aligner, chain, refseq, read, true, k).unwrap(),
+            extend_seed(aligner, chain, refseq, read, true, k, (0, read.len())).unwrap(),


TODO This is getting a bit hard to read.

extend_seed is back to it's original signature now. Also I think we can remove this function now, maybe a PR for it?

marcelm · 2026-11-05T14:01:33Z

    read: &Read,
    use_ssw: bool,
    k: usize,
+    query_interval: (usize, usize),


Now this uses two terms for the same thing: read and query. I’m fine with using query for both, that is, renaming read to query.

I changed the range toread_interval, I think it just makes more sense. Then query is the slice of the read that actually gets aligned in extend_seed_in_interval.

marcelm · 2026-11-05T14:22:23Z

+    n_chains: usize,
+}
+
+fn admit_chain(


This function needs a little bit of explanation (docstring). In particular

what it modifies in intervals

what it returns

I gave admit_chain and QueryInterval documentation, i think perhaps admit_chain should also be renamed? But I'm not sure how.

marcelm · 2026-11-12T13:36:51Z

+    /// Do not search for supplementary alignments
+    #[arg(long = "no-supplementary", help_heading = "Supplementary")]
+    no_supplementary: bool,
+
+    /// Minimum score for a supplementary alignment to be reported
+    #[arg(long = "sm", default_value_t = MappingParameters::default().min_supplementary_score, value_name = "N", help_heading = "Supplementary")]
+    min_supplementary_score: u32,
+
+    /// Maximum gap between the primary and a supplementary, on the read or the reference
+    #[arg(long = "ss", default_value_t = MappingParameters::default().max_segment_gap, value_name = "N", help_heading = "Supplementary")]
+    max_segment_gap: usize,
+
+    /// Longest deletion or tandem duplication a supplementary alignment may describe
+    #[arg(long = "sd", default_value_t = MappingParameters::default().max_gap_distance, value_name = "N", help_heading = "Supplementary")]
+    max_gap_distance: usize,
+
+    /// Maximum reference gap between the primary and an insertion candidate
+    #[arg(long = "si", default_value_t = MappingParameters::default().max_insertion_ref_gap, value_name = "N", help_heading = "Supplementary")]
+    max_insertion_ref_gap: usize,
+
+    /// Fraction of the primary's read span a supplementary candidate may overlap
+    #[arg(long = "so", default_value_t = MappingParameters::default().max_supplementary_overlap, value_name = "FRACTION", help_heading = "Supplementary")]
+    max_supplementary_overlap: f32,


We need to discuss which of these parameters we actually expose, but I’m fine with doing it in a follow-up PR.

Yes I think so too, a follow-up PR.

Co-authored-by: Marcel Martin <marcel.martin@scilifelab.se>

Use the existing alignment terminology for the alignment-with-mapq type and the supplementary search, and give the type an sa_entry() method.

Capping the number of chains kept per query interval at 20 made paired-end mapping regress in benchmarks. With 1000 the paired-end results no longer regress, and there is no visible runtime cost. Is-new-baseline: yes

NicolasBuchin marked this pull request as ready for review November 5, 2026 10:17

NicolasBuchin force-pushed the supplementary branch 2 times, most recently from 82d82b8 to 885fc08 Compare November 9, 2026 15:03

ksahlin added the high priority Extra attention is needed label Nov 10, 2026

NicolasBuchin force-pushed the supplementary branch from 5d351bc to 213c9ce Compare November 10, 2026 10:51

NicolasBuchin added 9 commits November 11, 2026 11:53

Report supplementary search statistics

5e3f650

Recover part of a chain sharing a parent with a better chain.

b447c3c

A chain that ran into an anchor an earlier chain had already taken used to be dropped whole. It is now cut and scored on its own unique anchors and kept if that still clears the floor. Is-new-baseline: yes

Supplementary candidates overlap tolerance with primary on the read

f0dfb24

Is-new-baseline: yes

Document supplementary alignments

7e8fe50

Add --no-supplementary to turn off the supplementary search

5734cda

Keep extend_seed's signature and add extend_seed_in_interval

c323508

NicolasBuchin force-pushed the supplementary branch from 213c9ce to c323508 Compare November 11, 2026 10:57

marcelm reviewed Nov 12, 2026

View reviewed changes

NicolasBuchin and others added 6 commits November 13, 2026 08:04

Apply review suggestions

572e145

Co-authored-by: Marcel Martin <marcel.martin@scilifelab.se>

Rename Segment to AlignmentWithMapq

7e4d2b5

Use the existing alignment terminology for the alignment-with-mapq type and the supplementary search, and give the type an sa_entry() method.

Clarify documentation of the supplementary search and chain filtering

03531de

Rename query_interval to read_interval in extend_seed_in_interval

abcf8ec

Use alignment terminology in duplication splitting

227643a

Raise the default --vn to 1000

6bc74de

Capping the number of chains kept per query interval at 20 made paired-end mapping regress in benchmarks. With 1000 the paired-end results no longer regress, and there is no visible runtime cost. Is-new-baseline: yes

NicolasBuchin force-pushed the supplementary branch from b7d87b1 to 6bc74de Compare November 13, 2026 11:25

Conversation

NicolasBuchin commented Nov 4, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

marcelm left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

This branch has not been deployed

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

NicolasBuchin commented Nov 4, 2026 •

edited

Loading