fastx/

qual.rs

//! Phred quality scores: decoding, statistics and trimming.
//!
//! FASTQ stores quality as printable ASCII. Sanger/Illumina 1.8+ files use an
//! offset of 33 ([`QualityEncoding::Phred33`]); old Illumina 1.3–1.7 files use
//! 64 ([`QualityEncoding::Phred64`]). Everything in this module takes the offset
//! explicitly so that neither is assumed silently.

use std::ops::Range;

/// ASCII offset of Sanger / Illumina 1.8+ quality strings.
pub const PHRED33: u8 = 33;
/// ASCII offset of Illumina 1.3–1.7 quality strings.
pub const PHRED64: u8 = 64;

/// The ASCII offset used by a FASTQ quality string.
///
/// Deliberately *not* `#[non_exhaustive]`: these two are the whole of the
/// offset-encoded world. Solexa's old scheme is not a third offset but a
/// different formula, so it could not join this enum without changing what
/// [`QualityEncoding::offset`] means.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Default)]
pub enum QualityEncoding {
    /// Sanger / Illumina 1.8+ / PacBio / Nanopore: `!` == Q0.
    #[default]
    Phred33,
    /// Illumina 1.3–1.7: `@` == Q0.
    Phred64,
}

impl QualityEncoding {
    /// The ASCII offset for this encoding.
    pub const fn offset(self) -> u8 {
        match self {
            QualityEncoding::Phred33 => PHRED33,
            QualityEncoding::Phred64 => PHRED64,
        }
    }

    /// Guess the encoding from observed quality characters.
    ///
    /// Returns `None` when the sample is compatible with both encodings, which
    /// is common for short high-quality samples — callers should then keep their
    /// current default rather than guessing.
    ///
    /// ```
    /// use fastx::qual::QualityEncoding;
    /// assert_eq!(QualityEncoding::detect(b"!!##$$"), Some(QualityEncoding::Phred33));
    /// assert_eq!(QualityEncoding::detect(b"hhhhhh"), Some(QualityEncoding::Phred64));
    /// assert_eq!(QualityEncoding::detect(b"IIIIII"), None);
    /// ```
    pub fn detect(quality: &[u8]) -> Option<QualityEncoding> {
        let mut min = u8::MAX;
        let mut max = 0u8;
        for &q in quality {
            min = min.min(q);
            max = max.max(q);
        }
        if quality.is_empty() {
            return None;
        }
        // Phred64 never emits characters below '@' (64); Phred33 rarely exceeds
        // 'J' (Q41) but long-read platforms can, so only the low end is decisive
        // for Phred33 and only the high end for Phred64.
        if min < PHRED64 {
            Some(QualityEncoding::Phred33)
        } else if max > b'J' {
            Some(QualityEncoding::Phred64)
        } else {
            None
        }
    }
}

/// Decode one quality character into a Phred score, saturating at 0.
#[inline]
pub fn score(ch: u8, offset: u8) -> u8 {
    ch.saturating_sub(offset)
}

/// Encode a Phred score as a quality character, clamped to printable ASCII.
#[inline]
pub fn encode(score: u8, offset: u8) -> u8 {
    offset.saturating_add(score.min(126 - offset))
}

/// Decode a whole quality string into Phred scores.
pub fn scores(quality: &[u8], offset: u8) -> Vec<u8> {
    quality.iter().map(|&q| score(q, offset)).collect()
}

/// Probability that a base with the given Phred score is wrong: `10^(-q/10)`.
///
/// ```
/// # use fastx::qual::error_probability;
/// assert!((error_probability(20) - 0.01).abs() < 1e-12);
/// assert!((error_probability(30) - 0.001).abs() < 1e-12);
/// ```
#[inline]
pub fn error_probability(score: u8) -> f64 {
    10f64.powf(-(score as f64) / 10.0)
}

/// Phred score corresponding to an error probability, clamped to `0..=93`.
pub fn probability_to_score(p: f64) -> u8 {
    if p <= 0.0 {
        return 93;
    }
    let q = -10.0 * p.log10();
    q.round().clamp(0.0, 93.0) as u8
}

/// Arithmetic mean of the Phred scores, or `None` for an empty string.
pub fn mean_score(quality: &[u8], offset: u8) -> Option<f64> {
    if quality.is_empty() {
        return None;
    }
    let sum: u64 = quality.iter().map(|&q| score(q, offset) as u64).sum();
    Some(sum as f64 / quality.len() as f64)
}

/// Sum of per-base error probabilities — the expected number of wrong bases.
///
/// This is the quantity `fastp`/`vsearch` filter on, and it is a better read
/// quality summary than the arithmetic mean of Phred scores.
pub fn expected_errors(quality: &[u8], offset: u8) -> f64 {
    quality
        .iter()
        .map(|&q| error_probability(score(q, offset)))
        .sum()
}

/// Mean quality expressed as a Phred score derived from the mean error rate.
pub fn mean_quality(quality: &[u8], offset: u8) -> Option<f64> {
    if quality.is_empty() {
        return None;
    }
    let mean_p = expected_errors(quality, offset) / quality.len() as f64;
    Some(-10.0 * mean_p.log10())
}

