fastx/
qual.rs

1//! Phred quality scores: decoding, statistics and trimming.
2//!
3//! FASTQ stores quality as printable ASCII. Sanger/Illumina 1.8+ files use an
4//! offset of 33 ([`QualityEncoding::Phred33`]); old Illumina 1.3–1.7 files use
5//! 64 ([`QualityEncoding::Phred64`]). Everything in this module takes the offset
6//! explicitly so that neither is assumed silently.
7
8use std::ops::Range;
9
10/// ASCII offset of Sanger / Illumina 1.8+ quality strings.
11pub const PHRED33: u8 = 33;
12/// ASCII offset of Illumina 1.3–1.7 quality strings.
13pub const PHRED64: u8 = 64;
14
15/// The ASCII offset used by a FASTQ quality string.
16#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Default)]
17pub enum QualityEncoding {
18    /// Sanger / Illumina 1.8+ / PacBio / Nanopore: `!` == Q0.
19    #[default]
20    Phred33,
21    /// Illumina 1.3–1.7: `@` == Q0.
22    Phred64,
23}
24
25impl QualityEncoding {
26    /// The ASCII offset for this encoding.
27    pub const fn offset(self) -> u8 {
28        match self {
29            QualityEncoding::Phred33 => PHRED33,
30            QualityEncoding::Phred64 => PHRED64,
31        }
32    }
33
34    /// Guess the encoding from observed quality characters.
35    ///
36    /// Returns `None` when the sample is compatible with both encodings, which
37    /// is common for short high-quality samples — callers should then keep their
38    /// current default rather than guessing.
39    ///
40    /// ```
41    /// use fastx::qual::QualityEncoding;
42    /// assert_eq!(QualityEncoding::detect(b"!!##$$"), Some(QualityEncoding::Phred33));
43    /// assert_eq!(QualityEncoding::detect(b"hhhhhh"), Some(QualityEncoding::Phred64));
44    /// assert_eq!(QualityEncoding::detect(b"IIIIII"), None);
45    /// ```
46    pub fn detect(quality: &[u8]) -> Option<QualityEncoding> {
47        let mut min = u8::MAX;
48        let mut max = 0u8;
49        for &q in quality {
50            min = min.min(q);
51            max = max.max(q);
52        }
53        if quality.is_empty() {
54            return None;
55        }
56        // Phred64 never emits characters below '@' (64); Phred33 rarely exceeds
57        // 'J' (Q41) but long-read platforms can, so only the low end is decisive
58        // for Phred33 and only the high end for Phred64.
59        if min < PHRED64 {
60            Some(QualityEncoding::Phred33)
61        } else if max > b'J' {
62            Some(QualityEncoding::Phred64)
63        } else {
64            None
65        }
66    }
67}
68
69/// Decode one quality character into a Phred score, saturating at 0.
70#[inline]
71pub fn score(ch: u8, offset: u8) -> u8 {
72    ch.saturating_sub(offset)
73}
74
75/// Encode a Phred score as a quality character, clamped to printable ASCII.
76#[inline]
77pub fn encode(score: u8, offset: u8) -> u8 {
78    offset.saturating_add(score.min(126 - offset))
79}
80
81/// Decode a whole quality string into Phred scores.
82pub fn scores(quality: &[u8], offset: u8) -> Vec<u8> {
83    quality.iter().map(|&q| score(q, offset)).collect()
84}
85
86/// Probability that a base with the given Phred score is wrong: `10^(-q/10)`.
87///
88/// ```
89/// # use fastx::qual::error_probability;
90/// assert!((error_probability(20) - 0.01).abs() < 1e-12);
91/// assert!((error_probability(30) - 0.001).abs() < 1e-12);
92/// ```
93#[inline]
94pub fn error_probability(score: u8) -> f64 {
95    10f64.powf(-(score as f64) / 10.0)
96}
97
98/// Phred score corresponding to an error probability, clamped to `0..=93`.
99pub fn probability_to_score(p: f64) -> u8 {
100    if p <= 0.0 {
101        return 93;
102    }
103    let q = -10.0 * p.log10();
104    q.round().clamp(0.0, 93.0) as u8
105}
106
107/// Arithmetic mean of the Phred scores, or `None` for an empty string.
108pub fn mean_score(quality: &[u8], offset: u8) -> Option<f64> {
109    if quality.is_empty() {
110        return None;
111    }
112    let sum: u64 = quality.iter().map(|&q| score(q, offset) as u64).sum();
113    Some(sum as f64 / quality.len() as f64)
114}
115
116/// Sum of per-base error probabilities — the expected number of wrong bases.
117///
118/// This is the quantity `fastp`/`vsearch` filter on, and it is a better read
119/// quality summary than the arithmetic mean of Phred scores.
120pub fn expected_errors(quality: &[u8], offset: u8) -> f64 {
121    quality
122        .iter()
123        .map(|&q| error_probability(score(q, offset)))
124        .sum()
125}
126
127/// Mean quality expressed as a Phred score derived from the mean error rate.
128pub fn mean_quality(quality: &[u8], offset: u8) -> Option<f64> {
129    if quality.is_empty() {
130        return None;
131    }
132    let mean_p = expected_errors(quality, offset) / quality.len() as f64;
133    Some(-10.0 * mean_p.log10())
134}
135
136/// Fraction of bases with a Phred score of at least `threshold` (e.g. Q30).
137pub fn fraction_at_least(quality: &[u8], offset: u8, threshold: u8) -> f64 {
138    if quality.is_empty() {
139        return 0.0;
140    }
141    let n = quality
142        .iter()
143        .filter(|&&q| score(q, offset) >= threshold)
144        .count();
145    n as f64 / quality.len() as f64
146}
147
148/// Trim low-quality bases from both ends, keeping the inner region.
149///
150/// Returns an empty range positioned at 0 when every base is below `min_score`.
151pub fn trim_ends(quality: &[u8], offset: u8, min_score: u8) -> Range<usize> {
152    let start = quality.iter().position(|&q| score(q, offset) >= min_score);
153    match start {
154        None => 0..0,
155        Some(start) => {
156            let end = quality
157                .iter()
158                .rposition(|&q| score(q, offset) >= min_score)
159                .unwrap()
160                + 1;
161            start..end
162        }
163    }
164}
165
166/// Mott-style trimming: the highest-scoring subsequence of `q - threshold`.
167///
168/// This is the algorithm used by `phred`/`seqtk trimfq` and it handles reads
169/// whose quality dips in the middle far better than end trimming.
170pub fn trim_mott(quality: &[u8], offset: u8, threshold: u8) -> Range<usize> {
171    let mut best = 0..0;
172    let mut best_score = 0i64;
173    let mut running = 0i64;
174    let mut start = 0usize;
175
176    for (i, &q) in quality.iter().enumerate() {
177        running += score(q, offset) as i64 - threshold as i64;
178        if running < 0 {
179            running = 0;
180            start = i + 1;
181        } else if running > best_score {
182            best_score = running;
183            best = start..i + 1;
184        }
185    }
186    best
187}
188
189/// Sliding-window trimming from the 3' end, in the spirit of
190/// `trimmomatic SLIDINGWINDOW`.
191///
192/// Scans left to right and cuts at the first position where the mean score of
193/// the following `window` bases drops below `min_mean`; the returned range ends
194/// at the start of that window. Reads shorter than the window are kept intact.
195pub fn trim_sliding_window(
196    quality: &[u8],
197    offset: u8,
198    window: usize,
199    min_mean: f64,
200) -> Range<usize> {
201    if window == 0 || quality.len() < window {
202        return 0..quality.len();
203    }
204    let threshold = min_mean * window as f64;
205    let mut sum: f64 = quality[..window]
206        .iter()
207        .map(|&q| score(q, offset) as f64)
208        .sum();
209    if sum < threshold {
210        return 0..0;
211    }
212    for i in window..quality.len() {
213        sum += score(quality[i], offset) as f64;
214        sum -= score(quality[i - window], offset) as f64;
215        if sum < threshold {
216            return 0..i + 1 - window;
217        }
218    }
219    0..quality.len()
220}
221
222#[cfg(test)]
223mod tests {
224    use super::*;
225
226    #[test]
227    fn scores_round_trip() {
228        for s in 0u8..=60 {
229            assert_eq!(score(encode(s, PHRED33), PHRED33), s);
230        }
231    }
232
233    #[test]
234    fn error_probabilities() {
235        assert_eq!(probability_to_score(0.001), 30);
236        assert_eq!(probability_to_score(1.0), 0);
237        assert!((expected_errors(b"!!!!", PHRED33) - 4.0).abs() < 1e-9);
238    }
239
240    #[test]
241    fn mean_quality_is_error_weighted() {
242        // One terrible base drags the error-weighted mean far below the
243        // arithmetic mean of the Phred scores.
244        let q = b"IIIII!"; // five Q40 and one Q0
245        let arithmetic = mean_score(q, PHRED33).unwrap();
246        let weighted = mean_quality(q, PHRED33).unwrap();
247        assert!(arithmetic > 33.0, "{arithmetic}");
248        assert!(weighted < 8.0, "{weighted}");
249    }
250
251    #[test]
252    fn trims_ends() {
253        assert_eq!(trim_ends(b"!!III!!", PHRED33, 30), 2..5);
254        assert_eq!(trim_ends(b"!!!!", PHRED33, 30), 0..0);
255        assert_eq!(trim_ends(b"IIII", PHRED33, 30), 0..4);
256    }
257
258    #[test]
259    fn mott_keeps_best_window() {
260        // Low quality at both ends, a good stretch in the middle.
261        let q = b"###IIIIIIII###";
262        assert_eq!(trim_mott(q, PHRED33, 20), 3..11);
263        assert_eq!(trim_mott(b"####", PHRED33, 20), 0..0);
264    }
265
266    #[test]
267    fn sliding_window_cuts_at_drop() {
268        // Eight Q40 bases then four Q2 bases: the first window that fails
269        // starts at index 7, so everything from there on is discarded.
270        let q = b"IIIIIIII####";
271        assert_eq!(trim_sliding_window(q, PHRED33, 4, 20.0), 0..7);
272        assert_eq!(trim_sliding_window(b"IIIIIIII", PHRED33, 4, 20.0), 0..8);
273        assert_eq!(trim_sliding_window(b"II", PHRED33, 4, 20.0), 0..2);
274        assert_eq!(trim_sliding_window(b"####IIII", PHRED33, 4, 20.0), 0..0);
275    }
276
277    #[test]
278    fn q30_fraction() {
279        assert!((fraction_at_least(b"IIII!!!!", PHRED33, 30) - 0.5).abs() < 1e-9);
280        assert_eq!(fraction_at_least(b"", PHRED33, 30), 0.0);
281    }
282}
fastx/qual.rs

fastx/
qual.rs
