Duplications detection and supplementary records

CHANGES.md

@@ -6,6 +6,13 @@
 * #661: Fix a bug in index construction that, for some references, made the
   randstrobe with the smallest hash value impossible to find.
   Pre-generated index files (`.sti` files) need to be re-generated.
+* #656: Report tandem duplications in single-end mode: an insertion of at least
+  `--dl` bases (default 50) is re-aligned against the reference on both sides of
+  the insertion point, and if it matches at least a fraction `--dm` (default
+  0.80) of the inserted bases, the read is output as one record per copy instead
+  of as an insertion. `--dw` (default 20) is the bandwidth used for the two arms
+  and `--dl 0` turns detection off. Strobealign supports multiple records per
+  read with an `SA:Z` tag.
 
 ## v0.19.0 (2026-10-08)
 

CONTRIBUTING.md

@@ -102,6 +102,7 @@ mapped reads.
 `X0`: Number of equally scored best alignments (greater than 1 for multimappers).
  For paired-end reads, the tag is output for both reads, but the value is
  identical and is the number of equally scored best alignment *pairs*.
+`du`: Number of insertions re-interpreted as tandem duplications.
 
 ## Hidden command-line options
 

src/aligner.rs

@@ -5,6 +5,7 @@ use std::cell::Cell;
 use crate::chainer::Anchor;
 use crate::cigar::{Cigar, CigarOperation};
 use crate::piecewisealigner::PiecewiseAligner;
+use crate::simdaligner::SplitReferenceAlignment;
 use crate::ssw::SswAligner;
 
 #[derive(Debug, Clone, Copy)]
@@ -29,6 +30,22 @@ impl Default for Scores {
     }
 }
 
