Fix unmapped secondary alignments

.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

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<SamRecord> {
+        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
     }
 }
@@ -1161,12 +1176,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<NamPair> {
+fn get_best_scoring_nam_pairs(nams1: &[Nam], nams2: &[Nam], mu: f32, sigma: f32) -> Vec<NamPair> {
     let mut nam_pairs = vec![];
     if nams1.is_empty() && nams2.is_empty() {
         return nam_pairs;
@@ -1258,14 +1268,14 @@ pub fn mapping_quality(nams: &[Nam]) -> u8 {
 }
 
 #[derive(Debug)]
-pub struct NamPair {
+struct NamPair {
     pub score: f32,
     pub nam1: Option<Nam>,
     pub nam2: Option<Nam>,
 }
 
 #[derive(Debug, Clone)]
-pub struct ScoredAlignmentPair {
+struct ScoredAlignmentPair {
     score: f64,
     alignment1: Option<Alignment>,
     alignment2: Option<Alignment>,

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

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)