fastars 0.1.0

//! Overlap-based error corrector implementation.
//!
//! This module implements the core overlap detection and base correction logic
//! for paired-end reads.

use super::{CorrectionConfig, CorrectionStats};

// ============================================================================
// Overlap Region
// ============================================================================

/// Represents an overlap region between R1 and R2.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub struct OverlapRegion {
    /// Start position in R1 where overlap begins.
    pub r1_start: usize,
    /// End position in R1 (exclusive) where overlap ends.
    pub r1_end: usize,
    /// Start position in R2 (reverse-complemented) where overlap begins.
    pub r2_start: usize,
    /// End position in R2 (reverse-complemented, exclusive) where overlap ends.
    pub r2_end: usize,
    /// Number of mismatches in the overlap region.
    pub mismatches: usize,
}

impl OverlapRegion {
    /// Get the length of the overlap region.
    #[inline]
    pub fn len(&self) -> usize {
        self.r1_end.saturating_sub(self.r1_start)
    }

    /// Check if the overlap is empty.
    #[inline]
    pub fn is_empty(&self) -> bool {
        self.r1_end <= self.r1_start
    }
}

// ============================================================================
// Complement Table
// ============================================================================

/// Lookup table for DNA base complement.
/// A <-> T, G <-> C, N -> N
const COMPLEMENT: [u8; 256] = {
    let mut table = [0u8; 256];
    let mut i = 0;
    while i < 256 {
        table[i] = i as u8;
        i += 1;
    }
    table[b'A' as usize] = b'T';
    table[b'T' as usize] = b'A';
    table[b'G' as usize] = b'C';
    table[b'C' as usize] = b'G';
    table[b'a' as usize] = b't';
    table[b't' as usize] = b'a';
    table[b'g' as usize] = b'c';
    table[b'c' as usize] = b'g';
    table[b'N' as usize] = b'N';
    table[b'n' as usize] = b'n';
    table
};

/// Get the complement of a DNA base.
#[inline]
pub(crate) fn complement_base(base: u8) -> u8 {
    COMPLEMENT[base as usize]
}

// ============================================================================
// Overlap Corrector
// ============================================================================

/// Overlap-based error corrector for paired-end reads.
///
/// This corrector finds the overlap region between the tail of R1 and the
/// head of the reverse-complemented R2, then corrects mismatches based on
/// quality scores.
#[derive(Debug, Clone)]
pub struct OverlapCorrector {
    config: CorrectionConfig,
}

/// Pre-allocated buffers for correction operations.
/// Reused across read pairs within a worker thread.
pub struct CorrectionBuffers {
    pub r1_seq: Vec<u8>,
    pub r1_qual: Vec<u8>,
    pub r2_seq: Vec<u8>,
    pub r2_qual: Vec<u8>,
    pub r2_seq_rc: Vec<u8>,
    pub r2_qual_rc: Vec<u8>,
}

impl Default for CorrectionBuffers {
    fn default() -> Self {
        Self::new()
    }
}

impl CorrectionBuffers {
    /// Create new buffers with capacity for typical Illumina reads (300bp).
    pub fn new() -> Self {
        Self {
            r1_seq: Vec::with_capacity(300),
            r1_qual: Vec::with_capacity(300),
            r2_seq: Vec::with_capacity(300),
            r2_qual: Vec::with_capacity(300),
            r2_seq_rc: Vec::with_capacity(300),
            r2_qual_rc: Vec::with_capacity(300),
        }
    }

    /// Clear all buffers while preserving capacity.
    #[inline]
    pub fn clear(&mut self) {
        self.r1_seq.clear();
        self.r1_qual.clear();
        self.r2_seq.clear();
        self.r2_qual.clear();
        self.r2_seq_rc.clear();
        self.r2_qual_rc.clear();
    }
}

impl OverlapCorrector {
    /// Create a new corrector with the given configuration.
    pub fn new(config: CorrectionConfig) -> Self {
        Self { config }
    }

    /// Compute the reverse complement of a sequence.
    ///
    /// Returns a new Vec<u8> with the reverse-complemented sequence.
    #[inline]
    pub fn reverse_complement(seq: &[u8]) -> Vec<u8> {
        seq.iter()
            .rev()
            .map(|&b| complement_base(b))
            .collect()
    }

