dnagent_domain/
feature_match.rs

1//! Sequence identity for the feature library: a feature's bases in its own direction,
2//! a strand-independent key, and exact occurrences on both strands (wrapping on circles).
3use crate::{Feature, SequenceRecord, Strand, Topology};
4
5/// Watson–Crick complement of an upper-case ACGT sequence, read 5′→3′; `None` for any other base.
6#[must_use]
7pub fn reverse_complement(sequence: &str) -> Option<String> {
8    sequence
9        .bytes()
10        .rev()
11        .map(|base| match base {
12            b'A' => Some('T'),
13            b'C' => Some('G'),
14            b'G' => Some('C'),
15            b'T' => Some('A'),
16            _ => None,
17        })
18        .collect()
19}
20
21/// The feature's bases 5′→3′ in its own direction: parts joined in source order, each
22/// wrapping through the origin if it does; reverse features are reverse-complemented.
23/// `None` if any base is not A, C, G or T.
24#[must_use]
25pub fn feature_sequence(record: &SequenceRecord, feature: &Feature) -> Option<String> {
26    let sequence = record.sequence().as_str();
27    let length = sequence.len();
28    let mut joined = String::new();
29    for part in feature.location().parts() {
30        let start = part.start().get();
31        let size = part.length().get();
32        let end = start + size;
33        if end <= length {
34            joined.push_str(&sequence[start..end]);
35        } else {
36            joined.push_str(&sequence[start..]);
37            joined.push_str(&sequence[..end - length]);
38        }
39    }
40    if !joined
41        .bytes()
42        .all(|b| matches!(b, b'A' | b'C' | b'G' | b'T'))
43    {
44        return None;
45    }
46    if feature.location().strand() == Strand::Reverse {
47        reverse_complement(&joined)
48    } else {
49        Some(joined)
50    }
51}
52
53/// Strand-independent identity: the lexicographically smaller of a sequence and its
54/// reverse complement. `None` for non-ACGT input.
55#[must_use]
56pub fn canonical(sequence: &str) -> Option<String> {
57    let reverse = reverse_complement(sequence)?;
58    Some(if reverse.as_str() < sequence {
59        reverse
60    } else {
61        sequence.to_owned()
62    })
63}
64
65/// One exact occurrence of a query on a molecule.
66#[derive(Debug, Clone, Copy, PartialEq, Eq)]
67pub struct Occurrence {
68    /// Start on the forward strand (zero-based); with `length`, may wrap on a circle.
69    pub start: usize,
70    pub length: usize,
71    /// `Forward`: the query reads 5′→3′ on the top strand; `Reverse`: on the bottom strand.
72    pub strand: Strand,
73}
74
75/// Every exact occurrence of `query` (ACGT, 5′→3′) on both strands, ordered by start then
76/// strand. A palindromic query is reported once, on the forward strand. On a circle,
77/// occurrences may run through the origin; a query longer than the molecule never matches.
78#[must_use]
79pub fn find_exact(sequence: &str, topology: Topology, query: &str) -> Vec<Occurrence> {
80    Searcher::new(sequence, topology, query.len()).find(query)
81}
82
83/// Repeated [`find_exact`] searches on one molecule, preparing the origin-wrapped text once.
84pub struct Searcher {
85    text: String,
86    length: usize,
87    circular: bool,
88}
89
90impl Searcher {
91    /// Prepare for queries up to `longest` bases.
92    #[must_use]
93    pub fn new(sequence: &str, topology: Topology, longest: usize) -> Self {
94        let circular = topology == Topology::Circular;
95        let length = sequence.len();
96        let wrap = if circular {
97            longest.saturating_sub(1).min(length)
98        } else {
99            0
100        };
101        Self {
102            text: format!("{sequence}{}", &sequence[..wrap]),
103            length,
104            circular,
105        }
106    }
107
108    /// See [`find_exact`]. Queries longer than the `longest` given to [`Searcher::new`]
109    /// miss origin-spanning occurrences, so size it for the longest query.
110    #[must_use]
111    pub fn find(&self, query: &str) -> Vec<Occurrence> {
112        let (n, m) = (self.length, query.len());
113        let Some(reverse) = reverse_complement(query) else {
114            return Vec::new();
115        };
116        if m == 0 || m > n {
117            return Vec::new();
118        }
119        let text = if self.circular {
120            &self.text[..(n + m - 1).min(self.text.len())]
121        } else {
122            &self.text[..n]
123        };
124        let mut found = Vec::new();
125        let mut scan = |needle: &str, strand: Strand| {
126            let mut from = 0;
127            while let Some(offset) = text[from..].find(needle) {
128                let start = from + offset;
129                if start >= n {
130                    break;
131                }
132                found.push(Occurrence {
133                    start,
134                    length: m,
135                    strand,
136                });
137                from = start + 1;
138            }
139        };
140        scan(query, Strand::Forward);
141        if reverse != query {
142            scan(&reverse, Strand::Reverse);
143        }
144        found.sort_by_key(|o| (o.start, o.strand == Strand::Reverse));
145        found
146    }
147}
148
149/// Fewest substitutions, insertions and deletions that align all of `short` to some
150/// stretch of `long` (end gaps in `long` are free), or `None` if more than `max`.
151#[must_use]
152pub fn semiglobal_edits(short: &[u8], long: &[u8], max: usize) -> Option<usize> {
153    let m = long.len();
154    let mut previous = vec![0usize; m + 1]; // aligning nothing costs nothing anywhere in `long`
155    let mut current = vec![0usize; m + 1];
156    for (i, &a) in short.iter().enumerate() {
157        current[0] = i + 1;
158        let mut best = current[0];
159        for j in 1..=m {
160            let substitution = previous[j - 1] + usize::from(a != long[j - 1]);
161            current[j] = substitution.min(previous[j] + 1).min(current[j - 1] + 1);
162            best = best.min(current[j]);
163        }
164        if best > max {
165            return None;
166        }
167        std::mem::swap(&mut previous, &mut current);
168    }
169    previous.iter().copied().min().filter(|&edits| edits <= max)
170}
171
172/// Standard-code translation of whole codons from the first base (stops as `*`).
173#[must_use]
174pub fn translate_standard(sequence: &str) -> String {
175    let code = crate::translation::GeneticCode::standard();
176    sequence
177        .as_bytes()
178        .as_chunks::<3>()
179        .0
180        .iter()
181        .map(|c| code.translate_codon(*c))
182        .collect()
183}
184
185#[cfg(test)]
186mod tests {
187    use super::*;
188
189    #[test]
190    fn both_strands_origin_and_palindromes() {
191        // GAATTC is palindromic: forward only. ACCG forward at 1; its reverse CGGT at 6.
192        let hits = find_exact("TACCGAACGGT", Topology::Linear, "ACCG");
193        assert_eq!(
194            hits,
195            vec![
196                Occurrence {
197                    start: 1,
198                    length: 4,
199                    strand: Strand::Forward
200                },
201                Occurrence {
202                    start: 7,
203                    length: 4,
204                    strand: Strand::Reverse
205                },
206            ]
207        );
208        assert_eq!(
209            find_exact("GAATTCGAATTC", Topology::Linear, "GAATTC").len(),
210            2
211        );
212        // Across the origin of a circle, not on a linear molecule.
213        assert_eq!(
214            find_exact("CGTTTTTAC", Topology::Circular, "ACCG"),
215            vec![Occurrence {
216                start: 7,
217                length: 4,
218                strand: Strand::Forward
219            }]
220        );
221        assert_eq!(find_exact("CGTTTTTAC", Topology::Linear, "ACCG"), []);
222        assert_eq!(find_exact("ACG", Topology::Circular, "ACGT"), []);
223    }
224
225    #[test]
226    fn semiglobal_edits_align_the_short_one_fully() {
227        assert_eq!(
228            semiglobal_edits(b"ACGT", b"TTACGTTT", 0),
229            Some(0),
230            "exact containment"
231        );
232        assert_eq!(
233            semiglobal_edits(b"ACCT", b"TTACGTTT", 1),
234            Some(1),
235            "one substitution"
236        );
237        assert_eq!(
238            semiglobal_edits(b"ACGGT", b"TTACGTTT", 1),
239            Some(1),
240            "one insertion in the short one"
241        );
242        assert_eq!(
243            semiglobal_edits(b"AGT", b"TTACGTTT", 1),
244            Some(1),
245            "one deletion"
246        );
247        assert_eq!(semiglobal_edits(b"GGGG", b"TTACGTTT", 2), None, "bounded");
248        assert_eq!(translate_standard("ATGAAATAGC"), "MK*");
249    }
250
251    #[test]
252    fn canonical_is_strand_independent() {
253        assert_eq!(canonical("TTGCA"), canonical("TGCAA"));
254        assert_eq!(canonical("ACGN"), None);
255        assert_eq!(reverse_complement("AACG").as_deref(), Some("CGTT"));
256    }
257}
dnagent_domain/feature_match.rs

dnagent_domain/
feature_match.rs
