dnagent_domain/
digest.rs

1//! Conservative complete restriction digests with explicit strand and end geometry.
2
3use crate::restriction::{self, Enzyme, OverhangPolarity, RestrictionError};
4use crate::{DnaSeq, Strand, Topology};
5use serde::Serialize;
6use thiserror::Error;
7
8#[derive(Debug, Clone, PartialEq, Eq, Error)]
9pub enum DigestError {
10    #[error(transparent)]
11    Restriction(#[from] RestrictionError),
12    #[error(
13        "complete digest refused: {enzyme} site at {start} has a cut outside the linear molecule; no partial products returned (use sites to inspect uncleavable sites)"
14    )]
15    OutOfBounds { enzyme: &'static str, start: usize },
16    #[error(
17        "complete digest refused: {enzyme} cuts at a linear terminus; terminal single-stranded products are not modelled and no partial products are returned"
18    )]
19    TerminalCut { enzyme: &'static str },
20    #[error(
21        "cannot digest overlapping/touching cleavage regions; products without a paired core are not modelled"
22    )]
23    OverlappingCuts,
24}
25
26#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
27pub struct DigestCut {
28    pub top_cut: usize,
29    pub bottom_cut: usize,
30    /// Enzymes producing the same physical cut pair are grouped, not cut twice.
31    pub enzymes: Vec<&'static str>,
32}
33
34#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
35pub struct FragmentEnd {
36    /// True only for an original linear end; input ends are assumed blunt.
37    pub original_terminus: bool,
38    pub polarity: OverhangPolarity,
39    pub protruding_strand: Option<Strand>,
40    /// Exposed oligo written 5′ to 3′ on its own strand; empty for blunt ends.
41    pub overhang_sequence: String,
42    pub enzymes: Vec<&'static str>,
43}
44
45#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
46pub struct FragmentStrand {
47    /// Start of the source interval along the input FORWARD coordinate axis.
48    /// For the bottom strand this is not the oligo's 5′ boundary.
49    pub source_start: usize,
50    pub length: usize,
51    pub sequence_5to3: String,
52}
53
54#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
55pub struct DigestFragment {
56    pub id: String,
57    pub topology: Topology,
58    pub top: FragmentStrand,
59    pub bottom: FragmentStrand,
60    pub paired_length: usize,
61    /// Null only for an uncut circular molecule.
62    pub left_end: Option<FragmentEnd>,
63    pub right_end: Option<FragmentEnd>,
64}
65
66#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
67pub struct Digest {
68    pub input_length: usize,
69    pub input_topology: Topology,
70    pub enzymes: Vec<Enzyme>,
71    pub cuts: Vec<DigestCut>,
72    pub fragments: Vec<DigestFragment>,
73    pub assumptions: Vec<&'static str>,
74}
75
76/// Internal unwrapped cut boundaries. Bottom may be negative on a circle.
77#[derive(Debug, Clone)]
78struct Cut {
79    top: i128,
80    bottom: i128,
81    enzymes: Vec<&'static str>,
82}
83
84/// Simulate complete cleavage. Ambiguous input, unavailable/terminal cuts and
85/// overlapping staggered regions are rejected rather than approximated.
86pub fn simulate_digest(
87    sequence: &DnaSeq,
88    topology: Topology,
89    names: &[String],
90) -> Result<Digest, DigestError> {
91    let scan = restriction::find_sites(sequence, topology, names)?;
92    let n = sequence.len();
93    let mut cuts: Vec<Cut> = Vec::new();
94    for site in &scan.sites {
95        let (Some(top), Some(bottom)) = (site.top_cut, site.bottom_cut) else {
96            return Err(DigestError::OutOfBounds {
97                enzyme: site.enzyme,
98                start: site.recognition.start().get(),
99            });
100        };
101        if topology == Topology::Linear && (top == 0 || top == n || bottom == 0 || bottom == n) {
102            return Err(DigestError::TerminalCut {
103                enzyme: site.enzyme,
104            });
105        }
106        // Deliberately do NOT reuse site.bottom_cut below: circular site cuts
107        // are modulo n, whereas partition needs a locally unwrapped cut pair.
108        // Reconstructing bottom = top + signed stagger preserves a 4-base end
109        // across the origin (including negative bottom coordinates for 3′ ends).
110        // Substituting the wrapped bottom would corrupt adjacent fragment windows.
111        let delta = match site.overhang_polarity {
112            OverhangPolarity::Blunt => 0,
113            OverhangPolarity::FivePrime => site.overhang_length as i128,
114            OverhangPolarity::ThreePrime => -(site.overhang_length as i128),
115        };
116        cuts.push(Cut {
117            top: top as i128,
118            bottom: top as i128 + delta,
119            enzymes: vec![site.enzyme],
120        });
121    }
122    let cuts = merge_cuts(cuts);
123    let fragments = if cuts.is_empty() {
124        vec![uncut_fragment(sequence, topology)]
125    } else {
126        partition(sequence, topology, &cuts)?
127    };
128    let normalise = |p: i128| {
129        if topology == Topology::Circular {
130            p.rem_euclid(n as i128) as usize
131        } else {
132            p as usize
133        }
134    };
135    Ok(Digest {
136        input_length: n,
137        input_topology: topology,
138        enzymes: scan.enzymes,
139        cuts: cuts
140            .into_iter()
141            .map(|cut| DigestCut {
142                top_cut: normalise(cut.top),
143                bottom_cut: normalise(cut.bottom),
144                enzymes: cut.enzymes,
145            })
146            .collect(),
147        fragments,
148        assumptions: vec![
149            "Complete cleavage at all selected recognition sites; double-stranded DNA assumed",
150            "Sequence-only prediction: methylation, star activity, reaction conditions, accessibility and flanking-DNA requirements not modelled",
151            "Original linear ends assumed blunt; terminal and overlapping cleavage regions rejected",
152            "Top and bottom sequences are both written 5-prime to 3-prime; source_start uses the stored forward coordinate axis",
153            "Fragments ordered along the stored top strand, not by size; no gel or experimental validation",
154            "Sequence-only fragments: feature annotations, primers and retained metadata are not propagated; source record is unchanged",
155        ],
156    })
157}
158
159fn merge_cuts(mut cuts: Vec<Cut>) -> Vec<Cut> {
160    cuts.sort_by_key(|c| (c.top, c.bottom));
161    let mut unique: Vec<Cut> = Vec::new();
162    for cut in cuts {
163        if let Some(previous) = unique
164            .last_mut()
165            .filter(|c| c.top == cut.top && c.bottom == cut.bottom)
166        {
167            previous.enzymes.extend(cut.enzymes);
168            previous.enzymes.sort_unstable();
169            previous.enzymes.dedup();
170        } else {
171            unique.push(cut);
172        }
173    }
174    unique
175}
176
177fn partition(
178    sequence: &DnaSeq,
179    topology: Topology,
180    cuts: &[Cut],
181) -> Result<Vec<DigestFragment>, DigestError> {
182    let n = sequence.len() as i128;
183    let mut boundaries = cuts.to_vec();
184    if topology == Topology::Linear {
185        boundaries.insert(
186            0,
187            Cut {
188                top: 0,
189                bottom: 0,
190                enzymes: vec![],
191            },
192        );
193        boundaries.push(Cut {
194            top: n,
195            bottom: n,
196            enzymes: vec![],
197        });
198    } else {
199        let first = &cuts[0];
200        boundaries.push(Cut {
201            top: first.top + n,
202            bottom: first.bottom + n,
203            enzymes: first.enzymes.clone(),
204        });
205    }
206    boundaries
207        .windows(2)
208        .enumerate()
209        .map(|(index, pair)| {
210            let (left, right) = (&pair[0], &pair[1]);
211            let paired_start = left.top.max(left.bottom);
212            let paired_end = right.top.min(right.bottom);
213            if paired_start >= paired_end {
214                return Err(DigestError::OverlappingCuts);
215            }
216            Ok(DigestFragment {
217                id: format!("fragment-{:04}", index + 1),
218                topology: Topology::Linear,
219                top: strand_interval(sequence, left.top, right.top, false),
220                bottom: strand_interval(sequence, left.bottom, right.bottom, true),
221                paired_length: (paired_end - paired_start) as usize,
222                left_end: Some(fragment_end(sequence, left, true)),
223                right_end: Some(fragment_end(sequence, right, false)),
224            })
225        })
226        .collect()
227}
228
229fn strand_interval(sequence: &DnaSeq, start: i128, end: i128, reverse: bool) -> FragmentStrand {
230    let bases = sequence.as_str().as_bytes();
231    let n = bases.len() as i128;
232    let interval: String = (start..end)
233        .map(|p| char::from(bases[p.rem_euclid(n) as usize]))
234        .collect();
235    FragmentStrand {
236        source_start: start.rem_euclid(n) as usize,
237        length: interval.len(),
238        sequence_5to3: if reverse {
239            reverse_complement(&interval)
240        } else {
241            interval
242        },
243    }
244}
245
246/// Internal ACGT-only helper for validated digest/ligation strand sequences.
247#[must_use]
248pub fn reverse_complement(sequence: &str) -> String {
249    sequence
250        .bytes()
251        .rev()
252        .map(|base| match base {
253            b'A' => 'T',
254            b'C' => 'G',
255            b'G' => 'C',
256            b'T' => 'A',
257            _ => unreachable!("restriction scan validated ACGT"),
258        })
259        .collect()
260}
261
262fn fragment_end(sequence: &DnaSeq, cut: &Cut, left: bool) -> FragmentEnd {
263    let delta = cut.bottom - cut.top;
264    let polarity = match delta.cmp(&0) {
265        std::cmp::Ordering::Equal => OverhangPolarity::Blunt,
266        std::cmp::Ordering::Greater => OverhangPolarity::FivePrime,
267        std::cmp::Ordering::Less => OverhangPolarity::ThreePrime,
268    };
269    let protruding_strand = if delta == 0 {
270        None
271    } else if (delta > 0) == left {
272        Some(Strand::Forward)
273    } else {
274        Some(Strand::Reverse)
275    };
276    let overhang_sequence = strand_interval(
277        sequence,
278        cut.top.min(cut.bottom),
279        cut.top.max(cut.bottom),
280        protruding_strand == Some(Strand::Reverse),
281    )
282    .sequence_5to3;
283    FragmentEnd {
284        original_terminus: cut.enzymes.is_empty(),
285        polarity,
286        protruding_strand,
287        overhang_sequence,
288        enzymes: cut.enzymes.clone(),
289    }
290}
291
292fn uncut_fragment(sequence: &DnaSeq, topology: Topology) -> DigestFragment {
293    let natural = FragmentEnd {
294        original_terminus: true,
295        polarity: OverhangPolarity::Blunt,
296        protruding_strand: None,
297        overhang_sequence: String::new(),
298        enzymes: vec![],
299    };
300    DigestFragment {
301        id: "fragment-0001".into(),
302        topology,
303        top: strand_interval(sequence, 0, sequence.len() as i128, false),
304        bottom: strand_interval(sequence, 0, sequence.len() as i128, true),
305        paired_length: sequence.len(),
306        left_end: (topology == Topology::Linear).then(|| natural.clone()),
307        right_end: (topology == Topology::Linear).then_some(natural),
308    }
309}
310
311#[cfg(test)]
312mod tests {
313    use super::*;
314    fn digest(seq: &str, topology: Topology, names: &[&str]) -> Result<Digest, DigestError> {
315        simulate_digest(
316            &DnaSeq::new(seq).unwrap(),
317            topology,
318            &names.iter().map(|n| (*n).to_owned()).collect::<Vec<_>>(),
319        )
320    }
321    fn conserved(result: &Digest) {
322        assert_eq!(
323            result.fragments.iter().map(|f| f.top.length).sum::<usize>(),
324            result.input_length
325        );
326        assert_eq!(
327            result
328                .fragments
329                .iter()
330                .map(|f| f.bottom.length)
331                .sum::<usize>(),
332            result.input_length
333        );
334        for f in &result.fragments {
335            assert_eq!(f.top.length, f.top.sequence_5to3.len());
336            assert_eq!(f.bottom.length, f.bottom.sequence_5to3.len());
337            assert!(f.paired_length > 0);
338        }
339    }
340    #[test]
341    fn ecori_linear_has_asymmetric_strands_and_complementary_ends() {
342        let d = digest("AAAGAATTCTTT", Topology::Linear, &["EcoRI"]).unwrap();
343        conserved(&d);
344        assert_eq!(d.fragments.len(), 2);
345        assert_eq!(d.fragments[0].top.sequence_5to3, "AAAG");
346        assert_eq!(d.fragments[0].bottom.sequence_5to3, "AATTCTTT");
347        assert_eq!(d.fragments[1].top.sequence_5to3, "AATTCTTT");
348        assert_eq!(d.fragments[1].bottom.sequence_5to3, "AAAG");
349        let right = d.fragments[0].right_end.as_ref().unwrap();
350        let left = d.fragments[1].left_end.as_ref().unwrap();
351        assert_eq!(right.protruding_strand, Some(Strand::Reverse));
352        assert_eq!(left.protruding_strand, Some(Strand::Forward));
353        assert_eq!(right.overhang_sequence, "AATT");
354        assert_eq!(left.overhang_sequence, "AATT");
355        assert_eq!(d.fragments[0].paired_length, 4);
356    }
357    #[test]
358    fn blunt_and_three_prime_geometry() {
359        for (enzyme, seq, polarity) in [
360            ("EcoRV", "AAAGATATCTTT", OverhangPolarity::Blunt),
361            ("KpnI", "AAAGGTACCTTT", OverhangPolarity::ThreePrime),
362        ] {
363            let d = digest(seq, Topology::Linear, &[enzyme]).unwrap();
364            conserved(&d);
365            let end = d.fragments[0].right_end.as_ref().unwrap();
366            assert_eq!(end.polarity, polarity);
367            if enzyme == "KpnI" {
368                assert_eq!(end.protruding_strand, Some(Strand::Forward));
369                assert_eq!(end.overhang_sequence, "GTAC");
370            } else {
371                assert_eq!(end.protruding_strand, None);
372                assert_eq!(d.fragments[0].top.length, 6);
373            }
374        }
375    }
376    #[test]
377    fn nonpalindromic_type_iis_overhang_is_oriented_per_fragment() {
378        let d = digest("AAAGGTCTCAACGATTT", Topology::Linear, &["BsaI"]).unwrap();
379        conserved(&d);
380        assert_eq!(
381            d.fragments[0].right_end.as_ref().unwrap().overhang_sequence,
382            "TCGT"
383        );
384        assert_eq!(
385            d.fragments[1].left_end.as_ref().unwrap().overhang_sequence,
386            "ACGA"
387        );
388        let reverse = digest(
389            &reverse_complement("AAAGGTCTCAACGATTT"),
390            Topology::Linear,
391            &["BsaI"],
392        )
393        .unwrap();
394        conserved(&reverse);
395        assert_eq!(
396            reverse.fragments[0]
397                .right_end
398                .as_ref()
399                .unwrap()
400                .overhang_sequence,
401            "ACGA"
402        );
403    }
404    #[test]
405    fn circular_single_cut_wraps_both_strands() {
406        let d = digest("TTCGAA", Topology::Circular, &["EcoRI"]).unwrap();
407        conserved(&d);
408        assert_eq!(d.fragments.len(), 1);
409        assert_eq!(d.fragments[0].topology, Topology::Linear);
410        assert_eq!(d.fragments[0].top.sequence_5to3, "AATTCG");
411        assert_eq!(d.fragments[0].bottom.sequence_5to3, "AATTCG");
412        assert_eq!(d.fragments[0].paired_length, 2);
413        let kpn = digest("ACCAAAGGT", Topology::Circular, &["KpnI"]).unwrap();
414        conserved(&kpn);
415        assert_eq!(
416            kpn.fragments[0].left_end.as_ref().unwrap().polarity,
417            OverhangPolarity::ThreePrime
418        );
419    }
420    #[test]
421    fn two_cut_circle_keeps_bottom_unwrapped_across_origin() {
422        for (prefix, suffix, enzyme, expected_cuts, expected_lengths) in [
423            (
424                "TTC",
425                "GAA",
426                "EcoRI",
427                [(14, 14), (26, 2)],
428                [(12, 16), (16, 12)],
429            ),
430            (
431                "ACC",
432                "GGT",
433                "KpnI",
434                [(2, 26), (14, 14)],
435                [(12, 16), (16, 12)],
436            ),
437        ] {
438            let seq = format!("{prefix}{}GATATC{}{suffix}", "A".repeat(8), "A".repeat(8));
439            let d = digest(&seq, Topology::Circular, &[enzyme, "EcoRV"]).unwrap();
440            conserved(&d);
441            assert_eq!(d.input_length, 28);
442            assert_eq!(
443                d.cuts
444                    .iter()
445                    .map(|c| (c.top_cut, c.bottom_cut))
446                    .collect::<Vec<_>>(),
447                expected_cuts
448            );
449            assert_eq!(
450                d.fragments
451                    .iter()
452                    .map(|f| (f.top.length, f.bottom.length))
453                    .collect::<Vec<_>>(),
454                expected_lengths
455            );
456            assert!(d.fragments.iter().all(|f| f.paired_length == 12));
457            // Independent reconstruction of each entire source strand in forward
458            // coordinates, starting at the first reported boundary on that strand.
459            for bottom in [false, true] {
460                let start = if bottom {
461                    d.cuts[0].bottom_cut
462                } else {
463                    d.cuts[0].top_cut
464                };
465                let reconstructed: String = d
466                    .fragments
467                    .iter()
468                    .map(|f| {
469                        if bottom {
470                            reverse_complement(&f.bottom.sequence_5to3)
471                        } else {
472                            f.top.sequence_5to3.clone()
473                        }
474                    })
475                    .collect();
476                assert_eq!(reconstructed, format!("{}{}", &seq[start..], &seq[..start]));
477            }
478        }
479    }
480
481    #[test]
482    fn no_cut_retains_topology_and_multiple_cuts_conserve_strands() {
483        for topology in [Topology::Linear, Topology::Circular] {
484            let no_cut = digest("ACGTACGTACGT", topology, &["EcoRI"]).unwrap();
485            conserved(&no_cut);
486            assert_eq!(no_cut.cuts, []);
487            assert_eq!(no_cut.fragments[0].topology, topology);
488            assert_eq!(
489                no_cut.fragments[0].left_end.is_none(),
490                topology == Topology::Circular
491            );
492            let d = digest(
493                "AAAGAATTCAAAAAAGGATCCTTT",
494                topology,
495                &["BamHI", "EcoRI", "ecori"],
496            )
497            .unwrap();
498            conserved(&d);
499            assert_eq!(d.cuts.len(), 2);
500            assert_eq!(
501                d.fragments.len(),
502                if topology == Topology::Circular { 2 } else { 3 }
503            );
504        }
505    }
506    #[test]
507    fn convergent_type_iis_sites_merge_only_identical_physical_cuts() {
508        let d = digest("GGTCTCAAAAAAGAGACG", Topology::Linear, &["BsaI", "BsmBI"]).unwrap();
509        conserved(&d);
510        assert_eq!(d.cuts.len(), 1);
511        assert_eq!(d.cuts[0].enzymes, ["BsaI", "BsmBI"]);
512        assert_eq!(d.fragments.len(), 2);
513        for seq in ["GGTCTCAAAAAGAGACC", "GGTCTCAAAAAAAAAAGAGACC"] {
514            assert!(matches!(
515                digest(seq, Topology::Linear, &["BsaI"]),
516                Err(DigestError::OverlappingCuts)
517            ));
518        }
519    }
520
521    #[test]
522    fn unsafe_geometries_are_rejected_and_duplicate_pairs_merged() {
523        assert!(matches!(
524            digest("GGTCTC", Topology::Linear, &["BsaI"]),
525            Err(DigestError::OutOfBounds { .. })
526        ));
527        assert!(matches!(
528            digest("AAAAAGAGACCAAA", Topology::Linear, &["BsaI"]),
529            Err(DigestError::TerminalCut { .. })
530        ));
531        let seq = DnaSeq::new("ACGTACGTACGT").unwrap();
532        let cuts = vec![
533            Cut {
534                top: 2,
535                bottom: 6,
536                enzymes: vec!["EcoRI"],
537            },
538            Cut {
539                top: 5,
540                bottom: 9,
541                enzymes: vec!["BamHI"],
542            },
543        ];
544        assert!(matches!(
545            partition(&seq, Topology::Linear, &cuts),
546            Err(DigestError::OverlappingCuts)
547        ));
548        let seam = vec![
549            Cut {
550                top: 1,
551                bottom: 5,
552                enzymes: vec!["EcoRI"],
553            },
554            Cut {
555                top: 9,
556                bottom: 13,
557                enzymes: vec!["BamHI"],
558            },
559        ];
560        assert!(matches!(
561            partition(&seq, Topology::Circular, &seam),
562            Err(DigestError::OverlappingCuts)
563        ));
564        let same = merge_cuts(vec![cuts[0].clone(), cuts[0].clone()]);
565        assert_eq!(same.len(), 1);
566        assert_eq!(same[0].enzymes, ["EcoRI"]);
567    }
568}
dnagent_domain/digest.rs

dnagent_domain/
digest.rs
