dnagent_domain/
compatibility.rs

1//! Sequence-only restriction-end compatibility; no ligation product prediction.
2use crate::Strand;
3use crate::digest::{Digest, FragmentEnd};
4use crate::restriction::OverhangPolarity;
5use serde::Serialize;
6use thiserror::Error;
7
8/// Output/resource guard: unordered distinct pairs scale as n*(n-1)/2,
9/// so 128 endpoints cap the matrix at 8,128 comparisons (not a biological limit).
10pub const MAX_ENDPOINTS: usize = 128;
11
12#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize)]
13#[serde(rename_all = "snake_case")]
14pub enum EndSide {
15    Left,
16    Right,
17}
18
19#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize)]
20#[serde(rename_all = "snake_case")]
21pub enum CompatibilityReason {
22    BluntEnds,
23    ComplementaryOverhangs,
24    PolarityMismatch,
25    OverhangLengthMismatch,
26    OverhangSequenceMismatch,
27}
28
29#[derive(Debug, Clone, PartialEq, Eq, Error)]
30pub enum CompatibilityError {
31    #[error("invalid fragment end: {0}")]
32    InvalidEnd(&'static str),
33    #[error("invalid digest endpoint structure: {0}")]
34    InvalidDigest(&'static str),
35    #[error("expected one or two digests")]
36    InputCount,
37    #[error(
38        "{count} endpoints exceeds the explicit comparison limit of {MAX_ENDPOINTS}; select fewer enzymes"
39    )]
40    TooManyEnds { count: usize },
41}
42
43#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
44pub struct EndCompatibility {
45    pub compatible: bool,
46    pub reason: CompatibilityReason,
47    /// Relative to the first fragment, which stays in its stored orientation.
48    /// This is a proposed orientation, not an applied sequence transformation.
49    pub second_fragment_orientation: Strand,
50    /// Position along the first fragment's fixed forward axis.
51    pub second_placement: Placement,
52}
53
54#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize)]
55#[serde(rename_all = "snake_case")]
56pub enum Placement {
57    BeforeFirst,
58    AfterFirst,
59}
60
61/// Compare physical oligos (each written 5′→3′). Enzyme names are provenance,
62/// not compatibility criteria. Same-side ends require reversing fragment two.
63pub fn compare_ends(
64    first: &FragmentEnd,
65    first_side: EndSide,
66    second: &FragmentEnd,
67    second_side: EndSide,
68) -> Result<EndCompatibility, CompatibilityError> {
69    validate_end(first, first_side)?;
70    validate_end(second, second_side)?;
71    let reason = if first.polarity != second.polarity {
72        CompatibilityReason::PolarityMismatch
73    } else if first.polarity == OverhangPolarity::Blunt {
74        CompatibilityReason::BluntEnds
75    } else if first.overhang_sequence.len() != second.overhang_sequence.len() {
76        CompatibilityReason::OverhangLengthMismatch
77    } else if first.overhang_sequence.bytes().eq(second
78        .overhang_sequence
79        .bytes()
80        .rev()
81        .map(complement))
82    {
83        CompatibilityReason::ComplementaryOverhangs
84    } else {
85        CompatibilityReason::OverhangSequenceMismatch
86    };
87    Ok(EndCompatibility {
88        compatible: matches!(
89            reason,
90            CompatibilityReason::BluntEnds | CompatibilityReason::ComplementaryOverhangs
91        ),
92        reason,
93        second_placement: if first_side == EndSide::Left {
94            Placement::BeforeFirst
95        } else {
96            Placement::AfterFirst
97        },
98        second_fragment_orientation: if first_side == second_side {
99            Strand::Reverse
100        } else {
101            Strand::Forward
102        },
103    })
104}
105
106fn complement(base: u8) -> u8 {
107    match base {
108        b'A' => b'T',
109        b'T' => b'A',
110        b'C' => b'G',
111        b'G' => b'C',
112        _ => unreachable!("end validated"),
113    }
114}
115
116fn validate_end(end: &FragmentEnd, side: EndSide) -> Result<(), CompatibilityError> {
117    if end.polarity == OverhangPolarity::Blunt {
118        if !end.overhang_sequence.is_empty() || end.protruding_strand.is_some() {
119            return Err(CompatibilityError::InvalidEnd(
120                "blunt ends cannot have a protruding oligo",
121            ));
122        }
123    } else {
124        if end.original_terminus {
125            return Err(CompatibilityError::InvalidEnd(
126                "original termini are assumed blunt in this model",
127            ));
128        }
129        if end.overhang_sequence.is_empty()
130            || !end
131                .overhang_sequence
132                .bytes()
133                .all(|b| matches!(b, b'A' | b'C' | b'G' | b'T'))
134        {
135            return Err(CompatibilityError::InvalidEnd(
136                "cohesive ends require a nonempty uppercase ACGT oligo",
137            ));
138        }
139        let top = (end.polarity == OverhangPolarity::FivePrime) == (side == EndSide::Left);
140        if end.protruding_strand
141            != Some(if top {
142                Strand::Forward
143            } else {
144                Strand::Reverse
145            })
146        {
147            return Err(CompatibilityError::InvalidEnd(
148                "strand, polarity and fragment side disagree",
149            ));
150        }
151    }
152    Ok(())
153}
154
155#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
156pub struct Endpoint {
157    pub id: String,
158    /// One-based input number; fragment IDs are local to this input's digest.
159    pub input: usize,
160    pub fragment_id: String,
161    pub side: EndSide,
162    pub end: FragmentEnd,
163}
164
165#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
166pub struct EndPair {
167    pub first: String,
168    pub second: String,
169    /// True only for opposite ends of the same fragment, not a product prediction.
170    pub same_fragment: bool,
171    pub assessment: EndCompatibility,
172}
173
174#[derive(Debug, Clone, PartialEq, Eq, Serialize)]
175pub struct CompatibilityReport {
176    pub endpoints: Vec<Endpoint>,
177    pub pairs: Vec<EndPair>,
178    pub assumptions: Vec<&'static str>,
179}
180
181/// All unordered pairs of distinct physical ends in one or two digests.
182/// Uncut circles contribute no ends. Reuse of the same end/copies is not inferred.
183pub fn compatible_ends(digests: &[Digest]) -> Result<CompatibilityReport, CompatibilityError> {
184    if !(1..=2).contains(&digests.len()) {
185        return Err(CompatibilityError::InputCount);
186    }
187    let mut endpoints = Vec::new();
188    for (index, digest) in digests.iter().enumerate() {
189        let mut ids = std::collections::HashSet::new();
190        for fragment in &digest.fragments {
191            if fragment.id.is_empty() || !ids.insert(&fragment.id) {
192                return Err(CompatibilityError::InvalidDigest(
193                    "empty or duplicate fragment ID",
194                ));
195            }
196            let has_ends = fragment.topology == crate::Topology::Linear;
197            if fragment.left_end.is_some() != has_ends || fragment.right_end.is_some() != has_ends {
198                return Err(CompatibilityError::InvalidDigest(
199                    "fragment topology and free ends disagree",
200                ));
201            }
202            for (side, end) in [
203                (EndSide::Left, &fragment.left_end),
204                (EndSide::Right, &fragment.right_end),
205            ] {
206                if let Some(end) = end {
207                    validate_end(end, side)?;
208                    let label = if side == EndSide::Left {
209                        "left"
210                    } else {
211                        "right"
212                    };
213                    endpoints.push(Endpoint {
214                        id: format!("input-{}:{}:{label}", index + 1, fragment.id),
215                        input: index + 1,
216                        fragment_id: fragment.id.clone(),
217                        side,
218                        end: end.clone(),
219                    });
220                }
221            }
222        }
223    }
224    if endpoints.len() > MAX_ENDPOINTS {
225        return Err(CompatibilityError::TooManyEnds {
226            count: endpoints.len(),
227        });
228    }
229    let mut pairs = Vec::new();
230    for (i, first) in endpoints.iter().enumerate() {
231        for second in &endpoints[i + 1..] {
232            pairs.push(EndPair {
233                first: first.id.clone(),
234                second: second.id.clone(),
235                same_fragment: first.input == second.input
236                    && first.fragment_id == second.fragment_id,
237                assessment: compare_ends(&first.end, first.side, &second.end, second.side)?,
238            });
239        }
240    }
241    Ok(CompatibilityReport {
242        endpoints,
243        pairs,
244        assumptions: vec![
245            "Sequence-only end compatibility, not ligation efficiency, an assembled product or experimental validation",
246            "Matching polarity and reverse-complementary 5-prime-to-3-prime oligos required; blunt ends are sequence-compatible",
247            "First fragment orientation is fixed; the reported second orientation is a proposal, not an applied transformation",
248            "Phosphorylation, damage, ligase, reaction conditions and circularisation geometry are not modelled; original linear ends are assumed blunt",
249            "All unordered distinct-end pairs, including within-input and same-fragment pairs; no end reuse, molecule copies or global assembly inference",
250            "Uncut circles have no free ends; identical enzyme names do not imply compatibility and different names do not imply incompatibility",
251        ],
252    })
253}
254
255#[cfg(test)]
256mod tests {
257    use super::*;
258    fn end(seq: &str, polarity: OverhangPolarity, side: EndSide) -> FragmentEnd {
259        FragmentEnd {
260            original_terminus: false,
261            polarity,
262            protruding_strand: if polarity == OverhangPolarity::Blunt {
263                None
264            } else {
265                Some(
266                    if (polarity == OverhangPolarity::FivePrime) == (side == EndSide::Left) {
267                        Strand::Forward
268                    } else {
269                        Strand::Reverse
270                    },
271                )
272            },
273            overhang_sequence: seq.into(),
274            enzymes: vec!["BsaI"],
275        }
276    }
277    #[test]
278    fn reasons_and_orientation_are_explicit() {
279        use OverhangPolarity::{Blunt, FivePrime, ThreePrime};
280        let first = end("ACGA", FivePrime, EndSide::Right);
281        for (seq, polarity, expected) in [
282            (
283                "TCGT",
284                FivePrime,
285                CompatibilityReason::ComplementaryOverhangs,
286            ),
287            (
288                "ACGA",
289                FivePrime,
290                CompatibilityReason::OverhangSequenceMismatch,
291            ),
292            (
293                "TCG",
294                FivePrime,
295                CompatibilityReason::OverhangLengthMismatch,
296            ),
297            ("TCGT", ThreePrime, CompatibilityReason::PolarityMismatch),
298            ("", Blunt, CompatibilityReason::PolarityMismatch),
299        ] {
300            for side in [EndSide::Left, EndSide::Right] {
301                let result =
302                    compare_ends(&first, EndSide::Right, &end(seq, polarity, side), side).unwrap();
303                assert_eq!(result.reason, expected);
304                assert_eq!(
305                    result.compatible,
306                    expected == CompatibilityReason::ComplementaryOverhangs
307                );
308                assert_eq!(
309                    result.second_fragment_orientation,
310                    if side == EndSide::Right {
311                        Strand::Reverse
312                    } else {
313                        Strand::Forward
314                    }
315                );
316            }
317        }
318        assert!(
319            compare_ends(
320                &end("", Blunt, EndSide::Left),
321                EndSide::Left,
322                &end("", Blunt, EndSide::Right),
323                EndSide::Right
324            )
325            .unwrap()
326            .compatible
327        );
328    }
329    #[test]
330    fn three_prime_and_cross_enzyme_compatibility() {
331        for polarity in [OverhangPolarity::FivePrime, OverhangPolarity::ThreePrime] {
332            let a = end("ACGA", polarity, EndSide::Left);
333            let mut b = end("TCGT", polarity, EndSide::Right);
334            b.enzymes = vec!["BsmBI"];
335            assert!(
336                compare_ends(&a, EndSide::Left, &b, EndSide::Right)
337                    .unwrap()
338                    .compatible
339            );
340        }
341    }
342    #[test]
343    fn malformed_public_end_values_are_rejected() {
344        let good = end("ACGA", OverhangPolarity::FivePrime, EndSide::Left);
345        for bad in [
346            end("N", OverhangPolarity::FivePrime, EndSide::Left),
347            end("", OverhangPolarity::FivePrime, EndSide::Left),
348            end("A", OverhangPolarity::Blunt, EndSide::Left),
349            end("ACGA", OverhangPolarity::FivePrime, EndSide::Right),
350        ] {
351            assert!(compare_ends(&bad, EndSide::Left, &good, EndSide::Left).is_err());
352        }
353    }
354    #[test]
355    fn circle_self_closure_and_no_free_ends() {
356        let seq = crate::DnaSeq::new("AAAGAATTCTTT").unwrap();
357        let digest =
358            crate::digest::simulate_digest(&seq, crate::Topology::Circular, &["EcoRI".into()])
359                .unwrap();
360        let report = compatible_ends(&[digest]).unwrap();
361        assert_eq!(report.pairs.len(), 1);
362        assert!(report.pairs[0].same_fragment && report.pairs[0].assessment.compatible);
363        let uncut =
364            crate::digest::simulate_digest(&seq, crate::Topology::Circular, &["BamHI".into()])
365                .unwrap();
366        assert_eq!(compatible_ends(&[uncut]).unwrap().endpoints, []);
367    }
368    #[test]
369    fn public_digest_endpoint_structure_is_checked() {
370        assert!(matches!(
371            compatible_ends(&[]),
372            Err(CompatibilityError::InputCount)
373        ));
374        let seq = crate::DnaSeq::new("AAAGAATTCTTT").unwrap();
375        let mut digest =
376            crate::digest::simulate_digest(&seq, crate::Topology::Linear, &["EcoRI".into()])
377                .unwrap();
378        digest.fragments[1].id = digest.fragments[0].id.clone();
379        assert!(matches!(
380            compatible_ends(&[digest.clone()]),
381            Err(CompatibilityError::InvalidDigest(_))
382        ));
383        digest.fragments[1].id = "unique".into();
384        digest.fragments[0].topology = crate::Topology::Circular;
385        assert!(matches!(
386            compatible_ends(&[digest]),
387            Err(CompatibilityError::InvalidDigest(_))
388        ));
389        let mut original = end("ACGA", OverhangPolarity::FivePrime, EndSide::Left);
390        original.original_terminus = true;
391        assert!(compare_ends(&original, EndSide::Left, &original, EndSide::Left).is_err());
392    }
393
394    #[test]
395    fn bounded_and_deterministic_pair_matrix() {
396        let seq = crate::DnaSeq::new("AAAGAATTCTTT").unwrap();
397        let mut digest =
398            crate::digest::simulate_digest(&seq, crate::Topology::Linear, &["EcoRI".into()])
399                .unwrap();
400        let report = compatible_ends(&[digest.clone(), digest.clone()]).unwrap();
401        assert_eq!(report.endpoints.len(), 8);
402        assert_eq!(report.pairs.len(), 28);
403        assert_eq!(
404            report,
405            compatible_ends(&[digest.clone(), digest.clone()]).unwrap()
406        );
407        let fragment = digest.fragments[0].clone();
408        digest.fragments = (0..MAX_ENDPOINTS)
409            .map(|i| {
410                let mut f = fragment.clone();
411                f.id = format!("f{i}");
412                f
413            })
414            .collect();
415        assert!(matches!(
416            compatible_ends(&[digest]),
417            Err(CompatibilityError::TooManyEnds { .. })
418        ));
419    }
420}
dnagent_domain/compatibility.rs

dnagent_domain/
compatibility.rs