    /// Reverse a quality string (for the reverse complement).
    #[inline]
    pub fn reverse_quality(qual: &[u8]) -> Vec<u8> {
        qual.iter().rev().copied().collect()
    }

    /// Find the overlap between R1 and R2.
    ///
    /// This method searches for the best alignment between the tail of R1
    /// and the head of R2 (reverse-complemented). The overlap must be at least
    /// `min_overlap` bases and satisfy BOTH:
    /// - At most `diff_limit` mismatches (absolute count)
    /// - At most `diff_percent_limit`% mismatches (percentage)
    ///
    /// # Arguments
    /// * `r1_seq` - R1 sequence
    /// * `r2_seq_rc` - R2 sequence (already reverse-complemented)
    ///
    /// # Returns
    /// `Some(OverlapRegion)` if a valid overlap is found, `None` otherwise.
    pub fn find_overlap(&self, r1_seq: &[u8], r2_seq_rc: &[u8]) -> Option<OverlapRegion> {
        let r1_len = r1_seq.len();
        let r2_len = r2_seq_rc.len();

        if r1_len < self.config.min_overlap || r2_len < self.config.min_overlap {
            return None;
        }

        let mut best_overlap: Option<OverlapRegion> = None;
        let mut best_score: i64 = i64::MIN;

        // Try different overlap lengths, from minimum to maximum possible
        let max_overlap = r1_len.min(r2_len);
        let min_overlap = self.config.min_overlap;

        for overlap_len in min_overlap..=max_overlap {
            // R1 tail aligns with R2 head
            let r1_start = r1_len - overlap_len;
            let r2_start = 0;

            let r1_region = &r1_seq[r1_start..];
            let r2_region = &r2_seq_rc[r2_start..overlap_len];

            // Count mismatches
            let mismatches = count_mismatches(r1_region, r2_region);

            // Calculate mismatch percentage
            let mismatch_percent = (mismatches as f64 / overlap_len as f64) * 100.0;

            // Check both absolute and percentage limits
            if mismatches <= self.config.diff_limit && mismatch_percent <= self.config.diff_percent_limit {
                // Score: prefer longer overlaps with fewer mismatches
                // Score = overlap_len - 5 * mismatches
                let score = overlap_len as i64 - 5 * mismatches as i64;

                if score > best_score {
                    best_score = score;
                    best_overlap = Some(OverlapRegion {
                        r1_start,
                        r1_end: r1_len,
                        r2_start,
                        r2_end: overlap_len,
                        mismatches,
                    });
                }
            }
        }

        best_overlap
    }

