From 506c374927b993e6b430bd71833d19951e17d3da Mon Sep 17 00:00:00 2001 From: NicolasBuchin Date: Sat, 10 Oct 2026 17:23:19 +0200 Subject: [PATCH 01/15] Report a supplementary alignment when a read crosses a structural variant 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 --- src/duplications.rs | 38 +++++++++++- src/lib.rs | 1 + src/main.rs | 22 +++++++ src/modes/align.rs | 75 +++++++++++++++++------ src/supplementary.rs | 139 +++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 256 insertions(+), 19 deletions(-) create mode 100644 src/supplementary.rs diff --git a/src/duplications.rs b/src/duplications.rs index 0b82589f1..7358791aa 100644 --- a/src/duplications.rs +++ b/src/duplications.rs @@ -1,12 +1,48 @@ +use std::cmp::Reverse; + use crate::aligner::{Aligner, Scores}; use crate::cigar::{Cigar, CigarOperation, OpLen}; use crate::details::AlignmentDetails; -use crate::modes::align::Alignment; +use crate::modes::align::{Alignment, MappingParameters, Segment}; use crate::packed_seq::PackedSeqSlice; use crate::read::Read; use crate::refseq::RefSequence; use crate::simdaligner::SplitReferenceAlignment; +/// Each segment holding a duplication split into one per copy, the copies carrying the mapping +/// quality of the segment they came out of. +pub fn split_duplications( + aligner: &Aligner, + refseq: &RefSequence, + read: &Read, + segments: Vec, + parameters: &MappingParameters, + details: &mut AlignmentDetails, +) -> Vec { + let mut resolved = Vec::new(); + for Segment { alignment, mapq } in segments { + let mut split = reinterpret_insertions( + aligner, + alignment, + refseq, + read, + parameters.min_duplication_length, + parameters.min_duplication_matches, + parameters.duplication_band, + details, + ); + + split.sort_by_key(|segment| Reverse(segment.score)); + resolved.extend( + split + .into_iter() + .map(|alignment| Segment { alignment, mapq }), + ); + } + + resolved +} + /// Cuts the alignment into one segment per copy wherever an insertion of at least `min_length` /// bases turns out to be a tandem duplication. A `min_length` of zero turns this off. Always /// returns at least the alignment it was given. diff --git a/src/lib.rs b/src/lib.rs index d35cb535c..d81f7f58a 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -22,3 +22,4 @@ pub mod seeding; pub mod shuffle; pub mod simdaligner; pub mod ssw; +pub mod supplementary; diff --git a/src/main.rs b/src/main.rs index 2999c3ff1..c7adcda73 100644 --- a/src/main.rs +++ b/src/main.rs @@ -299,6 +299,24 @@ struct Args { #[arg(long = "bd", default_value_t = 1036, value_name = "N", help_heading = "Alignment")] score_drop: u32, + // Supplementary parameters + + /// 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, + + /// 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, + /// Minimum insertion length that triggers tandem duplication detection (single-end only); 0 to disable #[arg(long = "dl", default_value_t = MappingParameters::default().min_duplication_length, value_name = "N", help_heading = "Duplications")] min_duplication_length: usize, @@ -584,6 +602,10 @@ fn run() -> Result<(), CliError> { args.max_secondary as usize }, secondary_threshold: args.secondary_threshold, + max_segment_gap: args.max_segment_gap, + max_gap_distance: args.max_gap_distance, + min_supplementary_score: args.min_supplementary_score, + max_supplementary_overlap: args.max_supplementary_overlap, max_tries: args.max_tries, dropoff_threshold: args.dropoff_threshold, rescue_distance: args.rescue_distance, diff --git a/src/modes/align.rs b/src/modes/align.rs index b88b6573b..468588bb8 100644 --- a/src/modes/align.rs +++ b/src/modes/align.rs @@ -18,7 +18,7 @@ use crate::chain::Chain; use crate::chainer::Anchor; use crate::cigar::{Cigar, CigarOperation}; use crate::details::{AlignmentDetails, ChainDetails, Details}; -use crate::duplications::reinterpret_insertions; +use crate::duplications::split_duplications; use crate::hit::HitsDetails; use crate::index::StrobemerIndex; use crate::insertsize::InsertSizeDistribution; @@ -34,6 +34,7 @@ use crate::read::Read; use crate::refseq::{ContigPosition, RefSequence}; use crate::revcomp::reverse_complement; use crate::seeding::SeedingParameters; +use crate::supplementary::get_supplementary_segment; const MAX_PAIR_CHAINS: usize = 1000; @@ -41,6 +42,10 @@ const MAX_PAIR_CHAINS: usize = 1000; pub struct MappingParameters { pub max_secondary: usize, pub secondary_threshold: f32, + pub max_segment_gap: usize, + pub max_gap_distance: usize, + pub min_supplementary_score: u32, + pub max_supplementary_overlap: f32, pub dropoff_threshold: f32, pub rescue_distance: usize, pub max_tries: usize, @@ -57,6 +62,10 @@ impl Default for MappingParameters { MappingParameters { max_secondary: 0, secondary_threshold: 0.95, + max_segment_gap: 50, + max_gap_distance: 100_000, + min_supplementary_score: 80, + max_supplementary_overlap: 0.5, dropoff_threshold: 0.5, rescue_distance: 100, max_tries: 20, @@ -126,9 +135,14 @@ impl Alignment { } } +/// An alignment with its mapping quality. +pub struct Segment { + pub alignment: Alignment, + pub mapq: u8, +} + fn build_sa_tags( - segments: &[Alignment], - mapq: u8, + segments: &[Segment], refseq: &RefSequence, cigar_eqx: bool, ) -> Vec> { @@ -138,7 +152,7 @@ fn build_sa_tags( let mut entries = Vec::with_capacity(segments.len()); for segment in segments { - entries.push(segment.sa_entry(refseq, mapq, cigar_eqx)); + entries.push(segment.alignment.sa_entry(refseq, segment.mapq, cigar_eqx)); } let mut tags = Vec::with_capacity(segments.len()); @@ -156,7 +170,7 @@ fn build_sa_tags( } impl AlignmentDetails { - fn add_alignment(&mut self, alignment: &Alignment) { + pub fn add_alignment(&mut self, alignment: &Alignment) { self.tried_alignment += 1; self.gapped += alignment.gapped as usize; self.collisions += alignment.collisions; @@ -541,20 +555,31 @@ pub fn align_single_end_read( }; let mapq = (60 * (best_score - second_best_score)).div_ceil(best_score) as u8; - let mut segments = reinterpret_insertions( + let supplementary = get_supplementary_segment( aligner, - best_alignment, refseq, &read, - mapping_parameters.min_duplication_length, - mapping_parameters.min_duplication_matches, - mapping_parameters.duplication_band, + index.k(), + &best_alignment, + chains, + mapping_parameters, &mut details.align, ); - // The best scoring segment becomes the primary and the rest supplementary - segments.sort_by_key(|segment| Reverse(segment.score)); - let mut sa_tags = build_sa_tags(&segments, mapq, refseq, sam_output.cigar_eqx); + let mut reported = vec![Segment { + alignment: best_alignment, + mapq, + }]; + reported.extend(supplementary); + let segments = split_duplications( + aligner, + refseq, + &read, + reported, + mapping_parameters, + &mut details.align, + ); + let mut sa_tags = build_sa_tags(&segments, refseq, sam_output.cigar_eqx); for i in 0..segments.len() { let status = if i == 0 { @@ -563,10 +588,10 @@ pub fn align_single_end_read( AlignmentStatus::Supplementary }; sam_records.push(sam_output.make_mapped_record( - &segments[i], + &segments[i].alignment, refseq, record, - mapq, + segments[i].mapq, status, details.clone(), sa_tags[i].take(), @@ -649,7 +674,7 @@ fn trace_log( /// Extend a chain so that it covers the entire read and return the resulting /// alignment. -fn extend_seed( +pub fn extend_seed( aligner: &Aligner, chain: &mut Chain, refseq: &RefSequence, @@ -1631,11 +1656,25 @@ mod tests { collisions: 0, }; - let tags = build_sa_tags(&[first.clone(), second], 60, &refseq, true); + let segments = [ + Segment { + alignment: first.clone(), + mapq: 60, + }, + Segment { + alignment: second, + mapq: 60, + }, + ]; + let tags = build_sa_tags(&segments, &refseq, true); assert_eq!(tags[0].as_deref(), Some("chr1,11,-,40S70=,60,0;")); assert_eq!(tags[1].as_deref(), Some("chr1,101,+,40=70S,60,0;")); - assert_eq!(build_sa_tags(&[first], 60, &refseq, true), vec![None]); + let only = [Segment { + alignment: first, + mapq: 60, + }]; + assert_eq!(build_sa_tags(&only, &refseq, true), vec![None]); } #[test] diff --git a/src/supplementary.rs b/src/supplementary.rs new file mode 100644 index 000000000..9656a6f96 --- /dev/null +++ b/src/supplementary.rs @@ -0,0 +1,139 @@ +use crate::aligner::Aligner; +use crate::chain::Chain; +use crate::details::AlignmentDetails; +use crate::modes::align::{Alignment, MappingParameters, Segment, extend_seed}; +use crate::read::Read; +use crate::refseq::RefSequence; + +/// Extend the chains that fit the pattern of a large variation with the +/// primary and return the best one. Returns `None` if no chain fits or clears +/// the score floor. +pub fn get_supplementary_segment( + aligner: &Aligner, + refseq: &RefSequence, + read: &Read, + k: usize, + primary: &Alignment, + chains: &mut [Chain], + parameters: &MappingParameters, + details: &mut AlignmentDetails, +) -> Option { + let mut candidates: Vec<&mut Chain> = chains + .iter_mut() + .filter(|chain| { + fits_variant_pattern( + primary, + refseq, + chain, + read.len(), + parameters.max_gap_distance as i64, + parameters.max_supplementary_overlap, + ) + }) + .collect(); + + let dropoff = candidates.first()?.score * parameters.dropoff_threshold; + let max_gap = parameters.max_segment_gap as i64; + let mut best_alignment = None; + let mut best_score = 0; + let mut second_best_score = 0; + + for chain in candidates.iter_mut().take(parameters.max_tries) { + if chain.score < dropoff { + break; + } + let Some(alignment) = extend_seed(aligner, chain, refseq, read, parameters.use_ssw, k) + else { + continue; + }; + details.add_alignment(&alignment); + + let (read_gap, ref_gap) = segment_gaps(primary, &alignment, read.len()); + if read_gap.abs() > max_gap && ref_gap.abs() > max_gap { + continue; + } + + if alignment.score > best_score { + second_best_score = best_score; + best_score = alignment.score; + best_alignment = Some(alignment); + } else if alignment.score > second_best_score { + second_best_score = alignment.score; + } + } + + if best_score < parameters.min_supplementary_score { + return None; + } + + let alignment = best_alignment?; + let mapq = (60 * best_score.saturating_sub(second_best_score)).div_ceil(best_score) as u8; + + Some(Segment { alignment, mapq }) +} + +/// Whether a chain and the primary are the two sides of a deletion, insertion +/// or duplication. +fn fits_variant_pattern( + primary: &Alignment, + refseq: &RefSequence, + chain: &Chain, + read_len: usize, + max_gap_distance: i64, + overlap_fraction: f32, +) -> bool { + let (contig_id, contig_start) = refseq.unflatten(chain.ref_start); + let contig_start = contig_start.0 as usize; + if primary.is_revcomp != chain.is_revcomp || primary.contig_id != contig_id { + return false; + } + + let (read_gap, ref_gap) = if chain.query_start >= primary.soft_clip_left { + ( + chain.query_start as i64 - (read_len - primary.soft_clip_right) as i64, + contig_start as i64 - (primary.contig_start + primary.length) as i64, + ) + } else { + ( + primary.soft_clip_left as i64 - chain.query_end as i64, + primary.contig_start as i64 - (contig_start + chain.ref_span()) as i64, + ) + }; + + let primary_span = (read_len - primary.soft_clip_left - primary.soft_clip_right) as i64; + let overlap_tolerance = (primary_span as f32 * overlap_fraction) as i64; + // A candidate has to bring read bases the primary does not have, which is also what + // keeps the chain the primary was built from out of its own search + if read_gap < 0 { + return false; + } + + // An insertion: the contig carries straight on while the read jumps the + // inserted bases + if (0..=overlap_tolerance).contains(&ref_gap) { + return true; + } + + // A deletion/duplication: the read carries straight on while the contig + // jumps bases + if read_gap <= overlap_tolerance && (-max_gap_distance..max_gap_distance).contains(&ref_gap) { + return true; + } + + false +} + +/// How far apart two alignments end up on the read and on the reference. +fn segment_gaps(primary: &Alignment, candidate: &Alignment, read_len: usize) -> (i64, i64) { + if candidate.soft_clip_left >= primary.soft_clip_left { + ( + candidate.soft_clip_left as i64 - (read_len - primary.soft_clip_right) as i64, + candidate.contig_start as i64 - (primary.contig_start + primary.length) as i64, + ) + } else { + ( + primary.soft_clip_left as i64 - (read_len - candidate.soft_clip_right) as i64, + primary.contig_start as i64 - (candidate.contig_start + candidate.length) as i64, + ) + } +} From 77b14ad791f5b5540721e218dfbe7b821fdb3c60 Mon Sep 17 00:00:00 2001 From: NicolasBuchin Date: Sat, 10 Oct 2026 17:23:19 +0200 Subject: [PATCH 02/15] Keep the best chains over each read interval --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 --- src/chainer.rs | 79 ++++++++++++++++++++++++++++++++++++++++++----- src/main.rs | 10 ++++++ tests/phix.pe.sam | 6 ++-- tests/phix.se.sam | 6 ++-- 4 files changed, 88 insertions(+), 13 deletions(-) diff --git a/src/chainer.rs b/src/chainer.rs index 770021096..509c53db7 100644 --- a/src/chainer.rs +++ b/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, 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: 20, max_ref_gap: None, matches_weight: 0.01, max_diagonal_ratio: 10.0, @@ -452,6 +456,54 @@ fn hits_to_anchors(hits: &Vec, index: &StrobemerIndex) -> Vec { anchors } +struct QueryInterval { + start: usize, + end: usize, + best_score: f32, + n_chains: usize, +} + +fn admit_chain( + intervals: &mut Vec, + 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 +513,8 @@ impl ChainingResult { chains: &mut Vec, ) { 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,8 +526,7 @@ 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()); + let mut intervals = Vec::new(); for (i, score) in candidates { if used[i] { continue; @@ -506,6 +558,20 @@ impl ChainingResult { 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; + } + // Walk the chain again to collect its anchors (from last to first) // into a Vec of exact size let mut chain_anchors = Vec::with_capacity(n_chain_anchors); @@ -515,13 +581,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, diff --git a/src/main.rs b/src/main.rs index c7adcda73..b9423e8e2 100644 --- a/src/main.rs +++ b/src/main.rs @@ -249,6 +249,14 @@ struct Args { #[arg(long = "vp", default_value_t = ChainingParameters::default().valid_score_threshold, help_heading = "Collinear chaining")] valid_score_threshold: f32, + /// Score at which a chain is kept whatever the best chain of the read scores + #[arg(long = "vs", default_value_t = ChainingParameters::default().min_chain_score, help_heading = "Collinear chaining")] + min_chain_score: f32, + + /// Maximum number of chains kept over one query interval + #[arg(long = "vn", default_value_t = ChainingParameters::default().max_chains_per_interval, help_heading = "Collinear chaining")] + max_chains_per_interval: usize, + /// Maximum allowed distance between anchors on the reference [default: same as read length] #[arg(long = "sg", help_heading = "Collinear chaining")] max_ref_gap: Option, @@ -622,6 +630,8 @@ fn run() -> Result<(), CliError> { diag_diff_penalty: args.diag_diff_penalty, gap_length_penalty: args.gap_length_penalty, valid_score_threshold: args.valid_score_threshold, + min_chain_score: args.min_chain_score, + max_chains_per_interval: args.max_chains_per_interval, max_ref_gap: args.max_ref_gap, matches_weight: args.matches_weight, max_diagonal_ratio: args.max_diagonal_ratio, diff --git a/tests/phix.pe.sam b/tests/phix.pe.sam index b4b2e4458..e8476211c 100644 --- a/tests/phix.pe.sam +++ b/tests/phix.pe.sam @@ -23,7 +23,7 @@ SRR1377138.10 99 NC_001422.1 1 60 26S275= = 58 358 NAAATGATTGGCGTATCCAACCTGCAGAG SRR1377138.10 147 NC_001422.1 58 60 300=1X = 1 -358 GAAAAATTATCTTGATAAAGCAGGAATTACTACTGCTTGTTTACGAATTAAATCGAAGTGGACTGCTGGCGGAAAATGAGAAAATTCGACCTATCCTTGCGCAGCTCGAGAAGCTCTTACTTTGCGACCTTTCGCCATCAACTAACGATTCTGTCAAAAACTGACGCGTTGGATGAGGAGAAGTGGCTTAATATGCTTGGCACGTTCGTCAAGGACTGGTTTAGATATGAGTCACATTTTGTTCATGGTAGAGATTCTCTTGTTGACATTTTAAAAGAGCGTGGATTACTATCTGAGTCCN ###FBDA<5)2?>FFDFFECC49FFD?3<@FFFFEDEC?3E?C?0EFFFFFFF>@5FECFFGGGFGD:FGGCDGFGGGGGG>GGGGGGGDGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGEGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCCB8# NM:i:1 AS:i:612 na:i:1 al:i:1 ga:i:0 X0:i:1 du:i:0 mr:i:0 RG:Z:1 SRR1377138.11 99 NC_001422.1 5117 60 1X268=32S = 5220 270 NAAGCTGTTCAGAATCAGAATGAGCCGCAACTTCGGGATGAAAATGCTCACAATGACAAATCTGTCCACGGAGTGCTTAATCCAACTTACCAAGCTGGGTTACGACGCGACGCCGTTCAACCAGATATTGAAGCAGAACGCAAAAAGAGAGATGAGATTGAGGCTGGGAAAAGTTACTGTAGCCGACGTTTTGGCGGCGCAACCTGTGACGACAAATCTGCTCAAATTTATGCGCGCTTCGATAAAAATGATTGGCGTATCCAACCTGCGGAGTTTTATCGCTTCCATGACGCAGAAGTTA #8ACCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGDGGGGF8FFGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGFGGGCFGGGGGGGGGGEE5EGGGGGGGGGFGGGGFDGEGDGGDGGGE7FFGGG>*?*:<>9BDECDFGF<*9@F?*7B?;F;F;B:AFFFFFF7*(03:0:AAAFFF:18?BA:82><9B>4<2BAA2A<<:?:4,.(B?C?.37;)EE<;)C?6(((;FFFFFEE?7)?9/);FDFB>A=FDAFFFFFFFDEDDCCFAGFFGGGGGGF:GGGFGGGGGFGGGGFGCGGEECDEFCCGGFGGGGFGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGFFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCCA8# NM:i:0 AS:i:344 na:i:2 al:i:1 ga:i:1 X0:i:1 du:i:0 mr:i:0 RG:Z:1 -SRR1377138.12 99 NC_001422.1 1 60 57S238=1X5= = 28 328 NAATCTGCTCAAATTTATGCGCGCTTCGATAAAAATGATTGGCGTATCCAACCTGCAGAGTTTTATCGCTTCCATGACGCAGAAGTTAACACTTTCGGATATTTCTGATGAGTCGAAAAATTATCTTGATAAAGCAGGAATTACTACTGCTTGTTTACGAATTAAATCGAAGTGGACTGCTGGCGGAAAATGAGAAAATTCGACCTATCCTTGCGCAGCTCGAGAAGCTCTTACTTTGCGACCTTTCGCCATCAACTAACGATTCTGTCAAAAACTGACGCGTTGGATGAGGAGATGTGGC #8BCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGDCGGGGGGGGGGFGGGGGGGFGFGGFGFGFFFFFFFFFFFAFFFFFFFFFFFA3.1.48?A<:?FFFFFFBDFFFF?1GGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCCB8# NM:i:2 AS:i:602 na:i:1 al:i:1 ga:i:0 X0:i:1 du:i:0 mr:i:0 RG:Z:1 SRR1377138.13 99 NC_001422.1 5006 60 1X300= = 5049 344 NTATGTGGCTAAATACGTTAACAAAAAGTCAGATATGGACCTTGCTGCTAAAGGTCTAGGAGCTAAAGAATGGAACAACTCACTAAAAACCAAGCTGTCGCTACTTCCCAAGAAGCTGTTCAGAATCAGAATGAGCCGCAACTTCGGGATGAAAATGCTCACAATGACAAATCTGTCCACGGAGTGCTTAATCCAACTTACCAAGCTGGGTTACGACGCGACGCCGTTCAACCAGATATTGAAGCAGAACGCAAAAAGAGAGATGAGATTGAGGCTGGGAAAAGTTACTGTAGCCGACGTT #8BCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFCFGGFGFGFFFFFFFFFFFDFFFFB:>FFFFF2.1.4::?FFFFF?BAAFFFFFFDFF>>9>9? NM:i:1 AS:i:612 na:i:1 al:i:1 ga:i:0 X0:i:1 du:i:0 mr:i:0 RG:Z:1 SRR1377138.13 147 NC_001422.1 5049 60 300=1X = 5006 -344 GCTGCTAAAGGTCTAGGAGCTAAAGAATGGAACAACTCACTAAAAACCAAGCTGTCGCTACTTCCCAAGAAGCTGTTCAGAATCAGAATGAGCCGCAACTTCGGGATGAAAATGCTCACAATGACAAATCTGTCCACGGAGTGCTTAATCCAACTTACCAAGCTGGGTTACGACGCGACGCCGTTCAACCAGATATTGAAGCAGAACGCAAAAAGAGAGATGAGATTGAGGCTGGGAAAAGTTACTGTAGCCGACGTTTTGGCGGCGCAACCTGTGACGACAAATCTGCTCAAATTTATGN #FFFFBFFFBD>A<::64>AAFFFB><?FFBABBAFFFFFFBAFFFFFFFFFFFFFE;@F@FFFFFFFD@FE?FFFFFFFFFFFFFEBFFFFFFFFGFFGGGDGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGDGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCGGCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCCC8# NM:i:1 AS:i:612 na:i:1 al:i:1 ga:i:0 X0:i:1 du:i:0 mr:i:0 RG:Z:1 @@ -59,8 +59,8 @@ SRR1377138.28 99 NC_001422.1 3016 60 1X116=1X183= = 3273 558 NCTCCGGTGGCATTCAAGG SRR1377138.28 147 NC_001422.1 3273 60 52=1X247=1X = 3016 -558 AAGGATACTCGTGATTATCTTGCTGCTGCATTTCCTGAGCTTAATGCTTGGGGGCGTGCTGGTGCTGATGCTTCCTCTGCTGGTATGGTTGACGCCGGATTTGAGAATCAAAAAGAGCTTACTAAAATGCAACTGGACAATCAGAAAGAGATTGCCGAGATGCAAAATGAGACTCAAAAAGAGATTGCTGGCATTCAGTCGGCGACTTCACGCCAGAATACGAAAGACCAGGTATATGCACAAAATGAGATGCTTGCTTATCAACAGAAGGAGTCTACTGCTCGCGTTGCGTCTATTATGN #####F?9FFFD<4<6,22:(FFFDA((FFFFEAFFFEAFDDFFFFEEFAFFFFFCDGGGGGGGFDDFGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGEFFCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGFGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCCA8# NM:i:2 AS:i:602 na:i:1 al:i:1 ga:i:0 X0:i:1 du:i:0 mr:i:0 RG:Z:1 SRR1377138.29 99 NC_001422.1 155 60 1X300= = 250 396 NGCGCAGCTCGAGAAGCTCTTACTTTGCGACCTTTCGCCATCAACTAACGATTCTGTCAAAAACTGACGCGTTGGATGAGGAGAAGTGGCTTAATATGCTTGGCACGTTCGTCAAGGACTGGTTTAGATATGAGTCACATTTTGTTCATGGTAGAGATTCTCTTGTTGACATTTTAAAAGAGCGTGGATTACTATCTGAGTCCGATGCTGTTCAACCACTAATAGGTAAGAAATCATGAGTCAAGTTACTGAACAATCCGTACGTTTCCAGACCGCTTTGGCCTCTATTAAGCTCATTCAG #8BCCGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGFGGGGGGGGGGGGFGGGGGGGDGGGGGGGGGGFFGFGGGGGGGGGGGGGGGGGEGGGGGFGFGGFGGFGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFGFACFGGFGEFGGGGFEGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGFGGGGGGG6C<>=4A:AEDFFFFFFFFFFFFFGFFGGGGGGGGGGGGGGGGGGGGGGGGGGCFGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGFGGGGFCGEGGGGGGGGGGD9FGGGGGGGGGFGFFGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGEGGFCCB8# NM:i:3 AS:i:592 na:i:1 al:i:1 ga:i:0 X0:i:1 du:i:0 mr:i:0 RG:Z:1 -SRR1377138.30 81 NC_001422.1 1 60 97S203=1X = 5176 5386 TTACTGTAGCCGACGTTTTGGCGGCGCAACCTGTGACGACAAATCTGCTCAAATTTATGCGCGCTTCGATAAAAATGATTGGCGTATCCAACCTGCAGAGTTTTATCGCTTCCATGACGCAGAAGTTAACACTTTCGGATATTTCTGATGAGTCGAAAAATTATCTTGATAAAGCAGGAATTACTACTGCTTGTTTACGAATTAAATCGAAGTGGACTGCTGGCGGAAAATGAGAAAATTCGACCTATCCTTGCGCAGCTCGAGAAGCTCTTACTTTGCGACCTTTCGCCATCAACTAACN FFFBFF?4(>7,.6FF?;96FFFB>??450.3:FFFFFBB,>EGFGGGGFGCEGCFEEEGDEGGGFGGGCF4CFGGGGC3*85CFF4:3?@FF6*.9*96<4:>?F24><?FB9AAFBFF2796;??BFFF99:A96-)))46?F>(:0<>4>5<2))-5>BF7,.6FF?;96FFFB>??450.3:FFFFFBB,>EGFGGGGFGCEGCFEEEGDEGGGFGGGCF4CFGGGGC3*85CFF4:3?@FF6*.9*96<4:>?F24><?FB9AAFBFF2796;??BFFF99:A96-)))46?F>(:0<>4>5<2))-5>BFBFBFFFAF?FF9B? NM:i:1 AS:i:612 na:i:1 al:i:1 ga:i:0 X0:i:1 du:i:0 mr:i:0 RG:Z:1 SRR1377138.31 147 NC_001422.1 4401 60 300=1X = 4320 -382 GCTGAACGCCCTCTTAAGGATATTCGCGATGAGTATAATTACCCCAAAAAGAAAGGTATTAAGGATGAGTGTTCAAGATTGCTGGAGGCCTCCACTATGAAATCGCGTAGAGGCTTTGCTATTCAGCGTTTGATGAATGCAATGCGACAGGCTCATGCTGATGGTTGGTTTATCGTTTTTGACACTCTCACGTTGGCTGACGACCGATTAGAGGCGTTTTATGATAATCCCAATGCTTTGCGTGACTATTTTCGTGATATTGGTCGTATGGTTCTTGCTGCCGAGGGTCGCAAGGCTAATN ################AA4)19>B913<96.))6<6<976;6B>9EGGGGGGGGGFGFEF9GEGGGGGGGGGDGGGFGGGGGGGGGGGGGGGFGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCCB8# NM:i:1 AS:i:612 na:i:1 al:i:1 ga:i:0 X0:i:1 du:i:0 mr:i:0 RG:Z:1 SRR1377138.32 99 NC_001422.1 1433 60 1X268=1X31= = 1487 355 NCCGAGGCTAACCCTAATGAGCTTAATCAAGATGATGCTCGTTATGGTTTCCGTTGCTGCCATCTCAAAAACATTTGGACTGCTCCGCTTCCTCCTGAGACTGAGCTTTCTCGCCAAATGACGACTTCTACCACATCTATTGACATTATGGGTCTGCAAGCTGCTTATGCTAATTTGCATACTGACCAAGAACGTGATTACTTCATGCAGCGTTACCATGATGTTATTTCTTCATTTGGAGGTAAAACCTCTTATGACGCTGACAACCGCCCTTTACTTGTCATGCGCTCTAATCTCTGGG #8BCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGFGGFGFGFFFFFFFFFFFFFFFFFFFFFFFF*.:6=;CFFFFFEFFEF9;(3:AFFB)66>?: NM:i:2 AS:i:602 na:i:1 al:i:1 ga:i:0 X0:i:1 du:i:0 mr:i:0 RG:Z:1 diff --git a/tests/phix.se.sam b/tests/phix.se.sam index f091aab66..9cc645182 100644 --- a/tests/phix.se.sam +++ b/tests/phix.se.sam @@ -12,7 +12,8 @@ SRR1377138.8 16 NC_001422.1 1254 60 124=1X132=1X42=1X * 0 0 TCATGAAGGATGGTGTTAAT SRR1377138.9 16 NC_001422.1 4731 60 300=1X * 0 0 TGTGTGCCTGAGTATGGTACAGCTAATGGCCGTCTTCATTTCCATGCGGTGCACTTTATGCGGACACTTCCTACAGGTAGCGTTGACCCTAATTTTGGTCGTCGGGTACGCAATCGCCGCCAGTTAAATAGCTTGCAAAATACGTGGCCTTATGGTTACAGTATGCCCATCGCAGTTCGCTACACGCAGGACGCTTTTTCACGTTCTGGTTGGTTGTGGCCTGTTGATGCTAAAGGTGAGCCGCTTAAAGCTACCAGTTATATGGCTGTTGGTTTCTATGTGGCTAAATACGTTAACAAAN ?BFF?<:?*?*:<>9BDECDFGF<*9@F?*7B?;F;F;B:AFFFFFF7*(03:0:AAAFFF:18?BA:FFFFF2.1.4::?FFFFF?BAAFFFFFFDFF>>9>9? NM:i:1 AS:i:612 na:i:1 al:i:1 ga:i:0 X0:i:0 du:i:0 RG:Z:1 SRR1377138.14 0 NC_001422.1 66 60 1X300= * 0 0 NTCTTGATAAAGCAGGAATTACTACTGCTTGTTTACGAATTAAATCGAAGTGGACTGCTGGCGGAAAATGAGAAAATTCGACCTATCCTTGCGCAGCTCGAGAAGCTCTTACTTTGCGACCTTTCGCCATCAACTAACGATTCTGTCAAAAACTGACGCGTTGGATGAGGAGAAGTGGCTTAATATGCTTGGCACGTTCGTCAAGGACTGGTTTAGATATGAGTCACATTTTGTTCATGGTAGAGATTCTCTTGTTGACATTTTAAAAGAGCGTGGATTACTATCTGAGTCCGATGCTGTT #8@CCGGGGGGFFGGGGGGGGGGGGGFGGGGGGGGGGGGGFGGEEFGGGGGGGGGGGGGFGDCFEFGGGGGFGGGGGGGGGGGGGGGGGGGG7FGGGGGGGGGGGGGGGGGGGDFGGGGGGGFGGEGGGGGFGGGGGGGGGGGGGGGGGFGGGGGGGGGGDGGGGGGGGDGGGGGGGFEGGGGGGGFGGFFGFGGGGGGGGGFGGGFG7FGGFGDGGGGGFFCGGGGGGGGGGGGGFGFGGC>FGCF=A@FFFFFFFFFFEFFFFFFE<=/04(34:501:AB>4 NM:i:1 AS:i:612 na:i:1 al:i:1 ga:i:0 X0:i:0 du:i:0 RG:Z:1 SRR1377138.15 0 NC_001422.1 1907 60 1X185=115S * 0 0 NTTAACGCTAAAGGTGCTTTGACTTATACCGATATTGCTGGCGACCCTGTTTTGTATGGCAACTTGCCGCCGCGTGAAATTTCTATGAAGGATGTTTTCCGTTCTGGTGATTCGTCTAAGAAGTTTAAGATTGCTGAGGGTCAGTGGTATCGTTATGCGCCTTCGTATGTTTCTCCTGCTTATCACAAGAACGGAAAACATCCTTCATAGAAATTTCACGCAGATCGGAAGAGCGGTTCAGCAGGAATGCCGAGACCGATCTCGTATGCCGTCTTCTGCTTGAAAAAAAAAAAAAAAGATA #8BCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGEGFGGFGFGFFFFFFFFFFFFF>73?FFFFFFF3-048:BDFFAFFFFFF06>BFF######### NM:i:1 AS:i:372 na:i:1 al:i:1 ga:i:1 X0:i:0 du:i:0 RG:Z:1 @@ -30,7 +31,8 @@ SRR1377138.26 16 NC_001422.1 1642 60 300=1X * 0 0 GCGTTACCATGATGTTATTTCTTCATTTGG SRR1377138.27 16 NC_001422.1 4289 60 31=1X268=1X * 0 0 TGTGCTCCCCCAACTTGATATTAATAACACTTTAGACCACCGCCCCGAAGGGGACGAAAAATGGTTTTTAGAGAACGAGAAGACGGTTACGCAGTTTTGCCGCAAGCTGGCTGCTGAACGCCCTCTTAAGGATATTCGCGATGAGTATAATTACCCCAAAAAGAAAGGTATTAAGGATGAGTGTTCAAGATTGCTGGAGGCCTCCACTATGAAATCGCGTAGAGGCTTTGCTATTCAGCGTTTGATGAATGCAATGCGACAGGCTCATGCTGATGGTTGGTTTATCGTTTTTGACACTCTN 0FFF>>96BBAE9FFFF@FFFDFFFFGFGFGGDGFGGDGGGGGGGGGFGGGGEGGGGGGGGGGGGGGGGEGDFDGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFGFGGGFGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGEGGGGGGGGGGFGGGGGGGGGGGGGGGGF?GGGGGGGGGGGGCCA8# NM:i:2 AS:i:602 na:i:1 al:i:1 ga:i:0 X0:i:0 du:i:0 RG:Z:1 SRR1377138.28 0 NC_001422.1 3016 60 1X116=1X183= * 0 0 NCTCCGGTGGCATTCAAGGTGATGTGCTTGCTACCGATAACAATACTGTAGGCATGGGTGATGCTGGTATTAAATCTGCCATTCAAGGCTCTAATGTTCCTAACCCTGATGAGGCCGTCCCTAGTTTTGTTTCTGGTGCTATGGCTAAAGCTGGTAAAGGACTTCTTGAAGGTACGTTGCAGGCTGGCACTTCTGCCGTTTCTGATAAGTTGCTTGATTTGGTTGGACTTGGTGGCAAGTCTGCCGCTGATAAAGGAAAGGATACTCGTGATTATCTTGCTGCTGCATTTCCTGAGCTTAA #8ACCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGCFGGGGGGGGGGGGGGFFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGEGGGGGGGGGGFGGGGGGGFGGFGGGGGDFGGFGGGGGFGFGGFGFGFFFFFFFFFFF:FEF@@F@=FEF@8/;?9))9CCEFFFFEFFEFF<7,.6FF?;96FFFB>??450.3:FFFFFB7,.6FF?;96FFFB>??450.3:FFFFFB7,.6FF?;96FFFB>??450.3:FFFFFBBFBFFFAF?FF9B? NM:i:1 AS:i:612 na:i:1 al:i:1 ga:i:0 X0:i:0 du:i:0 RG:Z:1 SRR1377138.32 0 NC_001422.1 1433 60 1X268=1X31= * 0 0 NCCGAGGCTAACCCTAATGAGCTTAATCAAGATGATGCTCGTTATGGTTTCCGTTGCTGCCATCTCAAAAACATTTGGACTGCTCCGCTTCCTCCTGAGACTGAGCTTTCTCGCCAAATGACGACTTCTACCACATCTATTGACATTATGGGTCTGCAAGCTGCTTATGCTAATTTGCATACTGACCAAGAACGTGATTACTTCATGCAGCGTTACCATGATGTTATTTCTTCATTTGGAGGTAAAACCTCTTATGACGCTGACAACCGCCCTTTACTTGTCATGCGCTCTAATCTCTGGG #8BCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGFGGFGFGFFFFFFFFFFFFFFFFFFFFFFFF*.:6=;CFFFFFEFFEF9;(3:AFFB)66>?: NM:i:2 AS:i:602 na:i:1 al:i:1 ga:i:0 X0:i:0 du:i:0 RG:Z:1 SRR1377138.33 0 NC_001422.1 3817 60 1X268=1X31= * 0 0 NTATTTTTCATGGTATTGATAAAGCTGTTGCCGATACTTGGAACAATTTCTGGAAAGACGGTAAAGCTGATGGTATTGGCTCTAATTTGTCTAGGAAATAACCGTCAGGATTGACACCCTCCCAATTGTATGTTTTCATGCCTCCAAATCTTGGAGGCTTTTTTATGGTTCGTTCTTATTACCCTTCTGAATGTCACGCTGATTATTTTGACTTTGAGCGTATCGAGGCTCTTAAACCTGCTATTGAGGCTTGTGGCATTTCTACTCTTCCTCAATCCCCAATGCTTGGCTTCCATAAGCA #8BCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGFEFGGGGGGGFGGFGGGGGGFFFGGFGFGFFFAFFFFFFBDFEFFEFFFFFFC,)217@EFFFDF<@EFFFF?FFFDE####### NM:i:2 AS:i:602 na:i:1 al:i:1 ga:i:0 X0:i:0 du:i:0 RG:Z:1 From 5e3f650802953464e181db80ce2f745c670b44b0 Mon Sep 17 00:00:00 2001 From: NicolasBuchin Date: Tue, 27 Oct 2026 17:38:41 +0100 Subject: [PATCH 03/15] Report supplementary search statistics --- src/details.rs | 18 ++++++++++++++++++ src/main.rs | 25 +++++++++++++++++++++++++ src/modes/align.rs | 5 +++++ src/supplementary.rs | 3 +++ 4 files changed, 51 insertions(+) diff --git a/src/details.rs b/src/details.rs index 74a7f24f7..c46df3bff 100644 --- a/src/details.rs +++ b/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 for AlignmentDetails { @@ -84,7 +97,12 @@ impl ops::AddAssign 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; } } diff --git a/src/main.rs b/src/main.rs index b9423e8e2..0e02150d4 100644 --- a/src/main.rs +++ b/src/main.rs @@ -990,6 +990,21 @@ fn print_statistics( "Insertions reinterpreted as duplications: {:12}", details.align.duplications ); + debug!( + "Supplementary candidates: {:12} Per read: {:7.1}", + details.align.supplementary_candidates, + details.align.supplementary_candidates as f64 / details.chain.n_reads as f64 + ); + debug!( + "Supplementary candidates extended: {:12} Per read: {:7.1}", + details.align.supplementary_extended, + details.align.supplementary_extended as f64 / details.chain.n_reads as f64 + ); + debug!( + "Supplementary alignments reported: {:12} Per read: {:7.1}", + details.align.supplementary, + details.align.supplementary as f64 / details.chain.n_reads as f64 + ); debug!(""); debug!("## Mapping runtime breakdown"); @@ -1034,6 +1049,16 @@ fn print_statistics( details.align.time_extend, details.align.time_extend * 100. / total ); + debug!( + " primary chains: {:8.2} s {:5.1} %", + details.align.time_extend_primary, + details.align.time_extend_primary * 100. / total + ); + debug!( + " supplementary search: {:8.2} s {:5.1} %", + details.align.time_extend_supplementary, + details.align.time_extend_supplementary * 100. / total + ); debug!( "Total: {:8.2} s {:5.1} %", total, 100. diff --git a/src/modes/align.rs b/src/modes/align.rs index 468588bb8..638cf0edc 100644 --- a/src/modes/align.rs +++ b/src/modes/align.rs @@ -488,6 +488,7 @@ pub fn align_single_end_read( let mut best_alignment = None; let read = Read::new(&record.sequence); + let primary_timer = Instant::now(); for (tries, chain) in chains .iter_mut() .take(mapping_parameters.max_tries) @@ -547,6 +548,8 @@ pub fn align_single_end_read( alignments.push(alignment); } } + details.align.time_extend_primary += primary_timer.elapsed().as_secs_f64(); + let Some(best_alignment) = best_alignment else { return ( vec![sam_output.make_unmapped_record(record, details.clone())], @@ -555,6 +558,7 @@ pub fn align_single_end_read( }; let mapq = (60 * (best_score - second_best_score)).div_ceil(best_score) as u8; + let supplementary_timer = Instant::now(); let supplementary = get_supplementary_segment( aligner, refseq, @@ -565,6 +569,7 @@ pub fn align_single_end_read( mapping_parameters, &mut details.align, ); + details.align.time_extend_supplementary += supplementary_timer.elapsed().as_secs_f64(); let mut reported = vec![Segment { alignment: best_alignment, diff --git a/src/supplementary.rs b/src/supplementary.rs index 9656a6f96..0a9fb6bb0 100644 --- a/src/supplementary.rs +++ b/src/supplementary.rs @@ -31,6 +31,7 @@ pub fn get_supplementary_segment( ) }) .collect(); + details.supplementary_candidates += candidates.len(); let dropoff = candidates.first()?.score * parameters.dropoff_threshold; let max_gap = parameters.max_segment_gap as i64; @@ -46,6 +47,7 @@ pub fn get_supplementary_segment( else { continue; }; + details.supplementary_extended += 1; details.add_alignment(&alignment); let (read_gap, ref_gap) = segment_gaps(primary, &alignment, read.len()); @@ -68,6 +70,7 @@ pub fn get_supplementary_segment( let alignment = best_alignment?; let mapq = (60 * best_score.saturating_sub(second_best_score)).div_ceil(best_score) as u8; + details.supplementary += 1; Some(Segment { alignment, mapq }) } From 1b7df9eb194e9765c7042cf98d7aeaefe539b15a Mon Sep 17 00:00:00 2001 From: NicolasBuchin Date: Sat, 10 Oct 2026 17:39:45 +0200 Subject: [PATCH 04/15] Extend a supplementary over the interval the primary does not cover 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 --- src/modes/align.rs | 72 ++++++++++++++++++++++++++++++++++---------- src/supplementary.rs | 19 ++++++++++-- tests/phix.se.sam | 8 ++--- 3 files changed, 77 insertions(+), 22 deletions(-) diff --git a/src/modes/align.rs b/src/modes/align.rs index 638cf0edc..43ca63b9a 100644 --- a/src/modes/align.rs +++ b/src/modes/align.rs @@ -509,6 +509,7 @@ pub fn align_single_end_read( &read, mapping_parameters.use_ssw, index.k(), + (0, read.len()), ) else { continue; }; @@ -652,13 +653,13 @@ fn trace_log( 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(), alignment.clone(), ) } else { ( alignment.clone(), - extend_seed(aligner, chain, refseq, read, false, k).unwrap(), + extend_seed(aligner, chain, refseq, read, false, k, (0, read.len())).unwrap(), ) }; // manually adds the soft clips @@ -686,20 +687,26 @@ pub fn extend_seed( read: &Read, use_ssw: bool, k: usize, + query_interval: (usize, usize), ) -> Option { - let query = if chain.is_revcomp { + let full = if chain.is_revcomp { read.rc() } else { read.seq() }; + let (start, end) = query_interval; + debug_assert!(start <= chain.query_start && chain.query_end <= end,); + let query = &full[start..end]; + let chain_start = chain.query_start.saturating_sub(start); + let chain_end = chain.query_end.saturating_sub(start).min(query.len()); let (contig_id, contig_start) = refseq.unflatten(chain.ref_start); let contig_start = contig_start.0 as usize; let contig_end = contig_start + chain.ref_span(); let contig = refseq.contig(contig_id); assert!(contig_end <= contig.len()); - let projected_start = contig_start.saturating_sub(chain.query_start); - let projected_end = min(contig_end + query.len() - chain.query_end, contig.len()); + let projected_start = contig_start.saturating_sub(chain_start); + let projected_end = min(contig_end + (query.len() - chain_end), contig.len()); // TODO ugly let mut info = AlignmentInfo { @@ -752,11 +759,10 @@ pub fn extend_seed( .iter() .filter_map(|a| { let ref_start = a.ref_start() - decode_start - refseq.starts.0[contig_id]; + let query_start = a.query_start() - start; // Skip anchors resulting from hash collisions - if query[a.query_start()..a.query_start() + k] - == segment[ref_start..ref_start + k] - { - Some(Anchor::new(ref_start, a.query_start())) + if query[query_start..query_start + k] == segment[ref_start..ref_start + k] { + Some(Anchor::new(ref_start, query_start)) } else { collisions += 1; @@ -772,8 +778,8 @@ pub fn extend_seed( Some(Alignment { cigar: mem::take(&mut info.cigar), edit_distance: info.edit_distance, - soft_clip_left: info.query_start, - soft_clip_right: query.len() - info.query_end, + soft_clip_left: start + info.query_start, + soft_clip_right: full.len() - start - info.query_end, score: info.score, contig_id, contig_start: result_start, @@ -1003,8 +1009,24 @@ fn extend_paired_seeds( let mut ch_max1 = chains[0][0].clone(); let mut ch_max2 = chains[1][0].clone(); - let alignment1 = extend_seed(aligner, &mut ch_max1, refseq, read1, use_ssw, k); - let alignment2 = extend_seed(aligner, &mut ch_max2, refseq, read2, use_ssw, k); + let alignment1 = extend_seed( + aligner, + &mut ch_max1, + refseq, + read1, + use_ssw, + k, + (0, read1.len()), + ); + let alignment2 = extend_seed( + aligner, + &mut ch_max2, + refseq, + read2, + use_ssw, + k, + (0, read2.len()), + ); if let (Some(alignment1), Some(alignment2)) = (alignment1, alignment2) { details[0].add_alignment(&alignment1); details[1].add_alignment(&alignment2); @@ -1028,7 +1050,15 @@ fn extend_paired_seeds( // the paired-end read as two single-end reads. let mut a_indv_max = [None, None]; for i in 0..2 { - a_indv_max[i] = extend_seed(aligner, &mut chains[i][0], refseq, reads[i], use_ssw, k); + a_indv_max[i] = extend_seed( + aligner, + &mut chains[i][0], + refseq, + reads[i], + use_ssw, + k, + (0, reads[i].len()), + ); details[i].tried_alignment += 1; details[i].gapped += a_indv_max[i].as_ref().map_or(0, |a| a.gapped as usize); alignment_cache[i].insert(chains[i][0].id, a_indv_max[i].clone()); @@ -1056,7 +1086,15 @@ fn extend_paired_seeds( let alignment; if let Some(mut this_chain) = chainsp[i].clone() { if let Entry::Vacant(e) = alignment_cache[i].entry(this_chain.id) { - alignment = extend_seed(aligner, &mut this_chain, refseq, reads[i], use_ssw, k); + alignment = extend_seed( + aligner, + &mut this_chain, + refseq, + reads[i], + use_ssw, + k, + (0, reads[i].len()), + ); details[i].tried_alignment += 1; details[i].gapped += alignment.as_ref().map_or(0, |a| a.gapped as usize); e.insert(alignment.clone()); @@ -1154,7 +1192,9 @@ fn rescue_read( if score_dropoff1 < dropoff { break; } - if let Some(alignment) = extend_seed(aligner, chain, refseq, read1, use_ssw, k) { + if let Some(alignment) = + extend_seed(aligner, chain, refseq, read1, use_ssw, k, (0, read1.len())) + { details[0].add_alignment(&alignment); alignments1.push(alignment); diff --git a/src/supplementary.rs b/src/supplementary.rs index 0a9fb6bb0..4dc8e187d 100644 --- a/src/supplementary.rs +++ b/src/supplementary.rs @@ -39,12 +39,27 @@ pub fn get_supplementary_segment( let mut best_score = 0; let mut second_best_score = 0; + let primary_start = primary.soft_clip_left; + let primary_end = read.len() - primary.soft_clip_right; for chain in candidates.iter_mut().take(parameters.max_tries) { if chain.score < dropoff { break; } - let Some(alignment) = extend_seed(aligner, chain, refseq, read, parameters.use_ssw, k) - else { + + let interval = if chain.query_start >= primary_start { + (primary_end, read.len()) + } else { + (0, primary_start) + }; + let Some(alignment) = extend_seed( + aligner, + chain, + refseq, + read, + parameters.use_ssw, + k, + interval, + ) else { continue; }; details.supplementary_extended += 1; diff --git a/tests/phix.se.sam b/tests/phix.se.sam index 9cc645182..40aba8fcc 100644 --- a/tests/phix.se.sam +++ b/tests/phix.se.sam @@ -12,8 +12,8 @@ SRR1377138.8 16 NC_001422.1 1254 60 124=1X132=1X42=1X * 0 0 TCATGAAGGATGGTGTTAAT SRR1377138.9 16 NC_001422.1 4731 60 300=1X * 0 0 TGTGTGCCTGAGTATGGTACAGCTAATGGCCGTCTTCATTTCCATGCGGTGCACTTTATGCGGACACTTCCTACAGGTAGCGTTGACCCTAATTTTGGTCGTCGGGTACGCAATCGCCGCCAGTTAAATAGCTTGCAAAATACGTGGCCTTATGGTTACAGTATGCCCATCGCAGTTCGCTACACGCAGGACGCTTTTTCACGTTCTGGTTGGTTGTGGCCTGTTGATGCTAAAGGTGAGCCGCTTAAAGCTACCAGTTATATGGCTGTTGGTTTCTATGTGGCTAAATACGTTAACAAAN ?BFF?<:?*?*:<>9BDECDFGF<*9@F?*7B?;F;F;B:AFFFFFF7*(03:0:AAAFFF:18?BA:FFFFF2.1.4::?FFFFF?BAAFFFFFFDFF>>9>9? NM:i:1 AS:i:612 na:i:1 al:i:1 ga:i:0 X0:i:0 du:i:0 RG:Z:1 SRR1377138.14 0 NC_001422.1 66 60 1X300= * 0 0 NTCTTGATAAAGCAGGAATTACTACTGCTTGTTTACGAATTAAATCGAAGTGGACTGCTGGCGGAAAATGAGAAAATTCGACCTATCCTTGCGCAGCTCGAGAAGCTCTTACTTTGCGACCTTTCGCCATCAACTAACGATTCTGTCAAAAACTGACGCGTTGGATGAGGAGAAGTGGCTTAATATGCTTGGCACGTTCGTCAAGGACTGGTTTAGATATGAGTCACATTTTGTTCATGGTAGAGATTCTCTTGTTGACATTTTAAAAGAGCGTGGATTACTATCTGAGTCCGATGCTGTT #8@CCGGGGGGFFGGGGGGGGGGGGGFGGGGGGGGGGGGGFGGEEFGGGGGGGGGGGGGFGDCFEFGGGGGFGGGGGGGGGGGGGGGGGGGG7FGGGGGGGGGGGGGGGGGGGDFGGGGGGGFGGEGGGGGFGGGGGGGGGGGGGGGGGFGGGGGGGGGGDGGGGGGGGDGGGGGGGFEGGGGGGGFGGFFGFGGGGGGGGGFGGGFG7FGGFGDGGGGGFFCGGGGGGGGGGGGGFGFGGC>FGCF=A@FFFFFFFFFFEFFFFFFE<=/04(34:501:AB>4 NM:i:1 AS:i:612 na:i:1 al:i:1 ga:i:0 X0:i:0 du:i:0 RG:Z:1 SRR1377138.15 0 NC_001422.1 1907 60 1X185=115S * 0 0 NTTAACGCTAAAGGTGCTTTGACTTATACCGATATTGCTGGCGACCCTGTTTTGTATGGCAACTTGCCGCCGCGTGAAATTTCTATGAAGGATGTTTTCCGTTCTGGTGATTCGTCTAAGAAGTTTAAGATTGCTGAGGGTCAGTGGTATCGTTATGCGCCTTCGTATGTTTCTCCTGCTTATCACAAGAACGGAAAACATCCTTCATAGAAATTTCACGCAGATCGGAAGAGCGGTTCAGCAGGAATGCCGAGACCGATCTCGTATGCCGTCTTCTGCTTGAAAAAAAAAAAAAAAGATA #8BCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGEGFGGFGFGFFFFFFFFFFFFF>73?FFFFFFF3-048:BDFFAFFFFFF06>BFF######### NM:i:1 AS:i:372 na:i:1 al:i:1 ga:i:1 X0:i:0 du:i:0 RG:Z:1 @@ -31,8 +31,8 @@ SRR1377138.26 16 NC_001422.1 1642 60 300=1X * 0 0 GCGTTACCATGATGTTATTTCTTCATTTGG SRR1377138.27 16 NC_001422.1 4289 60 31=1X268=1X * 0 0 TGTGCTCCCCCAACTTGATATTAATAACACTTTAGACCACCGCCCCGAAGGGGACGAAAAATGGTTTTTAGAGAACGAGAAGACGGTTACGCAGTTTTGCCGCAAGCTGGCTGCTGAACGCCCTCTTAAGGATATTCGCGATGAGTATAATTACCCCAAAAAGAAAGGTATTAAGGATGAGTGTTCAAGATTGCTGGAGGCCTCCACTATGAAATCGCGTAGAGGCTTTGCTATTCAGCGTTTGATGAATGCAATGCGACAGGCTCATGCTGATGGTTGGTTTATCGTTTTTGACACTCTN 0FFF>>96BBAE9FFFF@FFFDFFFFGFGFGGDGFGGDGGGGGGGGGFGGGGEGGGGGGGGGGGGGGGGEGDFDGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFFGFGGGFGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGEGGGGGGGGGGFGGGGGGGGGGGGGGGGF?GGGGGGGGGGGGCCA8# NM:i:2 AS:i:602 na:i:1 al:i:1 ga:i:0 X0:i:0 du:i:0 RG:Z:1 SRR1377138.28 0 NC_001422.1 3016 60 1X116=1X183= * 0 0 NCTCCGGTGGCATTCAAGGTGATGTGCTTGCTACCGATAACAATACTGTAGGCATGGGTGATGCTGGTATTAAATCTGCCATTCAAGGCTCTAATGTTCCTAACCCTGATGAGGCCGTCCCTAGTTTTGTTTCTGGTGCTATGGCTAAAGCTGGTAAAGGACTTCTTGAAGGTACGTTGCAGGCTGGCACTTCTGCCGTTTCTGATAAGTTGCTTGATTTGGTTGGACTTGGTGGCAAGTCTGCCGCTGATAAAGGAAAGGATACTCGTGATTATCTTGCTGCTGCATTTCCTGAGCTTAA #8ACCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGCFGGGGGGGGGGGGGGFFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGEGGGGGGGGGGFGGGGGGGFGGFGGGGGDFGGFGGGGGFGFGGFGFGFFFFFFFFFFF:FEF@@F@=FEF@8/;?9))9CCEFFFFEFFEFF<7,.6FF?;96FFFB>??450.3:FFFFFB7,.6FF?;96FFFB>??450.3:FFFFFB7,.6FF?;96FFFB>??450.3:FFFFFB7,.6FF?;96FFFB>??450.3:FFFFFBBFBFFFAF?FF9B? NM:i:1 AS:i:612 na:i:1 al:i:1 ga:i:0 X0:i:0 du:i:0 RG:Z:1 SRR1377138.32 0 NC_001422.1 1433 60 1X268=1X31= * 0 0 NCCGAGGCTAACCCTAATGAGCTTAATCAAGATGATGCTCGTTATGGTTTCCGTTGCTGCCATCTCAAAAACATTTGGACTGCTCCGCTTCCTCCTGAGACTGAGCTTTCTCGCCAAATGACGACTTCTACCACATCTATTGACATTATGGGTCTGCAAGCTGCTTATGCTAATTTGCATACTGACCAAGAACGTGATTACTTCATGCAGCGTTACCATGATGTTATTTCTTCATTTGGAGGTAAAACCTCTTATGACGCTGACAACCGCCCTTTACTTGTCATGCGCTCTAATCTCTGGG #8BCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGFGGFGFGFFFFFFFFFFFFFFFFFFFFFFFF*.:6=;CFFFFFEFFEF9;(3:AFFB)66>?: NM:i:2 AS:i:602 na:i:1 al:i:1 ga:i:0 X0:i:0 du:i:0 RG:Z:1 SRR1377138.33 0 NC_001422.1 3817 60 1X268=1X31= * 0 0 NTATTTTTCATGGTATTGATAAAGCTGTTGCCGATACTTGGAACAATTTCTGGAAAGACGGTAAAGCTGATGGTATTGGCTCTAATTTGTCTAGGAAATAACCGTCAGGATTGACACCCTCCCAATTGTATGTTTTCATGCCTCCAAATCTTGGAGGCTTTTTTATGGTTCGTTCTTATTACCCTTCTGAATGTCACGCTGATTATTTTGACTTTGAGCGTATCGAGGCTCTTAAACCTGCTATTGAGGCTTGTGGCATTTCTACTCTTCCTCAATCCCCAATGCTTGGCTTCCATAAGCA #8BCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGGGGGGGGGGGGGGFEFGGGGGGGFGGFGGGGGGFFFGGFGFGFFFAFFFFFFBDFEFFEFFFFFFC,)217@EFFFDF<@EFFFF?FFFDE####### NM:i:2 AS:i:602 na:i:1 al:i:1 ga:i:0 X0:i:0 du:i:0 RG:Z:1 From b447c3c0b6038af1d18e1c1db35722441ea23c9e Mon Sep 17 00:00:00 2001 From: NicolasBuchin Date: Sat, 10 Oct 2026 19:14:35 +0200 Subject: [PATCH 05/15] Recover part of a chain sharing a parent with a better chain. 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 --- src/chainer.rs | 56 ++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 52 insertions(+), 4 deletions(-) diff --git a/src/chainer.rs b/src/chainer.rs index 509c53db7..4f164de25 100644 --- a/src/chainer.rs +++ b/src/chainer.rs @@ -527,10 +527,11 @@ impl ChainingResult { let mut used = vec![false; n]; let mut intervals = Vec::new(); - for (i, score) in candidates { + for (i, mut score) in candidates { if used[i] { continue; } + used[i] = true; let mut j = i; let mut overlaps = false; @@ -540,11 +541,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; @@ -555,7 +557,10 @@ impl ChainingResult { } if overlaps { - continue; + score = self.dp[i] - self.dp[j] + k as f32; + if score < valid_score { + continue; + } } let first = &self.anchors[j]; @@ -663,6 +668,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()); From f0dfb24d790cf6900d764f1952a988c845a35b2a Mon Sep 17 00:00:00 2001 From: NicolasBuchin Date: Fri, 30 Oct 2026 08:44:32 +0100 Subject: [PATCH 06/15] Supplementary candidates overlap tolerance with primary on the read Is-new-baseline: yes --- src/main.rs | 5 +++++ src/modes/align.rs | 9 +++++++-- src/supplementary.rs | 8 ++++---- 3 files changed, 16 insertions(+), 6 deletions(-) diff --git a/src/main.rs b/src/main.rs index 0e02150d4..a6637dc0a 100644 --- a/src/main.rs +++ b/src/main.rs @@ -321,6 +321,10 @@ struct Args { #[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, @@ -612,6 +616,7 @@ fn run() -> Result<(), CliError> { secondary_threshold: args.secondary_threshold, max_segment_gap: args.max_segment_gap, max_gap_distance: args.max_gap_distance, + max_insertion_ref_gap: args.max_insertion_ref_gap, min_supplementary_score: args.min_supplementary_score, max_supplementary_overlap: args.max_supplementary_overlap, max_tries: args.max_tries, diff --git a/src/modes/align.rs b/src/modes/align.rs index 43ca63b9a..35ff0a476 100644 --- a/src/modes/align.rs +++ b/src/modes/align.rs @@ -44,6 +44,7 @@ pub struct MappingParameters { pub secondary_threshold: f32, pub max_segment_gap: usize, pub max_gap_distance: usize, + pub max_insertion_ref_gap: usize, pub min_supplementary_score: u32, pub max_supplementary_overlap: f32, pub dropoff_threshold: f32, @@ -64,6 +65,7 @@ impl Default for MappingParameters { secondary_threshold: 0.95, max_segment_gap: 50, max_gap_distance: 100_000, + max_insertion_ref_gap: 200, min_supplementary_score: 80, max_supplementary_overlap: 0.5, dropoff_threshold: 0.5, @@ -695,7 +697,6 @@ pub fn extend_seed( read.seq() }; let (start, end) = query_interval; - debug_assert!(start <= chain.query_start && chain.query_end <= end,); let query = &full[start..end]; let chain_start = chain.query_start.saturating_sub(start); let chain_end = chain.query_end.saturating_sub(start).min(query.len()); @@ -759,7 +760,11 @@ pub fn extend_seed( .iter() .filter_map(|a| { let ref_start = a.ref_start() - decode_start - refseq.starts.0[contig_id]; - let query_start = a.query_start() - start; + // Remove anchors outside the interval + let query_start = a.query_start().checked_sub(start)?; + if query_start + k > query.len() { + return None; + } // Skip anchors resulting from hash collisions if query[query_start..query_start + k] == segment[ref_start..ref_start + k] { Some(Anchor::new(ref_start, query_start)) diff --git a/src/supplementary.rs b/src/supplementary.rs index 4dc8e187d..95a85b484 100644 --- a/src/supplementary.rs +++ b/src/supplementary.rs @@ -27,6 +27,7 @@ pub fn get_supplementary_segment( chain, read.len(), parameters.max_gap_distance as i64, + parameters.max_insertion_ref_gap as i64, parameters.max_supplementary_overlap, ) }) @@ -98,6 +99,7 @@ fn fits_variant_pattern( chain: &Chain, read_len: usize, max_gap_distance: i64, + max_insertion_ref_gap: i64, overlap_fraction: f32, ) -> bool { let (contig_id, contig_start) = refseq.unflatten(chain.ref_start); @@ -120,15 +122,13 @@ fn fits_variant_pattern( let primary_span = (read_len - primary.soft_clip_left - primary.soft_clip_right) as i64; let overlap_tolerance = (primary_span as f32 * overlap_fraction) as i64; - // A candidate has to bring read bases the primary does not have, which is also what - // keeps the chain the primary was built from out of its own search - if read_gap < 0 { + if read_gap < -overlap_tolerance { return false; } // An insertion: the contig carries straight on while the read jumps the // inserted bases - if (0..=overlap_tolerance).contains(&ref_gap) { + if (0..=max_insertion_ref_gap).contains(&ref_gap) { return true; } From 7e8fe50110232a753fb230d357155d96934ddace Mon Sep 17 00:00:00 2001 From: NicolasBuchin Date: Thu, 5 Nov 2026 15:55:03 +0100 Subject: [PATCH 07/15] Document supplementary alignments --- CHANGES.md | 5 +++++ README.md | 31 +++++++++++++++++++++++++++++++ 2 files changed, 36 insertions(+) diff --git a/CHANGES.md b/CHANGES.md index 3197565a4..0fb087e43 100644 --- a/CHANGES.md +++ b/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) diff --git a/README.md b/README.md index 099047170..7b2029e7f 100644 --- a/README.md +++ b/README.md @@ -281,6 +281,37 @@ 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 (20 by +default) apply per interval. A chain is kept if it scores at least `--vs` (40 +by default) or `--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 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, and duplication detection runs on it the same way it runs on the +primary. + +The supplementary alignment search runs for single-end reads only. + ### Duplication detection Strobealign attempts to reinterpret insertions inside an alignment as tandem From 5734cdad6039d191da2cbdb3bda181fb871ce9b0 Mon Sep 17 00:00:00 2001 From: NicolasBuchin Date: Mon, 9 Nov 2026 16:03:47 +0100 Subject: [PATCH 08/15] Add --no-supplementary to turn off the supplementary search --- README.md | 4 +++- src/main.rs | 5 +++++ src/modes/align.rs | 31 +++++++++++++++++++------------ 3 files changed, 27 insertions(+), 13 deletions(-) diff --git a/README.md b/README.md index 7b2029e7f..cef23fbbe 100644 --- a/README.md +++ b/README.md @@ -310,7 +310,9 @@ supplementary gets its own mapping quality, computed from the second best candidate, and duplication detection runs on it the same way it runs on the primary. -The supplementary alignment search runs for single-end reads only. +The supplementary alignment search runs for single-end reads only. Use +`--no-supplementary` to turn it off, which leaves the primary alignment +unchanged. ### Duplication detection diff --git a/src/main.rs b/src/main.rs index a6637dc0a..b341697c4 100644 --- a/src/main.rs +++ b/src/main.rs @@ -309,6 +309,10 @@ struct Args { // Supplementary parameters + /// 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, @@ -619,6 +623,7 @@ fn run() -> Result<(), CliError> { max_insertion_ref_gap: args.max_insertion_ref_gap, min_supplementary_score: args.min_supplementary_score, max_supplementary_overlap: args.max_supplementary_overlap, + supplementary: !args.no_supplementary, max_tries: args.max_tries, dropoff_threshold: args.dropoff_threshold, rescue_distance: args.rescue_distance, diff --git a/src/modes/align.rs b/src/modes/align.rs index 35ff0a476..6fc5250b4 100644 --- a/src/modes/align.rs +++ b/src/modes/align.rs @@ -47,6 +47,7 @@ pub struct MappingParameters { pub max_insertion_ref_gap: usize, pub min_supplementary_score: u32, pub max_supplementary_overlap: f32, + pub supplementary: bool, pub dropoff_threshold: f32, pub rescue_distance: usize, pub max_tries: usize, @@ -68,6 +69,7 @@ impl Default for MappingParameters { max_insertion_ref_gap: 200, min_supplementary_score: 80, max_supplementary_overlap: 0.5, + supplementary: true, dropoff_threshold: 0.5, rescue_distance: 100, max_tries: 20, @@ -561,18 +563,23 @@ pub fn align_single_end_read( }; let mapq = (60 * (best_score - second_best_score)).div_ceil(best_score) as u8; - let supplementary_timer = Instant::now(); - let supplementary = get_supplementary_segment( - aligner, - refseq, - &read, - index.k(), - &best_alignment, - chains, - mapping_parameters, - &mut details.align, - ); - details.align.time_extend_supplementary += supplementary_timer.elapsed().as_secs_f64(); + let supplementary = if mapping_parameters.supplementary { + let supplementary_timer = Instant::now(); + let supplementary = get_supplementary_segment( + aligner, + refseq, + &read, + index.k(), + &best_alignment, + chains, + mapping_parameters, + &mut details.align, + ); + details.align.time_extend_supplementary += supplementary_timer.elapsed().as_secs_f64(); + supplementary + } else { + None + }; let mut reported = vec![Segment { alignment: best_alignment, From c3235089197e78a3f6642883166950cf764d5133 Mon Sep 17 00:00:00 2001 From: NicolasBuchin Date: Wed, 11 Nov 2026 11:57:04 +0100 Subject: [PATCH 09/15] Keep extend_seed's signature and add extend_seed_in_interval --- src/modes/align.rs | 70 ++++++++++++++++---------------------------- src/supplementary.rs | 8 ++--- 2 files changed, 29 insertions(+), 49 deletions(-) diff --git a/src/modes/align.rs b/src/modes/align.rs index 6fc5250b4..f6b8d1f35 100644 --- a/src/modes/align.rs +++ b/src/modes/align.rs @@ -6,6 +6,7 @@ use std::cmp::{Reverse, min}; use std::collections::hash_map::Entry; use std::collections::{HashMap, HashSet}; use std::mem; +use std::ops::Range; use std::time::Instant; use fastrand::Rng; @@ -513,7 +514,6 @@ pub fn align_single_end_read( &read, mapping_parameters.use_ssw, index.k(), - (0, read.len()), ) else { continue; }; @@ -662,13 +662,13 @@ fn trace_log( let (mut ssw, mut pw) = if !mapping_parameters.use_ssw { ( - extend_seed(aligner, chain, refseq, read, true, k, (0, read.len())).unwrap(), + extend_seed(aligner, chain, refseq, read, true, k).unwrap(), alignment.clone(), ) } else { ( alignment.clone(), - extend_seed(aligner, chain, refseq, read, false, k, (0, read.len())).unwrap(), + extend_seed(aligner, chain, refseq, read, false, k).unwrap(), ) }; // manually adds the soft clips @@ -696,15 +696,29 @@ pub fn extend_seed( read: &Read, use_ssw: bool, k: usize, - query_interval: (usize, usize), +) -> Option { + extend_seed_in_interval(aligner, chain, refseq, read, use_ssw, k, 0..read.len()) +} + +/// Extend a chain over the given interval of the read (in the chain's +/// orientation) and return the resulting alignment. Read bases outside the +/// interval are soft clipped. +pub fn extend_seed_in_interval( + aligner: &Aligner, + chain: &mut Chain, + refseq: &RefSequence, + read: &Read, + use_ssw: bool, + k: usize, + query_interval: Range, ) -> Option { let full = if chain.is_revcomp { read.rc() } else { read.seq() }; - let (start, end) = query_interval; - let query = &full[start..end]; + let start = query_interval.start; + let query = &full[query_interval]; let chain_start = chain.query_start.saturating_sub(start); let chain_end = chain.query_end.saturating_sub(start).min(query.len()); let (contig_id, contig_start) = refseq.unflatten(chain.ref_start); @@ -1021,24 +1035,8 @@ fn extend_paired_seeds( let mut ch_max1 = chains[0][0].clone(); let mut ch_max2 = chains[1][0].clone(); - let alignment1 = extend_seed( - aligner, - &mut ch_max1, - refseq, - read1, - use_ssw, - k, - (0, read1.len()), - ); - let alignment2 = extend_seed( - aligner, - &mut ch_max2, - refseq, - read2, - use_ssw, - k, - (0, read2.len()), - ); + let alignment1 = extend_seed(aligner, &mut ch_max1, refseq, read1, use_ssw, k); + let alignment2 = extend_seed(aligner, &mut ch_max2, refseq, read2, use_ssw, k); if let (Some(alignment1), Some(alignment2)) = (alignment1, alignment2) { details[0].add_alignment(&alignment1); details[1].add_alignment(&alignment2); @@ -1062,15 +1060,7 @@ fn extend_paired_seeds( // the paired-end read as two single-end reads. let mut a_indv_max = [None, None]; for i in 0..2 { - a_indv_max[i] = extend_seed( - aligner, - &mut chains[i][0], - refseq, - reads[i], - use_ssw, - k, - (0, reads[i].len()), - ); + a_indv_max[i] = extend_seed(aligner, &mut chains[i][0], refseq, reads[i], use_ssw, k); details[i].tried_alignment += 1; details[i].gapped += a_indv_max[i].as_ref().map_or(0, |a| a.gapped as usize); alignment_cache[i].insert(chains[i][0].id, a_indv_max[i].clone()); @@ -1098,15 +1088,7 @@ fn extend_paired_seeds( let alignment; if let Some(mut this_chain) = chainsp[i].clone() { if let Entry::Vacant(e) = alignment_cache[i].entry(this_chain.id) { - alignment = extend_seed( - aligner, - &mut this_chain, - refseq, - reads[i], - use_ssw, - k, - (0, reads[i].len()), - ); + alignment = extend_seed(aligner, &mut this_chain, refseq, reads[i], use_ssw, k); details[i].tried_alignment += 1; details[i].gapped += alignment.as_ref().map_or(0, |a| a.gapped as usize); e.insert(alignment.clone()); @@ -1204,9 +1186,7 @@ fn rescue_read( if score_dropoff1 < dropoff { break; } - if let Some(alignment) = - extend_seed(aligner, chain, refseq, read1, use_ssw, k, (0, read1.len())) - { + if let Some(alignment) = extend_seed(aligner, chain, refseq, read1, use_ssw, k) { details[0].add_alignment(&alignment); alignments1.push(alignment); diff --git a/src/supplementary.rs b/src/supplementary.rs index 95a85b484..d67bd7c64 100644 --- a/src/supplementary.rs +++ b/src/supplementary.rs @@ -1,7 +1,7 @@ use crate::aligner::Aligner; use crate::chain::Chain; use crate::details::AlignmentDetails; -use crate::modes::align::{Alignment, MappingParameters, Segment, extend_seed}; +use crate::modes::align::{Alignment, MappingParameters, Segment, extend_seed_in_interval}; use crate::read::Read; use crate::refseq::RefSequence; @@ -48,11 +48,11 @@ pub fn get_supplementary_segment( } let interval = if chain.query_start >= primary_start { - (primary_end, read.len()) + primary_end..read.len() } else { - (0, primary_start) + 0..primary_start }; - let Some(alignment) = extend_seed( + let Some(alignment) = extend_seed_in_interval( aligner, chain, refseq, From 572e145aa0635598cff3cdc9dbe41e170b2a600c Mon Sep 17 00:00:00 2001 From: NicolasBuchin Date: Fri, 13 Nov 2026 07:54:15 +0100 Subject: [PATCH 10/15] Apply review suggestions Co-authored-by: Marcel Martin --- README.md | 22 ++++++++++++---------- src/duplications.rs | 2 +- src/supplementary.rs | 2 +- 3 files changed, 14 insertions(+), 12 deletions(-) diff --git a/README.md b/README.md index cef23fbbe..1a7b90d7e 100644 --- a/README.md +++ b/README.md @@ -289,14 +289,14 @@ the read the best chain misses is kept as well. ### Supplementary alignments -Strobealign searches for a 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 +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. @@ -307,8 +307,10 @@ 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, and duplication detection runs on it the same way it runs on the -primary. +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, which leaves the primary alignment diff --git a/src/duplications.rs b/src/duplications.rs index 7358791aa..614a41991 100644 --- a/src/duplications.rs +++ b/src/duplications.rs @@ -32,7 +32,7 @@ pub fn split_duplications( details, ); - split.sort_by_key(|segment| Reverse(segment.score)); + split.sort_by_key(|alignment| Reverse(alignment.score)); resolved.extend( split .into_iter() diff --git a/src/supplementary.rs b/src/supplementary.rs index d67bd7c64..9da56ef60 100644 --- a/src/supplementary.rs +++ b/src/supplementary.rs @@ -91,7 +91,7 @@ pub fn get_supplementary_segment( Some(Segment { alignment, mapq }) } -/// Whether a chain and the primary are the two sides of a deletion, insertion +/// Whether a chain and the primary alignment are the two sides of a deletion, insertion /// or duplication. fn fits_variant_pattern( primary: &Alignment, From 7e4d2b52833fd4c9e0219204a9fe1875db13a662 Mon Sep 17 00:00:00 2001 From: NicolasBuchin Date: Fri, 13 Nov 2026 08:05:10 +0100 Subject: [PATCH 11/15] Rename Segment to AlignmentWithMapq Use the existing alignment terminology for the alignment-with-mapq type and the supplementary search, and give the type an sa_entry() method. --- src/duplications.rs | 10 ++++----- src/main.rs | 6 ++--- src/modes/align.rs | 52 ++++++++++++++++++++++++-------------------- src/supplementary.rs | 16 ++++++++------ 4 files changed, 46 insertions(+), 38 deletions(-) diff --git a/src/duplications.rs b/src/duplications.rs index 614a41991..7869c8119 100644 --- a/src/duplications.rs +++ b/src/duplications.rs @@ -3,7 +3,7 @@ use std::cmp::Reverse; use crate::aligner::{Aligner, Scores}; use crate::cigar::{Cigar, CigarOperation, OpLen}; use crate::details::AlignmentDetails; -use crate::modes::align::{Alignment, MappingParameters, Segment}; +use crate::modes::align::{Alignment, AlignmentWithMapq, MappingParameters}; use crate::packed_seq::PackedSeqSlice; use crate::read::Read; use crate::refseq::RefSequence; @@ -15,12 +15,12 @@ pub fn split_duplications( aligner: &Aligner, refseq: &RefSequence, read: &Read, - segments: Vec, + alignments: Vec, parameters: &MappingParameters, details: &mut AlignmentDetails, -) -> Vec { +) -> Vec { let mut resolved = Vec::new(); - for Segment { alignment, mapq } in segments { + for AlignmentWithMapq { alignment, mapq } in alignments { let mut split = reinterpret_insertions( aligner, alignment, @@ -36,7 +36,7 @@ pub fn split_duplications( resolved.extend( split .into_iter() - .map(|alignment| Segment { alignment, mapq }), + .map(|alignment| AlignmentWithMapq { alignment, mapq }), ); } diff --git a/src/main.rs b/src/main.rs index b341697c4..3898fc188 100644 --- a/src/main.rs +++ b/src/main.rs @@ -318,8 +318,8 @@ struct Args { 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, + #[arg(long = "ss", default_value_t = MappingParameters::default().max_supplementary_gap, value_name = "N", help_heading = "Supplementary")] + max_supplementary_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")] @@ -618,7 +618,7 @@ fn run() -> Result<(), CliError> { args.max_secondary as usize }, secondary_threshold: args.secondary_threshold, - max_segment_gap: args.max_segment_gap, + max_supplementary_gap: args.max_supplementary_gap, max_gap_distance: args.max_gap_distance, max_insertion_ref_gap: args.max_insertion_ref_gap, min_supplementary_score: args.min_supplementary_score, diff --git a/src/modes/align.rs b/src/modes/align.rs index f6b8d1f35..938e2881e 100644 --- a/src/modes/align.rs +++ b/src/modes/align.rs @@ -35,7 +35,7 @@ use crate::read::Read; use crate::refseq::{ContigPosition, RefSequence}; use crate::revcomp::reverse_complement; use crate::seeding::SeedingParameters; -use crate::supplementary::get_supplementary_segment; +use crate::supplementary::get_supplementary_alignment; const MAX_PAIR_CHAINS: usize = 1000; @@ -43,7 +43,7 @@ const MAX_PAIR_CHAINS: usize = 1000; pub struct MappingParameters { pub max_secondary: usize, pub secondary_threshold: f32, - pub max_segment_gap: usize, + pub max_supplementary_gap: usize, pub max_gap_distance: usize, pub max_insertion_ref_gap: usize, pub min_supplementary_score: u32, @@ -65,7 +65,7 @@ impl Default for MappingParameters { MappingParameters { max_secondary: 0, secondary_threshold: 0.95, - max_segment_gap: 50, + max_supplementary_gap: 50, max_gap_distance: 100_000, max_insertion_ref_gap: 200, min_supplementary_score: 80, @@ -141,26 +141,32 @@ impl Alignment { } /// An alignment with its mapping quality. -pub struct Segment { +pub struct AlignmentWithMapq { pub alignment: Alignment, pub mapq: u8, } +impl AlignmentWithMapq { + fn sa_entry(&self, refseq: &RefSequence, cigar_eqx: bool) -> String { + self.alignment.sa_entry(refseq, self.mapq, cigar_eqx) + } +} + fn build_sa_tags( - segments: &[Segment], + alignments: &[AlignmentWithMapq], refseq: &RefSequence, cigar_eqx: bool, ) -> Vec> { - if segments.len() < 2 { - return vec![None; segments.len()]; + if alignments.len() < 2 { + return vec![None; alignments.len()]; } - let mut entries = Vec::with_capacity(segments.len()); - for segment in segments { - entries.push(segment.alignment.sa_entry(refseq, segment.mapq, cigar_eqx)); + let mut entries = Vec::with_capacity(alignments.len()); + for alignment in alignments { + entries.push(alignment.sa_entry(refseq, cigar_eqx)); } - let mut tags = Vec::with_capacity(segments.len()); + let mut tags = Vec::with_capacity(alignments.len()); for i in 0..entries.len() { let mut tag = String::new(); for (j, entry) in entries.iter().enumerate() { @@ -565,7 +571,7 @@ pub fn align_single_end_read( let supplementary = if mapping_parameters.supplementary { let supplementary_timer = Instant::now(); - let supplementary = get_supplementary_segment( + let supplementary = get_supplementary_alignment( aligner, refseq, &read, @@ -581,12 +587,12 @@ pub fn align_single_end_read( None }; - let mut reported = vec![Segment { + let mut reported = vec![AlignmentWithMapq { alignment: best_alignment, mapq, }]; reported.extend(supplementary); - let segments = split_duplications( + let reported = split_duplications( aligner, refseq, &read, @@ -594,19 +600,19 @@ pub fn align_single_end_read( mapping_parameters, &mut details.align, ); - let mut sa_tags = build_sa_tags(&segments, refseq, sam_output.cigar_eqx); + let mut sa_tags = build_sa_tags(&reported, refseq, sam_output.cigar_eqx); - for i in 0..segments.len() { + for i in 0..reported.len() { let status = if i == 0 { AlignmentStatus::Primary } else { AlignmentStatus::Supplementary }; sam_records.push(sam_output.make_mapped_record( - &segments[i].alignment, + &reported[i].alignment, refseq, record, - segments[i].mapq, + reported[i].mapq, status, details.clone(), sa_tags[i].take(), @@ -1693,21 +1699,21 @@ mod tests { collisions: 0, }; - let segments = [ - Segment { + let alignments = [ + AlignmentWithMapq { alignment: first.clone(), mapq: 60, }, - Segment { + AlignmentWithMapq { alignment: second, mapq: 60, }, ]; - let tags = build_sa_tags(&segments, &refseq, true); + let tags = build_sa_tags(&alignments, &refseq, true); assert_eq!(tags[0].as_deref(), Some("chr1,11,-,40S70=,60,0;")); assert_eq!(tags[1].as_deref(), Some("chr1,101,+,40=70S,60,0;")); - let only = [Segment { + let only = [AlignmentWithMapq { alignment: first, mapq: 60, }]; diff --git a/src/supplementary.rs b/src/supplementary.rs index 9da56ef60..dd96884c1 100644 --- a/src/supplementary.rs +++ b/src/supplementary.rs @@ -1,14 +1,16 @@ use crate::aligner::Aligner; use crate::chain::Chain; use crate::details::AlignmentDetails; -use crate::modes::align::{Alignment, MappingParameters, Segment, extend_seed_in_interval}; +use crate::modes::align::{ + Alignment, AlignmentWithMapq, MappingParameters, extend_seed_in_interval, +}; use crate::read::Read; use crate::refseq::RefSequence; /// Extend the chains that fit the pattern of a large variation with the /// primary and return the best one. Returns `None` if no chain fits or clears /// the score floor. -pub fn get_supplementary_segment( +pub fn get_supplementary_alignment( aligner: &Aligner, refseq: &RefSequence, read: &Read, @@ -17,7 +19,7 @@ pub fn get_supplementary_segment( chains: &mut [Chain], parameters: &MappingParameters, details: &mut AlignmentDetails, -) -> Option { +) -> Option { let mut candidates: Vec<&mut Chain> = chains .iter_mut() .filter(|chain| { @@ -35,7 +37,7 @@ pub fn get_supplementary_segment( details.supplementary_candidates += candidates.len(); let dropoff = candidates.first()?.score * parameters.dropoff_threshold; - let max_gap = parameters.max_segment_gap as i64; + let max_gap = parameters.max_supplementary_gap as i64; let mut best_alignment = None; let mut best_score = 0; let mut second_best_score = 0; @@ -66,7 +68,7 @@ pub fn get_supplementary_segment( details.supplementary_extended += 1; details.add_alignment(&alignment); - let (read_gap, ref_gap) = segment_gaps(primary, &alignment, read.len()); + let (read_gap, ref_gap) = alignment_gaps(primary, &alignment, read.len()); if read_gap.abs() > max_gap && ref_gap.abs() > max_gap { continue; } @@ -88,7 +90,7 @@ pub fn get_supplementary_segment( let mapq = (60 * best_score.saturating_sub(second_best_score)).div_ceil(best_score) as u8; details.supplementary += 1; - Some(Segment { alignment, mapq }) + Some(AlignmentWithMapq { alignment, mapq }) } /// Whether a chain and the primary alignment are the two sides of a deletion, insertion @@ -142,7 +144,7 @@ fn fits_variant_pattern( } /// How far apart two alignments end up on the read and on the reference. -fn segment_gaps(primary: &Alignment, candidate: &Alignment, read_len: usize) -> (i64, i64) { +fn alignment_gaps(primary: &Alignment, candidate: &Alignment, read_len: usize) -> (i64, i64) { if candidate.soft_clip_left >= primary.soft_clip_left { ( candidate.soft_clip_left as i64 - (read_len - primary.soft_clip_right) as i64, From 03531de5bc79e17b98fc91bc0f97a6fbb64e8c01 Mon Sep 17 00:00:00 2001 From: NicolasBuchin Date: Fri, 13 Nov 2026 08:19:54 +0100 Subject: [PATCH 12/15] Clarify documentation of the supplementary search and chain filtering --- README.md | 9 ++++----- src/chainer.rs | 16 ++++++++++++++++ src/duplications.rs | 11 +++++++++-- src/main.rs | 6 +++--- src/supplementary.rs | 7 ++++--- 5 files changed, 36 insertions(+), 13 deletions(-) diff --git a/README.md b/README.md index 1a7b90d7e..676104d70 100644 --- a/README.md +++ b/README.md @@ -283,9 +283,9 @@ 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 (20 by -default) apply per interval. A chain is kept if it scores at least `--vs` (40 -by default) or `--vp` of the best chain found, so a chain that covers a part of -the read the best chain misses is kept as well. +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 @@ -313,8 +313,7 @@ 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, which leaves the primary alignment -unchanged. +`--no-supplementary` to turn it off. ### Duplication detection diff --git a/src/chainer.rs b/src/chainer.rs index 4f164de25..ab2f02c68 100644 --- a/src/chainer.rs +++ b/src/chainer.rs @@ -456,6 +456,8 @@ fn hits_to_anchors(hits: &Vec, index: &StrobemerIndex) -> Vec { 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, @@ -463,6 +465,20 @@ struct QueryInterval { 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, start: usize, diff --git a/src/duplications.rs b/src/duplications.rs index 7869c8119..7a4369ca1 100644 --- a/src/duplications.rs +++ b/src/duplications.rs @@ -9,8 +9,13 @@ use crate::read::Read; use crate::refseq::RefSequence; use crate::simdaligner::SplitReferenceAlignment; -/// Each segment holding a duplication split into one per copy, the copies carrying the mapping -/// quality of the segment they came out of. +/// Splits each alignment that contains a tandem duplication into one alignment per copy. Each +/// copy keeps the mapping quality of the alignment it came from. +/// +/// The returned alignments keep the order of the input alignments: the copies of the first +/// input alignment come first, then those of the second, and so on. The copies of one input +/// alignment are sorted by score, best first. With the primary alignment passed first, the +/// first returned alignment is therefore the best copy of the primary. pub fn split_duplications( aligner: &Aligner, refseq: &RefSequence, @@ -32,6 +37,8 @@ pub fn split_duplications( details, ); + // Sort only within one split, so that the copies of the primary alignment stay ahead of + // those of the supplementary alignment split.sort_by_key(|alignment| Reverse(alignment.score)); resolved.extend( split diff --git a/src/main.rs b/src/main.rs index 3898fc188..dd5d03b66 100644 --- a/src/main.rs +++ b/src/main.rs @@ -249,7 +249,7 @@ struct Args { #[arg(long = "vp", default_value_t = ChainingParameters::default().valid_score_threshold, help_heading = "Collinear chaining")] valid_score_threshold: f32, - /// Score at which a chain is kept whatever the best chain of the read scores + /// Chains scoring at least this are kept even when below the --vp fraction of the best chain #[arg(long = "vs", default_value_t = ChainingParameters::default().min_chain_score, help_heading = "Collinear chaining")] min_chain_score: f32, @@ -1060,12 +1060,12 @@ fn print_statistics( details.align.time_extend * 100. / total ); debug!( - " primary chains: {:8.2} s {:5.1} %", + " primary extension: {:8.2} s {:5.1} %", details.align.time_extend_primary, details.align.time_extend_primary * 100. / total ); debug!( - " supplementary search: {:8.2} s {:5.1} %", + " supplementary extension: {:8.2} s {:5.1} %", details.align.time_extend_supplementary, details.align.time_extend_supplementary * 100. / total ); diff --git a/src/supplementary.rs b/src/supplementary.rs index dd96884c1..de96ff313 100644 --- a/src/supplementary.rs +++ b/src/supplementary.rs @@ -8,8 +8,8 @@ use crate::read::Read; use crate::refseq::RefSequence; /// Extend the chains that fit the pattern of a large variation with the -/// primary and return the best one. Returns `None` if no chain fits or clears -/// the score floor. +/// primary alignment and return the best one. Returns `None` if no chain fits, +/// or if the best extension scores below `min_supplementary_score`. pub fn get_supplementary_alignment( aligner: &Aligner, refseq: &RefSequence, @@ -143,7 +143,8 @@ fn fits_variant_pattern( false } -/// How far apart two alignments end up on the read and on the reference. +/// Returns the gaps between two alignments on the read and on the reference, +/// as (read_gap, ref_gap). A negative gap means the alignments overlap. fn alignment_gaps(primary: &Alignment, candidate: &Alignment, read_len: usize) -> (i64, i64) { if candidate.soft_clip_left >= primary.soft_clip_left { ( From abcf8ecebe3f2d8961fef899497aa76957579ab2 Mon Sep 17 00:00:00 2001 From: NicolasBuchin Date: Fri, 13 Nov 2026 08:20:04 +0100 Subject: [PATCH 13/15] Rename query_interval to read_interval in extend_seed_in_interval --- src/modes/align.rs | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/modes/align.rs b/src/modes/align.rs index 938e2881e..e659b5480 100644 --- a/src/modes/align.rs +++ b/src/modes/align.rs @@ -716,15 +716,15 @@ pub fn extend_seed_in_interval( read: &Read, use_ssw: bool, k: usize, - query_interval: Range, + read_interval: Range, ) -> Option { let full = if chain.is_revcomp { read.rc() } else { read.seq() }; - let start = query_interval.start; - let query = &full[query_interval]; + let start = read_interval.start; + let query = &full[read_interval]; let chain_start = chain.query_start.saturating_sub(start); let chain_end = chain.query_end.saturating_sub(start).min(query.len()); let (contig_id, contig_start) = refseq.unflatten(chain.ref_start); From 227643ab2b41964e2d4dee0e041ae5ae103c739d Mon Sep 17 00:00:00 2001 From: NicolasBuchin Date: Fri, 13 Nov 2026 12:19:07 +0100 Subject: [PATCH 14/15] Use alignment terminology in duplication splitting --- src/duplications.rs | 50 ++++++++++++++++++++++----------------------- 1 file changed, 25 insertions(+), 25 deletions(-) diff --git a/src/duplications.rs b/src/duplications.rs index 7a4369ca1..62cd31151 100644 --- a/src/duplications.rs +++ b/src/duplications.rs @@ -50,7 +50,7 @@ pub fn split_duplications( resolved } -/// Cuts the alignment into one segment per copy wherever an insertion of at least `min_length` +/// Cuts the alignment into one alignment per copy wherever an insertion of at least `min_length` /// bases turns out to be a tandem duplication. A `min_length` of zero turns this off. Always /// returns at least the alignment it was given. pub fn reinterpret_insertions( @@ -73,10 +73,10 @@ pub fn reinterpret_insertions( }; let contig = refseq.contig(alignment.contig_id); - let mut segments = Vec::new(); - let mut segment_cigar = Cigar::new(); - let mut segment_start = alignment.contig_start; - let mut segment_clip_left = alignment.soft_clip_left; + let mut alignments = Vec::new(); + let mut current_cigar = Cigar::new(); + let mut current_start = alignment.contig_start; + let mut current_clip_left = alignment.soft_clip_left; let mut ref_pos = alignment.contig_start; let mut query_pos = alignment.soft_clip_left; @@ -94,52 +94,52 @@ pub fn reinterpret_insertions( details, ) { - // The head arm carries the segment being built on past the insertion point, and + // The head arm carries the alignment being built on past the insertion point, and // the tail arm opens the next one at the second copy. - let mut upstream = segment_cigar.clone(); + let mut upstream = current_cigar.clone(); upstream.extend(&split.head.cigar); - if let Some(segment) = build_segment( + if let Some(head) = build_alignment( &alignment, upstream, - segment_start, - segment_clip_left, + current_start, + current_clip_left, query.len(), &aligner.scores, ) { - segments.push(segment); - segment_cigar = split.tail.cigar; - segment_start = copy_start; - segment_clip_left = query_pos + split.tail.query_start; + alignments.push(head); + current_cigar = split.tail.cigar; + current_start = copy_start; + current_clip_left = query_pos + split.tail.query_start; query_pos += len; continue; } } - segment_cigar.push(op, len); + current_cigar.push(op, len); let (query_len, ref_len) = query_ref_len(op, len); query_pos += query_len; ref_pos += ref_len; } - if segments.is_empty() { + if alignments.is_empty() { return vec![alignment]; } - let Some(last) = build_segment( + let Some(last) = build_alignment( &alignment, - segment_cigar, - segment_start, - segment_clip_left, + current_cigar, + current_start, + current_clip_left, query.len(), &aligner.scores, ) else { return vec![alignment]; }; - segments.push(last); - details.duplications += segments.len() - 1; + alignments.push(last); + details.duplications += alignments.len() - 1; - segments + alignments } /// Aligns the bases inserted at `ref_pos` against both references surrounding the @@ -179,10 +179,10 @@ fn find_duplication( Some((left_start + split.tail.ref_start, split)) } -/// Makes an alignment for a segment of a duplication split. It needs its own +/// Makes an alignment for one copy of a duplication split. It needs its own /// cigar, reference start, soft clips, score and edit distance, the rest comes /// from the original alignment. Returns None if it scores zero or less. -fn build_segment( +fn build_alignment( source: &Alignment, cigar: Cigar, ref_start: usize, From 6bc74de20f64f5dfb16a1f72d6e794e0dbdd7532 Mon Sep 17 00:00:00 2001 From: NicolasBuchin Date: Fri, 13 Nov 2026 12:21:47 +0100 Subject: [PATCH 15/15] Raise the default --vn to 1000 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 --- README.md | 2 +- src/chainer.rs | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 676104d70..679cd5f45 100644 --- a/README.md +++ b/README.md @@ -282,7 +282,7 @@ controls the chaining look-back window heuristic, `--gd` sets the diagonal gap c 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 (20 by +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. diff --git a/src/chainer.rs b/src/chainer.rs index ab2f02c68..05fc4d787 100644 --- a/src/chainer.rs +++ b/src/chainer.rs @@ -69,7 +69,7 @@ impl Default for ChainingParameters { gap_length_penalty: 0.05, valid_score_threshold: 0.7, min_chain_score: 40.0, - max_chains_per_interval: 20, + max_chains_per_interval: 1000, max_ref_gap: None, matches_weight: 0.01, max_diagonal_ratio: 10.0,