fastx/
stats.rs

1//! Summary statistics over a set of records — the `seqkit stats` equivalent.
2
3use std::fmt;
4
5use crate::qual::{self, PHRED33};
6use crate::record::Sequence;
7use crate::seq::{self, BaseCounts};
8
9/// Accumulates length, composition and quality statistics.
10///
11/// ```
12/// use fastx::{FastxReader, SeqStats};
13///
14/// let data = b">a\nACGTACGTAC\n>b\nGGCC\n";
15/// let mut stats = SeqStats::default();
16/// FastxReader::new(&data[..]).for_each_record(|r| { stats.push(r); Ok(()) })?;
17///
18/// assert_eq!(stats.count, 2);
19/// assert_eq!(stats.total_length, 14);
20/// assert_eq!(stats.min_length, Some(4));
21/// assert_eq!(stats.max_length, Some(10));
22/// assert_eq!(stats.n50(), Some(10));
23/// # Ok::<(), fastx::Error>(())
24/// ```
25#[derive(Debug, Clone, Default)]
26pub struct SeqStats {
27    /// Number of records seen.
28    pub count: u64,
29    /// Sum of all sequence lengths.
30    pub total_length: u64,
31    /// Shortest record, if any.
32    pub min_length: Option<u64>,
33    /// Longest record, if any.
34    pub max_length: Option<u64>,
35    /// Base composition across all records.
36    pub bases: BaseCounts,
37    /// Every observed length, kept so that N50 and the median can be computed.
38    lengths: Vec<u64>,
39    /// Counts per Phred score, from which every quality figure is derived.
40    quality: QualityHistogram,
41}
42
43/// Number of distinct Phred scores tracked: 0..=93 covers all printable ASCII.
44pub const PHRED_SCORES: usize = 94;
45
46/// Highest Phred score the histogram can hold.
47const MAX_SCORE: u8 = (PHRED_SCORES - 1) as u8;
48
49/// Counts of each Phred score, indexed by the score itself.
50///
51/// Accumulating a histogram keeps the hot loop to one increment per base — no
52/// floating point at all — and every quality statistic falls out of it
53/// afterwards at a cost that does not depend on the size of the input.
54#[derive(Debug, Clone)]
55struct QualityHistogram([u64; PHRED_SCORES]);
56
57impl Default for QualityHistogram {
58    fn default() -> Self {
59        QualityHistogram([0; PHRED_SCORES])
60    }
61}
62
63impl SeqStats {
64    /// An empty accumulator.
65    pub fn new() -> SeqStats {
66        SeqStats::default()
67    }
68
69    /// Fold one record in.
70    pub fn push(&mut self, record: &Sequence) {
71        let len = record.len() as u64;
72        self.count += 1;
73        self.total_length += len;
74        self.min_length = Some(self.min_length.map_or(len, |m| m.min(len)));
75        self.max_length = Some(self.max_length.map_or(len, |m| m.max(len)));
76        self.bases.merge(&record.base_counts());
77        self.lengths.push(len);
78        if let Some(quality) = &record.quality {
79            let histogram = &mut self.quality.0;
80            for &c in quality {
81                // `min` also keeps the index provably in range, which removes
82                // the bounds check from the loop.
83                histogram[qual::score(c, PHRED33).min(MAX_SCORE) as usize] += 1;
84            }
85        }
86    }
87
88    /// Merge another accumulator, for parallel or per-file aggregation.
89    pub fn merge(&mut self, other: &SeqStats) {
90        self.count += other.count;
91        self.total_length += other.total_length;
92        self.min_length = min_option(self.min_length, other.min_length);
93        self.max_length = max_option(self.max_length, other.max_length);
94        self.bases.merge(&other.bases);
95        self.lengths.extend_from_slice(&other.lengths);
96        for (slot, count) in self.quality.0.iter_mut().zip(other.quality.0.iter()) {
97            *slot += count;
98        }
99    }
100
101    /// True when no records have been seen.
102    pub fn is_empty(&self) -> bool {
103        self.count == 0
104    }
105
106    /// Mean sequence length.
107    pub fn mean_length(&self) -> Option<f64> {
108        if self.count == 0 {
109            None
110        } else {
111            Some(self.total_length as f64 / self.count as f64)
112        }
113    }
114
115    /// Median sequence length (lower median for an even count).
116    pub fn median_length(&self) -> Option<u64> {
117        if self.lengths.is_empty() {
118            return None;
119        }
120        let mut lengths = self.lengths.clone();
121        lengths.sort_unstable();
122        Some(lengths[(lengths.len() - 1) / 2])
123    }
124
125    /// N50: contigs of at least this length cover half the total.
126    pub fn n50(&self) -> Option<u64> {
127        seq::n50(&mut self.lengths.clone())
128    }
129
130    /// N90.
131    pub fn n90(&self) -> Option<u64> {
132        seq::nx(&mut self.lengths.clone(), 0.9)
133    }
134
135    /// L50: the number of contigs that make up the N50.
136    pub fn l50(&self) -> Option<u64> {
137        let n50 = self.n50()?;
138        let mut lengths = self.lengths.clone();
139        lengths.sort_unstable_by(|a, b| b.cmp(a));
140        Some(lengths.iter().take_while(|&&l| l >= n50).count() as u64)
141    }
142
143    /// GC fraction over unambiguous bases.
144    pub fn gc_content(&self) -> Option<f64> {
145        self.bases.gc_content()
146    }
147
148    /// Number of quality characters seen across all records.
149    pub fn quality_bases(&self) -> u64 {
150        self.quality.0.iter().sum()
151    }
152
153    /// Counts per Phred score, indexed by score — the same data FastQC plots.
154    ///
155    /// ```
156    /// # use fastx::{SeqStats, Sequence};
157    /// let mut stats = SeqStats::new();
158    /// stats.push(&Sequence::fastq("r", b"ACGT", b"IIII")?);
159    /// assert_eq!(stats.quality_histogram()[40], 4); // 'I' is Q40
160    /// # Ok::<(), fastx::Error>(())
161    /// ```
162    pub fn quality_histogram(&self) -> &[u64; PHRED_SCORES] {
163        &self.quality.0
164    }
165
166    /// Expected number of wrong bases across all records.
167    pub fn expected_errors(&self) -> f64 {
168        self.quality
169            .0
170            .iter()
171            .enumerate()
172            .map(|(score, &count)| count as f64 * qual::error_probability(score as u8))
173            .sum()
174    }
175
176    /// Mean quality as a Phred score, derived from the mean error rate.
177    pub fn mean_quality(&self) -> Option<f64> {
178        let bases = self.quality_bases();
179        if bases == 0 {
180            return None;
181        }
182        let mean_p = self.expected_errors() / bases as f64;
183        Some(-10.0 * mean_p.log10())
184    }
185
186    /// Fraction of bases at Q20 or better.
187    pub fn q20_fraction(&self) -> Option<f64> {
188        self.fraction_at_least(20)
189    }
190
191    /// Fraction of bases at Q30 or better.
192    pub fn q30_fraction(&self) -> Option<f64> {
193        self.fraction_at_least(30)
194    }
195
196    /// Fraction of bases whose Phred score is at least `score`.
197    pub fn fraction_at_least(&self, score: u8) -> Option<f64> {
198        let bases = self.quality_bases();
199        if bases == 0 {
200            return None;
201        }
202        let at_least: u64 = self.quality.0[(score as usize).min(PHRED_SCORES)..]
203            .iter()
204            .sum();
205        Some(at_least as f64 / bases as f64)
206    }
207
208    /// All observed lengths, in the order the records were seen.
209    pub fn lengths(&self) -> &[u64] {
210        &self.lengths
211    }
212
213    /// Render the statistics as a JSON object, for pipelines that parse output.
214    ///
215    /// The object is on one line, so a run over several files is valid
216    /// line-delimited JSON and streams straight into `jq`. Pipe through `jq .`
217    /// if you want it laid out.
218    ///
219    /// Absent figures — quality for FASTA, N50 after
220    /// [`SeqStats::forget_lengths`] — come out as `null` rather than being
221    /// omitted, so the shape of the object never changes.
222    ///
223    /// ```
224    /// # use fastx::{SeqStats, Sequence};
225    /// let mut stats = SeqStats::new();
226    /// stats.push(&Sequence::fasta("a", b"ACGT"));
227    /// let json = stats.to_json();
228    /// assert!(json.starts_with('{') && json.ends_with('}'));
229    /// assert!(!json.contains('\n'), "must stay on one line");
230    /// assert!(json.contains("\"records\":1"));
231    /// assert!(json.contains("\"mean_quality\":null"));
232    /// ```
233    pub fn to_json(&self) -> String {
234        fn number(value: Option<f64>, decimals: usize) -> String {
235            match value {
236                // JSON has no NaN or Infinity, so anything not finite is null.
237                Some(v) if v.is_finite() => format!("{v:.decimals$}"),
238                _ => "null".to_string(),
239            }
240        }
241        fn integer(value: Option<u64>) -> String {
242            value.map_or_else(|| "null".to_string(), |v| v.to_string())
243        }
244
245        let fields = [
246            ("records".to_string(), self.count.to_string()),
247            ("total_length".to_string(), self.total_length.to_string()),
248            ("min_length".to_string(), integer(self.min_length)),
249            ("max_length".to_string(), integer(self.max_length)),
250            ("mean_length".to_string(), number(self.mean_length(), 2)),
251            ("median_length".to_string(), integer(self.median_length())),
252            ("n50".to_string(), integer(self.n50())),
253            ("n90".to_string(), integer(self.n90())),
254            ("l50".to_string(), integer(self.l50())),
255            ("gc_content".to_string(), number(self.gc_content(), 6)),
256            ("a".to_string(), self.bases.a.to_string()),
257            ("c".to_string(), self.bases.c.to_string()),
258            ("g".to_string(), self.bases.g.to_string()),
259            ("t".to_string(), self.bases.t.to_string()),
260            ("ambiguous".to_string(), self.bases.n.to_string()),
261            ("other".to_string(), self.bases.other.to_string()),
262            (
263                "quality_bases".to_string(),
264                self.quality_bases().to_string(),
265            ),
266            ("mean_quality".to_string(), number(self.mean_quality(), 4)),
267            ("q20_fraction".to_string(), number(self.q20_fraction(), 6)),
268            ("q30_fraction".to_string(), number(self.q30_fraction(), 6)),
269        ];
270        let body = fields
271            .iter()
272            .map(|(key, value)| format!("\"{key}\":{value}"))
273            .collect::<Vec<_>>()
274            .join(",");
275        format!("{{{body}}}")
276    }
277
278    /// Drop the retained per-record lengths to cap memory on huge inputs.
279    ///
280    /// After this call [`SeqStats::n50`], [`SeqStats::median_length`] and
281    /// [`SeqStats::l50`] return `None`, but counts and means stay correct.
282    pub fn forget_lengths(&mut self) {
283        self.lengths = Vec::new();
284    }
285}
286
287fn min_option(a: Option<u64>, b: Option<u64>) -> Option<u64> {
288    match (a, b) {
289        (Some(a), Some(b)) => Some(a.min(b)),
290        (a, b) => a.or(b),
291    }
292}
293
294fn max_option(a: Option<u64>, b: Option<u64>) -> Option<u64> {
295    match (a, b) {
296        (Some(a), Some(b)) => Some(a.max(b)),
297        (a, b) => a.or(b),
298    }
299}
300
301impl fmt::Display for SeqStats {
302    /// A block of aligned `label value` lines, without a trailing newline.
303    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
304        let mut lines: Vec<String> = Vec::with_capacity(14);
305        let mut push = |label: &str, value: String| lines.push(format!("{label:<12} {value}"));
306
307        push("records", self.count.to_string());
308        push("total bases", self.total_length.to_string());
309        if let Some(v) = self.min_length {
310            push("min length", v.to_string());
311        }
312        if let Some(v) = self.max_length {
313            push("max length", v.to_string());
314        }
315        if let Some(v) = self.mean_length() {
316            push("avg length", format!("{v:.1}"));
317        }
318        if let Some(v) = self.median_length() {
319            push("median len", v.to_string());
320        }
321        if let Some(v) = self.n50() {
322            push("N50", v.to_string());
323        }
324        if let Some(v) = self.n90() {
325            push("N90", v.to_string());
326        }
327        if let Some(v) = self.l50() {
328            push("L50", v.to_string());
329        }
330        if let Some(v) = self.gc_content() {
331            push("GC%", format!("{:.2}", v * 100.0));
332        }
333        push("N bases", self.bases.n.to_string());
334        if let Some(v) = self.mean_quality() {
335            push("avg quality", format!("Q{v:.1}"));
336        }
337        if let Some(v) = self.q20_fraction() {
338            push("Q20%", format!("{:.2}", v * 100.0));
339        }
340        if let Some(v) = self.q30_fraction() {
341            push("Q30%", format!("{:.2}", v * 100.0));
342        }
343        f.write_str(&lines.join("\n"))
344    }
345}
346
347#[cfg(test)]
348mod tests {
349    use super::*;
350
351    fn stats_of(records: &[Sequence]) -> SeqStats {
352        let mut stats = SeqStats::new();
353        for record in records {
354            stats.push(record);
355        }
356        stats
357    }
358
359    #[test]
360    fn empty_stats_have_no_summaries() {
361        let stats = SeqStats::new();
362        assert!(stats.is_empty());
363        assert_eq!(stats.mean_length(), None);
364        assert_eq!(stats.n50(), None);
365        assert_eq!(stats.l50(), None);
366        assert_eq!(stats.gc_content(), None);
367        assert_eq!(stats.mean_quality(), None);
368    }
369
370    #[test]
371    fn length_statistics() {
372        let stats = stats_of(&[
373            Sequence::fasta("a", b"A".repeat(50)),
374            Sequence::fasta("b", b"C".repeat(30)),
375            Sequence::fasta("c", b"G".repeat(15)),
376            Sequence::fasta("d", b"T".repeat(5)),
377        ]);
378        assert_eq!(stats.count, 4);
379        assert_eq!(stats.total_length, 100);
380        assert_eq!(stats.min_length, Some(5));
381        assert_eq!(stats.max_length, Some(50));
382        assert_eq!(stats.mean_length(), Some(25.0));
383        assert_eq!(stats.median_length(), Some(15));
384        assert_eq!(stats.n50(), Some(50));
385        assert_eq!(stats.n90(), Some(15));
386        assert_eq!(stats.l50(), Some(1));
387        assert_eq!(stats.gc_content(), Some(0.45));
388    }
389
390    #[test]
391    fn quality_statistics() {
392        let stats = stats_of(&[
393            Sequence::fastq("a", b"ACGT", b"IIII").unwrap(),
394            Sequence::fastq("b", b"ACGT", b"!!!!").unwrap(),
395        ]);
396        assert_eq!(stats.q30_fraction(), Some(0.5));
397        assert_eq!(stats.q20_fraction(), Some(0.5));
398        let mean = stats.mean_quality().unwrap();
399        assert!(mean > 2.0 && mean < 4.0, "{mean}");
400    }
401
402    #[test]
403    fn quality_histogram_is_the_source_of_truth() {
404        let stats = stats_of(&[
405            Sequence::fastq("a", b"ACGT", b"IIII").unwrap(), // Q40
406            Sequence::fastq("b", b"AC", b"!5").unwrap(),     // Q0 and Q20
407        ]);
408        let histogram = stats.quality_histogram();
409        assert_eq!(histogram[40], 4);
410        assert_eq!(histogram[20], 1);
411        assert_eq!(histogram[0], 1);
412        assert_eq!(stats.quality_bases(), 6);
413
414        // Derived figures must agree with computing them the direct way.
415        let direct: f64 = qual::expected_errors(b"IIII!5", PHRED33);
416        assert!((stats.expected_errors() - direct).abs() < 1e-12);
417        assert_eq!(stats.fraction_at_least(0), Some(1.0));
418        assert_eq!(stats.fraction_at_least(20), Some(5.0 / 6.0));
419        assert_eq!(stats.fraction_at_least(40), Some(4.0 / 6.0));
420        assert_eq!(stats.fraction_at_least(93), Some(0.0));
421        // Scores beyond the histogram are simply never reached.
422        assert_eq!(stats.fraction_at_least(200), Some(0.0));
423    }
424
425    #[test]
426    fn quality_scores_are_clamped_not_wrapped() {
427        // A byte above '~' cannot appear in valid FASTQ, but it must not panic
428        // or corrupt neighbouring buckets if it does.
429        let mut record = Sequence::fastq("r", b"AC", b"II").unwrap();
430        record.quality = Some(vec![255, 33]);
431        let stats = stats_of(&[record]);
432        assert_eq!(stats.quality_histogram()[93], 1);
433        assert_eq!(stats.quality_histogram()[0], 1);
434        assert_eq!(stats.quality_bases(), 2);
435    }
436
437    #[test]
438    fn merging_matches_sequential() {
439        let records: Vec<Sequence> = (1..20)
440            .map(|i| Sequence::fasta(format!("s{i}"), b"ACGT".repeat(i)))
441            .collect();
442        let sequential = stats_of(&records);
443        let (left, right) = records.split_at(7);
444        let mut merged = stats_of(left);
445        merged.merge(&stats_of(right));
446
447        assert_eq!(merged.count, sequential.count);
448        assert_eq!(merged.total_length, sequential.total_length);
449        assert_eq!(merged.min_length, sequential.min_length);
450        assert_eq!(merged.max_length, sequential.max_length);
451        assert_eq!(merged.n50(), sequential.n50());
452        assert_eq!(merged.gc_content(), sequential.gc_content());
453    }
454
455    #[test]
456    fn forgetting_lengths_keeps_counts() {
457        let mut stats = stats_of(&[Sequence::fasta("a", b"ACGT")]);
458        stats.forget_lengths();
459        assert_eq!(stats.count, 1);
460        assert_eq!(stats.total_length, 4);
461        assert_eq!(stats.n50(), None);
462    }
463
464    #[test]
465    fn display_is_readable() {
466        let text = stats_of(&[Sequence::fastq("a", b"ACGT", b"IIII").unwrap()]).to_string();
467        assert!(text.contains("records      1"));
468        assert!(text.contains("GC%"));
469        assert!(text.contains("Q30%"));
470    }
471}
fastx/stats.rs

fastx/
stats.rs
