#![allow(non_snake_case)]
use na_seq::{
Nucleotide::{self, A, C, G, T},
calc_gc,
};
use crate::primer::{IonConcentrations, MIN_PRIMER_LEN};
const R: f32 = 1.987;
fn _dH_dS_tmm(nts: (Nucleotide, Nucleotide)) -> Option<(f32, f32)> {
match nts {
(A, A) => Some((-7.6, -21.3)),
(A, T) => Some((-7.2, -20.4)),
(T, A) => Some((-7.2, -21.3)),
(C, A) => Some((-8.5, -22.7)),
(C, G) => Some((-10.6, -27.2)),
(G, A) => Some((-8.2, -22.2)),
(G, C) => Some((-9.8, -24.4)),
(G, T) => Some((-8.4, -22.4)),
(G, G) => Some((-8.0, -19.9)),
_ => None,
}
}
fn _dH_dS_imm(nts: (Nucleotide, Nucleotide)) -> Option<(f32, f32)> {
match nts {
(A, A) => Some((-7.6, -21.3)),
(A, T) => Some((-7.2, -20.4)),
(T, A) => Some((-7.2, -21.3)),
(C, A) => Some((-8.5, -22.7)),
(C, G) => Some((-10.6, -27.2)),
(G, A) => Some((-8.2, -22.2)),
(G, C) => Some((-9.8, -24.4)),
(G, T) => Some((-8.4, -22.4)),
(G, G) => Some((-8.0, -19.9)),
_ => None,
}
}
fn _dH_dS_de(nts: (Nucleotide, Nucleotide)) -> Option<(f32, f32)> {
match nts {
(A, A) => Some((0.2, 2.3)),
(A, T) => Some((-7.2, -20.4)),
(T, A) => Some((-7.2, -21.3)),
(C, A) => Some((-8.5, -22.7)),
(C, G) => Some((-10.6, -27.2)),
(G, A) => Some((-8.2, -22.2)),
(G, C) => Some((-9.8, -24.4)),
(G, T) => Some((-8.4, -22.4)),
(G, G) => Some((-8.0, -19.9)),
_ => None,
}
}
fn dH_dS_neighbors(neighbors: (Nucleotide, Nucleotide)) -> (f32, f32) {
match neighbors {
(A, A) | (T, T) => (-7.6, -21.3),
(A, T) => (-7.2, -20.4),
(T, A) => (-7.2, -21.3),
(C, A) | (T, G) => (-8.5, -22.7),
(G, T) | (A, C) => (-8.4, -22.4),
(C, T) | (A, G) => (-7.8, -21.0),
(G, A) | (T, C) => (-8.2, -22.2),
(C, G) => (-10.6, -27.2),
(G, C) => (-9.8, -24.4),
(G, G) | (C, C) => (-8.0, -19.9),
}
}
fn salt_correction(seq: &[Nucleotide], ion: &IonConcentrations) -> Option<f32> {
let method = 5;
let tris = 0.;
if (5..=7).contains(&method) && seq.is_empty() {
return None;
}
let mut corr = 0.0;
if method == 0 {
return Some(corr);
}
let mut mon = ion.monovalent + tris / 2.0;
let mon_molar = mon * 1e-3;
let mg_molar = ion.divalent * 1e-3;
if (ion.monovalent > 0.0 || ion.divalent > 0.0 || tris > 0.0 || ion.dntp > 0.0)
&& method != 7
&& ion.dntp < ion.divalent
{
mon += 120.0 * (ion.divalent - ion.dntp).sqrt();
}
if (1..=6).contains(&method) && mon_molar == 0.0 {
return None;
}
match method {
1 => corr = 16.6 * mon_molar.log10(),
2 => corr = 16.6 * (mon_molar / (1.0 + 0.7 * mon_molar)).log10(),
3 => corr = 12.5 * mon_molar.log10(),
4 => corr = 11.7 * mon_molar.log10(),
5 => {
corr = 0.368 * (seq.len() as f32 - 1.0) * mon_molar.ln();
}
6 => {
let gc_fraction = calc_gc(seq);
corr = ((4.29 * gc_fraction - 3.95) * 1e-5 * mon_molar.ln())
+ 9.40e-6 * mon_molar.ln().powi(2);
}
7 => {
let (mut a, b, c, mut d, e, f, mut g) = (3.92, -0.911, 6.26, 1.42, -48.2, 52.5, 8.31);
let mut mg_free = mg_molar;
if ion.dntp > 0.0 {
let dntps_molar = ion.dntp * 1e-3;
let ka = 3e4;
mg_free = (-(ka * dntps_molar - ka * mg_free + 1.0)
+ ((ka * dntps_molar - ka * mg_free + 1.0).powi(2) + 4.0 * ka * mg_free)
.sqrt())
/ (2.0 * ka);
}
if mon > 0.0 {
let r = (mg_free).sqrt() / mon_molar;
if r < 0.22 {
let gc_fraction = calc_gc(seq);
corr = (4.29 * gc_fraction - 3.95) * 1e-5 * mon_molar.ln()
+ 9.40e-6 * mon_molar.ln().powi(2);
return Some(corr);
} else if r < 6.0 {
a = 3.92 * (0.843 - 0.352 * mon_molar.sqrt() * mon_molar.ln());
d = 1.42
* (1.279 - 4.03e-3 * mon_molar.ln() - 8.03e-3 * mon_molar.ln().powi(2));
g = 8.31 * (0.486 - 0.258 * mon_molar.ln() + 5.25e-3 * mon_molar.ln().powi(3));
}
let gc_fraction = calc_gc(seq);
corr = (a
+ b * mg_free.ln()
+ gc_fraction * (c + d * mg_free.ln())
+ (1.0 / (2.0 * (seq.len() as f32 - 1.0)))
* (e + f * mg_free.ln() + g * mg_free.ln().powi(2)))
* 1e-5;
}
}
_ => return None,
}
Some(corr)
}
pub fn calc_tm(seq: &[Nucleotide], ion_concentrations: &IonConcentrations) -> Option<f32> {
if seq.len() < MIN_PRIMER_LEN {
return None;
}
let mut dH = 0.2;
let mut dS = -5.7;
if calc_gc(seq) < 0.001 {
dH += 2.2;
dS += 6.9;
}
{
let term_pair = vec![seq[0], seq[seq.len() - 1]];
let mut at_term_count = 0;
for nt in term_pair {
if nt == A || nt == T {
at_term_count += 1;
}
}
dH += 2.2 * at_term_count as f32;
dS += 6.9 * at_term_count as f32;
}
for (i, nt) in seq.iter().enumerate() {
if i + 1 >= seq.len() {
break;
}
let neighbors = (*nt, seq[i + 1]);
let (dH_nn, dS_nn) = dH_dS_neighbors(neighbors);
dH += dH_nn;
dS += dS_nn;
}
let C_T = ion_concentrations.primer * 1.0e-9;
if let Some(sc) = salt_correction(seq, ion_concentrations) {
dS += sc;
} else {
eprintln!("Error calculating salt correction.");
}
let result = (1_000. * dH) / (dS + R * (C_T / 2.).ln()) - 273.15;
Some(result)
}