+impl Scores {
+    pub fn score(&self, cigar: &Cigar) -> i32 {
+        cigar
+            .iter()
+            .map(|(op, len)| match op {
+                CigarOperation::Eq => self.match_ as i32 * len as i32,
+                CigarOperation::X => -(self.mismatch as i32) * len as i32,
+                CigarOperation::Insertion | CigarOperation::Deletion => {
+                    -(self.gap_open as i32 + (len as i32 - 1) * self.gap_extend as i32)
+                }
+                _ => unreachable!("scoring a CIGAR that is not =/X/I/D: {len}{op}"),
+            })
+            .sum()
+    }
+}
+
 #[derive(Debug)]
 pub struct AlignmentInfo {
     pub cigar: Cigar,
@@ -174,6 +191,21 @@ impl Aligner {
             .expect("this Aligner was built for SSW extension only")
             .extend_piecewise(query, refseq, chain, padding)
     }
+
+    pub fn realign_insertion(
+        &self,
+        inserted: &[u8],
+        left_reference: &[u8],
+        right_reference: &[u8],
+        band: usize,
+    ) -> Option<SplitReferenceAlignment> {
+        Some(self.piecewise_aligner.as_ref()?.realign_insertion(
+            inserted,
+            left_reference,
+            right_reference,
+            band,
+        ))
+    }
 }
 
 pub fn hamming_distance(s: &[u8], t: &[u8]) -> Option<u32> {
@@ -308,6 +340,21 @@ mod test {
     use crate::aligner::{
         Aligner, Scores, hamming_align, hamming_align_global, highest_scoring_segment,
     };
+    use crate::cigar::Cigar;
+
+    #[test]
+    fn score_cigar() {
+        let scores = Scores {
+            match_: 2,
+            mismatch: 8,
+            gap_open: 12,
+            gap_extend: 1,
+            end_bonus: 10,
+        };
+        let cigar: Cigar = "10X100=52D5=".parse().unwrap();
+
+        assert_eq!(scores.score(&cigar), 67);
+    }
 
     #[test]
     fn ssw_align_no_result() {

src/cigar.rs

@@ -155,6 +155,20 @@ impl Cigar {
         dist
     }
 
+    /// Bases matched exactly, zero for a CIGAR written with `M`.
+    pub fn matches(&self) -> usize {
+        self.ops
+            .iter()
+            .filter(|op_len| op_len.op == CigarOperation::Eq)
+            .map(|op_len| op_len.len)
+            .sum()
+    }
+
+    /// Iterate over the operations as `(operation, length)` pairs.
+    pub fn iter(&self) -> impl ExactSizeIterator<Item = (CigarOperation, usize)> {
+        self.ops.iter().map(|oplen| (oplen.op, oplen.len))
+    }
+
     /// Return a new Cigar that uses = and X operations instead of M
     pub fn with_eqx(&self, query: &[u8], refseq: &[u8]) -> Cigar {
         let mut cigar = Cigar::new();
@@ -270,6 +284,20 @@ mod test {
         assert!(!cigar.is_empty());
     }
 
+    #[test]
+    fn iter_yields_operations_and_lengths() {
+        let cigar = Cigar::from_str("10=2X3D5=").unwrap();
+        assert_eq!(
+            cigar.iter().collect::<Vec<_>>(),
+            vec![
+                (CigarOperation::Eq, 10),
+                (CigarOperation::X, 2),
+                (CigarOperation::Deletion, 3),
+                (CigarOperation::Eq, 5),
+            ]
+        );
+    }
+
     #[test]
     fn construct() {
         let cigar = Cigar::new();

src/details.rs

@@ -65,6 +65,9 @@ pub struct AlignmentDetails {
     /// No. of best alignments with same score
     pub best_alignments: usize,
 
+    /// No. of insertions re-interpreted as tandem duplications
+    pub duplications: usize,
+
     pub time_extend: f64,
 }
 
@@ -76,6 +79,7 @@ impl ops::AddAssign<AlignmentDetails> for AlignmentDetails {
         self.anchors += rhs.anchors;
         self.collisions += rhs.collisions;
         self.best_alignments += rhs.best_alignments;
+        self.duplications += rhs.duplications;
         self.time_extend += rhs.time_extend;
     }
 }

src/duplications.rs

@@ -0,0 +1,186 @@
+use crate::aligner::{Aligner, Scores};
+use crate::cigar::{Cigar, CigarOperation};
+use crate::modes::align::Alignment;
+use crate::packed_seq::PackedSeqSlice;
+use crate::read::Read;
+use crate::refseq::RefSequence;
+use crate::simdaligner::SplitReferenceAlignment;
+
+/// 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.
+pub fn reinterpret_insertions(
+    aligner: &Aligner,
+    alignment: Alignment,
+    refseq: &RefSequence,
+    read: &Read,
+    min_length: usize,
+    min_matches: f32,
+    band: usize,
+) -> Vec<Alignment> {
+    if min_length == 0 {
+        return vec![alignment];
+    }
+    let query = if alignment.is_revcomp {
+        read.rc()
+    } else {
+        read.seq()
+    };
+    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 ref_pos = alignment.contig_start;
+    let mut query_pos = alignment.soft_clip_left;
+
+    for (op, len) in alignment.cigar.iter() {
+        if op == CigarOperation::Insertion
+            && len >= min_length
+            && let Some((copy_start, split)) = find_duplication(
+                aligner,
+                &contig,
+                &query[query_pos..query_pos + len],
+                ref_pos,
+                min_length,
+                min_matches,
+                band,
+            )
+        {
+            // The right arm carries the segment being built on past the insertion point, and
+            // the left arm opens the next one at the second copy.
+            let mut upstream = segment_cigar.clone();
+            upstream.extend(&split.right.cigar);
+
+            if let Some(segment) = clip_segment(
+                &alignment,
+                upstream,
+                segment_start,
+                segment_clip_left,
+                query.len(),
+                &aligner.scores,
+            ) {
+                segments.push(segment);
+                segment_cigar = split.left.cigar;
+                segment_start = copy_start;
+                segment_clip_left = query_pos + split.left.query_start;
+                query_pos += len;
+                continue;
+            }
+        }
+
+        segment_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() {
+        return vec![alignment];
+    }
+
+    let Some(last) = clip_segment(
+        &alignment,
+        segment_cigar,
+        segment_start,
+        segment_clip_left,
+        query.len(),
+        &aligner.scores,
+    ) else {
+        return vec![alignment];
+    };
+    segments.push(last);
+
+    segments
+}
+
+/// Aligns the bases inserted at `ref_pos` against both references surrounding the
+/// insertion region and returns a duplication pattern: 2 alignments bridging the
+/// query region by jumping from one reference to the other. None when the
+/// duplication pattern is not good enough to report.
+fn find_duplication(
+    aligner: &Aligner,
+    contig: &PackedSeqSlice,
+    insertion: &[u8],
+    ref_pos: usize,
+    min_length: usize,
+    min_matches: f32,
+    band: usize,
+) -> Option<(usize, SplitReferenceAlignment)> {
+    let len = insertion.len();
+    let padding = len + band;
+    let left_start = ref_pos.saturating_sub(len + padding);
+    let right_end = (ref_pos + len + padding).min(contig.len());
+    let left_reference = contig.decode(left_start, ref_pos);
+    let right_reference = contig.decode(ref_pos, right_end);
+
+    let split = aligner.realign_insertion(insertion, &left_reference, &right_reference, band)?;
+
+    let matched = split.left.cigar.matches() + split.right.cigar.matches();
+    let copy_len =
+        (split.right.ref_end - split.right.ref_start) + (split.left.ref_end - split.left.ref_start);
+
+    // The duplication must match most of the inserted sequence.
+    if (matched as f32) < min_matches * len as f32 || copy_len < min_length {
+        return None;
+    }
+
+    Some((left_start + split.left.ref_start, split))
+}
+
+/// Builds the alignment of one segment from the cigar it covers, soft-clipping the read on
+/// both sides and scoring the segment on its own. None if it scores zero or less.
+fn clip_segment(
+    source: &Alignment,
+    cigar: Cigar,
+    ref_start: usize,
+    soft_clip_left: usize,
+    read_len: usize,
+    scores: &Scores,
+) -> Option<Alignment> {
+    let (query_span, ref_span) = query_ref_span(&cigar);
+    let soft_clip_right = read_len - (soft_clip_left + query_span);
+    let unclipped_ends = (soft_clip_left == 0) as u32 + (soft_clip_right == 0) as u32;
+    let score = scores.score(&cigar) + (unclipped_ends * scores.end_bonus) as i32;
+    if score <= 0 {
+        return None;
+    }
+
+    Some(Alignment {
+        contig_id: source.contig_id,
+        is_revcomp: source.is_revcomp,
+        contig_start: ref_start,
+        length: ref_span,
+        soft_clip_left,
+        soft_clip_right,
+        edit_distance: cigar.edit_distance(),
+        score: score as u32,
+        cigar,
+        gapped: source.gapped,
+        anchors: source.anchors,
+        collisions: source.collisions,
+    })
+}
+
+fn query_ref_span(cigar: &Cigar) -> (usize, usize) {
+    let mut query = 0;
+    let mut reference = 0;
+
+    for (op, len) in cigar.iter() {
+        let (q, r) = query_ref_len(op, len);
+        query += q;
+        reference += r;
+    }
+
+    (query, reference)
+}
+
+fn query_ref_len(op: CigarOperation, len: usize) -> (usize, usize) {
+    match op {
+        CigarOperation::Match | CigarOperation::Eq | CigarOperation::X => (len, len),
+        CigarOperation::Insertion | CigarOperation::Softclip => (len, 0),
+        CigarOperation::Deletion | CigarOperation::Skip => (0, len),
+        CigarOperation::Hardclip | CigarOperation::Pad => (0, 0),
+    }
+}

src/io/sam.rs

@@ -35,6 +35,8 @@ pub struct SamRecord {
     pub alignment_score: Option<u32>,
     pub details: Option<Details>,
     pub rg_id: Option<String>,
+    /// Lists the alignments for a read split in multiple records
+    pub sa_tag: Option<String>,
     /// Copied from FASTQ header
     pub extra: Option<String>,
 }
@@ -116,11 +118,12 @@ impl Display for SamRecord {
         if let Some(details) = &self.details {
             write!(
                 f,
-                "\tna:i:{}\tal:i:{}\tga:i:{}\tX0:i:{}",
+                "\tna:i:{}\tal:i:{}\tga:i:{}\tX0:i:{}\tdu:i:{}",
                 details.chain.n_chains,
                 details.align.tried_alignment,
                 details.align.gapped,
                 details.align.best_alignments,
+                details.align.duplications,
             )?;
             if self.flags & PAIRED != 0 {
                 write!(f, "\tmr:i:{}", details.align.mate_rescue)?;
@@ -129,6 +132,9 @@ impl Display for SamRecord {
         if let Some(rg_id) = &self.rg_id {
             write!(f, "\tRG:Z:{}", rg_id)?;
         }
+        if let Some(sa_tag) = &self.sa_tag {
+            write!(f, "\tSA:Z:{}", sa_tag)?;
+        }
         if let Some(extra) = &self.extra {
             write!(f, "\t{}", extra)?;
         }

src/lib.rs

@@ -3,6 +3,7 @@ pub mod chain;
 pub mod chainer;
 pub mod cigar;
 pub mod details;
+pub mod duplications;
 pub mod hit;
 pub mod index;
 pub mod indexer;