fastqc_rust/modules/
per_sequence_quality.rs

1// Per Sequence Quality Scores module
2// Corresponds to Modules/PerSequenceQualityScores.java
3
4use std::io;
5
6use crate::config::{Limits, LimitsExt};
7use crate::modules::QCModule;
8use crate::report::charts::line_graph::{render_line_graph, LineGraphData};
9use crate::report::charts::CHART_WIDTH;
10use crate::sequence::Sequence;
11use crate::utils::format::java_format_double;
12use crate::utils::phred;
13
14/// Maximum raw ASCII average quality score we can track.
15/// Quality chars are in range 0-127; the average of any sequence of such chars
16/// fits in the same range. 128 slots covers all possible values.
17const MAX_QUALITY_SCORE: usize = 128;
18
19pub struct PerSequenceQualityScores {
20    // JAVA COMPAT: Indexed by raw ASCII average quality (before offset subtraction),
21    // computed using integer arithmetic: sum of quality chars / length.
22    // Using a fixed array instead of HashMap eliminates hashing on every read.
23    average_score_counts: [u64; MAX_QUALITY_SCORE],
24    has_data: bool,
25    lowest_char: u16,
26    // Set by QCModule::set_phred_encoding; see the trait docs.
27    known_encoding: Option<phred::PhredEncoding>,
28    limits: Limits,
29}
30
31impl PerSequenceQualityScores {
32    pub fn new(limits: &Limits) -> Self {
33        PerSequenceQualityScores {
34            average_score_counts: [0u64; MAX_QUALITY_SCORE],
35            has_data: false,
36            lowest_char: phred::NO_QUALITY_SEEN,
37            known_encoding: None,
38            limits: limits.clone(),
39        }
40    }
41
42    fn calculate(&self) -> Option<DistributionData> {
43        if !self.has_data {
44            return None;
45        }
46
47        let encoding = phred::resolve(self.known_encoding, self.lowest_char)
48            .unwrap_or(phred::PhredEncoding::SANGER);
49
50        // Find the range of scores with non-zero counts
51        let mut range_start: Option<usize> = None;
52        let mut range_end: usize = 0;
53        for i in 0..MAX_QUALITY_SCORE {
54            if self.average_score_counts[i] > 0 {
55                if range_start.is_none() {
56                    range_start = Some(i);
57                }
58                range_end = i;
59            }
60        }
61
62        let range_start = range_start? as i32;
63        let range_end = range_end as i32;
64
65        // Distribution runs from lowest to highest raw score
66        let len = (1 + range_end - range_start) as usize;
67
68        let mut quality_distribution = vec![0.0f64; len];
69        let mut x_categories = Vec::with_capacity(len);
70
71        // Build distribution and x_categories arrays in parallel
72        for (i, qd) in quality_distribution.iter_mut().enumerate() {
73            x_categories.push(range_start + i as i32 - encoding.offset as i32);
74            let key = (range_start + i as i32) as usize;
75            *qd = self.average_score_counts[key] as f64;
76        }
77
78        // Find most frequent score
79        let mut max_count = 0.0;
80        let mut most_frequent_score = 0;
81        // index needed for both quality_distribution and x_categories
82        for (&qd, &xc) in quality_distribution.iter().zip(x_categories.iter()) {
83            if qd > max_count {
84                max_count = qd;
85                most_frequent_score = xc;
86            }
87        }
88
89        Some(DistributionData {
90            quality_distribution,
91            x_categories,
92            most_frequent_score,
93        })
94    }
95}
96
97impl PerSequenceQualityScores {
98    fn build_chart_svg(&self) -> Option<String> {
99        let data = self.calculate()?;
100
101        let max_count = data
102            .quality_distribution
103            .iter()
104            .cloned()
105            .fold(0.0_f64, f64::max);
106        // Java passes raw max to LineGraph (no ceil rounding)
107        let max_y = max_count;
108
109        let x_categories: Vec<String> =
110            data.x_categories.iter().map(|v| format!("{}", v)).collect();
111
112        Some(render_line_graph(&LineGraphData {
113            width: CHART_WIDTH,
114            data: vec![data.quality_distribution],
115            min_y: 0.0,
116            max_y,
117            // These labels match the Java constructor call
118            x_label: "Mean Sequence Quality (Phred Score)".to_string(),
119            series_names: vec!["Average Quality per read".to_string()],
120            x_categories,
121            title: "Quality score distribution over all sequences".to_string(),
122        }))
123    }
124}
125
126impl QCModule for PerSequenceQualityScores {
127    fn cost_hint(&self) -> u32 {
128        1
129    }
130
131    fn process_sequence(&mut self, sequence: &Sequence) {
132        let qual = &sequence.quality;
133
134        // JAVA COMPAT: Average quality computed using integer arithmetic on raw ASCII values.
135        // sum of quality chars (as int) / length (integer division), stored as raw value.
136        let mut average_quality: i32 = 0;
137
138        for &q in qual.iter() {
139            self.lowest_char = self.lowest_char.min(q as u16);
140            average_quality += q as i32;
141        }
142
143        if !qual.is_empty() {
144            // JAVA COMPAT: Integer division truncates towards zero, matching Java's `/`
145            average_quality /= qual.len() as i32;
146
147            // Clamp rather than index blindly. Java throws
148            // ArrayIndexOutOfBoundsException and kills the run; a file whose
149            // quality line holds bytes at or above MAX_QUALITY_SCORE (a
150            // mis-encoded or non-ASCII line) would otherwise panic here, and
151            // under the parallel pipeline take a worker thread with it. Same
152            // treatment as QualityCount::add_value: the value is wrong for that
153            // read, but the rest of the file still gets analysed. Valid input
154            // never reaches this, so byte-identical output is unaffected.
155            let slot = (average_quality as usize).min(MAX_QUALITY_SCORE - 1);
156            if slot != average_quality as usize {
157                static WARNED: crate::progress::OncePerRun = crate::progress::OncePerRun::new();
158                WARNED.log(|| {
159                    format!(
160                        "Warning: mean quality {} exceeds maximum {}; clamping",
161                        average_quality,
162                        MAX_QUALITY_SCORE - 1
163                    )
164                });
165            }
166            self.average_score_counts[slot] += 1;
167            self.has_data = true;
168        }
169    }
170
171    fn set_phred_encoding(&mut self, encoding: phred::PhredEncoding) {
172        self.known_encoding = Some(encoding);
173    }
174
175    fn name(&self) -> &str {
176        "Per sequence quality scores"
177    }
178
179    fn description(&self) -> &str {
180        "Shows the distribution of average quality scores for whole sequences"
181    }
182
183    fn reset(&mut self) {
184        self.average_score_counts = [0u64; MAX_QUALITY_SCORE];
185        self.has_data = false;
186        self.lowest_char = phred::NO_QUALITY_SEEN;
187    }
188
189    fn raises_error(&self) -> bool {
190        let error_threshold = self.limits.threshold("quality_sequence\terror", 20.0);
191        // Error if most frequent quality score <= threshold
192        self.calculate()
193            .is_some_and(|data| (data.most_frequent_score as f64) <= error_threshold)
194    }
195
196    fn raises_warning(&self) -> bool {
197        let warn_threshold = self.limits.threshold("quality_sequence\twarn", 27.0);
198        self.calculate()
199            .is_some_and(|data| (data.most_frequent_score as f64) <= warn_threshold)
200    }
201
202    fn ignore_filtered_sequences(&self) -> bool {
203        true
204    }
205
206    fn ignore_in_report(&self) -> bool {
207        self.limits.is_ignored("quality_sequence") || !self.has_data
208    }
209
210    fn write_text_report(&self, writer: &mut dyn io::Write) -> io::Result<()> {
211        let data = match self.calculate() {
212            Some(d) => d,
213            None => return Ok(()),
214        };
215
216        // Header format matches Java's makeReport
217        writeln!(writer, "#Quality\tCount")?;
218
219        for i in 0..data.x_categories.len() {
220            writeln!(
221                writer,
222                "{}\t{}",
223                data.x_categories[i],
224                java_format_double(data.quality_distribution[i]),
225            )?;
226        }
227
228        Ok(())
229    }
230
231    // Image filename matches Java's "per_sequence_quality.png" in Images/
232    fn chart_image_name(&self) -> Option<&str> {
233        Some("per_sequence_quality")
234    }
235    fn chart_alt_text(&self) -> Option<&str> {
236        Some("Per Sequence quality graph")
237    }
238    fn generate_chart_svg(&self) -> Option<String> {
239        self.build_chart_svg()
240    }
241}
242
243struct DistributionData {
244    quality_distribution: Vec<f64>,
245    x_categories: Vec<i32>,
246    most_frequent_score: i32,
247}
fastqc_rust/modules/per_sequence_quality.rs

fastqc_rust/modules/
per_sequence_quality.rs