    /// Correct a read pair based on overlap.
    ///
    /// # Deprecated
    /// This method allocates new Vecs for each call. For hot paths, use
    /// [`correct_pair_into`] with a [`CorrectionBuffers`] instance for
    /// zero-allocation correction.
    ///
    /// This method:
    /// 1. Finds the overlap between R1 and R2
    /// 2. For each mismatched position, chooses the base with higher quality
    /// 3. Updates both reads with the corrected base and adjusted quality
    ///
    /// # Arguments
    /// * `r1_seq` - R1 sequence (mutable for in-place correction)
    /// * `r1_qual` - R1 quality scores (mutable for in-place correction)
    /// * `r2_seq` - R2 sequence (mutable for in-place correction)
    /// * `r2_qual` - R2 quality scores (mutable for in-place correction)
    ///
    /// # Returns
    /// `Some(...)` if correction was applied, `None` if correction is disabled or no overlap found.
    /// When `Some`, returns tuple of (corrected_r1_seq, corrected_r1_qual, corrected_r2_seq, corrected_r2_qual, stats).
    #[deprecated(since = "0.4.0", note = "Use correct_pair_into() with CorrectionBuffers for better performance")]
    #[allow(clippy::type_complexity)]
    pub fn correct_pair(
        &self,
        r1_seq: &[u8],
        r1_qual: &[u8],
        r2_seq: &[u8],
        r2_qual: &[u8],
    ) -> Option<(Vec<u8>, Vec<u8>, Vec<u8>, Vec<u8>, CorrectionStats)> {
        if !self.config.enabled {
            return None;
        }

        // Quick length check before any allocation
        let r1_len = r1_seq.len();
        let r2_len = r2_seq.len();
        if r1_len < self.config.min_overlap || r2_len < self.config.min_overlap {
            return None;
        }

        // NOW compute RC (only if we might have overlap)
        let r2_seq_rc = Self::reverse_complement(r2_seq);
        let r2_qual_rc = Self::reverse_quality(r2_qual);

        // Find overlap
        let overlap = self.find_overlap(r1_seq, &r2_seq_rc)?;

        let mut stats = CorrectionStats::new();
        stats.pairs_processed = 1;
        stats.pairs_with_overlap = 1;

        // Create mutable copies for correction
        let mut corrected_r1_seq = r1_seq.to_vec();
        let mut corrected_r1_qual = r1_qual.to_vec();
        let mut corrected_r2_seq = r2_seq.to_vec();
        let mut corrected_r2_qual = r2_qual.to_vec();

        // Apply corrections in the overlap region
        let overlap_len = overlap.len();
        for i in 0..overlap_len {
            let r1_pos = overlap.r1_start + i;
            let r2_pos_rc = overlap.r2_start + i;
            // Convert back to original R2 position (un-reverse)
            let r2_pos = r2_seq.len() - 1 - r2_pos_rc;

            let r1_base = r1_seq[r1_pos];
            let r2_base_rc = r2_seq_rc[r2_pos_rc];
            let r1_q = r1_qual[r1_pos];
            let r2_q_rc = r2_qual_rc[r2_pos_rc];

            if r1_base != r2_base_rc {
                // Bases differ - choose the one with higher quality
                if r1_q >= r2_q_rc {
                    // R1 has higher or equal quality - correct R2
                    // The corrected base for R2 is the complement of R1's base
                    // (since R2 is in opposite orientation)
                    corrected_r2_seq[r2_pos] = complement_base(r1_base);
                    // Adjust quality: use the higher quality, minus the difference
                    let new_qual = calculate_corrected_quality(r1_q, r2_q_rc);
                    corrected_r2_qual[r2_pos] = new_qual;
                    stats.bases_corrected_r2 += 1;
                } else {
                    // R2 has higher quality - correct R1
                    corrected_r1_seq[r1_pos] = r2_base_rc;
                    let new_qual = calculate_corrected_quality(r2_q_rc, r1_q);
                    corrected_r1_qual[r1_pos] = new_qual;
                    stats.bases_corrected_r1 += 1;
                }
                stats.bases_corrected += 1;
            }
        }

        if stats.bases_corrected > 0 {
            stats.pairs_corrected = 1;
        }

        Some((
            corrected_r1_seq,
            corrected_r1_qual,
            corrected_r2_seq,
            corrected_r2_qual,
            stats,
        ))
    }

    /// Apply corrections in-place to the output buffers.
    fn apply_corrections_in_place(
        &self,
        r1_seq_out: &mut [u8],
        r1_qual_out: &mut [u8],
        r2_seq_out: &mut [u8],
        r2_qual_out: &mut [u8],
        overlap: &OverlapRegion,
        r2_seq_rc: &[u8],
        r2_qual_rc: &[u8],
    ) -> CorrectionStats {
        let mut stats = CorrectionStats::new();
        stats.pairs_processed = 1;
        stats.pairs_with_overlap = 1;

        let r2_seq_len = r2_seq_out.len();
        let overlap_len = overlap.len();

        for i in 0..overlap_len {
            let r1_pos = overlap.r1_start + i;
            let r2_pos_rc = overlap.r2_start + i;
            let r2_pos = r2_seq_len - 1 - r2_pos_rc;

            let r1_base = r1_seq_out[r1_pos];
            let r2_base_rc = r2_seq_rc[r2_pos_rc];
            let r1_q = r1_qual_out[r1_pos];
            let r2_q_rc = r2_qual_rc[r2_pos_rc];

            if r1_base != r2_base_rc {
                if r1_q >= r2_q_rc {
                    r2_seq_out[r2_pos] = complement_base(r1_base);
                    let new_qual = calculate_corrected_quality(r1_q, r2_q_rc);
                    r2_qual_out[r2_pos] = new_qual;
                    stats.bases_corrected_r2 += 1;
                } else {
                    r1_seq_out[r1_pos] = r2_base_rc;
                    let new_qual = calculate_corrected_quality(r2_q_rc, r1_q);
                    r1_qual_out[r1_pos] = new_qual;
                    stats.bases_corrected_r1 += 1;
                }
                stats.bases_corrected += 1;
            }
        }

        if stats.bases_corrected > 0 {
            stats.pairs_corrected = 1;
        }

        stats
    }

