dnagent_domain/
existing_overlaps.rs

1//! Assembly of declared linear fragment views with explicit exact overlaps.
2use crate::fragment_annotations::FeatureMapping;
3use crate::gibson::{
4    CoreSelection, GibsonError, GibsonJunction, MAX_BASES, core, core_annotations,
5    unique_duplex_site, validate_templates,
6};
7use crate::{SequenceRecord, Topology};
8use serde::Serialize;
9
10#[derive(Debug, Serialize)]
11pub struct ExistingComponent {
12    pub selection: CoreSelection,
13    pub sequence_5to3: String,
14    pub product_start: usize,
15    pub wraps_origin: bool,
16    pub annotations: Vec<FeatureMapping>,
17}
18
19#[derive(Debug, Serialize)]
20pub struct ExistingAssembly {
21    pub inputs: Vec<SequenceRecord>,
22    pub topology: Topology,
23    pub product_sequence_5to3: String,
24    pub components: Vec<ExistingComponent>,
25    pub junctions: Vec<GibsonJunction>,
26    pub assumptions: Vec<&'static str>,
27}
28
29pub fn validate_request(
30    count: usize,
31    selected: &[CoreSelection],
32    topology: Topology,
33    overlaps: &[usize],
34) -> Result<(), GibsonError> {
35    if !(1..=16).contains(&count) || selected.is_empty() || selected.len() > 128 {
36        return Err(GibsonError::Invalid(
37            "expected 1..=16 sources and 1..=128 existing fragments",
38        ));
39    }
40    let junction_count = selected.len() - usize::from(topology == Topology::Linear);
41    if overlaps.len() != junction_count || overlaps.iter().any(|n| !(20..=60).contains(n)) {
42        return Err(GibsonError::Invalid(
43            "supply exactly one 20..=60-base overlap per junction, including circular closure",
44        ));
45    }
46    let mut total = 0usize;
47    for (i, s) in selected.iter().enumerate() {
48        let incoming = if i > 0 {
49            overlaps[i - 1]
50        } else if topology == Topology::Circular {
51            overlaps[selected.len() - 1]
52        } else {
53            0
54        };
55        let outgoing = overlaps.get(i).copied().unwrap_or(0);
56        if s.input == 0 || s.input > count || s.length <= incoming + outgoing {
57            return Err(GibsonError::Invalid(
58                "invalid source index or fragment without a non-overlapping interior",
59            ));
60        }
61        total = total
62            .checked_add(s.length)
63            .filter(|n| *n <= MAX_BASES)
64            .ok_or(GibsonError::Invalid(
65                "total supplied fragment bases exceed 1,000,000",
66            ))?;
67    }
68    Ok(())
69}
70
71pub fn assemble(
72    records: &[SequenceRecord],
73    selected: &[CoreSelection],
74    topology: Topology,
75    overlaps: &[usize],
76) -> Result<ExistingAssembly, GibsonError> {
77    validate_request(records.len(), selected, topology, overlaps)?;
78    validate_templates(records)?;
79    let sequences = selected
80        .iter()
81        .map(|s| core(&records[s.input - 1], s))
82        .collect::<Result<Vec<_>, _>>()?;
83    let mut product = sequences[0].clone();
84    let mut starts = vec![0];
85    for (i, next) in sequences.iter().enumerate().skip(1) {
86        let overlap = overlaps[i - 1];
87        check_overlap(&product, next, overlap, i)?;
88        starts.push(product.len() - overlap);
89        product.push_str(&next[overlap..]);
90    }
91    if topology == Topology::Circular {
92        let overlap = overlaps[overlaps.len() - 1];
93        check_overlap(&product, &sequences[0], overlap, selected.len())?;
94        product.truncate(product.len() - overlap);
95    }
96    let mut junctions = Vec::new();
97    for (i, &length) in overlaps.iter().enumerate() {
98        let next = (i + 1) % selected.len();
99        let motif = &sequences[next][..length];
100        if !unique_duplex_site(&product, motif, topology) {
101            return Err(GibsonError::Component {
102                component: i + 1,
103                reason: "existing overlap is repeated or reverse-complement ambiguous in the merged product",
104            });
105        }
106        junctions.push(GibsonJunction {
107            after_component: i + 1,
108            before_component: next + 1,
109            closure: next == 0,
110            product_start: starts[next],
111            overlap_sequence_5to3: motif.into(),
112            added_by: None,
113        });
114    }
115    let mut components = Vec::new();
116    for (i, (s, sequence)) in selected.iter().zip(sequences).enumerate() {
117        // Conservation includes both source associations at merged homology bases.
118        if !sequence
119            .bytes()
120            .enumerate()
121            .all(|(j, b)| product.as_bytes()[(starts[i] + j) % product.len()] == b)
122        {
123            return Err(GibsonError::Invalid(
124                "existing-fragment source-to-product conservation failed",
125            ));
126        }
127        components.push(ExistingComponent {
128            selection: s.clone(),
129            product_start: starts[i],
130            wraps_origin: starts[i] + sequence.len() > product.len(),
131            annotations: core_annotations(&records[s.input - 1], s, &sequence)?,
132            sequence_5to3: sequence,
133        });
134    }
135    Ok(ExistingAssembly {
136        inputs: records.to_vec(),
137        topology,
138        product_sequence_5to3: product,
139        components,
140        junctions,
141        assumptions: vec![
142            "Selections declare existing linear fragment views, even when source coordinate topology is circular; no PCR or digestion preparation is inferred",
143            "Only explicitly declared exact overlaps are merged; no ordering, orientation, polishing or overlap-length search is performed",
144            "Each overlap must be unique across both orientations of the merged product; partial/mismatched homology and reaction kinetics are not assessed",
145            "Component-local annotations map by product_start plus local coordinate, modulo product length for circles; both source associations at shared homology are retained",
146            "Merged homology has one sequence copy but may have multiple source annotations; genes, translations and biological function are not reconstructed",
147            "Assumes ideal fully paired assembly; no experimental validation, yield prediction, raw XML or opaque metadata export",
148        ],
149    })
150}
151
152fn check_overlap(
153    first: &str,
154    second: &str,
155    length: usize,
156    component: usize,
157) -> Result<(), GibsonError> {
158    if first.len() < length
159        || second.len() < length
160        || first[first.len() - length..] != second[..length]
161    {
162        return Err(GibsonError::Component {
163            component,
164            reason: "declared existing suffix/prefix overlap does not match",
165        });
166    }
167    Ok(())
168}
169
170#[cfg(test)]
171mod tests {
172    use super::*;
173    use crate::{DnaSeq, ligation::Orientation};
174    fn source() -> SequenceRecord {
175        let mut state = 1_234_567_u64;
176        let sequence: String = (0..240)
177            .map(|_| {
178                state ^= state << 13;
179                state ^= state >> 7;
180                state ^= state << 17;
181                char::from(b"ACGT"[(state % 4) as usize])
182            })
183            .collect();
184        SequenceRecord::new(
185            "synthetic",
186            DnaSeq::new(sequence).unwrap(),
187            Topology::Circular,
188            vec![],
189            vec![],
190        )
191        .unwrap()
192    }
193    fn selections() -> Vec<CoreSelection> {
194        vec![
195            CoreSelection {
196                input: 1,
197                start: 0,
198                length: 140,
199                orientation: Orientation::Forward,
200            },
201            CoreSelection {
202                input: 1,
203                start: 115,
204                length: 150,
205                orientation: Orientation::Forward,
206            },
207        ]
208    }
209    #[test]
210    fn merges_existing_homology_and_preserves_circular_source_associations() {
211        let r = source();
212        let out = assemble(
213            std::slice::from_ref(&r),
214            &selections(),
215            Topology::Circular,
216            &[25, 25],
217        )
218        .unwrap();
219        assert_eq!(out.product_sequence_5to3, r.sequence().as_str());
220        assert_eq!(out.components[1].product_start, 115);
221        assert!(out.components[1].wraps_origin);
222        assert_eq!(out.junctions[1].product_start, 0);
223        let mut reverse = selections();
224        reverse.reverse();
225        for s in &mut reverse {
226            s.orientation = Orientation::Reverse;
227        }
228        let out = assemble(
229            std::slice::from_ref(&r),
230            &reverse,
231            Topology::Circular,
232            &[25, 25],
233        )
234        .unwrap();
235        let expected = crate::digest::reverse_complement(r.sequence().as_str());
236        assert_eq!(out.product_sequence_5to3.len(), 240);
237        assert!(expected.repeat(2).contains(&out.product_sequence_5to3));
238    }
239    #[test]
240    fn single_linear_view_needs_no_overlap_or_pcr_length() {
241        let r = source();
242        let selected = [CoreSelection {
243            input: 1,
244            start: 0,
245            length: 1,
246            orientation: Orientation::Forward,
247        }];
248        let out = assemble(std::slice::from_ref(&r), &selected, Topology::Linear, &[]).unwrap();
249        assert_eq!(out.product_sequence_5to3, &r.sequence().as_str()[..1]);
250        assert_eq!(out.junctions.len(), 0, "no junction without an overlap");
251        assert!(!out.components[0].wraps_origin);
252    }
253
254    #[test]
255    fn rejects_mismatch_bad_counts_exhausted_interiors_and_repeated_homology() {
256        assert!(assemble(&[source()], &selections(), Topology::Circular, &[24, 25]).is_err());
257        assert!(validate_request(1, &selections(), Topology::Circular, &[25]).is_err());
258        assert!(validate_request(0, &selections(), Topology::Circular, &[25, 25]).is_err());
259        let mut s = selections();
260        s[0].length = 50;
261        assert!(validate_request(1, &s, Topology::Circular, &[25, 25]).is_err());
262        s[0].length = usize::MAX;
263        assert!(validate_request(1, &s, Topology::Circular, &[25, 25]).is_err());
264        let repeated = SequenceRecord::new(
265            "synthetic",
266            DnaSeq::new("A".repeat(240)).unwrap(),
267            Topology::Circular,
268            vec![],
269            vec![],
270        )
271        .unwrap();
272        assert!(assemble(&[repeated], &selections(), Topology::Circular, &[25, 25]).is_err());
273    }
274}
dnagent_domain/existing_overlaps.rs

dnagent_domain/
existing_overlaps.rs
