dnagent_domain/
primer_optimisation.rs

1//! Deterministic bounded optimisation: NN annealing Tm and sequence-only structure screens.
2//! Numerical parameter provenance and scope: docs/primer-optimisation.md.
3use crate::digest::reverse_complement;
4use crate::gibson::{
5    CoreSelection, GibsonError, GibsonReport, Preparation, PrimerCandidate, core,
6    design_with_lengths, primer, unique_duplex_site, validate_request, validate_templates,
7};
8use crate::{SequenceRecord, Topology};
9use serde::{Deserialize, Serialize};
10
11#[derive(Debug, Clone, Serialize, Deserialize)]
12#[serde(deny_unknown_fields)]
13pub struct Solution {
14    pub sodium_mm: f64,
15    pub potassium_mm: f64,
16    pub tris_mm: f64,
17    pub magnesium_mm: f64,
18    pub dntp_mm: f64,
19    pub primer_nm: f64,
20}
21impl Solution {
22    fn equivalent_sodium(&self) -> f64 {
23        self.sodium_mm
24            + self.potassium_mm
25            + self.tris_mm / 2.0
26            + 120.0 * (self.magnesium_mm - self.dntp_mm).max(0.0).sqrt()
27    }
28    fn validate(&self) -> bool {
29        [
30            (self.sodium_mm, 500.0),
31            (self.potassium_mm, 500.0),
32            (self.tris_mm, 200.0),
33            (self.magnesium_mm, 10.0),
34            (self.dntp_mm, 10.0),
35        ]
36        .iter()
37        .all(|(v, max)| v.is_finite() && (0.0..=*max).contains(v))
38            && self.primer_nm.is_finite()
39            && (1.0..=5000.0).contains(&self.primer_nm)
40            && self.equivalent_sodium() > 0.0
41    }
42}
43
44#[derive(Debug, Clone, Serialize, Deserialize)]
45#[serde(deny_unknown_fields)]
46pub struct PrimerConstraints {
47    pub min_length: usize,
48    pub max_length: usize,
49    pub min_tm_c: f64,
50    pub max_tm_c: f64,
51    pub target_tm_c: f64,
52    pub max_pair_tm_difference_c: f64,
53    pub min_gc_fraction: f64,
54    pub max_gc_fraction: f64,
55    pub max_hairpin_stem: usize,
56    pub max_dimer_run: usize,
57    pub max_three_prime_run: usize,
58    pub solution: Solution,
59}
60impl PrimerConstraints {
61    pub fn validate(&self) -> Result<(), GibsonError> {
62        if !(18..=40).contains(&self.min_length)
63            || !(self.min_length..=40).contains(&self.max_length)
64            || !self.solution.validate()
65            || ![
66                self.min_tm_c,
67                self.max_tm_c,
68                self.target_tm_c,
69                self.max_pair_tm_difference_c,
70                self.min_gc_fraction,
71                self.max_gc_fraction,
72            ]
73            .iter()
74            .all(|x| x.is_finite())
75            || !(0.0..=100.0).contains(&self.min_tm_c)
76            || !(self.min_tm_c..=100.0).contains(&self.max_tm_c)
77            || !(self.min_tm_c..=self.max_tm_c).contains(&self.target_tm_c)
78            || !(0.0..=20.0).contains(&self.max_pair_tm_difference_c)
79            || !(0.0..=1.0).contains(&self.min_gc_fraction)
80            || !(self.min_gc_fraction..=1.0).contains(&self.max_gc_fraction)
81            || self.max_hairpin_stem > 20
82            || self.max_dimer_run > 40
83            || self.max_three_prime_run > 20
84        {
85            return Err(GibsonError::Invalid(
86                "invalid primer optimisation constraints or solution concentrations",
87            ));
88        }
89        Ok(())
90    }
91}
92
93#[derive(Debug, Clone, Copy, PartialEq, Eq, Serialize)]
94pub struct DimerScreen {
95    pub longest_run: usize,
96    pub longest_three_prime_run: usize,
97}
98#[derive(Debug, Clone, Serialize)]
99pub struct AssessedPrimer {
100    pub primer: PrimerCandidate,
101    pub annealing_tm_c: f64,
102    pub hairpin_stem: usize,
103    pub self_dimer: DimerScreen,
104}
105#[derive(Debug, Serialize)]
106pub struct OptimisedPair {
107    /// 1-based component this pair amplifies; provided fragments have no pair.
108    pub component: usize,
109    pub forward: AssessedPrimer,
110    pub reverse: AssessedPrimer,
111    pub heterodimer: DimerScreen,
112    pub score: f64,
113    pub candidate_pairs_examined: usize,
114    pub feasible_pairs: usize,
115}
116#[derive(Debug, Serialize)]
117pub struct OptimisedGibson {
118    pub design: GibsonReport,
119    pub constraints: PrimerConstraints,
120    pub pairs: Vec<OptimisedPair>,
121    pub assumptions: Vec<&'static str>,
122}
123
124pub fn optimise(
125    records: &[SequenceRecord],
126    selected: &[CoreSelection],
127    preparation: &[Preparation],
128    topology: Topology,
129    overlap: usize,
130    constraints: &PrimerConstraints,
131) -> Result<OptimisedGibson, GibsonError> {
132    constraints.validate()?;
133    validate_request(records.len(), selected, overlap, constraints.min_length)?;
134    validate_templates(records)?;
135    let cores = selected
136        .iter()
137        .map(|s| core(&records[s.input - 1], s))
138        .collect::<Result<Vec<_>, _>>()?;
139    // Tails follow the same rule as the design: the upstream fragment's reverse primer
140    // carries the overlap, unless that fragment is provided and so cannot be extended,
141    // in which case the downstream fragment's forward primer does.
142    let joined = |i: usize| i + 1 < selected.len() || topology == Topology::Circular;
143    let mut forward_tails = vec![String::new(); selected.len()];
144    let mut reverse_tails = vec![String::new(); selected.len()];
145    for i in 0..selected.len() {
146        if !joined(i) {
147            continue;
148        }
149        let next = (i + 1) % selected.len();
150        match (preparation[i], preparation[next]) {
151            (Preparation::Provided, Preparation::Pcr) => {
152                forward_tails[next].clone_from(&cores[i][cores[i].len() - overlap..].to_owned());
153            }
154            (Preparation::Provided, Preparation::Provided) => {
155                return Err(GibsonError::Component {
156                    component: i + 1,
157                    reason: "both fragments at this junction are provided, so no primer can add the overlap; give them existing overlaps instead",
158                });
159            }
160            (Preparation::Pcr, _) => {
161                reverse_tails[i] = reverse_complement(&cores[next][..overlap]);
162            }
163        }
164    }
165    let mut pairs = Vec::new();
166    let mut lengths = Vec::with_capacity(selected.len());
167    for (i, (selection, sequence)) in selected.iter().zip(&cores).enumerate() {
168        if preparation[i] == Preparation::Provided {
169            // Not amplified: no primers to choose, and no annealing lengths to report.
170            lengths.push((0, 0));
171            continue;
172        }
173        let template = &records[selection.input - 1];
174        let forward = candidates(sequence, &forward_tails[i], false, template, constraints)?;
175        let reverse = candidates(sequence, &reverse_tails[i], true, template, constraints)?;
176        let pair = choose_pair(&forward, &reverse, sequence.len(), constraints).ok_or(GibsonError::Component { component: i + 1, reason: "no primer pair satisfies exact-site, Tm, GC and sequence-structure constraints; no fallback design returned" })?;
177        lengths.push((
178            pair.forward.primer.annealing_sequence_5to3.len(),
179            pair.reverse.primer.annealing_sequence_5to3.len(),
180        ));
181        pairs.push(OptimisedPair {
182            component: i + 1,
183            ..pair
184        });
185    }
186    let design = design_with_lengths(records, selected, preparation, topology, overlap, &lengths)?;
187    Ok(OptimisedGibson {
188        design,
189        constraints: constraints.clone(),
190        pairs,
191        assumptions: vec![
192            "Annealing Tm: Allawi/SantaLucia 1997 DNA NN parameters, SantaLucia 1998 entropy salt correction, von Ahsen 2001 sodium-equivalent Mg/dNTP approximation; primer concentration in excess, complementary template concentration neglected",
193            "Tm describes the template-matching annealing segment, not the full tailed oligo; no recommended PCR annealing temperature is inferred",
194            "Hairpin and dimer checks use full oligos but only contiguous Watson-Crick sequence complementarity, not structure energies, bulges, mismatches or temperature-dependent folding",
195            "Hairpin stems require at least a three-base unpaired loop; dimer screens report all contiguous runs and runs reaching either oligo's 3-prime end",
196            "Each component primer pair is optimised separately for separate PCR reactions; cross-component multiplex primer interactions are not evaluated",
197            "Among feasible pairs minimise summed absolute deviations from target Tm; ties prefer shorter total annealing length, then shorter forward length, then reverse length",
198            "Exact annealing-site uniqueness is tested on the full source in both orientations; off-template and mismatched binding are not evaluated",
199            "Thermodynamic estimates and sequence screens do not establish ordering readiness or experimental validation; bounds are implementation constraints, not calibration guarantees",
200        ],
201    })
202}
203
204fn candidates(
205    sequence: &str,
206    tail: &str,
207    reverse: bool,
208    template: &SequenceRecord,
209    c: &PrimerConstraints,
210) -> Result<Vec<AssessedPrimer>, GibsonError> {
211    let mut result = Vec::new();
212    for length in c.min_length..=c.max_length.min(sequence.len()) {
213        let anneal = if reverse {
214            reverse_complement(&sequence[sequence.len() - length..])
215        } else {
216            sequence[..length].into()
217        };
218        if !unique_duplex_site(template.sequence().as_str(), &anneal, template.topology()) {
219            continue;
220        }
221        let p = primer(&anneal, tail);
222        let tm = annealing_tm(&anneal, &c.solution)?;
223        let gc = f64::from(u32::try_from(p.annealing_gc_bases).expect("bounded oligo"))
224            / f64::from(u32::try_from(length).expect("bounded oligo"));
225        let hairpin = hairpin_stem(&p.sequence_5to3);
226        let dimer = dimer_screen(&p.sequence_5to3, &p.sequence_5to3);
227        if (c.min_tm_c..=c.max_tm_c).contains(&tm)
228            && (c.min_gc_fraction..=c.max_gc_fraction).contains(&gc)
229            && hairpin <= c.max_hairpin_stem
230            && passes(dimer, c)
231        {
232            result.push(AssessedPrimer {
233                primer: p,
234                annealing_tm_c: tm,
235                hairpin_stem: hairpin,
236                self_dimer: dimer,
237            });
238        }
239    }
240    Ok(result)
241}
242fn passes(d: DimerScreen, c: &PrimerConstraints) -> bool {
243    d.longest_run <= c.max_dimer_run && d.longest_three_prime_run <= c.max_three_prime_run
244}
245
246fn choose_pair(
247    forward: &[AssessedPrimer],
248    reverse: &[AssessedPrimer],
249    length: usize,
250    c: &PrimerConstraints,
251) -> Option<OptimisedPair> {
252    let mut best: Option<OptimisedPair> = None;
253    let mut examined = 0;
254    let mut feasible = 0;
255    for f in forward {
256        for r in reverse {
257            examined += 1;
258            let fl = f.primer.annealing_sequence_5to3.len();
259            let rl = r.primer.annealing_sequence_5to3.len();
260            let dimer = dimer_screen(&f.primer.sequence_5to3, &r.primer.sequence_5to3);
261            if fl + rl > length
262                || (f.annealing_tm_c - r.annealing_tm_c).abs() > c.max_pair_tm_difference_c
263                || !passes(dimer, c)
264            {
265                continue;
266            }
267            feasible += 1;
268            let score =
269                (f.annealing_tm_c - c.target_tm_c).abs() + (r.annealing_tm_c - c.target_tm_c).abs();
270            let key = (score, fl + rl, fl, rl);
271            if best.as_ref().is_none_or(|b| {
272                let bf = b.forward.primer.annealing_sequence_5to3.len();
273                let br = b.reverse.primer.annealing_sequence_5to3.len();
274                key < (b.score, bf + br, bf, br)
275            }) {
276                best = Some(OptimisedPair {
277                    component: 0, // set by the caller, which knows the component
278                    forward: f.clone(),
279                    reverse: r.clone(),
280                    heterodimer: dimer,
281                    score,
282                    candidate_pairs_examined: 0,
283                    feasible_pairs: 0,
284                });
285            }
286        }
287    }
288    best.map(|mut b| {
289        b.candidate_pairs_examined = examined;
290        b.feasible_pairs = feasible;
291        b
292    })
293}
294
295/// Exact-match, non-self-complementary annealing duplex model. Concentrations are explicit.
296pub fn annealing_tm(sequence: &str, solution: &Solution) -> Result<f64, GibsonError> {
297    if !solution.validate()
298        || !(18..=40).contains(&sequence.len())
299        || !sequence.bytes().all(|b| b"ACGT".contains(&b))
300        || reverse_complement(sequence) == sequence
301    {
302        return Err(GibsonError::Invalid(
303            "Tm requires 18..=40 unambiguous non-self-complementary bases and valid solution conditions",
304        ));
305    }
306    let (mut h, mut s) = (0.0, 0.0);
307    for b in [
308        sequence.as_bytes()[0],
309        sequence.as_bytes()[sequence.len() - 1],
310    ] {
311        let (dh, ds) = if matches!(b, b'A' | b'T') {
312            (2.3, 4.1)
313        } else {
314            (0.1, -2.8)
315        };
316        h += dh;
317        s += ds;
318    }
319    for pair in sequence.as_bytes().windows(2) {
320        let (dh, ds) = match pair {
321            b"AA" | b"TT" => (-7.9, -22.2),
322            b"AT" => (-7.2, -20.4),
323            b"TA" => (-7.2, -21.3),
324            b"CA" | b"TG" => (-8.5, -22.7),
325            b"GT" | b"AC" => (-8.4, -22.4),
326            b"CT" | b"AG" => (-7.8, -21.0),
327            b"GA" | b"TC" => (-8.2, -22.2),
328            b"CG" => (-10.6, -27.2),
329            b"GC" => (-9.8, -24.4),
330            b"GG" | b"CC" => (-8.0, -19.9),
331            _ => unreachable!("validated ACGT"),
332        };
333        h += dh;
334        s += ds;
335    }
336    s += 0.368
337        * f64::from(u32::try_from(sequence.len() - 1).expect("bounded oligo"))
338        * (solution.equivalent_sodium() * 1e-3).ln();
339    Ok(1000.0 * h / (s + 1.987 * (solution.primer_nm * 1e-9).ln()) - 273.15)
340}
341
342// Private screens only receive bounded, validated ACGT oligos.
343pub(crate) fn dimer_screen(a: &str, b: &str) -> DimerScreen {
344    let b = reverse_complement(b);
345    let mut previous = vec![0; b.len() + 1];
346    let mut result = DimerScreen {
347        longest_run: 0,
348        longest_three_prime_run: 0,
349    };
350    for (i, base) in a.bytes().enumerate() {
351        let mut row = vec![0; b.len() + 1];
352        for (j, other) in b.bytes().enumerate() {
353            if base == other {
354                let run = previous[j] + 1;
355                row[j + 1] = run;
356                result.longest_run = result.longest_run.max(run);
357                if i + 1 == a.len() || j + 1 == run {
358                    result.longest_three_prime_run = result.longest_three_prime_run.max(run);
359                }
360            }
361        }
362        previous = row;
363    }
364    result
365}
366pub(crate) fn hairpin_stem(s: &str) -> usize {
367    let rc = reverse_complement(s);
368    let mut best = 0;
369    for left in 0..s.len() {
370        for right in left + 4..s.len() {
371            let mut k = 0;
372            while left + k + 3 < right - k
373                && s.as_bytes()[left + k] == rc.as_bytes()[s.len() - 1 - right + k]
374            {
375                k += 1;
376            }
377            best = best.max(k);
378        }
379    }
380    best
381}
382
383#[cfg(test)]
384mod tests {
385    use super::*;
386    fn solution() -> Solution {
387        Solution {
388            sodium_mm: 50.0,
389            potassium_mm: 0.0,
390            tris_mm: 0.0,
391            magnesium_mm: 1.5,
392            dntp_mm: 0.2,
393            primer_nm: 250.0,
394        }
395    }
396    #[test]
397    fn nn_tm_matches_pinned_biopython_reference_and_rc() {
398        let seq = "ACGTTGCAATGCCGTAGCTAGC";
399        let tm = annealing_tm(seq, &solution()).unwrap();
400        assert!((tm - 68.544_083_450_990_97).abs() < 1e-10);
401        assert!((tm - annealing_tm(&reverse_complement(seq), &solution()).unwrap()).abs() < 1e-10);
402        assert!(annealing_tm("ACGT", &solution()).is_err());
403        assert!(annealing_tm("ACGTTGCAATGCCGTAGCTAGN", &solution()).is_err());
404        let mut invalid = solution();
405        invalid.primer_nm = f64::NAN;
406        assert!(annealing_tm(seq, &invalid).is_err());
407        invalid = solution();
408        invalid.sodium_mm = 0.0;
409        invalid.magnesium_mm = 0.0;
410        assert!(annealing_tm(seq, &invalid).is_err());
411    }
412    #[test]
413    fn a_longer_primer_can_resolve_a_repeated_minimum_length_site() {
414        use crate::{DnaSeq, ligation::Orientation};
415        let mut state = 1_234_567_u64;
416        let mut bases: Vec<_> = (0..180)
417            .map(|_| {
418                state ^= state << 13;
419                state ^= state >> 7;
420                state ^= state << 17;
421                b"ACGT"[(state % 4) as usize]
422            })
423            .collect();
424        bases.copy_within(0..18, 120);
425        bases[138] = if bases[18] == b'A' { b'C' } else { b'A' };
426        let record = SequenceRecord::new(
427            "synthetic",
428            DnaSeq::new(String::from_utf8(bases).unwrap()).unwrap(),
429            Topology::Linear,
430            vec![],
431            vec![],
432        )
433        .unwrap();
434        let c = PrimerConstraints {
435            min_length: 18,
436            max_length: 19,
437            min_tm_c: 0.0,
438            max_tm_c: 100.0,
439            target_tm_c: 60.0,
440            max_pair_tm_difference_c: 20.0,
441            min_gc_fraction: 0.0,
442            max_gc_fraction: 1.0,
443            max_hairpin_stem: 20,
444            max_dimer_run: 40,
445            max_three_prime_run: 20,
446            solution: solution(),
447        };
448        let selected = [CoreSelection {
449            input: 1,
450            start: 0,
451            length: 90,
452            orientation: Orientation::Forward,
453        }];
454        assert!(!unique_duplex_site(
455            record.sequence().as_str(),
456            &record.sequence().as_str()[..18],
457            Topology::Linear
458        ));
459        let result = optimise(
460            &[record],
461            &selected,
462            &[Preparation::Pcr],
463            Topology::Linear,
464            25,
465            &c,
466        )
467        .unwrap();
468        assert_eq!(
469            result.pairs[0].forward.primer.annealing_sequence_5to3.len(),
470            19
471        );
472    }
473
474    #[test]
475    fn structure_screens_distinguish_internal_and_three_prime_runs() {
476        assert_eq!(
477            dimer_screen("AAAAA", "TTTTT"),
478            DimerScreen {
479                longest_run: 5,
480                longest_three_prime_run: 5
481            }
482        );
483        assert_eq!(
484            dimer_screen("GAAAAC", "ATTTTG"),
485            DimerScreen {
486                longest_run: 4,
487                longest_three_prime_run: 1
488            }
489        );
490        assert_eq!(
491            dimer_screen("GAAAAC", "ATTTTG"),
492            dimer_screen("ATTTTG", "GAAAAC")
493        );
494        assert_eq!(hairpin_stem("GCGAAACGC"), 3);
495        assert_eq!(hairpin_stem("GCGCGC"), 1);
496        assert_eq!(hairpin_stem("AAAAAAAAAA"), 0);
497    }
498}
dnagent_domain/primer_optimisation.rs

dnagent_domain/
primer_optimisation.rs