    /// Correct pair using pre-allocated buffers (zero per-call allocation).
    ///
    /// Returns `Some(stats)` if correction was applied, `None` if no overlap.
    pub fn correct_pair_into(
        &self,
        r1_seq: &[u8],
        r1_qual: &[u8],
        r2_seq: &[u8],
        r2_qual: &[u8],
        buffers: &mut CorrectionBuffers,
    ) -> Option<CorrectionStats> {
        if !self.config.enabled {
            return None;
        }

        // Quick length check before any work
        if r1_seq.len() < self.config.min_overlap || r2_seq.len() < self.config.min_overlap {
            return None;
        }

        // Clear buffers (preserves capacity)
        buffers.clear();

        // Compute RC into pre-allocated buffer
        buffers.r2_seq_rc.extend(r2_seq.iter().rev().map(|&b| complement_base(b)));
        buffers.r2_qual_rc.extend(r2_qual.iter().rev().copied());

        // Find overlap
        let overlap = self.find_overlap(r1_seq, &buffers.r2_seq_rc)?;

        // Copy input to output buffers
        buffers.r1_seq.extend_from_slice(r1_seq);
        buffers.r1_qual.extend_from_slice(r1_qual);
        buffers.r2_seq.extend_from_slice(r2_seq);
        buffers.r2_qual.extend_from_slice(r2_qual);

        // Apply corrections directly to buffers
        let stats = self.apply_corrections_in_place(
            &mut buffers.r1_seq,
            &mut buffers.r1_qual,
            &mut buffers.r2_seq,
            &mut buffers.r2_qual,
            &overlap,
            &buffers.r2_seq_rc,
            &buffers.r2_qual_rc,
        );

        Some(stats)
    }
}

/// Count mismatches between two sequences of equal length.
#[inline]
fn count_mismatches(seq1: &[u8], seq2: &[u8]) -> usize {
    seq1.iter()
        .zip(seq2.iter())
        .filter(|(&a, &b)| a != b)
        .count()
}

/// Calculate the quality score for a corrected base.
///
/// The corrected quality is based on the difference between the two quality scores.
/// If quality scores are similar, we're less confident in the correction.
/// If they differ significantly, we're more confident.
///
/// Formula: new_qual = max_qual - min(33, abs(q1 - q2) * 2)
/// This ensures the quality is never worse than the lower of the two,
/// but reflects uncertainty when qualities are close.
#[inline]
fn calculate_corrected_quality(winner_q: u8, loser_q: u8) -> u8 {
    // Quality scores are Phred+33 encoded
    let winner_phred = winner_q.saturating_sub(33);
    let loser_phred = loser_q.saturating_sub(33);

    // Calculate the new Phred score based on the quality difference
    // If difference is large, we're confident -> use winner's quality
    // If difference is small, we're less confident -> reduce quality
    let diff = winner_phred.saturating_sub(loser_phred);

    // If qualities are very close (diff < 3), reduce quality
    // Otherwise, use the winner's quality
    let new_phred = if diff < 3 {
        // Average minus a penalty for uncertainty
        ((winner_phred as u16 + loser_phred as u16) / 2) as u8
    } else {
        // Use winner's quality, slightly reduced
        winner_phred.saturating_sub(2)
    };

    // Convert back to Phred+33, ensuring minimum of 33 (Q0)
    new_phred.saturating_add(33).max(33)
}

#[cfg(test)]
mod tests {
    use super::*;

    fn make_qual(scores: &[u8]) -> Vec<u8> {
        scores.iter().map(|&s| s + 33).collect()
    }

