New pairing of chains for extension mode

src/chainer.rs

@@ -392,7 +392,6 @@ fn extract_chains_from_dp(
         }
 
         let mut j = i;
-        let mut c = 1;
         let mut overlaps = false;
         let mut chain_anchors = vec![anchors[i]];
 
@@ -407,7 +406,6 @@ fn extract_chains_from_dp(
             }
             chain_anchors.push(anchors[j]);
             used[j] = true;
-            c += 1;
 
             matching_bases += ref_coverage.saturating_sub(anchors[j].ref_start).min(k);
             ref_coverage = anchors[j].ref_start;
@@ -426,10 +424,9 @@ fn extract_chains_from_dp(
             query_end: last.query_start + k,
             ref_start: first.ref_start,
             ref_end: last.ref_start + k,
-            n_matches: c,
             matching_bases,
             ref_id: last.ref_id,
-            score: score + c as f32 * chaining_parameters.matches_weight,
+            score: score + chain_anchors.len() as f32 * chaining_parameters.matches_weight,
             is_revcomp,
             anchors: chain_anchors,
         });

src/details.rs

@@ -53,6 +53,9 @@ pub struct Details {
     /// No. of times rescue by local alignment was attempted
     pub mate_rescue: usize,
 
+    /// No. of best alignments found using mate rescue
+    pub best_rescued: usize,
+
     /// No. of computed alignments (get_alignment or rescue_mate)
     pub tried_alignment: usize,
 
@@ -70,6 +73,7 @@ impl ops::AddAssign<Details> for Details {
         self.nam += rhs.nam;
         self.inconsistent_nams += rhs.inconsistent_nams;
         self.mate_rescue += rhs.mate_rescue;
+        self.best_rescued += rhs.best_rescued;
         self.tried_alignment += rhs.tried_alignment;
         self.gapped += rhs.gapped;
         self.best_alignments += rhs.best_alignments;

src/lib.rs

@@ -11,9 +11,11 @@ pub mod mapper;
 pub mod math;
 pub mod mcsstrategy;
 pub mod nam;
+pub mod pairing;
 pub mod partition;
 pub mod piecewisealigner;
 pub mod read;
 pub mod revcomp;
 pub mod seeding;
+pub mod shuffle;
 pub mod ssw;

src/main.rs

@@ -800,6 +800,10 @@ fn run() -> Result<(), CliError> {
     debug!("Total mapping sites tried: {}", details.tried_alignment);
     debug!("Inconsistent NAM ends: {}", details.inconsistent_nams);
     debug!("Mates rescued by alignment: {}", details.mate_rescue);
+    debug!(
+        "Best alignments with rescued pairs: {}",
+        details.best_rescued
+    );
 
     info!("Total time mapping: {:.2} s", timer.elapsed().as_secs_f64());
     //info!("Total time reading read-file(s): {:.2} s", );

src/maponly.rs

@@ -1,16 +1,16 @@
-use fastrand::Rng;
-
 use crate::chainer::Chainer;
-use crate::details::{Details, NamDetails};
+use crate::details::Details;
 use crate::index::StrobemerIndex;
 use crate::insertsize::InsertSizeDistribution;
 use crate::io::fasta::RefSequence;
 use crate::io::paf::PafRecord;
 use crate::io::record::{End, SequenceRecord};
 use crate::mapper::mapping_quality;
-use crate::math::normal_pdf;
 use crate::mcsstrategy::McsStrategy;
 use crate::nam::{Nam, get_nams_by_chaining, sort_nams};
+use crate::pairing::{PairedNams, get_paired_nams};
+use crate::shuffle::shuffle_best;
+use fastrand::Rng;
 
 /// Map a single-end read to the reference and return PAF records
 ///
@@ -132,20 +132,21 @@ pub fn map_paired_end_read(
         return (vec![], nam_details1.into());
     }
 
-    let nam_pairs = get_nam_pairs(
+    let mut paired_nams = get_paired_nams(
         &mut nams1,
         &mut nams2,
         insert_size_distribution.mu,
         insert_size_distribution.sigma,
         &nam_details1,
         &nam_details2,
     );
+    shuffle_best(&mut paired_nams, |p| p.score, rng);
 
     nam_details1.time_sort_nams = sort_nams(&mut nams1, rng);
     nam_details2.time_sort_nams = sort_nams(&mut nams2, rng);
 
     let mut records = vec![];
-    match get_best_paired_mapping_location(&nam_pairs, &nams1, &nams2, insert_size_distribution) {
+    match get_best_paired_mapping_location(&paired_nams, &nams1, &nams2, insert_size_distribution) {
         MappedNams::Individual(best1, best2) => {
             if let Some(nam1) = best1 {
                 records.push(paf_record_from_nam(
@@ -215,19 +216,20 @@ pub fn abundances_paired_end_read(
         return;
     }
 
-    let nam_pairs = get_nam_pairs(
+    let mut paired_nams = get_paired_nams(
         &mut nams1,
         &mut nams2,
         insert_size_distribution.mu,
         insert_size_distribution.sigma,
         &nam_details1,
         &nam_details2,
     );
+    shuffle_best(&mut paired_nams, |p| p.score, rng);
 
     sort_nams(&mut nams1, rng);
     sort_nams(&mut nams2, rng);
 
-    match get_best_paired_mapping_location(&nam_pairs, &nams1, &nams2, insert_size_distribution) {
+    match get_best_paired_mapping_location(&paired_nams, &nams1, &nams2, insert_size_distribution) {
         MappedNams::Individual(_, _) => {
             for (nams, read_len) in [(&nams1, r1.sequence.len()), (&nams2, r2.sequence.len())] {
                 let n_best = nams
@@ -241,17 +243,17 @@ pub fn abundances_paired_end_read(
             }
         }
         MappedNams::Pair(_, _, joint_score) => {
-            let n_best = nam_pairs
+            let n_best = paired_nams
                 .iter()
                 .take_while(|nam_pair| nam_pair.score == joint_score)
                 .count();
             let weight_r1 = r1.sequence.len() as f64 / n_best as f64;
             let weight_r2 = r2.sequence.len() as f64 / n_best as f64;
-            for NamPair {
+            for PairedNams {
                 nam1,
                 nam2,
                 score: _,
-            } in &nam_pairs[..n_best]
+            } in &paired_nams[..n_best]
             {
                 abundances[nam1.ref_id] += weight_r1;
                 abundances[nam2.ref_id] += weight_r2;
@@ -268,14 +270,14 @@ enum MappedNams<'a> {
 }
 
 /// Choose between:
-/// - the best proper pair of mappings
-/// - the best individual mappings
+/// - the best individual NAMs
+/// - the best proper pair of NAMs
 ///
 /// Also updates the insert size distribution using confident pairs.
 ///
 /// For paired-end mapping and abundance estimation modes only
 fn get_best_paired_mapping_location<'a>(
-    nam_pairs: &'a [NamPair],
+    paired_nams: &'a [PairedNams],
     nams1: &'a [Nam],
     nams2: &'a [Nam],
     insert_size_distribution: &mut InsertSizeDistribution,
@@ -289,7 +291,7 @@ fn get_best_paired_mapping_location<'a>(
 
     // Prefer a proper pair only if it beats a penalized individual mapping.
     // Divisor 2 is penalty for being mapped individually
-    if let Some(NamPair { nam1, nam2, score }) = nam_pairs.first()
+    if let Some(PairedNams { nam1, nam2, score }) = paired_nams.first()
         && *score >= individual_score / 2.0
     {
         // Update insert size using confident proper pairs.
@@ -302,157 +304,3 @@ fn get_best_paired_mapping_location<'a>(
         MappedNams::Individual(best_nam1, best_nam2)
     }
 }
-
-#[derive(Debug)]
-pub struct NamPair {
-    pub nam1: Nam,
-    pub nam2: Nam,
-    pub score: f64,
-}
-
-/// Build all plausible forward/revcomp mapping pairings
-fn get_nam_pairs(
-    nams1: &mut [Nam],
-    nams2: &mut [Nam],
-    mu: f32,
-    sigma: f32,
-    details1: &NamDetails,
-    details2: &NamDetails,
-) -> Vec<NamPair> {
-    let mut nam_pairs = vec![];
-    if nams1.is_empty() || nams2.is_empty() {
-        return nam_pairs;
-    }
-
-    let (fwd1, rev1): (&mut [Nam], &mut [Nam]) =
-        split_nams_by_orientation_checked(nams1, details1.both_orientations);
-    let (fwd2, rev2): (&mut [Nam], &mut [Nam]) =
-        split_nams_by_orientation_checked(nams2, details2.both_orientations);
-
-    if !fwd1.is_empty() && !rev2.is_empty() {
-        fwd1.sort_unstable_by_key(|nam| (nam.ref_id, nam.projected_ref_start()));
-        rev2.sort_unstable_by_key(|nam| (nam.ref_id, nam.projected_ref_start()));
-        nam_pairs.extend(find_pairs(fwd1, rev2, mu, sigma, false));
-    }
-    if !fwd2.is_empty() && !rev1.is_empty() {
-        fwd2.sort_unstable_by_key(|nam| (nam.ref_id, nam.projected_ref_start()));
-        rev1.sort_unstable_by_key(|nam| (nam.ref_id, nam.projected_ref_start()));
-        nam_pairs.extend(find_pairs(fwd2, rev1, mu, sigma, true));
-    }
-
-    nam_pairs.sort_unstable_by(|a, b| b.score.total_cmp(&a.score));
-    nam_pairs
-}
-
-/// Split nams into (forward, revcomp),
-/// if only 1 orientation exists, returns it and a empty slice
-fn split_nams_by_orientation_checked(nams: &mut [Nam], both: bool) -> (&mut [Nam], &mut [Nam]) {
-    if both {
-        split_nams_by_orientation(nams)
-    } else if nams[0].is_revcomp {
-        (&mut [], nams)
-    } else {
-        (nams, &mut [])
-    }
-}
-
-/// In-place partition of NAMs by orientation:
-/// forward on the left, revcomp on the right.
-/// Returns two slices separating (forward, revcomp)
-fn split_nams_by_orientation(nams: &mut [Nam]) -> (&mut [Nam], &mut [Nam]) {
-    let mut left = 0;
-    let mut right = nams.len();
-
-    while left < right {
-        if nams[left].is_revcomp {
-            right -= 1;
-            nams.swap(left, right);
-        } else {
-            left += 1;
-        }
-    }
-
-    nams.split_at_mut(left)
-}
-
-/// Find most forward/revcomp pairs using a two-pointer scan.
-/// Assumes both slices are sorted by (ref_id, projected_ref_start).
-fn find_pairs(fwd: &[Nam], rev: &[Nam], mu: f32, sigma: f32, swap_order: bool) -> Vec<NamPair> {
-    let mut out = Vec::new();
-    let max_dist = (mu + 10.0 * sigma).ceil() as usize; // distance cutoff from insert size distribution
-    let mut rev_ptr = 0;
-    let mut last_paired = None;
-
-    for f in fwd {
-        // Advance revcomp pointer to the first possible candidate
-        while rev_ptr < rev.len()
-            && (rev[rev_ptr].ref_id < f.ref_id
-                || rev[rev_ptr].ref_id == f.ref_id
-                    && rev[rev_ptr].projected_ref_start() < f.projected_ref_start())
-        {
-            rev_ptr += 1;
-        }
-        if rev_ptr == rev.len() {
-            break;
-        }
-        if rev[rev_ptr].ref_id > f.ref_id {
-            continue;
-        }
-
-        // Scan window of revcomp nams within distance limit.
-        let mut best = None;
-        let mut i = rev_ptr;
-        while i < rev.len()
-            && rev[i].ref_id == f.ref_id
-            && (rev[i].projected_ref_start() - f.projected_ref_start()) <= max_dist
-        {
-            let r = &rev[i];
-            // The pairing score gets a bonus based on the reference distance of the two chosen nam
-            // paired and from our current knowledge of the reference distance distribution
-            let x = f.ref_start.abs_diff(r.ref_start);
-            let score = f.score as f64
-                + r.score as f64
-                + 0.001f64.max((normal_pdf(x as f32, mu, sigma) + 1.0).ln() as f64);
-
-            if best.is_none_or(|(_, highest_score)| score > highest_score) {
-                best = Some((i, score));
-            }
-            i += 1;
-        }
-
-        // Highest scoring candidate
-        let Some((best_id, score)) = best else {
-            continue;
-        };
-        let r = &rev[best_id];
-
-        // If the same revcomp nam was paired previously, keep only the better scoring pair.
-        if let Some((last_id, prev_score)) = last_paired
-            && last_id == best_id
-        {
-            if score <= prev_score {
-                continue; // keep the previous better pair
-            }
-            out.pop(); // replace it 
-        }
-
-        out.push(if swap_order {
-            NamPair {
-                nam1: r.clone(),
-                nam2: f.clone(),
-                score,
-            }
-        } else {
-            NamPair {
-                nam1: f.clone(),
-                nam2: r.clone(),
-                score,
-            }
-        });
-
-        last_paired = Some((best_id, score));
-        rev_ptr = best_id;
-    }
-
-    out
-}