fastqc_rust/modules/
adapter_content.rs

1// Adapter Content module
2// Corresponds to Modules/AdapterContent.java
3
4use std::io;
5
6use aho_corasick::{packed, Span};
7use memchr::memmem;
8
9use crate::config::{Limits, LimitsExt};
10use crate::modules::QCModule;
11use crate::report::charts::line_graph::{render_line_graph, LineGraphData};
12use crate::report::charts::scaled_chart_width;
13use crate::sequence::Sequence;
14use crate::utils::base_group::BaseGroup;
15use crate::utils::format::java_format_double;
16
17/// A single adapter to search for in sequences.
18struct Adapter {
19    name: String,
20    /// Reads whose first match starts at each position; `cumulative_positions`
21    /// gives Java's running totals.
22    first_hits: Vec<u64>,
23}
24
25impl Adapter {
26    fn new(name: &str) -> Self {
27        Adapter {
28            name: name.to_string(),
29            first_hits: vec![0; 1],
30        }
31    }
32
33    fn record_hit(&mut self, position: usize) {
34        // Like Java, a match past the currently tracked positions is dropped,
35        // not counted once a longer read extends them.
36        if let Some(count) = self.first_hits.get_mut(position) {
37            *count += 1;
38        }
39    }
40
41    /// Java's `positions`: reads with a match at or before each position.
42    fn cumulative_positions(&self) -> Vec<u64> {
43        self.first_hits
44            .iter()
45            .scan(0u64, |total, &hits| {
46                *total += hits;
47                Some(*total)
48            })
49            .collect()
50    }
51}
52
53/// Finds each adapter's first match in a read. One SIMD multi-pattern pass
54/// (aho-corasick's Teddy) where possible, else a substring search per adapter.
55struct AdapterSearch {
56    finders: Vec<memmem::Finder<'static>>,
57    teddy: Option<packed::Searcher>,
58}
59
60impl AdapterSearch {
61    fn new(adapters: &[(String, String)]) -> Self {
62        let seqs: Vec<&[u8]> = adapters.iter().map(|(_, s)| s.as_bytes()).collect();
63        AdapterSearch {
64            finders: seqs
65                .iter()
66                .map(|s| memmem::Finder::new(s).into_owned())
67                .collect(),
68            teddy: Self::build_teddy(&seqs),
69        }
70    }
71
72    /// `None` if Teddy can't be built for these patterns, or if one adapter is
73    /// a prefix of another: leftmost-first reports one match per start position.
74    fn build_teddy(seqs: &[&[u8]]) -> Option<packed::Searcher> {
75        if seqs.is_empty() || seqs.len() > 64 || seqs.iter().any(|s| s.is_empty()) {
76            return None;
77        }
78        for (i, a) in seqs.iter().enumerate() {
79            for (j, b) in seqs.iter().enumerate() {
80                if i != j && b.starts_with(a) {
81                    return None;
82                }
83            }
84        }
85        packed::Config::new()
86            .match_kind(packed::MatchKind::LeftmostFirst)
87            .builder()
88            .extend(seqs)
89            .build()
90    }
91
92    /// Calls `on_hit(adapter_index, start)` once for the leftmost match of
93    /// each adapter in `seq`.
94    fn first_matches(&self, seq: &[u8], mut on_hit: impl FnMut(usize, usize)) {
95        let Some(teddy) = &self.teddy else {
96            for (a, finder) in self.finders.iter().enumerate() {
97                if let Some(start) = finder.find(seq) {
98                    on_hit(a, start);
99                }
100            }
101            return;
102        };
103
104        let mut remaining: u64 = u64::MAX >> (64 - self.finders.len());
105        let mut at = 0;
106        while remaining != 0 && at < seq.len() {
107            let Some(m) = teddy.find_in(seq, Span::from(at..seq.len())) else {
108                return;
109            };
110            let a = m.pattern().as_usize();
111            if remaining & (1 << a) == 0 {
112                // A repeat (typically inside a polyA/polyG run) would otherwise
113                // cost one restart per base. No adapter is a prefix of another,
114                // so nothing unseen starts before here: search the rest singly.
115                break;
116            }
117            remaining &= !(1 << a);
118            on_hit(a, m.start());
119            // One past the start, not the end, so overlapping matches of other
120            // adapters aren't skipped.
121            at = m.start() + 1;
122        }
123        let rest = &seq[at.min(seq.len())..];
124        while remaining != 0 {
125            let a = remaining.trailing_zeros() as usize;
126            remaining &= remaining - 1;
127            if let Some(start) = self.finders[a].find(rest) {
128                on_hit(a, at + start);
129            }
130        }
131    }
132}
133
134pub struct AdapterContent {
135    adapters: Vec<Adapter>,
136    search: AdapterSearch,
137    longest_sequence: usize,
138    longest_adapter: usize,
139    total_count: u64,
140    limits: Limits,
141    nogroup: bool,
142    expgroup: bool,
143    // Lazily computed
144    computed: Option<ComputedEnrichment>,
145}
146
147struct ComputedEnrichment {
148    enrichments: Vec<Vec<f64>>,
149    x_labels: Vec<String>,
150}
151
152impl AdapterContent {
153    pub fn new(
154        limits: &Limits,
155        adapter_entries: &[(String, String)],
156        nogroup: bool,
157        expgroup: bool,
158    ) -> Self {
159        let mut longest_adapter = 0;
160        let mut adapters = Vec::with_capacity(adapter_entries.len());
161
162        for (name, seq) in adapter_entries {
163            if seq.len() > longest_adapter {
164                longest_adapter = seq.len();
165            }
166            adapters.push(Adapter::new(name));
167        }
168
169        if adapter_entries
170            .iter()
171            .any(|(_, seq)| seq.len() != longest_adapter)
172        {
173            static WARNED: crate::progress::OncePerRun = crate::progress::OncePerRun::new();
174            WARNED.log(|| "[Warning] You are using adapter sequences with different lengths. Matches will only be reported up to the position where the longest adapter could match. Matches to shorter adapters at the end of sequences will not be recorded.".to_string());
175        }
176
177        AdapterContent {
178            adapters,
179            search: AdapterSearch::new(adapter_entries),
180            longest_sequence: 0,
181            longest_adapter,
182            total_count: 0,
183            limits: limits.clone(),
184            nogroup,
185            expgroup,
186            computed: None,
187        }
188    }
189
190    /// Replicates calculateEnrichment() from AdapterContent.java.
191    fn calculate_enrichment(&mut self) {
192        if self.computed.is_some() {
193            return;
194        }
195
196        let all_positions: Vec<Vec<u64>> = self
197            .adapters
198            .iter()
199            .map(Adapter::cumulative_positions)
200            .collect();
201        let max_length = all_positions.iter().map(Vec::len).max().unwrap_or(0);
202
203        // Group positions using BaseGroup
204        let groups = BaseGroup::make_base_groups(max_length, self.nogroup, self.expgroup);
205
206        let x_labels: Vec<String> = groups.iter().map(|g| g.label()).collect();
207
208        let mut enrichments = vec![vec![0.0f64; groups.len()]; self.adapters.len()];
209
210        for (a, positions) in all_positions.iter().enumerate() {
211            for (g, group) in groups.iter().enumerate() {
212                // lowerCount() is 1-based in Java, we use 0-based internally
213                // Java: p=groups[g].lowerCount()-1; p<groups[g].upperCount()
214                let lower = group.lower_count; // already 0-based
215                let upper = group.upper_count; // already 0-based, inclusive
216
217                for p in lower..=upper {
218                    if p < positions.len() {
219                        enrichments[a][g] +=
220                            (positions[p] as f64 * 100.0) / self.total_count as f64;
221                    }
222                }
223
224                // Average over the group width
225                enrichments[a][g] /= (upper - lower + 1) as f64;
226            }
227        }
228
229        self.computed = Some(ComputedEnrichment {
230            enrichments,
231            x_labels,
232        });
233    }
234
235    /// No read was longer than the longest adapter, so Java skips the analysis
236    /// (no table or chart) and just warns.
237    fn reads_too_short(&self) -> bool {
238        self.longest_adapter > self.longest_sequence
239    }
240
241    /// Derive adapter names from the adapters Vec (avoids storing a redundant copy).
242    fn adapter_names(&self) -> Vec<String> {
243        self.adapters.iter().map(|a| a.name.clone()).collect()
244    }
245
246    fn ensure_calculated(&self) -> &ComputedEnrichment {
247        static DEFAULT: ComputedEnrichment = ComputedEnrichment {
248            enrichments: Vec::new(),
249            x_labels: Vec::new(),
250        };
251        self.computed.as_ref().unwrap_or(&DEFAULT)
252    }
253}
254
255impl AdapterContent {
256    fn build_chart_svg(&self) -> String {
257        let computed = self.ensure_calculated();
258
259        // Matches Java's `new LineGraph(enrichments, 0, 100, "Position in read (bp)", labels, xLabels, "% Adapter")`
260        render_line_graph(&LineGraphData {
261            width: scaled_chart_width(computed.x_labels.len()),
262            data: computed.enrichments.clone(),
263            min_y: 0.0,
264            max_y: 100.0,
265            x_label: "Position in read (bp)".to_string(),
266            series_names: self.adapter_names(),
267            x_categories: computed.x_labels.clone(),
268            title: "% Adapter".to_string(),
269        })
270    }
271}
272
273impl QCModule for AdapterContent {
274    fn cost_hint(&self) -> u32 {
275        9
276    }
277
278    fn process_sequence(&mut self, sequence: &Sequence) {
279        self.computed = None;
280        self.total_count += 1;
281
282        let seq_len = sequence.sequence.len();
283
284        // Java only grows the arrays once a read is longer than the longest adapter.
285        if seq_len > self.longest_sequence && seq_len > self.longest_adapter {
286            self.longest_sequence = seq_len;
287            let new_len = (self.longest_sequence - self.longest_adapter) + 1;
288            for adapter in &mut self.adapters {
289                adapter.first_hits.resize(new_len, 0);
290            }
291        }
292
293        let adapters = &mut self.adapters;
294        self.search
295            .first_matches(&sequence.sequence, |a, start| adapters[a].record_hit(start));
296    }
297
298    fn finalize(&mut self) {
299        self.calculate_enrichment();
300    }
301
302    fn name(&self) -> &str {
303        "Adapter Content"
304    }
305
306    fn description(&self) -> &str {
307        "Searches for specific adapter sequences in a library"
308    }
309
310    fn reset(&mut self) {
311        self.total_count = 0;
312        self.longest_sequence = 0;
313        self.computed = None;
314        for adapter in &mut self.adapters {
315            adapter.first_hits = vec![0; 1];
316        }
317    }
318
319    fn raises_error(&self) -> bool {
320        let threshold = self.limits.threshold("adapter\terror", 10.0);
321        let computed = self.ensure_calculated();
322        computed
323            .enrichments
324            .iter()
325            .any(|enrichments| enrichments.iter().any(|&val| val > threshold))
326    }
327
328    fn raises_warning(&self) -> bool {
329        if self.reads_too_short() {
330            return true;
331        }
332
333        let threshold = self.limits.threshold("adapter\twarn", 5.0);
334        let computed = self.ensure_calculated();
335        computed
336            .enrichments
337            .iter()
338            .any(|enrichments| enrichments.iter().any(|&val| val > threshold))
339    }
340
341    fn ignore_filtered_sequences(&self) -> bool {
342        true
343    }
344
345    fn ignore_in_report(&self) -> bool {
346        self.limits.is_ignored("adapter")
347    }
348
349    fn write_html_report(&self, writer: &mut dyn io::Write, png: bool) -> io::Result<()> {
350        if self.reads_too_short() {
351            return write!(
352                writer,
353                "<p>Can't analyse adapters as read length is too short ({} vs {})</p>",
354                self.longest_adapter, self.longest_sequence
355            );
356        }
357        crate::report::html::write_chart(self, "Adapter graph", png, writer)
358    }
359
360    fn write_text_report(&self, writer: &mut dyn io::Write) -> io::Result<()> {
361        let computed = self.ensure_calculated();
362
363        if self.reads_too_short() {
364            return Ok(());
365        }
366
367        // Header line with Position tab and all adapter names
368        write!(writer, "#Position")?;
369        for adapter in &self.adapters {
370            write!(writer, "\t{}", adapter.name)?;
371        }
372        writeln!(writer)?;
373
374        // One row per base group, columns are Position then each adapter's enrichment
375        for (row, x_label) in computed.x_labels.iter().enumerate() {
376            write!(writer, "{}", x_label)?;
377            for a in 0..self.adapters.len() {
378                write!(
379                    writer,
380                    "\t{}",
381                    java_format_double(computed.enrichments[a][row])
382                )?;
383            }
384            writeln!(writer)?;
385        }
386
387        Ok(())
388    }
389
390    // Image filename matches Java's "adapter_content.png" in Images/
391    fn chart_image_name(&self) -> Option<&str> {
392        Some("adapter_content")
393    }
394    fn chart_alt_text(&self) -> Option<&str> {
395        Some("Adapter graph")
396    }
397    fn generate_chart_svg(&self) -> Option<String> {
398        if self.reads_too_short() {
399            return None;
400        }
401        Some(self.build_chart_svg())
402    }
403}
404
405#[cfg(test)]
406mod tests {
407    use super::*;
408
409    const DEFAULT_ADAPTERS: [&str; 6] = [
410        "AGATCGGAAGAG",
411        "TGGAATTCTCGG",
412        "GATCGTCGGACT",
413        "CTGTCTCTTATA",
414        "AAAAAAAAAAAA",
415        "GGGGGGGGGGGG",
416    ];
417
418    fn entries(seqs: &[&str]) -> Vec<(String, String)> {
419        seqs.iter()
420            .map(|s| (s.to_string(), s.to_string()))
421            .collect()
422    }
423
424    fn memmem_first(adapters: &[(String, String)], seq: &[u8]) -> Vec<Option<usize>> {
425        adapters
426            .iter()
427            .map(|(_, a)| memmem::find(seq, a.as_bytes()))
428            .collect()
429    }
430
431    fn search_first(search: &AdapterSearch, seq: &[u8]) -> Vec<Option<usize>> {
432        let mut found = vec![None; search.finders.len()];
433        search.first_matches(seq, |a, start| {
434            assert!(found[a].is_none(), "adapter {a} reported twice");
435            found[a] = Some(start);
436        });
437        found
438    }
439
440    /// The multi-pattern search must report exactly the first match that a
441    /// separate substring search per adapter would.
442    #[test]
443    fn test_adapter_search_matches_memmem() {
444        let adapters = entries(&DEFAULT_ADAPTERS);
445        let search = AdapterSearch::new(&adapters);
446        assert!(search.teddy.is_some(), "default adapters use Teddy");
447
448        let poly_g = vec![b'G'; 150];
449        let mut poly_a_tail = b"ACGT".repeat(500);
450        poly_a_tail.extend_from_slice(&[b'A'; 5000]);
451        poly_a_tail.extend_from_slice(b"AGATCGGAAGAG");
452        let mut seqs: Vec<Vec<u8>> = vec![
453            b"".to_vec(),
454            b"AGATCGGAAGA".to_vec(),
455            b"NNNAGATCGGAAGAGNNN".to_vec(),
456            // Adapter 0 twice, overlapping itself.
457            b"CCAGATCGGAAGAGATCGGAAGAGTT".to_vec(),
458            // PolyA overlaps the last base of adapter 3; restarting at a
459            // match's end would miss it.
460            b"CTGTCTCTTATAAAAAAAAAAA".to_vec(),
461            // Other adapters inside and after homopolymer runs.
462            b"GGGGGGGGGGGGGGGGGGGGAGATCGGAAGAGGGGGGGGGGGGGGGAAAAAAAAAAAAAA".to_vec(),
463            b"AAAAAAAAAAAAAAAAAAAAAAAAACTGTCTCTTATAAAAAAAAAAAAAAAATGGAATTCTCGG".to_vec(),
464            poly_g,
465            poly_a_tail,
466            b"GGGGGGGGGGGGAAAAAAAAAAAACTGTCTCTTATAGATCGTCGGACTTGGAATTCTCGGAGATCGGAAGAG".to_vec(),
467        ];
468        // Deterministic pseudo-random reads biased towards adapter fragments
469        // and homopolymer runs.
470        let mut state = 0x2545_f491_4f6c_dd1du64;
471        let mut next = || {
472            state ^= state << 13;
473            state ^= state >> 7;
474            state ^= state << 17;
475            state
476        };
477        for _ in 0..5000 {
478            let len = (next() % 300) as usize;
479            let mut seq = Vec::with_capacity(len + 40);
480            while seq.len() < len {
481                match next() % 16 {
482                    0 | 1 => seq.extend_from_slice(adapters[(next() % 6) as usize].1.as_bytes()),
483                    2 => {
484                        let run = (next() % 40) as usize;
485                        seq.extend(std::iter::repeat_n(b"AG"[(next() % 2) as usize], run));
486                    }
487                    _ => seq.push(b"ACGTN"[(next() % 5) as usize]),
488                }
489            }
490            seqs.push(seq);
491        }
492
493        for seq in &seqs {
494            assert_eq!(
495                search_first(&search, seq),
496                memmem_first(&adapters, seq),
497                "{}",
498                String::from_utf8_lossy(seq)
499            );
500        }
501    }
502
503    /// Leftmost-first cannot report two adapters starting at the same
504    /// position, so a prefix pair must fall back to per-adapter search.
505    #[test]
506    fn test_adapter_search_prefix_fallback() {
507        let prefix = entries(&["ACGTACGT", "ACGTACGTAA"]);
508        let search = AdapterSearch::new(&prefix);
509        assert!(search.teddy.is_none());
510        assert_eq!(
511            search_first(&search, b"TTACGTACGTAA"),
512            vec![Some(2), Some(2)]
513        );
514
515        assert!(AdapterSearch::new(&entries(&["ACGTACGT", "ACGTACGT"]))
516            .teddy
517            .is_none());
518        assert!(AdapterSearch::new(&entries(&["ACGTACGT", "CGTACGTA"]))
519            .teddy
520            .is_some());
521    }
522
523    /// Hits are tallied per first-match position and summed at the end; the
524    /// cumulative curve must equal incrementing every later position per hit.
525    #[test]
526    fn test_cumulative_positions() {
527        let mut adapter = Adapter::new("x");
528        adapter.first_hits = vec![0; 5];
529        for p in [0, 2, 2, 4, 9] {
530            adapter.record_hit(p);
531        }
532        assert_eq!(adapter.cumulative_positions(), vec![1, 1, 3, 3, 4]);
533    }
534}
fastqc_rust/modules/adapter_content.rs

fastqc_rust/modules/
adapter_content.rs