    #[test]
    fn test_complement_base() {
        assert_eq!(complement_base(b'A'), b'T');
        assert_eq!(complement_base(b'T'), b'A');
        assert_eq!(complement_base(b'G'), b'C');
        assert_eq!(complement_base(b'C'), b'G');
        assert_eq!(complement_base(b'N'), b'N');
    }

    #[test]
    fn test_reverse_complement() {
        let seq = b"ACGT";
        let rc = OverlapCorrector::reverse_complement(seq);
        assert_eq!(rc, b"ACGT"); // ACGT reversed is TGCA, complement is ACGT

        let seq2 = b"AAAAGGGG";
        let rc2 = OverlapCorrector::reverse_complement(seq2);
        assert_eq!(rc2, b"CCCCTTTT");
    }

    #[test]
    fn test_reverse_quality() {
        let qual = b"IIIHHH";
        let reversed = OverlapCorrector::reverse_quality(qual);
        assert_eq!(reversed, b"HHHIII");
    }

    #[test]
    fn test_count_mismatches() {
        assert_eq!(count_mismatches(b"ACGT", b"ACGT"), 0);
        assert_eq!(count_mismatches(b"ACGT", b"TCGT"), 1);
        assert_eq!(count_mismatches(b"ACGT", b"TGCA"), 4);
    }

    #[test]
    fn test_find_overlap_no_overlap() {
        let config = CorrectionConfig::new()
            .enabled()
            .with_min_overlap(10);
        let corrector = OverlapCorrector::new(config);

        // Sequences too short
        let r1 = b"ACGT";
        let r2_rc = b"TGCA";
        assert!(corrector.find_overlap(r1, r2_rc).is_none());
    }

    #[test]
    fn test_find_overlap_exact_match() {
        let config = CorrectionConfig::new()
            .enabled()
            .with_min_overlap(10)
            .with_diff_limit(5);
        let corrector = OverlapCorrector::new(config);

        // R1: AAAAAAAAAA CCCCCCCCCC
        // R2 (rc): CCCCCCCCCC GGGGGGGGGG
        // Overlap should be CCCCCCCCCC (10 bases)
        let r1 = b"AAAAAAAAAACCCCCCCCCC";
        let r2_rc = b"CCCCCCCCCCGGGGGGGGGG";

        let overlap = corrector.find_overlap(r1, r2_rc);
        assert!(overlap.is_some());
        let o = overlap.unwrap();
        assert_eq!(o.len(), 10);
        assert_eq!(o.mismatches, 0);
    }

    #[test]
    fn test_find_overlap_with_mismatches() {
        let config = CorrectionConfig::new()
            .enabled()
            .with_min_overlap(10)
            .with_diff_limit(2)
            .with_diff_percent_limit(25.0); // 2 mismatches in 10bp = 20%, need higher limit
        let corrector = OverlapCorrector::new(config);

        // R1: AAAAAAAAAA CCCCCCCCCC
        // R2 (rc): CCTCCCCTCC (2 mismatches in 10-base overlap)
        // R1 tail (10 bases): CCCCCCCCCC
        // R2 RC head (10 bases): CCTCCCCTCC
        // Mismatches at positions 2 and 7
        let r1 = b"AAAAAAAAAACCCCCCCCCC";
        let r2_rc = b"CCTCCCCTCCGGGGGGGGGG";

        let overlap = corrector.find_overlap(r1, r2_rc);
        assert!(overlap.is_some());
        let o = overlap.unwrap();
        assert_eq!(o.mismatches, 2);
    }

    #[test]
    fn test_find_overlap_too_many_mismatches() {
        let config = CorrectionConfig::new()
            .enabled()
            .with_min_overlap(10)
            .with_diff_limit(1);
        let corrector = OverlapCorrector::new(config);

        // 2 mismatches but limit is 1
        let r1 = b"AAAAAAAAAACCCCCCCCCC";
        let r2_rc = b"CCTCCCCTCCGGGGGGGGGG";

        let overlap = corrector.find_overlap(r1, r2_rc);
        assert!(overlap.is_none());
    }