/// Fraction of bases with a Phred score of at least `threshold` (e.g. Q30).
pub fn fraction_at_least(quality: &[u8], offset: u8, threshold: u8) -> f64 {
    if quality.is_empty() {
        return 0.0;
    }
    let n = quality
        .iter()
        .filter(|&&q| score(q, offset) >= threshold)
        .count();
    n as f64 / quality.len() as f64
}

/// Trim low-quality bases from both ends, keeping the inner region.
///
/// Returns an empty range positioned at 0 when every base is below `min_score`.
pub fn trim_ends(quality: &[u8], offset: u8, min_score: u8) -> Range<usize> {
    let start = quality.iter().position(|&q| score(q, offset) >= min_score);
    match start {
        None => 0..0,
        Some(start) => {
            let end = quality
                .iter()
                .rposition(|&q| score(q, offset) >= min_score)
                .unwrap()
                + 1;
            start..end
        }
    }
}

/// Mott-style trimming: the highest-scoring subsequence of `q - threshold`.
///
/// This is the algorithm used by `phred`/`seqtk trimfq` and it handles reads
/// whose quality dips in the middle far better than end trimming.
pub fn trim_mott(quality: &[u8], offset: u8, threshold: u8) -> Range<usize> {
    let mut best = 0..0;
    let mut best_score = 0i64;
    let mut running = 0i64;
    let mut start = 0usize;

    for (i, &q) in quality.iter().enumerate() {
        running += score(q, offset) as i64 - threshold as i64;
        if running < 0 {
            running = 0;
            start = i + 1;
        } else if running > best_score {
            best_score = running;
            best = start..i + 1;
        }
    }
    best
}

/// Sliding-window trimming from the 3' end, in the spirit of
/// `trimmomatic SLIDINGWINDOW`.
///
/// Scans left to right and cuts at the first position where the mean score of
/// the following `window` bases drops below `min_mean`; the returned range ends
/// at the start of that window. Reads shorter than the window are kept intact.
pub fn trim_sliding_window(
    quality: &[u8],
    offset: u8,
    window: usize,
    min_mean: f64,
) -> Range<usize> {
    if window == 0 || quality.len() < window {
        return 0..quality.len();
    }
    let threshold = min_mean * window as f64;
    let mut sum: f64 = quality[..window]
        .iter()
        .map(|&q| score(q, offset) as f64)
        .sum();
    if sum < threshold {
        return 0..0;
    }
    for i in window..quality.len() {
        sum += score(quality[i], offset) as f64;
        sum -= score(quality[i - window], offset) as f64;
        if sum < threshold {
            return 0..i + 1 - window;
        }
    }
    0..quality.len()
}

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

    #[test]
    fn scores_round_trip() {
        for s in 0u8..=60 {
            assert_eq!(score(encode(s, PHRED33), PHRED33), s);
        }
    }

    #[test]
    fn error_probabilities() {
        assert_eq!(probability_to_score(0.001), 30);
        assert_eq!(probability_to_score(1.0), 0);
        assert!((expected_errors(b"!!!!", PHRED33) - 4.0).abs() < 1e-9);
    }

    #[test]
    fn mean_quality_is_error_weighted() {
        // One terrible base drags the error-weighted mean far below the
        // arithmetic mean of the Phred scores.
        let q = b"IIIII!"; // five Q40 and one Q0
        let arithmetic = mean_score(q, PHRED33).unwrap();
        let weighted = mean_quality(q, PHRED33).unwrap();
        assert!(arithmetic > 33.0, "{arithmetic}");
        assert!(weighted < 8.0, "{weighted}");
    }

    #[test]
    fn trims_ends() {
        assert_eq!(trim_ends(b"!!III!!", PHRED33, 30), 2..5);
        assert_eq!(trim_ends(b"!!!!", PHRED33, 30), 0..0);
        assert_eq!(trim_ends(b"IIII", PHRED33, 30), 0..4);
    }

    #[test]
    fn mott_keeps_best_window() {
        // Low quality at both ends, a good stretch in the middle.
        let q = b"###IIIIIIII###";
        assert_eq!(trim_mott(q, PHRED33, 20), 3..11);
        assert_eq!(trim_mott(b"####", PHRED33, 20), 0..0);
    }

    #[test]
    fn sliding_window_cuts_at_drop() {
        // Eight Q40 bases then four Q2 bases: the first window that fails
        // starts at index 7, so everything from there on is discarded.
        let q = b"IIIIIIII####";
        assert_eq!(trim_sliding_window(q, PHRED33, 4, 20.0), 0..7);
        assert_eq!(trim_sliding_window(b"IIIIIIII", PHRED33, 4, 20.0), 0..8);
        assert_eq!(trim_sliding_window(b"II", PHRED33, 4, 20.0), 0..2);
        assert_eq!(trim_sliding_window(b"####IIII", PHRED33, 4, 20.0), 0..0);
    }

    #[test]
    fn q30_fraction() {
        assert!((fraction_at_least(b"IIII!!!!", PHRED33, 30) - 0.5).abs() < 1e-9);
        assert_eq!(fraction_at_least(b"", PHRED33, 30), 0.0);
    }
}