dnagent_domain/
ligation.rs

1//! Explicit complete restriction/ligation products; no implicit copies or end reuse.
2use crate::compatibility::{CompatibilityError, EndCompatibility, EndSide, compare_ends};
3use crate::digest::{DigestFragment, FragmentEnd, reverse_complement};
4use crate::fragment_annotations::{AnnotatedDigest, AnnotationError, annotated_digest};
5use crate::restriction::OverhangPolarity;
6use crate::{SequenceRecord, Strand, Topology};
7use serde::{Deserialize, Serialize};
8use std::collections::HashSet;
9use thiserror::Error;
10
11pub const MAX_INPUTS: usize = 16;
12pub const MAX_COMPONENTS: usize = 128;
13pub const MAX_STRAND_BASES: usize = 10_000_000;
14
15#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize, Deserialize)]
16#[serde(rename_all = "snake_case")]
17pub enum Orientation {
18    Forward,
19    Reverse,
20}
21
22#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)]
23#[serde(deny_unknown_fields)]
24pub struct FragmentSelection {
25    /// One-based digest-instance number, not an implicit molecule-copy count.
26    pub input: usize,
27    pub fragment_id: String,
28    pub orientation: Orientation,
29}
30
31#[derive(Debug, Error)]
32pub enum LigationError {
33    #[error("invalid ligation plan: {0}")]
34    InvalidPlan(&'static str),
35    #[error(transparent)]
36    Annotation(#[from] AnnotationError),
37    #[error(transparent)]
38    Compatibility(#[from] CompatibilityError),
39    #[error(
40        "incompatible junction after component {after} and before component {before}: {reason:?}"
41    )]
42    Incompatible {
43        after: usize,
44        before: usize,
45        reason: crate::compatibility::CompatibilityReason,
46    },
47    #[error("ligation duplex conservation check failed: {0}")]
48    Conservation(&'static str),
49}
50
51#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
52pub struct ProductSpan {
53    /// Zero-based offset on this product strand's own 5′→3′ sequence.
54    pub start: usize,
55    pub length: usize,
56    /// Which original digest strand supplies these bases and its annotations.
57    pub source_strand: Strand,
58}
59
60#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
61pub struct ComponentPlacement {
62    pub selection: FragmentSelection,
63    pub top: ProductSpan,
64    pub bottom: ProductSpan,
65}
66
67#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
68pub struct Junction {
69    pub after_component: usize,
70    pub before_component: usize,
71    pub closure: bool,
72    /// Nick boundary on each product strand's own 5′→3′ axis; circles modulo length.
73    pub top_boundary: usize,
74    pub bottom_boundary: usize,
75    pub assessment: EndCompatibility,
76}
77
78#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
79pub struct UnusedFragment {
80    pub input: usize,
81    pub fragment_id: String,
82}
83
84#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
85pub struct LigationProduct {
86    pub topology: Topology,
87    pub top_sequence_5to3: String,
88    pub bottom_sequence_5to3: String,
89    /// Start of bottom's forward-coordinate interval relative to top base zero.
90    /// Signed for linear products; normalised modulo product length for circles.
91    pub bottom_forward_start: i128,
92    pub paired_length: usize,
93    pub left_end: Option<FragmentEnd>,
94    pub right_end: Option<FragmentEnd>,
95    pub components: Vec<ComponentPlacement>,
96    pub junctions: Vec<Junction>,
97}
98
99#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
100pub struct LigationReport {
101    /// Original component annotations are linked through product placements.
102    pub inputs: Vec<AnnotatedDigest>,
103    pub product: LigationProduct,
104    pub unused_fragments: Vec<UnusedFragment>,
105    pub assumptions: Vec<&'static str>,
106}
107
108/// Validate bounded request shape before expensive imports or digestion.
109pub fn validate_request(
110    input_count: usize,
111    selected: &[FragmentSelection],
112) -> Result<(), LigationError> {
113    if !(1..=MAX_INPUTS).contains(&input_count) {
114        return Err(LigationError::InvalidPlan(
115            "expected 1..=16 input digest instances",
116        ));
117    }
118    if selected.is_empty() || selected.len() > MAX_COMPONENTS {
119        return Err(LigationError::InvalidPlan(
120            "expected 1..=128 selected components",
121        ));
122    }
123    let mut used = HashSet::new();
124    for selection in selected {
125        if selection.input == 0 || selection.input > input_count || selection.fragment_id.is_empty()
126        {
127            return Err(LigationError::InvalidPlan(
128                "invalid input number or empty fragment ID",
129            ));
130        }
131        if !used.insert((selection.input, &selection.fragment_id)) {
132            return Err(LigationError::InvalidPlan(
133                "fragment reused; declare a separate input digest instance for another copy",
134            ));
135        }
136    }
137    Ok(())
138}
139
140pub fn simulate_ligation(
141    inputs: &[(&SequenceRecord, &[String])],
142    selected: &[FragmentSelection],
143    topology: Topology,
144) -> Result<LigationReport, LigationError> {
145    validate_request(inputs.len(), selected)?;
146    // Accept validated source records, not caller-forged public digest structs.
147    let inputs = inputs
148        .iter()
149        .map(|(r, names)| annotated_digest(r, names))
150        .collect::<Result<Vec<_>, _>>()?;
151    let components = selected
152        .iter()
153        .map(|selection| {
154            let fragment = inputs[selection.input - 1]
155                .digest
156                .fragments
157                .iter()
158                .find(|f| f.id == selection.fragment_id)
159                .ok_or(LigationError::InvalidPlan(
160                    "selected fragment does not exist",
161                ))?;
162            orient(fragment, selection.orientation)
163        })
164        .collect::<Result<Vec<_>, _>>()?;
165    let product = assemble(&components, selected, topology)?;
166    let unused_fragments = inputs
167        .iter()
168        .enumerate()
169        .flat_map(|(index, input)| {
170            input
171                .digest
172                .fragments
173                .iter()
174                .filter(move |f| {
175                    !selected
176                        .iter()
177                        .any(|s| s.input == index + 1 && s.fragment_id == f.id)
178                })
179                .map(move |f| UnusedFragment {
180                    input: index + 1,
181                    fragment_id: f.id.clone(),
182                })
183        })
184        .collect();
185    Ok(LigationReport {
186        inputs,
187        product,
188        unused_fragments,
189        assumptions: vec![
190            "Explicit digest instances and component order; no implicit fragment copies or endpoint reuse",
191            "Assumes selected sequence-compatible junctions seal completely; phosphorylation, nicks, reaction conditions and experimental yield are not predicted",
192            "Reverse orientation exchanges whole top/bottom 5-prime-to-3-prime sequences and ends; exposed physical oligos are not reverse-complemented",
193            "Both selected strands are conserved without trimming, fill-in or overlap deduplication; unused fragments are reported separately",
194            "Component-local source annotations remain in inputs and map through product spans; gene repair, fusion, joined feature reconstruction and functional integrity are not inferred",
195            "For circular products bottom_forward_start records the phase between strand coordinate origins; bottom sequence is not silently rotated",
196            "Linear single-component plans are pass-through products with no ligation junctions",
197        ],
198    })
199}
200
201struct Oriented {
202    top: String,
203    bottom: String,
204    left: FragmentEnd,
205    right: FragmentEnd,
206}
207
208fn flip(strand: Strand) -> Strand {
209    match strand {
210        Strand::Forward => Strand::Reverse,
211        Strand::Reverse => Strand::Forward,
212        Strand::Unknown => Strand::Unknown,
213    }
214}
215
216fn orient(fragment: &DigestFragment, orientation: Orientation) -> Result<Oriented, LigationError> {
217    let (Some(left), Some(right)) = (&fragment.left_end, &fragment.right_end) else {
218        return Err(LigationError::InvalidPlan(
219            "uncut circular fragments have no free ends",
220        ));
221    };
222    if orientation == Orientation::Forward {
223        return Ok(Oriented {
224            top: fragment.top.sequence_5to3.clone(),
225            bottom: fragment.bottom.sequence_5to3.clone(),
226            left: left.clone(),
227            right: right.clone(),
228        });
229    }
230    let mut left = right.clone();
231    let mut right = fragment
232        .left_end
233        .clone()
234        .ok_or(LigationError::InvalidPlan("missing end"))?;
235    left.protruding_strand = left.protruding_strand.map(flip);
236    right.protruding_strand = right.protruding_strand.map(flip);
237    Ok(Oriented {
238        top: fragment.bottom.sequence_5to3.clone(),
239        bottom: fragment.top.sequence_5to3.clone(),
240        left,
241        right,
242    })
243}
244
245fn total_length(components: &[Oriented], bottom: bool) -> Result<usize, LigationError> {
246    components.iter().try_fold(0usize, |total, component| {
247        total
248            .checked_add(if bottom {
249                component.bottom.len()
250            } else {
251                component.top.len()
252            })
253            .filter(|n| *n <= MAX_STRAND_BASES)
254            .ok_or(LigationError::InvalidPlan(
255                "product exceeds 10,000,000 bases per strand",
256            ))
257    })
258}
259
260fn assemble(
261    oriented: &[Oriented],
262    selected: &[FragmentSelection],
263    topology: Topology,
264) -> Result<LigationProduct, LigationError> {
265    let nt = total_length(oriented, false)?;
266    let nb = total_length(oriented, true)?;
267    let mut components = Vec::new();
268    let mut junctions = Vec::new();
269    let (mut top, mut bottom) = (0, 0);
270    for (index, (component, selection)) in oriented.iter().zip(selected).enumerate() {
271        let source_strand = if selection.orientation == Orientation::Forward {
272            Strand::Forward
273        } else {
274            Strand::Reverse
275        };
276        components.push(ComponentPlacement {
277            selection: selection.clone(),
278            top: ProductSpan {
279                start: top,
280                length: component.top.len(),
281                source_strand,
282            },
283            bottom: ProductSpan {
284                start: nb - bottom - component.bottom.len(),
285                length: component.bottom.len(),
286                source_strand: flip(source_strand),
287            },
288        });
289        top += component.top.len();
290        bottom += component.bottom.len();
291        let closure = index + 1 == oriented.len();
292        if !closure || topology == Topology::Circular {
293            let next = (index + 1) % oriented.len();
294            let assessment = compare_ends(
295                &component.right,
296                EndSide::Right,
297                &oriented[next].left,
298                EndSide::Left,
299            )?;
300            if !assessment.compatible {
301                return Err(LigationError::Incompatible {
302                    after: index + 1,
303                    before: next + 1,
304                    reason: assessment.reason,
305                });
306            }
307            junctions.push(Junction {
308                after_component: index + 1,
309                before_component: next + 1,
310                closure,
311                top_boundary: if closure { 0 } else { top },
312                bottom_boundary: nb - bottom,
313                assessment,
314            });
315        }
316    }
317    let first = &oriented[0];
318    let last = &oriented[oriented.len() - 1];
319    let top_sequence_5to3: String = oriented.iter().map(|c| c.top.as_str()).collect();
320    let bottom_sequence_5to3: String = oriented.iter().rev().map(|c| c.bottom.as_str()).collect();
321    if top_sequence_5to3.len() != nt || bottom_sequence_5to3.len() != nb {
322        return Err(LigationError::Conservation(
323            "selected strand lengths changed",
324        ));
325    }
326    let delta = stagger(&first.left);
327    let (bottom_forward_start, paired_length) = verify_duplex(
328        &top_sequence_5to3,
329        &bottom_sequence_5to3,
330        delta,
331        stagger(&last.right),
332        topology,
333    )?;
334    Ok(LigationProduct {
335        topology,
336        top_sequence_5to3,
337        bottom_sequence_5to3,
338        bottom_forward_start,
339        paired_length,
340        left_end: (topology == Topology::Linear).then(|| first.left.clone()),
341        right_end: (topology == Topology::Linear).then(|| last.right.clone()),
342        components,
343        junctions,
344    })
345}
346
347fn stagger(end: &FragmentEnd) -> i128 {
348    let n = end.overhang_sequence.len() as i128;
349    match end.polarity {
350        OverhangPolarity::Blunt => 0,
351        OverhangPolarity::FivePrime => n,
352        OverhangPolarity::ThreePrime => -n,
353    }
354}
355
356fn verify_duplex(
357    top: &str,
358    bottom: &str,
359    delta: i128,
360    right_delta: i128,
361    topology: Topology,
362) -> Result<(i128, usize), LigationError> {
363    let bottom_forward = reverse_complement(bottom);
364    if topology == Topology::Circular {
365        if top.len() != bottom.len() {
366            return Err(LigationError::Conservation(
367                "unequal circular strand lengths",
368            ));
369        }
370        let offset = delta.rem_euclid(top.len() as i128) as usize;
371        if bottom_forward != format!("{}{}", &top[offset..], &top[..offset]) {
372            return Err(LigationError::Conservation(
373                "circular strand phase mismatch",
374            ));
375        }
376        return Ok((offset as i128, top.len()));
377    }
378    let end = delta + bottom.len() as i128;
379    if end - top.len() as i128 != right_delta {
380        return Err(LigationError::Conservation("terminal stagger mismatch"));
381    }
382    let lo = delta.max(0);
383    let hi = end.min(top.len() as i128);
384    if lo >= hi
385        || top[lo as usize..hi as usize]
386            != bottom_forward[(lo - delta) as usize..(hi - delta) as usize]
387    {
388        return Err(LigationError::Conservation("linear paired-core mismatch"));
389    }
390    Ok((delta, (hi - lo) as usize))
391}
392
393#[cfg(test)]
394mod tests {
395    use super::*;
396    use crate::DnaSeq;
397    fn record(seq: &str, topology: Topology) -> SequenceRecord {
398        SequenceRecord::new(
399            "synthetic",
400            DnaSeq::new(seq).unwrap(),
401            topology,
402            vec![],
403            vec![],
404        )
405        .unwrap()
406    }
407    fn selection(input: usize, number: usize, orientation: Orientation) -> FragmentSelection {
408        FragmentSelection {
409            input,
410            fragment_id: format!("fragment-{number:04}"),
411            orientation,
412        }
413    }
414    #[test]
415    fn rejoining_linear_digest_conserves_both_strands() {
416        for (seq, enzyme) in [
417            ("AAAGAATTCTTT", "EcoRI"),
418            ("AAAGGTACCTTT", "KpnI"),
419            ("AAAGATATCTTT", "EcoRV"),
420        ] {
421            let r = record(seq, Topology::Linear);
422            let names = vec![enzyme.into()];
423            let selected = [
424                selection(1, 1, Orientation::Forward),
425                selection(1, 2, Orientation::Forward),
426            ];
427            let report = simulate_ligation(&[(&r, &names)], &selected, Topology::Linear).unwrap();
428            assert_eq!(report.product.top_sequence_5to3, seq);
429            assert_eq!(report.product.bottom_sequence_5to3, reverse_complement(seq));
430            assert_eq!(report.product.junctions.len(), 1);
431            assert_eq!(report.product.paired_length, seq.len());
432            assert_eq!(report.unused_fragments, []);
433        }
434    }
435    #[test]
436    fn circle_closure_keeps_explicit_bottom_phase() {
437        for (seq, enzyme, phase) in [
438            ("AAAGAATTCTTT", "EcoRI", 4),
439            ("AAAGGTACCTTT", "KpnI", 8),
440            ("AAAGATATCTTT", "EcoRV", 0),
441        ] {
442            let r = record(seq, Topology::Circular);
443            let names = vec![enzyme.into()];
444            for orientation in [Orientation::Forward, Orientation::Reverse] {
445                let report = simulate_ligation(
446                    &[(&r, &names)],
447                    &[selection(1, 1, orientation)],
448                    Topology::Circular,
449                )
450                .unwrap();
451                assert_eq!(report.product.bottom_forward_start, phase);
452                assert!(report.product.left_end.is_none());
453                assert!(report.product.junctions[0].closure);
454                assert_eq!(report.product.paired_length, 12);
455            }
456        }
457    }
458    #[test]
459    fn reverse_order_and_orientation_reverse_the_whole_duplex() {
460        let r = record("AAAGAATTCTGC", Topology::Linear);
461        let names = vec!["EcoRI".into()];
462        let report = simulate_ligation(
463            &[(&r, &names)],
464            &[
465                selection(1, 2, Orientation::Reverse),
466                selection(1, 1, Orientation::Reverse),
467            ],
468            Topology::Linear,
469        )
470        .unwrap();
471        assert_eq!(report.product.top_sequence_5to3, "GCAGAATTCTTT");
472        assert_eq!(report.product.bottom_sequence_5to3, r.sequence().as_str());
473        assert_eq!(
474            report.product.components[0].top.source_strand,
475            Strand::Reverse
476        );
477    }
478    #[test]
479    fn pass_through_preserves_sticky_termini_and_paired_core() {
480        for (seq, enzyme) in [("AAAGAATTCTGC", "EcoRI"), ("AAAGGTACCTGC", "KpnI")] {
481            let r = record(seq, Topology::Linear);
482            let names = vec![enzyme.into()];
483            for number in [1, 2] {
484                for orientation in [Orientation::Forward, Orientation::Reverse] {
485                    let report = simulate_ligation(
486                        &[(&r, &names)],
487                        &[selection(1, number, orientation)],
488                        Topology::Linear,
489                    )
490                    .unwrap();
491                    let source = &report.inputs[0].digest.fragments[number - 1];
492                    assert_eq!(report.product.paired_length, source.paired_length);
493                    assert_eq!(report.product.junctions, []);
494                    assert_eq!(report.unused_fragments.len(), 1);
495                }
496            }
497        }
498    }
499
500    #[test]
501    fn resource_limits_are_checked_at_the_boundary() {
502        let s = selection(1, 1, Orientation::Forward);
503        assert!(validate_request(1, &vec![s; MAX_COMPONENTS + 1]).is_err());
504        // The length helper's threshold is also tested at the exact boundary.
505        let base = record("AAAA", Topology::Linear);
506        let names = vec!["EcoRI".into()];
507        let fragment = annotated_digest(&base, &names)
508            .unwrap()
509            .digest
510            .fragments
511            .remove(0);
512        let mut component = orient(&fragment, Orientation::Forward).unwrap();
513        component.top = "A".repeat(MAX_STRAND_BASES);
514        assert_eq!(
515            total_length(std::slice::from_ref(&component), false).unwrap(),
516            MAX_STRAND_BASES
517        );
518        component.top.push('A');
519        assert!(total_length(&[component], false).is_err());
520    }
521
522    #[test]
523    fn refuses_reuse_missing_fragments_and_uncut_circles() {
524        let r = record("AAAGAATTCTTT", Topology::Circular);
525        let names = vec!["BamHI".into()];
526        let s = selection(1, 1, Orientation::Forward);
527        assert!(
528            simulate_ligation(
529                &[(&r, &names)],
530                std::slice::from_ref(&s),
531                Topology::Circular
532            )
533            .is_err()
534        );
535        assert!(validate_request(1, &[s.clone(), s]).is_err());
536        assert!(
537            simulate_ligation(
538                &[(&r, &names)],
539                &[selection(1, 99, Orientation::Forward)],
540                Topology::Linear
541            )
542            .is_err()
543        );
544        assert!(validate_request(17, &[selection(1, 1, Orientation::Forward)]).is_err());
545    }
546    #[test]
547    fn incompatible_junctions_fail_and_distinct_instances_allow_explicit_copies() {
548        let a = record("AAAGAATTCTTT", Topology::Linear);
549        let b = record("AAAGGTACCTTT", Topology::Linear);
550        let ea = vec!["EcoRI".into()];
551        let eb = vec!["KpnI".into()];
552        let selected = [
553            selection(1, 1, Orientation::Forward),
554            selection(2, 2, Orientation::Forward),
555        ];
556        assert!(matches!(
557            simulate_ligation(&[(&a, &ea), (&b, &eb)], &selected, Topology::Linear),
558            Err(LigationError::Incompatible { .. })
559        ));
560        let report =
561            simulate_ligation(&[(&a, &ea), (&a, &ea)], &selected, Topology::Linear).unwrap();
562        assert_eq!(report.unused_fragments.len(), 2);
563        assert_eq!(report.product.top_sequence_5to3, a.sequence().as_str());
564    }
565}
dnagent_domain/ligation.rs

dnagent_domain/
ligation.rs
