From 196400276a0032a73050ffe886b2b51ef965948c Mon Sep 17 00:00:00 2001 From: Marcel Martin Date: Wed, 3 Jun 2026 16:49:06 +0200 Subject: [PATCH 1/4] Make fewer functions and structs public --- src/mapper.rs | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/src/mapper.rs b/src/mapper.rs index 84a6e342..4a93c8e6 100644 --- a/src/mapper.rs +++ b/src/mapper.rs @@ -1161,12 +1161,7 @@ pub fn is_proper_nam_pair(nam1: &Nam, nam2: &Nam, mu: f32, sigma: f32) -> bool { /// Find high-scoring NAMs and NAM pairs. Proper pairs are preferred, but also /// high-scoring NAMs that could not be paired up are returned (these are paired /// with None in the returned vector). -pub fn get_best_scoring_nam_pairs( - nams1: &[Nam], - nams2: &[Nam], - mu: f32, - sigma: f32, -) -> Vec { +fn get_best_scoring_nam_pairs(nams1: &[Nam], nams2: &[Nam], mu: f32, sigma: f32) -> Vec { let mut nam_pairs = vec![]; if nams1.is_empty() && nams2.is_empty() { return nam_pairs; @@ -1258,14 +1253,14 @@ pub fn mapping_quality(nams: &[Nam]) -> u8 { } #[derive(Debug)] -pub struct NamPair { +struct NamPair { pub score: f32, pub nam1: Option, pub nam2: Option, } #[derive(Debug, Clone)] -pub struct ScoredAlignmentPair { +struct ScoredAlignmentPair { score: f64, alignment1: Option, alignment2: Option, From 3f3f50e22f2b35b5a165181a31b7b7b9d2b51149 Mon Sep 17 00:00:00 2001 From: Marcel Martin Date: Wed, 3 Jun 2026 16:49:06 +0200 Subject: [PATCH 2/4] Skip secondary alignments in samdiff --- tests/samdiff.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/tests/samdiff.py b/tests/samdiff.py index 084d5490..dec0cf38 100755 --- a/tests/samdiff.py +++ b/tests/samdiff.py @@ -9,6 +9,7 @@ (RNAME, POS) only. Differences in the other SAM fields (flags, MAPQ, CIGAR, RNEXT, PNEXT, TLEN) are ignored. """ + import sys from argparse import ArgumentParser from contextlib import ExitStack @@ -35,6 +36,8 @@ def main(): with ExitStack() as stack: before = stack.enter_context(AlignmentFile(args.before)) after = stack.enter_context(AlignmentFile(args.after)) + before = skip_secondary(before) + after = skip_secondary(after) if has_truth: truth = stack.enter_context(AlignmentFile(args.truth)) else: @@ -164,6 +167,13 @@ def stat(description, value, should_be_zero: bool = True): sys.exit(1) +def skip_secondary(records): + for record in records: + if record.is_secondary: + continue + yield record + + def print_comparison(b: AlignedSegment, a: AlignedSegment): print(b.query_name) assert b.query_name == a.query_name, (b, a) From 7ca0599f905304745f2f03722a4e33afb1dcda77 Mon Sep 17 00:00:00 2001 From: Marcel Martin Date: Wed, 3 Jun 2026 16:49:06 +0200 Subject: [PATCH 3/4] Produce secondary alignments during baseline comparison These are ignored when comparing the outputs, but it allows us to test whether the BAM output is valid (according to picard ValidateSamFile). --- .github/workflows/rust.yml | 5 +++-- tests/compare-baseline.sh | 4 ++-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml index b62108ac..b13657fa 100644 --- a/.github/workflows/rust.yml +++ b/.github/workflows/rust.yml @@ -83,9 +83,10 @@ jobs: run: tests/compare-baseline.sh - name: Validate with Picard run: | - PicardCommandLine ValidateSamFile IGNORE=RECORD_MISSING_READ_GROUP IGNORE=MISSING_READ_GROUP I=head.bam + # IGNORE=QUALITY_NOT_STORED because Picard incorrectly warns about missing qualities for secondary alignments + PicardCommandLine ValidateSamFile IGNORE=RECORD_MISSING_READ_GROUP IGNORE=MISSING_READ_GROUP IGNORE=QUALITY_NOT_STORED I=head.bam - name: Compare to baseline (single-end) run: tests/compare-baseline.sh -s - name: Validate with Picard run: | - PicardCommandLine ValidateSamFile IGNORE=RECORD_MISSING_READ_GROUP IGNORE=MISSING_READ_GROUP I=head.bam + PicardCommandLine ValidateSamFile IGNORE=RECORD_MISSING_READ_GROUP IGNORE=MISSING_READ_GROUP IGNORE=QUALITY_NOT_STORED I=head.bam diff --git a/tests/compare-baseline.sh b/tests/compare-baseline.sh index 3879adb6..14844037 100755 --- a/tests/compare-baseline.sh +++ b/tests/compare-baseline.sh @@ -61,14 +61,14 @@ if ! test -f ${baseline_bam}; then mv ${srcdir}/target/release/strobealign ${baseline_binary} rm -rf "${srcdir}" fi - ${baseline_binary} -v -t ${threads} ${ref} ${reads[@]} | samtools view -o ${baseline_bam}.tmp.bam + ${baseline_binary} -N 2 -v -t ${threads} ${ref} ${reads[@]} | samtools view -o ${baseline_bam}.tmp.bam mv ${baseline_bam}.tmp.bam ${baseline_bam} fi # Build and run strobealign cargo build --release set -x -target/release/strobealign -v -t ${threads} ${ref} ${reads[@]} | samtools view -o head.bam +target/release/strobealign -N 2 -v -t ${threads} ${ref} ${reads[@]} | samtools view -o head.bam # Do the actual comparison tests/samdiff.py ${baseline_bam} head.bam From fdd5d6147073cf1108424b1ff2b03c7f3e7b95b7 Mon Sep 17 00:00:00 2001 From: Marcel Martin Date: Wed, 3 Jun 2026 16:49:06 +0200 Subject: [PATCH 4/4] Fix: Do not output unmapped secondary alignments Is-new-baseline: yes --- src/mapper.rs | 23 +++++++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) diff --git a/src/mapper.rs b/src/mapper.rs index 4a93c8e6..24412601 100644 --- a/src/mapper.rs +++ b/src/mapper.rs @@ -237,6 +237,11 @@ impl SamOutput { sam_records } + /// Create two SAM records, one for each end of a paired-end read. + /// At least one of the reads must have an `Alignment` + /// (use `make_unmapped_pair` otherwise). + /// No unmapped secondary alignments are returned (this is the only case + /// for which the returned vector can have one element). fn make_paired_records( &self, alignments: [Option<&Alignment>; 2], @@ -246,9 +251,10 @@ impl SamOutput { details: &[Details; 2], alignment_status: AlignmentStatus, pair_status: PairStatus, - ) -> [SamRecord; 2] { - // Create single-end records - let mut sam_records = [ + ) -> Vec { + assert!(alignments[0].is_some() || alignments[1].is_some()); + // Create single-end records ... + let mut sam_records = vec![ self.make_record( alignments[0], references, @@ -266,7 +272,7 @@ impl SamOutput { details[1].clone(), ), ]; - // Then make them paired + // ... then make them paired sam_records[0].flags |= READ1; sam_records[1].flags |= READ2; @@ -326,6 +332,15 @@ impl SamOutput { } } + // Avoid unmapped secondary alignments + if alignment_status == AlignmentStatus::Secondary { + if !sam_records[1].is_mapped() { + sam_records.pop(); + } else if !sam_records[0].is_mapped() { + sam_records.swap_remove(0); + } + } + sam_records } }