    #[test]
    fn test_correct_pair_disabled() {
        let config = CorrectionConfig::new().disabled();
        let corrector = OverlapCorrector::new(config);

        let r1_seq = b"ACGT";
        let r1_qual = make_qual(&[30, 30, 30, 30]);
        let r2_seq = b"ACGT";
        let r2_qual = make_qual(&[30, 30, 30, 30]);

        let result = corrector.correct_pair(r1_seq, &r1_qual, r2_seq, &r2_qual);
        assert!(result.is_none());
    }

    #[test]
    fn test_correct_pair_with_correction() {
        let config = CorrectionConfig::new()
            .enabled()
            .with_min_overlap(4)
            .with_diff_limit(2)
            .with_diff_percent_limit(20.0); // Allow up to 20% mismatch for the test
        let corrector = OverlapCorrector::new(config);

        // R1: ACGTACGT (8 bases)
        // R2: TTTTACGT (8 bases)
        // R2 reversed: TGCATTTT
        // R2 RC: ACGTAAAA
        //
        // R1 tail (4 bases): ACGT
        // R2 RC head (4 bases): ACGT -> perfect match at the end!
        //
        // Actually for a better test, let's create an overlap with a mismatch:
        // R1: AAAACCGT (8 bases) - ends with CCGT
        // R2: ACGGTTTT (8 bases)
        // R2 reversed: TTTTGGCA
        // R2 RC: AAAACCGT -> Wait, need to think more carefully
        //
        // Let me construct a simpler test case:
        // R1: AAAAAAAAACGT (12 bases)
        // R2 original: ACGAAAAAAAA (11 bases) - this should RC to TTTTTTTTTCGT
        // No wait, let me think again...
        //
        // For overlap to work:
        // R1 tail overlaps with R2 head (after R2 is RC'd)
        //
        // R1 = AAAAAAAAACGT (ends with ACGT)
        // R2 = ACGTTTTTTTTT (starts with ACGT in original orientation)
        // R2 reversed = TTTTTTTTTGCA
        // R2 RC = AAAAAAAAATGC... no
        //
        // Let me be more careful:
        // R2 = ACGTTTTTTTTT
        // Reverse of R2 = TTTTTTTTTGCA
        // Complement of reversed = AAAAAAAAACGT
        // So R2 RC = AAAAAAAAACGT
        //
        // Now R1 = AAAACGT ends in ACGT
        // R2 RC = AAAAAAAAACGT starts with AAAAAAAAAC
        // The overlap region would be where R1 tail matches R2 RC head
        //
        // Simpler approach - make them basically identical with one mismatch:
        let _r1_seq = b"AAAAAAAAAACCCCCCCCCC";
        let _r1_qual = make_qual(&[30; 20]);

        // For R2, we need its RC to match R1's tail
        // R1 tail (10 bases): CCCCCCCCCC
        // So R2 RC head should be: CCCCCCCCCC
        // Which means R2 original head (reversed) = GGGGGGGGGG complement = CCCCCCCCCC reversed = CCCCCCCCCC
        // So R2 = CCCCCCCCCC AAAAAAAAAA ... wait no
        //
        // R2 = some sequence
        // R2 reversed = rev
        // R2 RC = complement(rev)
        //
        // We want R2 RC to start with CCCCCCCCCC
        // complement(rev) starts with CCCCCCCCCC
        // rev starts with GGGGGGGGGG
        // So R2 ends with GGGGGGGGGG
        //
        // R2 = AAAAAAAAAGGGGGGGGGG (but with one base different for testing)
        // Let's put a mismatch: AAAAAAAAAGGGGTGGGGG (T instead of G at position 14)
        //
        // R2 reversed = GGGGTGGGGAAAAAAAAA... wait let me be precise
        // R2 = "AAAAAAAAAGGGGTGGGGG" (20 bases, G->T at position 13 from 0)
        // R2 reversed = "GGGGTGGGGAAAAAAAAAA"
        // R2 RC = complement of "GGGGTGGGGAAAAAAAAAA" = "CCCCACCCCTTTTTTTTT"
        //
        // R1 tail = "CCCCCCCCCC" (last 10 bases)
        // R2 RC head = "CCCCACCCCT" (first 10 bases)
        // Mismatch at position 4 (C vs A) and position 9 (C vs T)
        //
        // This has 2 mismatches. Let me simplify to 1:

        // R2 = "AAAAAAAAAGGGGGGGGG" + "G" with one G changed to T at a position
        // Let me just construct it directly:

        // R2 original ends with GGGGGGGGGG
        // R2 reversed starts with GGGGGGGGGG
        // R2 RC = complement starts with CCCCCCCCCC
        //
        // To get one mismatch, change one G in R2 to T:
        // R2 = AAAAAAAAAGGGGGGGGGT (last G is T)
        // R2 reversed = TGGGGGGGGAAAAAAAAAA... no wait
        //
        // Let me trace through step by step:
        // R2 = "AAAAAAAAAAGGGGGGGGT" (20 bases, ends in T)
        // R2 reversed = "TGGGGGGGGAAAAAAAAAA"
        // R2 RC = complement of reversed:
        //   T -> A
        //   G -> C
        // So R2 RC = "ACCCCCCCCTTTTTTTTT" (starts with A, then 8 C's, then T's)
        //
        // R1 tail (10 bases) = CCCCCCCCCC
        // R2 RC head (10 bases) = ACCCCCCCCT
        // Mismatches: position 0 (C vs A), position 9 (C vs T) = 2 mismatches

        // Let me try yet again with a cleaner setup:
        // R1 = AAAAAAAAAACCCCCCCCCC
        // R2 should be such that RC(R2) starts with CCCCCCCCCC (or close to it)
        //
        // If R2 RC = CCCCCCCCCC...
        // Then reverse(R2) = GGGGGGGGGG...
        // So R2 ends with GGGGGGGGGG

        // R2 = TTTTTTTTTTGGGGGGGGGG (perfect match)
        // R2 reversed = GGGGGGGGGGTTTTTTTTTT
        // R2 RC = CCCCCCCCCCAAAAAAAAAA

        // R1 = AAAAAAAAAACCCCCCCCCC (20 bases)
        // R1 tail (10 bases) = CCCCCCCCCC
        // R2 RC head (10 bases) = CCCCCCCCCC
        // Perfect match!

        // Now introduce one mismatch with different qualities:
        // Change one C in R1 to T with low quality:
        let mut r1_seq_with_error = b"AAAAAAAAAACCCCCCCCCC".to_vec();
        r1_seq_with_error[15] = b'T'; // Change one C to T
        let mut r1_qual_with_error = vec![30 + 33; 20]; // Q30 everywhere
        r1_qual_with_error[15] = 10 + 33; // Low quality at error position

        let r2_seq = b"TTTTTTTTTTGGGGGGGGGG";
        let r2_qual = make_qual(&[30; 20]); // Q30 everywhere

        let result = corrector.correct_pair(
            &r1_seq_with_error,
            &r1_qual_with_error,
            r2_seq,
            &r2_qual,
        );

        // The T at position 15 in R1 should be corrected to C (higher quality from R2)
        assert!(result.is_some());
        let (cr1_seq, _cr1_qual, _cr2_seq, _cr2_qual, stats) = result.unwrap();
        assert_eq!(stats.pairs_with_overlap, 1);
        assert_eq!(stats.bases_corrected, 1);
        assert_eq!(stats.bases_corrected_r1, 1);
        assert_eq!(cr1_seq[15], b'C'); // Should be corrected from T to C
    }

    #[test]
    fn test_overlap_region() {
        let region = OverlapRegion {
            r1_start: 10,
            r1_end: 20,
            r2_start: 0,
            r2_end: 10,
            mismatches: 2,
        };

        assert_eq!(region.len(), 10);
        assert!(!region.is_empty());
    }

    #[test]
    fn test_overlap_region_empty() {
        let region = OverlapRegion {
            r1_start: 10,
            r1_end: 10,
            r2_start: 0,
            r2_end: 0,
            mismatches: 0,
        };

        assert!(region.is_empty());
        assert_eq!(region.len(), 0);
    }

    #[test]
    fn test_calculate_corrected_quality() {
        // High difference - confident correction
        let q1 = calculate_corrected_quality(40 + 33, 10 + 33);
        assert!(q1 > 33); // Should be reasonable quality

        // Low difference - less confident
        let q2 = calculate_corrected_quality(30 + 33, 28 + 33);
        assert!(q2 > 33);
        // When qualities are close, result should be lower than when confident
    }
}