perbase_lib/position/
mate_fix.rs

1//! Check for mate fixing.
2//!
3//! The strategy here is as follows:
4//! - Always check first for the user-based filter status of each mate, if one fails, go with the other mate, if both fail, default to `a`
5//! - Try one of the [`MateResolutionStrategy`] variants
6//!
7//! ## MateResolutionStrategy::BaseQualMapQualFirstInPair
8//! If either is deletion / ref skip -> MAPQ -> first in pair
9//! Else BASEQ -> MAPQ -> first in pair
10//!
11//! ## MateResolutionStrategy::BaseQualMapQualIUPAC
12//! If either is deletion / ref skip -> MAPQ -> first in pair
13//! Else BASEQ -> MAPQ -> IAUPAC
14//!
15//! ## MateResolutionStrategy::BaseQualMapQualN
16//! If either is deletion / ref skip -> MAPQ -> first in pair
17//! Else BASEQ -> MAPQ -> N (if ambiguous)
18//!
19//! ## MateResolutionStrategy::MapQualBaseQualFirstInPair
20//! If either is deletion / ref skip -> MAPQ -> first in pair
21//! Else MAPQ -> BASEQ-> first in pair
22//!
23//! ## MateResolutionStrategy::MapQualBaseQualIUPAC
24//! If either is deletion / ref skip -> MAPQ -> first in pair
25//! Else MAPQ -> BASEQ -> IUPAC
26//!
27//! ## MateResolutionStrategy::MapQualBaseQualN
28//! If either is deletion / ref skip -> MAPQ -> first in pair
29//! Else MAPQ -> BASEQ -> N (if ambiguous)
30//!
31//! ## MateResolutionStrategy::IUPAC
32//! If either is deletion / ref skip -> MAPQ -> first in pair
33//! Else IUPAC
34//!
35//! ## MateResolutionStrategy::N
36//! If either is deletion / ref skip -> MAPQ -> first in pair
37//! Else N
38//!
39//! ## MateResolutionStrategy::Original
40//! MAPQ -> first in pair
41use crate::read_filter::ReadFilter;
42use rust_htslib::bam::{pileup::Alignment, record::Record};
43use strum::EnumString;
44
45use std::cmp::Ordering;
46
47/// IUPAC Bases
48#[derive(Copy, Clone, Debug)]
49pub enum Base {
50    /// A
51    A,
52    /// C
53    C,
54    /// T or U
55    T,
56    /// G
57    G,
58    /// N (any base)
59    N,
60    /// A or G
61    R,
62    /// C or T
63    Y,
64    /// G or C
65    S,
66    /// A or T
67    W,
68    /// G or T
69    K,
70    /// A or C
71    M,
72    /// C or G or T
73    B,
74    /// A or G or T
75    D,
76    /// A or C or T
77    H,
78    /// A or C or G
79    V,
80    //Gap, // Technically a gap but not needed for us
81}
82
83impl Base {
84    /// Maps uncertainty of bases
85    pub(crate) fn either_or<const DEFAULT_TO_N: bool>(a: Self, b: Self) -> Self {
86        match (a, b) {
87            (Self::A, Self::A) => Self::A,
88            (Self::C, Self::C) => Self::C,
89            (Self::G, Self::G) => Self::G,
90            (Self::T, Self::T) => Self::T,
91            (Self::A, Self::G) | (Self::G, Self::A) if !DEFAULT_TO_N => Self::R,
92            (Self::C, Self::T) | (Self::T, Self::C) if !DEFAULT_TO_N => Self::Y,
93            (Self::G, Self::C) | (Self::C, Self::G) if !DEFAULT_TO_N => Self::S,
94            (Self::A, Self::T) | (Self::T, Self::A) if !DEFAULT_TO_N => Self::W,
95            (Self::G, Self::T) | (Self::T, Self::G) if !DEFAULT_TO_N => Self::K,
96            (Self::A, Self::C) | (Self::C, Self::A) if !DEFAULT_TO_N => Self::M,
97            (_, _) => Self::N,
98        }
99    }
100}
101
102impl From<char> for Base {
103    /// Lossy conversion, we don't support reading IUPAC bases, only outputting them for now
104    ///
105    /// If we wanted to support them, then we get crazy combos of either_or
106    #[inline]
107    fn from(value: char) -> Self {
108        match value.to_ascii_uppercase() {
109            'A' => Self::A,
110            'C' => Self::C,
111            'T' | 'U' => Self::T,
112            'G' => Self::G,
113            _ => Self::N,
114        }
115    }
116}
117
118/// Result of mate resolution comparison containing ordering and optional base recommendation.
119///
120/// This struct represents the decision made when comparing two overlapping mate reads
121/// at the same genomic position. It contains both the ordering decision (which mate
122/// to prefer) and an optional base recommendation for ambiguous cases.
123#[derive(Debug, Copy, Clone)]
124pub struct MateResolution {
125    /// The ordering of mate a relative to mate b (Greater = choose a, Less = choose b)
126    pub(crate) ordering: Ordering,
127    /// Optional base recommendation for cases where both mates contribute to the final call
128    pub(crate) recommended_base: Option<Base>,
129}
130
131impl MateResolution {
132    pub(crate) fn new(ordering: Ordering, recommended_base: Option<Base>) -> Self {
133        Self {
134            ordering,
135            recommended_base,
136        }
137    }
138}
139
140/// Strategies for resolving which mate to use when overlapping mates disagree on a base call.
141///
142/// All strategies first check user-based read filters. If one mate fails filters, the other
143/// is chosen. If both fail, the first mate is chosen by default.
144///
145/// For reads that are deletions / ref skips or lack a base call, all strategies fall back to the Original
146/// strategy (MAPQ → first in pair).
147#[derive(Debug, Copy, Clone, EnumString)]
148#[strum(ascii_case_insensitive)]
149pub enum MateResolutionStrategy {
150    /// Base quality → MAPQ → first in pair
151    ///
152    /// Priority order:
153    /// 1. Higher base quality score wins
154    /// 2. If base qualities are equal, higher MAPQ wins  
155    /// 3. If both are equal, first mate in pair wins
156    BaseQualMapQualFirstInPair,
157
158    /// Base quality → MAPQ → IUPAC ambiguity code
159    ///
160    /// Priority order:
161    /// 1. Higher base quality score wins
162    /// 2. If base qualities are equal, higher MAPQ wins
163    /// 3. If both are equal, return IUPAC ambiguity code (e.g., A+G→R)
164    BaseQualMapQualIUPAC,
165
166    /// Base quality → MAPQ → N
167    ///
168    /// Priority order:
169    /// 1. Higher base quality score wins
170    /// 2. If base qualities are equal, higher MAPQ wins
171    /// 3. If both are equal, return N (unknown base), if ambiguous
172    BaseQualMapQualN,
173
174    /// MAPQ → base quality → first in pair
175    ///
176    /// Priority order:
177    /// 1. Higher MAPQ wins
178    /// 2. If MAPQ is equal, higher base quality wins
179    /// 3. If both are equal, first mate in pair wins
180    MapQualBaseQualFirstInPair,
181
182    /// MAPQ → base quality → IUPAC ambiguity code
183    ///
184    /// Priority order:
185    /// 1. Higher MAPQ wins
186    /// 2. If MAPQ is equal, higher base quality wins
187    /// 3. If both are equal, return IUPAC ambiguity code (e.g., A+G→R)
188    MapQualBaseQualIUPAC,
189
190    /// MAPQ → base quality → N
191    ///
192    /// Priority order:
193    /// 1. Higher MAPQ wins
194    /// 2. If MAPQ is equal, higher base quality wins
195    /// 3. If both are equal, return N (unknown base), if ambiguous
196    MapQualBaseQualN,
197
198    /// Always return IUPAC ambiguity code for different bases
199    ///
200    /// Returns the appropriate IUPAC code for the base combination:
201    /// - Same bases: return the base (A+A→A)
202    /// - Different bases: return IUPAC code (A+G→R, C+T→Y, etc.)
203    IUPAC,
204
205    /// Always return N for different bases
206    ///
207    /// Returns N for different bases, but identical bases return themselves:
208    /// - Same bases: return the base (A+A→A)
209    /// - Different bases: return N
210    N,
211
212    /// Original strategy: MAPQ → first in pair
213    ///
214    /// Priority order:
215    /// 1. Higher MAPQ wins
216    /// 2. If MAPQ is equal, first mate in pair wins
217    /// 3. If neither is marked as first in pair, choose first mate by default
218    Original,
219}
220
221impl MateResolutionStrategy {
222    pub(crate) fn cmp<F: ReadFilter>(
223        &self,
224        a: &(Alignment<'_>, Record),
225        b: &(Alignment<'_>, Record),
226        read_filter: &F,
227    ) -> MateResolution {
228        // Handle user-set filters.
229        let a_pass = read_filter.filter_read(&a.1, Some(&a.0));
230        let b_pass = read_filter.filter_read(&b.1, Some(&b.0));
231        if a_pass && !b_pass {
232            return MateResolution::new(Ordering::Greater, None);
233        } else if b_pass && !a_pass {
234            return MateResolution::new(Ordering::Less, None);
235        } else if !a_pass && !b_pass {
236            return MateResolution::new(Ordering::Greater, None);
237        }
238
239        match self {
240            MateResolutionStrategy::BaseQualMapQualFirstInPair => {
241                Self::base_qual_map_qual_first_in_pair(a, b)
242            }
243            MateResolutionStrategy::BaseQualMapQualIUPAC => Self::base_qual_map_qual_iupac(a, b),
244            MateResolutionStrategy::BaseQualMapQualN => Self::base_qual_map_qual_n(a, b),
245            MateResolutionStrategy::MapQualBaseQualFirstInPair => {
246                Self::map_qual_base_qual_first_in_pair(a, b)
247            }
248            MateResolutionStrategy::MapQualBaseQualIUPAC => Self::map_qual_base_qual_iupac(a, b),
249            MateResolutionStrategy::MapQualBaseQualN => Self::map_qual_base_qual_n(a, b),
250            MateResolutionStrategy::IUPAC => Self::resolve_base::<false>(a, b),
251            MateResolutionStrategy::N => Self::resolve_base::<true>(a, b),
252            MateResolutionStrategy::Original => Self::original(a, b),
253        }
254    }
255
256    pub(crate) fn resolve_base<const DEFAULT_TO_N: bool>(
257        a: &(Alignment<'_>, Record),
258        b: &(Alignment<'_>, Record),
259    ) -> MateResolution {
260        // First check that we have a base
261        let a_is_not_base = a.0.qpos().is_none();
262        let b_is_not_base = b.0.qpos().is_none();
263        if a_is_not_base || b_is_not_base {
264            return Self::original(a, b);
265        }
266
267        let a_base = Base::from(a.1.seq()[a.0.qpos().unwrap()] as char);
268        let b_base = Base::from(b.1.seq()[b.0.qpos().unwrap()] as char);
269        MateResolution::new(
270            Ordering::Greater,
271            Some(Base::either_or::<DEFAULT_TO_N>(a_base, b_base)),
272        )
273    }
274
275    pub(crate) fn base_qual_map_qual_n(
276        a: &(Alignment<'_>, Record),
277        b: &(Alignment<'_>, Record),
278    ) -> MateResolution {
279        // First check that we have a base
280        let a_is_not_base = a.0.qpos().is_none();
281        let b_is_not_base = b.0.qpos().is_none();
282        if a_is_not_base || b_is_not_base {
283            return Self::original(a, b);
284        }
285
286        // compare baseq -> mapq -> default to first in pair
287        let a_qual = a.1.qual()[a.0.qpos().unwrap()];
288        let b_qual = b.1.qual()[b.0.qpos().unwrap()];
289
290        match a_qual.cmp(&b_qual) {
291            Ordering::Greater => MateResolution::new(Ordering::Greater, None),
292            Ordering::Less => MateResolution::new(Ordering::Less, None),
293            Ordering::Equal => match a.1.mapq().cmp(&b.1.mapq()) {
294                Ordering::Greater => MateResolution::new(Ordering::Greater, None),
295                Ordering::Less => MateResolution::new(Ordering::Less, None),
296                Ordering::Equal => {
297                    let a_base = Base::from(a.1.seq()[a.0.qpos().unwrap()] as char);
298                    let b_base = Base::from(b.1.seq()[b.0.qpos().unwrap()] as char);
299                    let base = Base::either_or::<true>(a_base, b_base);
300                    MateResolution::new(Ordering::Greater, Some(base)) // how to pass this back
301                }
302            },
303        }
304    }
305
306    pub(crate) fn base_qual_map_qual_iupac(
307        a: &(Alignment<'_>, Record),
308        b: &(Alignment<'_>, Record),
309    ) -> MateResolution {
310        // First check that we have a base
311        let a_is_not_base = a.0.qpos().is_none();
312        let b_is_not_base = b.0.qpos().is_none();
313        if a_is_not_base || b_is_not_base {
314            return Self::original(a, b);
315        }
316
317        // compare baseq -> mapq -> default to first in pair
318        let a_qual = a.1.qual()[a.0.qpos().unwrap()];
319        let b_qual = b.1.qual()[b.0.qpos().unwrap()];
320
321        match a_qual.cmp(&b_qual) {
322            Ordering::Greater => MateResolution::new(Ordering::Greater, None),
323            Ordering::Less => MateResolution::new(Ordering::Less, None),
324            Ordering::Equal => match a.1.mapq().cmp(&b.1.mapq()) {
325                Ordering::Greater => MateResolution::new(Ordering::Greater, None),
326                Ordering::Less => MateResolution::new(Ordering::Less, None),
327                Ordering::Equal => {
328                    let a_base = Base::from(a.1.seq()[a.0.qpos().unwrap()] as char);
329                    let b_base = Base::from(b.1.seq()[b.0.qpos().unwrap()] as char);
330                    MateResolution::new(
331                        Ordering::Greater,
332                        Some(Base::either_or::<false>(a_base, b_base)),
333                    )
334                }
335            },
336        }
337    }
338
339    pub(crate) fn base_qual_map_qual_first_in_pair(
340        a: &(Alignment<'_>, Record),
341        b: &(Alignment<'_>, Record),
342    ) -> MateResolution {
343        // First check that we have a base
344        let a_is_not_base = a.0.qpos().is_none();
345        let b_is_not_base = b.0.qpos().is_none();
346        if a_is_not_base || b_is_not_base {
347            return Self::original(a, b);
348        }
349
350        // compare baseq -> mapq -> default to first in pair
351        let a_qual = a.1.qual()[a.0.qpos().unwrap()];
352        let b_qual = b.1.qual()[b.0.qpos().unwrap()];
353
354        match a_qual.cmp(&b_qual) {
355            Ordering::Greater => MateResolution::new(Ordering::Greater, None),
356            Ordering::Less => MateResolution::new(Ordering::Less, None),
357            Ordering::Equal => Self::original(a, b),
358        }
359    }
360
361    pub(crate) fn map_qual_base_qual_n(
362        a: &(Alignment<'_>, Record),
363        b: &(Alignment<'_>, Record),
364    ) -> MateResolution {
365        // First check that we have a base
366        let a_is_not_base = a.0.qpos().is_none();
367        let b_is_not_base = b.0.qpos().is_none();
368        if a_is_not_base || b_is_not_base {
369            return Self::original(a, b);
370        }
371
372        match a.1.mapq().cmp(&b.1.mapq()) {
373            Ordering::Greater => MateResolution::new(Ordering::Greater, None),
374            Ordering::Less => MateResolution::new(Ordering::Less, None),
375            Ordering::Equal => {
376                let a_qual = a.1.qual()[a.0.qpos().unwrap()];
377                let b_qual = b.1.qual()[b.0.qpos().unwrap()];
378
379                match a_qual.cmp(&b_qual) {
380                    Ordering::Greater => MateResolution::new(Ordering::Greater, None),
381                    Ordering::Less => MateResolution::new(Ordering::Less, None),
382                    Ordering::Equal => {
383                        let a_base = Base::from(a.1.seq()[a.0.qpos().unwrap()] as char);
384                        let b_base = Base::from(b.1.seq()[b.0.qpos().unwrap()] as char);
385                        let base = Base::either_or::<true>(a_base, b_base);
386                        MateResolution::new(Ordering::Greater, Some(base)) // how to pass this back
387                    }
388                }
389            }
390        }
391    }
392
393    pub(crate) fn map_qual_base_qual_iupac(
394        a: &(Alignment<'_>, Record),
395        b: &(Alignment<'_>, Record),
396    ) -> MateResolution {
397        // First check that we have a base
398        let a_is_not_base = a.0.qpos().is_none();
399        let b_is_not_base = b.0.qpos().is_none();
400        if a_is_not_base || b_is_not_base {
401            return Self::original(a, b);
402        }
403
404        match a.1.mapq().cmp(&b.1.mapq()) {
405            Ordering::Greater => MateResolution::new(Ordering::Greater, None),
406            Ordering::Less => MateResolution::new(Ordering::Less, None),
407            Ordering::Equal => {
408                let a_qual = a.1.qual()[a.0.qpos().unwrap()];
409                let b_qual = b.1.qual()[b.0.qpos().unwrap()];
410
411                match a_qual.cmp(&b_qual) {
412                    Ordering::Greater => MateResolution::new(Ordering::Greater, None),
413                    Ordering::Less => MateResolution::new(Ordering::Less, None),
414                    Ordering::Equal => {
415                        let a_base = Base::from(a.1.seq()[a.0.qpos().unwrap()] as char);
416                        let b_base = Base::from(b.1.seq()[b.0.qpos().unwrap()] as char);
417                        MateResolution::new(
418                            Ordering::Greater,
419                            Some(Base::either_or::<false>(a_base, b_base)),
420                        )
421                    }
422                }
423            }
424        }
425    }
426
427    pub(crate) fn map_qual_base_qual_first_in_pair(
428        a: &(Alignment<'_>, Record),
429        b: &(Alignment<'_>, Record),
430    ) -> MateResolution {
431        // First check that we have a base
432        let a_is_not_base = a.0.qpos().is_none();
433        let b_is_not_base = b.0.qpos().is_none();
434        if a_is_not_base || b_is_not_base {
435            return Self::original(a, b);
436        }
437
438        match a.1.mapq().cmp(&b.1.mapq()) {
439            Ordering::Greater => MateResolution::new(Ordering::Greater, None),
440            Ordering::Less => MateResolution::new(Ordering::Less, None),
441            Ordering::Equal => {
442                let a_qual = a.1.qual()[a.0.qpos().unwrap()];
443                let b_qual = b.1.qual()[b.0.qpos().unwrap()];
444
445                match a_qual.cmp(&b_qual) {
446                    Ordering::Greater => MateResolution::new(Ordering::Greater, None),
447                    Ordering::Less => MateResolution::new(Ordering::Less, None),
448                    Ordering::Equal => {
449                        if a.1.flags() & 64 != 0 {
450                            MateResolution::new(Ordering::Greater, None)
451                        } else if b.1.flags() & 64 != 0 {
452                            MateResolution::new(Ordering::Less, None)
453                        } else {
454                            // Default to `a` in the event that there is no first in pair for some reason
455                            MateResolution::new(Ordering::Greater, None)
456                        }
457                    }
458                }
459            }
460        }
461    }
462
463    /// Whichever has higher MAPQ, or if equal, whichever is first in pair
464    pub(crate) fn original(
465        a: &(Alignment<'_>, Record),
466        b: &(Alignment<'_>, Record),
467    ) -> MateResolution {
468        match a.1.mapq().cmp(&b.1.mapq()) {
469            Ordering::Greater => MateResolution::new(Ordering::Greater, None),
470            Ordering::Less => MateResolution::new(Ordering::Less, None),
471            Ordering::Equal => {
472                // Check if a is first in pair
473                if a.1.flags() & 64 != 0 {
474                    MateResolution::new(Ordering::Greater, None)
475                } else if b.1.flags() & 64 != 0 {
476                    MateResolution::new(Ordering::Less, None)
477                } else {
478                    // Default to `a` in the event that there is no first in pair for some reason
479                    MateResolution::new(Ordering::Greater, None)
480                }
481            }
482        }
483    }
484}
485
486#[cfg(test)]
487mod tests {
488    use super::*;
489    use rust_htslib::bam::{self, record::Record};
490
491    // Since we can't directly create Alignment objects, we need to work around this
492    // by testing the internal functions directly
493
494    // Base enum tests
495    #[test]
496    fn test_base_either_or_identical() {
497        assert!(matches!(
498            Base::either_or::<false>(Base::A, Base::A),
499            Base::A
500        ));
501        assert!(matches!(
502            Base::either_or::<false>(Base::C, Base::C),
503            Base::C
504        ));
505        assert!(matches!(
506            Base::either_or::<false>(Base::G, Base::G),
507            Base::G
508        ));
509        assert!(matches!(
510            Base::either_or::<false>(Base::T, Base::T),
511            Base::T
512        ));
513    }
514
515    #[test]
516    fn test_base_either_or_iupac() {
517        // Test IUPAC combinations
518        assert!(matches!(
519            Base::either_or::<false>(Base::A, Base::G),
520            Base::R
521        ));
522        assert!(matches!(
523            Base::either_or::<false>(Base::G, Base::A),
524            Base::R
525        ));
526
527        assert!(matches!(
528            Base::either_or::<false>(Base::C, Base::T),
529            Base::Y
530        ));
531        assert!(matches!(
532            Base::either_or::<false>(Base::T, Base::C),
533            Base::Y
534        ));
535
536        assert!(matches!(
537            Base::either_or::<false>(Base::G, Base::C),
538            Base::S
539        ));
540        assert!(matches!(
541            Base::either_or::<false>(Base::C, Base::G),
542            Base::S
543        ));
544
545        assert!(matches!(
546            Base::either_or::<false>(Base::A, Base::T),
547            Base::W
548        ));
549        assert!(matches!(
550            Base::either_or::<false>(Base::T, Base::A),
551            Base::W
552        ));
553
554        assert!(matches!(
555            Base::either_or::<false>(Base::G, Base::T),
556            Base::K
557        ));
558        assert!(matches!(
559            Base::either_or::<false>(Base::T, Base::G),
560            Base::K
561        ));
562
563        assert!(matches!(
564            Base::either_or::<false>(Base::A, Base::C),
565            Base::M
566        ));
567        assert!(matches!(
568            Base::either_or::<false>(Base::C, Base::A),
569            Base::M
570        ));
571    }
572
573    #[test]
574    fn test_base_either_or_default_to_n() {
575        // When DEFAULT_TO_N is true, all non-identical should return N
576        assert!(matches!(Base::either_or::<true>(Base::A, Base::G), Base::N));
577        assert!(matches!(Base::either_or::<true>(Base::C, Base::T), Base::N));
578        assert!(matches!(Base::either_or::<true>(Base::G, Base::C), Base::N));
579
580        // Identical bases should still return themselves
581        assert!(matches!(Base::either_or::<true>(Base::A, Base::A), Base::A));
582        assert!(matches!(Base::either_or::<true>(Base::C, Base::C), Base::C));
583    }
584
585    #[test]
586    fn test_base_from_char() {
587        // Uppercase
588        assert!(matches!(Base::from('A'), Base::A));
589        assert!(matches!(Base::from('C'), Base::C));
590        assert!(matches!(Base::from('T'), Base::T));
591        assert!(matches!(Base::from('G'), Base::G));
592
593        // Lowercase
594        assert!(matches!(Base::from('a'), Base::A));
595        assert!(matches!(Base::from('c'), Base::C));
596        assert!(matches!(Base::from('t'), Base::T));
597        assert!(matches!(Base::from('g'), Base::G));
598
599        // U -> T
600        assert!(matches!(Base::from('U'), Base::T));
601        assert!(matches!(Base::from('u'), Base::T));
602
603        // Unknown -> N
604        assert!(matches!(Base::from('X'), Base::N));
605        assert!(matches!(Base::from('?'), Base::N));
606        assert!(matches!(Base::from('1'), Base::N));
607    }
608
609    // Since we can't easily test the strategy functions that take Alignment parameters,
610    // we'll add integration tests with the actual BAM reading in pileup_position.rs
611
612    // Test that we can at least construct the types
613    #[test]
614    fn test_mate_resolution_construction() {
615        let res = MateResolution::new(Ordering::Greater, Some(Base::N));
616        assert_eq!(res.ordering, Ordering::Greater);
617        assert!(matches!(res.recommended_base, Some(Base::N)));
618    }
619
620    // Helper function to create test BAM data and test a strategy
621    #[allow(clippy::too_many_arguments)]
622    fn test_strategy_with_bam(
623        strategy: MateResolutionStrategy,
624        read1_seq: &[u8],
625        read1_qual: &[u8],
626        read1_mapq: u8,
627        read1_flags: u16,
628        read2_seq: &[u8],
629        read2_qual: &[u8],
630        read2_mapq: u8,
631        read2_flags: u16,
632    ) -> MateResolution {
633        use rust_htslib::bam::{HeaderView, IndexedReader, Read, Writer, index};
634        use tempfile::tempdir;
635
636        let tempdir = tempdir().unwrap();
637        let bam_path = tempdir.path().join("test.bam");
638
639        // Create header
640        let mut header = bam::header::Header::new();
641        let mut chr1 = bam::header::HeaderRecord::new(b"SQ");
642        chr1.push_tag(b"SN", &"chr1".to_owned());
643        chr1.push_tag(b"LN", &"100".to_owned());
644        header.push_record(&chr1);
645        let view = HeaderView::from_header(&header);
646
647        // Create overlapping mate pair
648        // Read 1: positions 10-19 (1-based: 11-20)
649        // Read 2: positions 15-24 (1-based: 16-25)
650        // They overlap at positions 15-19 (1-based: 16-20)
651        let records = vec![
652            Record::from_sam(
653                &view,
654                format!(
655                    "TESTPAIR\t{}\tchr1\t11\t{}\t10M\tchr1\t16\t30\t{}\t{}",
656                    read1_flags,
657                    read1_mapq,
658                    std::str::from_utf8(read1_seq).unwrap(),
659                    std::str::from_utf8(read1_qual).unwrap()
660                )
661                .as_bytes(),
662            )
663            .unwrap(),
664            Record::from_sam(
665                &view,
666                format!(
667                    "TESTPAIR\t{}\tchr1\t16\t{}\t10M\tchr1\t11\t30\t{}\t{}",
668                    read2_flags,
669                    read2_mapq,
670                    std::str::from_utf8(read2_seq).unwrap(),
671                    std::str::from_utf8(read2_qual).unwrap()
672                )
673                .as_bytes(),
674            )
675            .unwrap(),
676        ];
677
678        // Write BAM file
679        let mut writer = Writer::from_path(&bam_path, &header, bam::Format::Bam).unwrap();
680        for record in &records {
681            writer.write(record).unwrap();
682        }
683        drop(writer);
684
685        // Index the BAM file
686        index::build(&bam_path, None, index::Type::Bai, 1).unwrap();
687
688        // Read back and test strategy at overlapping position
689        let mut reader = IndexedReader::from_path(&bam_path).unwrap();
690        reader.fetch(("chr1", 15, 20)).unwrap(); // Fetch the overlapping region
691        let pileup_iter = reader.pileup();
692
693        let read_filter = crate::read_filter::DefaultReadFilter::new(0, 0, 0);
694
695        for pileup_result in pileup_iter {
696            let pileup = pileup_result.unwrap();
697            if pileup.pos() >= 15 && pileup.pos() < 20 {
698                // Check overlapping positions
699                let alns: Vec<_> = pileup
700                    .alignments()
701                    .map(|aln| {
702                        let rec = aln.record();
703                        (aln, rec)
704                    })
705                    .collect();
706
707                if alns.len() == 2 {
708                    return strategy.cmp(&alns[0], &alns[1], &read_filter);
709                }
710            }
711        }
712
713        panic!("Failed to find overlapping reads at test position");
714    }
715
716    #[test]
717    fn test_base_qual_map_qual_first_in_pair() {
718        // Test different base qualities - higher base quality should win
719        let result = test_strategy_with_bam(
720            MateResolutionStrategy::BaseQualMapQualFirstInPair,
721            b"AAAAAAAAAA",
722            b"##########",
723            30,
724            67, // First mate: A, qual 2, MAPQ 30, first in pair
725            b"CCCCCCCCCC",
726            b">>>>>>>>>>",
727            30,
728            147, // Second mate: C, qual 30, MAPQ 30, second in pair
729        );
730        assert_eq!(result.ordering, Ordering::Less); // Second mate wins due to higher base quality
731        assert!(result.recommended_base.is_none());
732
733        // Test equal base qualities, different MAPQ - higher MAPQ should win
734        let result = test_strategy_with_bam(
735            MateResolutionStrategy::BaseQualMapQualFirstInPair,
736            b"AAAAAAAAAA",
737            b">>>>>>>>>>",
738            40,
739            67, // First mate: A, qual 30, MAPQ 40, first in pair
740            b"CCCCCCCCCC",
741            b">>>>>>>>>>",
742            30,
743            147, // Second mate: C, qual 30, MAPQ 30, second in pair
744        );
745        assert_eq!(result.ordering, Ordering::Greater); // First mate wins due to higher MAPQ
746        assert!(result.recommended_base.is_none());
747
748        // Test equal base qualities and MAPQ - first in pair should win
749        let result = test_strategy_with_bam(
750            MateResolutionStrategy::BaseQualMapQualFirstInPair,
751            b"AAAAAAAAAA",
752            b">>>>>>>>>>",
753            30,
754            67, // First mate: A, qual 30, MAPQ 30, first in pair
755            b"CCCCCCCCCC",
756            b">>>>>>>>>>",
757            30,
758            147, // Second mate: C, qual 30, MAPQ 30, second in pair
759        );
760        assert_eq!(result.ordering, Ordering::Greater); // First mate wins due to being first in pair
761        assert!(result.recommended_base.is_none());
762    }
763
764    #[test]
765    fn test_base_qual_map_qual_iupac() {
766        // Test different base qualities - higher base quality should win
767        let result = test_strategy_with_bam(
768            MateResolutionStrategy::BaseQualMapQualIUPAC,
769            b"AAAAAAAAAA",
770            b"##########",
771            30,
772            67, // First mate: A, qual 2, MAPQ 30
773            b"CCCCCCCCCC",
774            b">>>>>>>>>>",
775            30,
776            147, // Second mate: C, qual 30, MAPQ 30
777        );
778        assert_eq!(result.ordering, Ordering::Less); // Second mate wins due to higher base quality
779        assert!(result.recommended_base.is_none());
780
781        // Test equal base qualities and MAPQ - should return IUPAC code
782        let result = test_strategy_with_bam(
783            MateResolutionStrategy::BaseQualMapQualIUPAC,
784            b"AAAAAAAAAA",
785            b">>>>>>>>>>",
786            30,
787            67, // First mate: A, qual 30, MAPQ 30
788            b"GGGGGGGGGG",
789            b">>>>>>>>>>",
790            30,
791            147, // Second mate: G, qual 30, MAPQ 30
792        );
793        assert_eq!(result.ordering, Ordering::Greater); // Chooses first by default
794        assert!(matches!(result.recommended_base, Some(Base::R))); // A + G = R
795    }
796
797    #[test]
798    fn test_base_qual_map_qual_n() {
799        // Test different base qualities - higher base quality should win
800        let result = test_strategy_with_bam(
801            MateResolutionStrategy::BaseQualMapQualN,
802            b"AAAAAAAAAA",
803            b"##########",
804            30,
805            67, // First mate: A, qual 2, MAPQ 30
806            b"CCCCCCCCCC",
807            b">>>>>>>>>>",
808            30,
809            147, // Second mate: C, qual 30, MAPQ 30
810        );
811        assert_eq!(result.ordering, Ordering::Less); // Second mate wins due to higher base quality
812        assert!(result.recommended_base.is_none());
813
814        // Test equal base qualities and MAPQ - should return N
815        let result = test_strategy_with_bam(
816            MateResolutionStrategy::BaseQualMapQualN,
817            b"AAAAAAAAAA",
818            b">>>>>>>>>>",
819            30,
820            67, // First mate: A, qual 30, MAPQ 30
821            b"CCCCCCCCCC",
822            b">>>>>>>>>>",
823            30,
824            147, // Second mate: C, qual 30, MAPQ 30
825        );
826        assert_eq!(result.ordering, Ordering::Greater); // Chooses first by default
827        assert!(matches!(result.recommended_base, Some(Base::N))); // Returns N
828    }
829
830    #[test]
831    fn test_map_qual_base_qual_first_in_pair() {
832        // Test different MAPQ - higher MAPQ should win
833        let result = test_strategy_with_bam(
834            MateResolutionStrategy::MapQualBaseQualFirstInPair,
835            b"AAAAAAAAAA",
836            b">>>>>>>>>>",
837            40,
838            67, // First mate: A, qual 30, MAPQ 40
839            b"CCCCCCCCCC",
840            b">>>>>>>>>>",
841            30,
842            147, // Second mate: C, qual 30, MAPQ 30
843        );
844        assert_eq!(result.ordering, Ordering::Greater); // First mate wins due to higher MAPQ
845        assert!(result.recommended_base.is_none());
846
847        // Test equal MAPQ, different base qualities - higher base quality should win
848        let result = test_strategy_with_bam(
849            MateResolutionStrategy::MapQualBaseQualFirstInPair,
850            b"AAAAAAAAAA",
851            b"##########",
852            30,
853            67, // First mate: A, qual 2, MAPQ 30
854            b"CCCCCCCCCC",
855            b">>>>>>>>>>",
856            30,
857            147, // Second mate: C, qual 30, MAPQ 30
858        );
859        assert_eq!(result.ordering, Ordering::Less); // Second mate wins due to higher base quality
860        assert!(result.recommended_base.is_none());
861
862        // Test equal MAPQ and base qualities - first in pair should win
863        let result = test_strategy_with_bam(
864            MateResolutionStrategy::MapQualBaseQualFirstInPair,
865            b"AAAAAAAAAA",
866            b">>>>>>>>>>",
867            30,
868            67, // First mate: A, qual 30, MAPQ 30, first in pair
869            b"CCCCCCCCCC",
870            b">>>>>>>>>>",
871            30,
872            147, // Second mate: C, qual 30, MAPQ 30, second in pair
873        );
874        assert_eq!(result.ordering, Ordering::Greater); // First mate wins due to being first in pair
875        assert!(result.recommended_base.is_none());
876    }
877
878    #[test]
879    fn test_map_qual_base_qual_iupac() {
880        // Test different MAPQ - higher MAPQ should win
881        let result = test_strategy_with_bam(
882            MateResolutionStrategy::MapQualBaseQualIUPAC,
883            b"AAAAAAAAAA",
884            b">>>>>>>>>>",
885            40,
886            67, // First mate: A, qual 30, MAPQ 40
887            b"CCCCCCCCCC",
888            b">>>>>>>>>>",
889            30,
890            147, // Second mate: C, qual 30, MAPQ 30
891        );
892        assert_eq!(result.ordering, Ordering::Greater); // First mate wins due to higher MAPQ
893        assert!(result.recommended_base.is_none());
894
895        // Test equal MAPQ and base qualities - should return IUPAC code
896        let result = test_strategy_with_bam(
897            MateResolutionStrategy::MapQualBaseQualIUPAC,
898            b"AAAAAAAAAA",
899            b">>>>>>>>>>",
900            30,
901            67, // First mate: A, qual 30, MAPQ 30
902            b"GGGGGGGGGG",
903            b">>>>>>>>>>",
904            30,
905            147, // Second mate: G, qual 30, MAPQ 30
906        );
907        assert_eq!(result.ordering, Ordering::Greater); // Chooses first by default
908        assert!(matches!(result.recommended_base, Some(Base::R))); // A + G = R
909    }
910
911    #[test]
912    fn test_map_qual_base_qual_n() {
913        // Test different MAPQ - higher MAPQ should win
914        let result = test_strategy_with_bam(
915            MateResolutionStrategy::MapQualBaseQualN,
916            b"AAAAAAAAAA",
917            b">>>>>>>>>>",
918            40,
919            67, // First mate: A, qual 30, MAPQ 40
920            b"CCCCCCCCCC",
921            b">>>>>>>>>>",
922            30,
923            147, // Second mate: C, qual 30, MAPQ 30
924        );
925        assert_eq!(result.ordering, Ordering::Greater); // First mate wins due to higher MAPQ
926        assert!(result.recommended_base.is_none());
927
928        // Test equal MAPQ and base qualities - should return N
929        let result = test_strategy_with_bam(
930            MateResolutionStrategy::MapQualBaseQualN,
931            b"AAAAAAAAAA",
932            b">>>>>>>>>>",
933            30,
934            67, // First mate: A, qual 30, MAPQ 30
935            b"CCCCCCCCCC",
936            b">>>>>>>>>>",
937            30,
938            147, // Second mate: C, qual 30, MAPQ 30
939        );
940        assert_eq!(result.ordering, Ordering::Greater); // Chooses first by default
941        assert!(matches!(result.recommended_base, Some(Base::N))); // Returns N
942    }
943
944    #[test]
945    fn test_iupac_strategy() {
946        // Test with different bases - should return IUPAC code
947        let result = test_strategy_with_bam(
948            MateResolutionStrategy::IUPAC,
949            b"AAAAAAAAAA",
950            b">>>>>>>>>>",
951            30,
952            67, // First mate: A
953            b"GGGGGGGGGG",
954            b">>>>>>>>>>",
955            40,
956            147, // Second mate: G
957        );
958        assert_eq!(result.ordering, Ordering::Greater); // Always chooses first read for ordering
959        assert!(matches!(result.recommended_base, Some(Base::R))); // A + G = R
960
961        // Test with same bases - should return that base
962        let result = test_strategy_with_bam(
963            MateResolutionStrategy::IUPAC,
964            b"AAAAAAAAAA",
965            b">>>>>>>>>>",
966            30,
967            67, // First mate: A
968            b"AAAAAAAAAA",
969            b">>>>>>>>>>",
970            40,
971            147, // Second mate: A
972        );
973        assert_eq!(result.ordering, Ordering::Greater); // Always chooses first read for ordering
974        assert!(matches!(result.recommended_base, Some(Base::A))); // A + A = A
975    }
976
977    #[test]
978    fn test_n_strategy() {
979        // Test with different bases - should always return N
980        let result = test_strategy_with_bam(
981            MateResolutionStrategy::N,
982            b"AAAAAAAAAA",
983            b">>>>>>>>>>",
984            30,
985            67, // First mate: A
986            b"CCCCCCCCCC",
987            b">>>>>>>>>>",
988            40,
989            147, // Second mate: C
990        );
991        assert_eq!(result.ordering, Ordering::Greater); // Always chooses first read for ordering
992        assert!(matches!(result.recommended_base, Some(Base::N))); // Always N
993
994        // Test with same bases - should return the actual base since they're identical
995        let result = test_strategy_with_bam(
996            MateResolutionStrategy::N,
997            b"AAAAAAAAAA",
998            b">>>>>>>>>>",
999            30,
1000            67, // First mate: A
1001            b"AAAAAAAAAA",
1002            b">>>>>>>>>>",
1003            40,
1004            147, // Second mate: A
1005        );
1006        assert_eq!(result.ordering, Ordering::Greater); // Always chooses first read for ordering
1007        assert!(matches!(result.recommended_base, Some(Base::A))); // Identical bases return themselves
1008    }
1009
1010    #[test]
1011    fn test_original_strategy() {
1012        // Test different MAPQ - higher MAPQ should win
1013        let result = test_strategy_with_bam(
1014            MateResolutionStrategy::Original,
1015            b"AAAAAAAAAA",
1016            b">>>>>>>>>>",
1017            40,
1018            67, // First mate: A, MAPQ 40, first in pair
1019            b"CCCCCCCCCC",
1020            b">>>>>>>>>>",
1021            30,
1022            147, // Second mate: C, MAPQ 30, second in pair
1023        );
1024        assert_eq!(result.ordering, Ordering::Greater); // First mate wins due to higher MAPQ
1025        assert!(result.recommended_base.is_none());
1026
1027        // Test equal MAPQ - first in pair should win
1028        let result = test_strategy_with_bam(
1029            MateResolutionStrategy::Original,
1030            b"AAAAAAAAAA",
1031            b">>>>>>>>>>",
1032            30,
1033            67, // First mate: A, MAPQ 30, first in pair
1034            b"CCCCCCCCCC",
1035            b">>>>>>>>>>",
1036            30,
1037            147, // Second mate: C, MAPQ 30, second in pair
1038        );
1039        assert_eq!(result.ordering, Ordering::Greater); // First mate wins due to being first in pair
1040        assert!(result.recommended_base.is_none());
1041
1042        // Test equal MAPQ, second read is first in pair
1043        let result = test_strategy_with_bam(
1044            MateResolutionStrategy::Original,
1045            b"AAAAAAAAAA",
1046            b">>>>>>>>>>",
1047            30,
1048            147, // First mate: A, MAPQ 30, second in pair
1049            b"CCCCCCCCCC",
1050            b">>>>>>>>>>",
1051            30,
1052            67, // Second mate: C, MAPQ 30, first in pair
1053        );
1054        assert_eq!(result.ordering, Ordering::Less); // Second mate wins due to being first in pair
1055        assert!(result.recommended_base.is_none());
1056    }
1057}
perbase_lib/position/mate_fix.rs

perbase_lib/position/
mate_fix.rs
