fastqc_rust/modules/
per_base_sequence_content.rs

1// Per Base Sequence Content module
2// Corresponds to Modules/PerBaseSequenceContent.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::scaled_chart_width;
10use crate::sequence::Sequence;
11use crate::utils::base_counts::{BASE_INDEX, IDX_A, IDX_C, IDX_G, IDX_T};
12use crate::utils::base_group::BaseGroup;
13use crate::utils::format::java_format_double;
14
15pub struct PerBaseSequenceContent {
16    /// Per-position base counts stored as [A, C, G, T] per position.
17    /// Using a single Vec of arrays gives better cache locality than
18    /// four separate Vecs, since all counts for a position are adjacent.
19    counts: Vec<[u64; 4]>,
20    nogroup: bool,
21    expgroup: bool,
22    limits: Limits,
23}
24
25impl PerBaseSequenceContent {
26    pub fn new(limits: &Limits, nogroup: bool, expgroup: bool) -> Self {
27        PerBaseSequenceContent {
28            counts: Vec::new(),
29            nogroup,
30            expgroup,
31            limits: limits.clone(),
32        }
33    }
34
35    fn calculate(&self) -> ContentData {
36        let groups = BaseGroup::make_base_groups(self.counts.len(), self.nogroup, self.expgroup);
37
38        let mut x_categories = Vec::with_capacity(groups.len());
39        let mut g_percent = vec![0.0f64; groups.len()];
40        let mut a_percent = vec![0.0f64; groups.len()];
41        let mut t_percent = vec![0.0f64; groups.len()];
42        let mut c_percent = vec![0.0f64; groups.len()];
43
44        for (i, group) in groups.iter().enumerate() {
45            x_categories.push(group.label());
46
47            let mut a_count: u64 = 0;
48            let mut c_count: u64 = 0;
49            let mut g_count: u64 = 0;
50            let mut t_count: u64 = 0;
51            let mut total: u64 = 0;
52
53            // Java iterates `for (int bp=groups[i].lowerCount()-1;bp<groups[i].upperCount();bp++)`
54            // which is 0-based lowerCount-1 to upperCount-1 inclusive. Our lower_count/upper_count
55            // are already 0-based.
56            for bp in group.lower_count..=group.upper_count {
57                let c = &self.counts[bp];
58                a_count += c[IDX_A];
59                c_count += c[IDX_C];
60                g_count += c[IDX_G];
61                t_count += c[IDX_T];
62                total += c[IDX_A] + c[IDX_C] + c[IDX_G] + c[IDX_T];
63            }
64
65            g_percent[i] = (g_count as f64 / total as f64) * 100.0;
66            a_percent[i] = (a_count as f64 / total as f64) * 100.0;
67            t_percent[i] = (t_count as f64 / total as f64) * 100.0;
68            c_percent[i] = (c_count as f64 / total as f64) * 100.0;
69        }
70
71        // percentages stored in order [T, C, A, G] matching Java's array layout
72        ContentData {
73            x_categories,
74            t_percent,
75            c_percent,
76            a_percent,
77            g_percent,
78        }
79    }
80}
81
82impl PerBaseSequenceContent {
83    fn build_chart_svg(&self) -> String {
84        let data = self.calculate();
85
86        // Series order in LineGraph is [%T, %C, %A, %G], matching Java's
87        // `new LineGraph(percentages, 0d, 100d, ..., new String[] {"%T","%C","%A","%G"}, ...)`
88        render_line_graph(&LineGraphData {
89            width: scaled_chart_width(data.x_categories.len()),
90            data: vec![
91                data.t_percent,
92                data.c_percent,
93                data.a_percent,
94                data.g_percent,
95            ],
96            min_y: 0.0,
97            max_y: 100.0,
98            x_label: "Position in read (bp)".to_string(),
99            series_names: vec![
100                "%T".to_string(),
101                "%C".to_string(),
102                "%A".to_string(),
103                "%G".to_string(),
104            ],
105            x_categories: data.x_categories,
106            title: "Sequence content across all bases".to_string(),
107        })
108    }
109}
110
111impl QCModule for PerBaseSequenceContent {
112    fn cost_hint(&self) -> u32 {
113        9
114    }
115
116    fn process_sequence(&mut self, sequence: &Sequence) {
117        let seq = &sequence.sequence;
118
119        // Grow array if needed
120        if self.counts.len() < seq.len() {
121            self.counts.resize(seq.len(), [0; 4]);
122        }
123
124        // Use lookup table for branchless base classification.
125        // Only A/C/G/T (indices 0-3) are counted; N and others (indices 4-5)
126        // are filtered by the single comparison `idx < 4`.
127        for (i, &b) in seq.iter().enumerate() {
128            let idx = BASE_INDEX[b as usize] as usize;
129            if idx < 4 {
130                self.counts[i][idx] += 1;
131            }
132        }
133    }
134
135    fn name(&self) -> &str {
136        "Per base sequence content"
137    }
138
139    fn description(&self) -> &str {
140        "Shows the relative amounts of each base at each position in a sequencing run"
141    }
142
143    fn reset(&mut self) {
144        self.counts.clear();
145    }
146
147    fn raises_error(&self) -> bool {
148        let error_threshold = self.limits.threshold("sequence\terror", 20.0);
149        let data = self.calculate();
150
151        // Check GC diff (C vs G) and AT diff (T vs A) per group
152        for i in 0..data.g_percent.len() {
153            let gc_diff = (data.c_percent[i] - data.g_percent[i]).abs();
154            let at_diff = (data.t_percent[i] - data.a_percent[i]).abs();
155
156            if gc_diff > error_threshold || at_diff > error_threshold {
157                return true;
158            }
159        }
160        false
161    }
162
163    fn raises_warning(&self) -> bool {
164        let warn_threshold = self.limits.threshold("sequence\twarn", 10.0);
165        let data = self.calculate();
166
167        for i in 0..data.g_percent.len() {
168            let gc_diff = (data.c_percent[i] - data.g_percent[i]).abs();
169            let at_diff = (data.t_percent[i] - data.a_percent[i]).abs();
170
171            if gc_diff > warn_threshold || at_diff > warn_threshold {
172                return true;
173            }
174        }
175        false
176    }
177
178    fn ignore_filtered_sequences(&self) -> bool {
179        true
180    }
181
182    fn ignore_in_report(&self) -> bool {
183        self.limits.is_ignored("sequence")
184    }
185
186    fn write_text_report(&self, writer: &mut dyn io::Write) -> io::Result<()> {
187        let data = self.calculate();
188
189        // Column order is G, A, T, C matching Java's makeReport
190        writeln!(writer, "#Base\tG\tA\tT\tC")?;
191
192        for i in 0..data.x_categories.len() {
193            // Java outputs percentages[3][i] (G), [2][i] (A), [0][i] (T), [1][i] (C)
194            writeln!(
195                writer,
196                "{}\t{}\t{}\t{}\t{}",
197                data.x_categories[i],
198                java_format_double(data.g_percent[i]),
199                java_format_double(data.a_percent[i]),
200                java_format_double(data.t_percent[i]),
201                java_format_double(data.c_percent[i]),
202            )?;
203        }
204
205        Ok(())
206    }
207
208    // Image filename matches Java's "per_base_sequence_content.png" in Images/
209    fn chart_image_name(&self) -> Option<&str> {
210        Some("per_base_sequence_content")
211    }
212    fn chart_alt_text(&self) -> Option<&str> {
213        Some("Per base sequence content")
214    }
215    fn generate_chart_svg(&self) -> Option<String> {
216        Some(self.build_chart_svg())
217    }
218}
219
220struct ContentData {
221    x_categories: Vec<String>,
222    t_percent: Vec<f64>,
223    c_percent: Vec<f64>,
224    a_percent: Vec<f64>,
225    g_percent: Vec<f64>,
226}
fastqc_rust/modules/per_base_sequence_content.rs

fastqc_rust/modules/
per_base_sequence_content.rs
