use super::{Domain, Strand};
use crate::Design;
use cadnano_format::{Cadnano, VStrand};
use std::collections::{HashMap, HashSet};
use std::io::Write;
use std::path::Path;
const SKIP_STEP: isize = 48;
const NO_HELIX: usize = std::usize::MAX;
pub enum GridType {
Other,
Square,
Honeycomb,
}
impl<StrandLabel, DomainLabel> Design<StrandLabel, DomainLabel> {
pub fn to_cadnano<W: Write>(
&mut self,
w: W,
try_do_skip: bool,
grid_type: GridType,
) -> Result<(), failure::Error> {
self.clean_up();
let mut xovers = if try_do_skip {
self.get_cross_overs()
} else {
Vec::new()
};
let mut vstrands = Vec::new();
let mut max_x = 0;
let mut min_x = 0;
let mut max_y = 0;
let mut min_helix_y = 0.;
let mut min_helix_z = 0.;
for s in self.strands.iter() {
for dom in s.domains.iter() {
if dom.helix < 0 {
continue;
}
max_x = max_x.max(dom.start).max(dom.end - 1);
min_x = min_x.min(dom.start).min(dom.end - 1);
max_y = max_y.max(dom.helix as usize);
let helix = &self.helices[dom.helix as usize];
if helix.origin.y < min_helix_y {
min_helix_y = helix.origin.y;
}
if helix.origin.z < min_helix_z {
min_helix_z = helix.origin.z;
}
}
}
if try_do_skip {
for s in &self.strands {
if let Some((h, x, b)) = s.first_nucl() {
xovers.push((h, x, b, h, x, b));
}
if let Some((h, x, b)) = s.last_nucl() {
xovers.push((h, x, b, h, x, b));
}
}
}
for i in 0..xovers.len() {
xovers[i].1 -= min_x;
xovers[i].4 -= min_x;
}
max_x -= min_x;
let skips = get_skip_points(&xovers, max_x);
let do_skip = try_do_skip && skips.len() > 0;
max_x = if do_skip {
convert_skip(max_x, &skips)
} else {
max_x
};
let max_x = match grid_type {
GridType::Honeycomb => ((1 + max_x / 21) * 21) as usize,
_ => ((1 + max_x / 32) * 32) as usize,
};
println!("max_x {}", max_x);
let mut father = make_group(self, max_y);
println!("father {:?}", father);
let graph = make_graph(self, max_y, &mut father);
println!("graph {:?}", graph);
let parity = color_graph(&graph, max_y, &mut father);
let mut odd_num = 1;
let mut even_num = 0;
let mut helices = Vec::new();
let mut vstrand_idx = Vec::new();
let mut last_parit = 0;
let mut rows = Vec::new();
let mut cols = Vec::new();
for (i, h) in self.helices.iter().enumerate().take(max_y + 1) {
let num;
let idx;
match grid_type {
GridType::Other => {
if !parity[i] {
num = odd_num;
odd_num += 2;
if last_parit == 1 {
even_num += 2;
}
last_parit = 1;
idx = 2 * (num as isize / 2);
} else {
num = even_num;
even_num += 2;
if last_parit == 0 {
odd_num += 2;
}
last_parit = 0;
idx = 2 * (num as isize / 2) + 1;
};
}
GridType::Square => {
let param = &self.parameters.expect("self.parameters");
let r = param.inter_helix_gap + param.helix_radius * 2.;
let row = (h.origin.y - min_helix_y).div_euclid(r) as isize;
let col = (h.origin.z - min_helix_z).div_euclid(r) as isize;
if row % 2 == col % 2 {
num = even_num;
even_num += 2;
} else {
num = odd_num;
odd_num += 2;
}
rows.push(row);
cols.push(col);
idx = i as isize;
}
GridType::Honeycomb => {
let param = &self.parameters.expect("self.parameters");
let r = param.helix_radius + param.inter_helix_gap / 2.;
let row = (h.origin.y - min_helix_y).div_euclid(3. * r) as isize;
let col = (h.origin.z - min_helix_z).div_euclid(3f64.sqrt() * r) as isize;
if row % 2 == col % 2 {
num = even_num;
even_num += 2;
} else {
num = odd_num;
odd_num += 2;
}
rows.push(row);
cols.push(col);
idx = i as isize;
}
}
println!("helix {} -> num {}", i, num);
helices.push(num as isize);
vstrand_idx.push(idx);
}
let max = if rows.len() == 0 {
odd_num.max(even_num)
} else {
rows.len() - 1
};
odd_num = 1;
even_num = 0;
for i in 0..(max + 1) {
let num;
let row;
let col;
if rows.len() == 0 {
if i % 2 == 0 {
num = odd_num;
odd_num += 2;
row = 2 * (i as isize / 2);
} else {
num = even_num;
even_num += 2;
row = 2 * (i as isize / 2) + 1;
};
col = 1;
} else {
num = helices[i] as usize;
row = rows[i];
col = cols[i];
}
vstrands.push(cadnano_format::VStrand {
col,
loop_: vec![0; max_x],
num: num as isize,
row,
scaf: vec![(-1, -1, -1, -1); max_x],
scaf_loop: vec![0; max_x],
skip: vec![0; max_x],
stap: vec![(-1, -1, -1, -1); max_x],
stap_loop: vec![0; max_x],
stap_colors: vec![],
})
}
for s in self.strands.iter() {
let mut prev: Option<(usize, usize)> = None;
let mut cyclic = s.cyclic;
if s.domains.len() > 1 {
let first_dom = s.domains.first().unwrap();
let last_dom = s.domains.last().unwrap();
if first_dom.helix == last_dom.helix && first_dom.right == last_dom.right {
if (first_dom.right && last_dom.end - 1 == last_dom.start)
|| (!first_dom.right && last_dom.start == first_dom.end - 1)
{
cyclic = true;
println!(
"cyclic, {} {} {}",
first_dom.helix, first_dom.start, first_dom.right
);
}
}
}
for (ndom, dom) in s.domains.iter().enumerate() {
let mut do_convert = true;
if dom.length() <= 0 {
continue;
}
let max_dom = s.domains.len();
let x_range = if dom.right {
dom.start..dom.end
} else {
(-(dom.end - 1))..(-(dom.start - 1))
};
for x_ in x_range {
let is_scaf = (helices[dom.helix as usize] % 2 == 0) ^ !dom.right;
let mut a_helix = -1;
let mut c_helix = -1;
let x = if dom.right {
(x_ - min_x) as usize
} else {
(-x_ - min_x) as usize
};
let mut result: (isize, isize, isize, isize) = (-1, -1, -1, -1);
if let Some((ref prev_helix, ref prev_x)) = prev {
result.0 = helices[*prev_helix] as isize;
a_helix = *prev_helix as isize;
result.1 = *prev_x as isize;
}
let last = if dom.right { dom.end - 1 } else { -dom.start };
if x_ != last {
result.2 = helices[dom.helix as usize] as isize;
c_helix = dom.helix;
result.3 = if dom.right {
(x + 1) as isize
} else {
(x - 1) as isize
};
} else {
if ndom + 1 < max_dom {
result.2 = helices[s.domains[ndom + 1].helix as usize] as isize;
c_helix = s.domains[ndom + 1].helix;
result.3 = if s.domains[ndom + 1].right {
s.domains[ndom + 1].start - min_x
} else {
s.domains[ndom + 1].end - 1 - min_x
}
} else if cyclic {
let self_pos = if do_skip {
convert_skip(x as isize, &skips)
} else {
x as isize
} as usize;
result.2 = if is_scaf {
vstrands[vstrand_idx[dom.helix as usize] as usize].scaf[self_pos].2
} else {
vstrands[vstrand_idx[dom.helix as usize] as usize].stap[self_pos].2
};
result.3 = if is_scaf {
vstrands[vstrand_idx[dom.helix as usize] as usize].scaf[self_pos].3
} else {
vstrands[vstrand_idx[dom.helix as usize] as usize].stap[self_pos].3
};
do_convert = false;
}
}
prev = Some((dom.helix as usize, x));
let self_pos = if do_skip {
convert_skip(x as isize, &skips)
} else {
x as isize
};
let result = if do_skip && do_convert {
convert_result(result, x as isize, helices[dom.helix as usize], &skips)
} else if do_skip {
let old = (result.2, result.3);
let convert =
convert_result(result, x as isize, helices[dom.helix as usize], &skips);
do_convert = true;
(convert.0, convert.1, old.0, old.1)
} else {
result
};
let (_, b, _, d) = result;
println!("result {:?}; self_pos {:?}", result, self_pos);
if is_scaf {
println!(
"changing {} scaf {} {:?}",
vstrand_idx[dom.helix as usize], self_pos, result
);
vstrands[vstrand_idx[dom.helix as usize] as usize].scaf
[self_pos as usize] = result;
if a_helix > -1 {
vstrands[vstrand_idx[a_helix as usize] as usize].scaf[b as usize].2 =
helices[dom.helix as usize];
vstrands[vstrand_idx[a_helix as usize] as usize].scaf[b as usize].3 =
self_pos;
}
if c_helix > -1 {
println!(
"changing {} scaf {} {} {}",
vstrand_idx[c_helix as usize],
d,
helices[dom.helix as usize],
self_pos
);
vstrands[vstrand_idx[c_helix as usize] as usize].scaf[d as usize].0 =
helices[dom.helix as usize];
vstrands[vstrand_idx[c_helix as usize] as usize].scaf[d as usize].1 =
self_pos;
}
} else {
vstrands[vstrand_idx[dom.helix as usize] as usize].stap
[self_pos as usize] = result;
if a_helix > -1 {
vstrands[vstrand_idx[a_helix as usize] as usize].stap[b as usize].2 =
helices[dom.helix as usize];
vstrands[vstrand_idx[a_helix as usize] as usize].stap[b as usize].3 =
self_pos;
}
if c_helix > -1 {
vstrands[vstrand_idx[c_helix as usize] as usize].stap[d as usize].0 =
helices[dom.helix as usize];
vstrands[vstrand_idx[c_helix as usize] as usize].stap[d as usize].1 =
self_pos;
}
}
}
}
}
for i in 0..vstrands.len() {
for j in &skips {
let j = *j as usize;
if vstrands[i].scaf[j].0 != -1
|| vstrands[i].scaf[j].2 != -1
|| vstrands[i].stap[j].0 != -1
|| vstrands[i].stap[j].2 != -1
{
vstrands[i].skip[j] = -1;
}
}
}
for st in self.strands.iter() {
if !st.domains.is_empty() {
let ref dom = st.domains[0];
let is_scaf = (helices[dom.helix as usize] % 2 == 0) ^ !dom.right;
if !is_scaf {
let color = st.color.unwrap_or_else(|| st.default_color());
let idx = if do_skip {
convert_skip(dom.start - min_x, &skips)
} else {
dom.start - min_x
};
vstrands[vstrand_idx[dom.helix as usize] as usize]
.stap_colors
.push((idx, color as isize));
}
}
}
let nano = cadnano_format::Cadnano {
name: "name".to_string(),
vstrands,
};
Ok(serde_json::to_writer(w, &nano).unwrap())
}
}
fn make_graph<S, T>(
design: &Design<S, T>,
max_y: usize,
father: &mut Vec<usize>,
) -> Vec<Vec<bool>> {
let mut ret = vec![vec![false; max_y + 1]; max_y + 1];
for s in &design.strands {
let mut group_sens: Vec<usize> = Vec::new();
let mut group_anti: Vec<usize> = Vec::new();
for d in &s.domains {
if d.right {
group_sens.push(d.helix as usize);
} else {
group_anti.push(d.helix as usize);
}
}
for i in group_sens.iter() {
for j in group_anti.iter() {
let repr_i = find(*i, father);
let repr_j = find(*j, father);
if repr_i == repr_j {
println!("not doable in cadnano");
}
ret[repr_i][repr_j] = true;
ret[repr_j][repr_i] = true;
}
}
}
ret
}
fn color_graph(graph: &Vec<Vec<bool>>, max_y: usize, father: &mut Vec<usize>) -> Vec<bool> {
let mut color = vec![false; max_y + 1];
let mut seen: Vec<bool> = (0..(max_y + 1)).map(|i| i != father[i]).collect();
for i in 0..(max_y + 1) {
if !seen[i] {
seen[i] = true;
let mut to_do: Vec<usize> = vec![i];
while to_do.len() > 0 {
let i = to_do.pop().unwrap();
let i = find(i, father);
for j in 0..(max_y + 1) {
let j = find(j, father);
if graph[i][j] && !seen[j] {
if seen[j] && color[j] == color[i] {
println!("not doable in cadnano");
}
seen[j] = true;
color[j] = !color[i];
to_do.push(j);
}
}
}
}
}
for i in 0..(max_y + 1) {
let repr_i = find(i, father);
color[i] = color[repr_i];
}
color
}
fn make_group<S, T>(design: &Design<S, T>, max_y: usize) -> Vec<usize> {
let mut father: Vec<usize> = (0..(max_y + 1)).map(|i| i).collect();
let mut rank: Vec<usize> = vec![0; max_y + 1];
for s in &design.strands {
let mut group_sens: Vec<usize> = Vec::new();
let mut group_anti: Vec<usize> = Vec::new();
for d in &s.domains {
if d.right {
group_sens.push(d.helix as usize);
} else {
group_anti.push(d.helix as usize);
}
}
for i in group_sens.iter() {
for j in group_sens.iter() {
union(*i, *j, &mut father, &mut rank);
}
}
for i in group_anti.iter() {
for j in group_anti.iter() {
union(*i, *j, &mut father, &mut rank);
}
}
}
father
}
fn union(i: usize, j: usize, father: &mut Vec<usize>, rank: &mut Vec<usize>) {
let i_root = find(i, father);
let j_root = find(j, father);
if i_root != j_root {
if rank[i_root] < rank[j_root] {
father[i_root] = j_root;
} else {
father[j_root] = i_root;
if rank[j_root] == rank[i_root] {
rank[i_root] += 1;
}
}
}
}
fn find(i: usize, father: &mut Vec<usize>) -> usize {
if father[i] != i {
father[i] = find(father[i], father);
}
father[i]
}
fn convert_skip(n: isize, skips: &Vec<isize>) -> isize {
let mut i = 0;
let mut ret = n;
while i < skips.len() && ret >= skips[i] {
ret += 1;
i += 1;
}
ret
}
fn convert_result(
t: (isize, isize, isize, isize),
n: isize,
h: isize,
skips: &Vec<isize>,
) -> (isize, isize, isize, isize) {
let (a, b, c, d) = t;
let f = |x, b| {
if x == -1 {
-1
} else if b {
let shift = convert_skip(n, skips) - n;
x + shift
} else {
let shift = convert_skip(x, skips) - x;
x + shift
}
};
(a, f(b, a == h), c, f(d, c == h))
}
fn get_skip_points(
xovers: &Vec<(isize, isize, bool, isize, isize, bool)>,
max_x: isize,
) -> Vec<isize> {
let mut ret = Vec::new();
let mut guess = SKIP_STEP;
while guess < max_x {
let skip = get_one_skip(xovers, guess, ret.len());
if skip == -1 {
return Vec::new();
}
ret.push(skip);
guess += SKIP_STEP + 1;
}
ret
}
fn get_one_skip(
xovers: &Vec<(isize, isize, bool, isize, isize, bool)>,
guess: isize,
nb_skips: usize,
) -> isize {
let delta = vec![0, -1, 1, -2, 2, -3, 3, -4, 4];
for d in delta {
let mut ok = true;
for (_, a, _, _, b, _) in xovers {
if *a + nb_skips as isize == guess + d || *b + nb_skips as isize == guess + d {
ok = false;
break;
}
}
if ok {
return guess + d;
}
}
-1
}
impl Design<(), ()> {
pub fn from_cadnano<P: AsRef<Path>>(file: P) -> Self {
let nano = Cadnano::from_file(file).unwrap();
let vstrands = nano.vstrands;
let mut seen: HashSet<(usize, usize, bool)> = HashSet::new();
let mut design: Design<(), ()> = Design::new();
let mut num_to_helix: HashMap<isize, usize> = HashMap::new();
let honneycomb = vstrands[0].scaf.len() % 21 == 0;
for (i, v) in vstrands.iter().enumerate() {
num_to_helix.insert(v.num, i);
if honneycomb {
design.add_honneycomb_helix(v.row, v.col);
} else {
design.add_adjusted_grid_helix(v.row, v.col);
}
}
num_to_helix.insert(-1, NO_HELIX);
for scaf in vec![false, true] {
for i in 0..vstrands.len() {
let v = &vstrands[i];
for j in 0..v.stap.len() {
let result = if scaf { v.scaf[j] } else { v.stap[j] };
if seen.insert((i, j, scaf)) && result != (-1, -1, -1, -1) {
println!("{}, {}, {}", scaf, i, j);
let end_5 = find_5_end(i, j, &vstrands, &num_to_helix, scaf);
println!("end: {:?}", end_5);
let strand = make_strand(end_5, &vstrands, &num_to_helix, &mut seen, scaf);
design.strands.push(strand);
}
}
}
}
design
}
}
fn find_5_end(
i: usize,
j: usize,
vstrands: &Vec<VStrand>,
num_to_helix: &HashMap<isize, usize>,
scaf: bool,
) -> (usize, usize, bool) {
let (mut a, mut b, mut c, mut d) = (i, j, 0, 0);
let mut cyclic = false;
while a != NO_HELIX {
let result = if scaf {
vstrands[a].scaf[b]
} else {
vstrands[a].stap[b]
};
c = a;
d = b;
a = num_to_helix[&result.0];
b = result.1 as usize;
if a == i && b == j {
cyclic = true;
a = NO_HELIX;
}
}
(c, d, cyclic)
}
fn make_strand(
end_5: (usize, usize, bool),
vstrands: &Vec<VStrand>,
num_to_helix: &HashMap<isize, usize>,
seen: &mut HashSet<(usize, usize, bool)>,
scaf: bool,
) -> Strand<(), ()> {
println!("making strand {:?}", end_5);
let cyclic = end_5.2;
let (mut i, mut j) = (end_5.0, end_5.1);
let mut ret = Strand {
domains: Vec::new(),
label: None,
sequence: None,
cyclic,
color: None,
};
let mut curent_dom = 0;
while curent_dom == 0 || i != end_5.0 || j != end_5.1 {
let curent_helix = i;
let curent_5 = j;
let mut curent_3 = j;
let mut once = false;
while i == curent_helix && (i != end_5.0 || j != end_5.1 || !once) {
once = true;
curent_3 = j;
println!("nucl {}, {}", i, j);
seen.insert((i, j, scaf));
let result = if scaf {
vstrands[i].scaf[j]
} else {
vstrands[i].stap[j]
};
println!("result {:?}", result);
i = num_to_helix[&result.2];
j = result.3 as usize;
}
println!("ready to build domain");
let right = curent_3 >= curent_5;
let start = if right {
substract_skips(curent_5, curent_helix, vstrands)
} else {
substract_skips(curent_3, curent_helix, vstrands)
};
let end = if right {
substract_skips(curent_3, curent_helix, vstrands)
} else {
substract_skips(curent_5, curent_helix, vstrands)
};
println!("pushing {} {} {} {}", curent_helix, start, end, right);
ret.domains.push(Domain {
helix: curent_helix as isize,
start,
end: end + 1,
right,
label: None,
sequence: None,
});
if i == NO_HELIX {
break;
}
curent_dom += 1;
}
if cyclic {
let dom0 = &ret.domains[0];
let last_dom = &ret.domains[ret.domains.len() - 1];
if last_dom.helix != dom0.helix {
let helix = dom0.helix;
let start = dom0.start;
let end = dom0.start + 1;
let right = dom0.right;
ret.domains.push(Domain {
helix,
start,
end,
right,
label: None,
sequence: None,
});
} else {
let len = ret.domains.len();
let right = dom0.right;
let start = dom0.start;
let end = dom0.end;
let last_dom = &mut ret.domains[len - 1];
if right {
last_dom.end = start + 1;
} else {
last_dom.start = end - 1;
}
}
}
ret
}
fn substract_skips(nucl: usize, helix: usize, vstrands: &Vec<VStrand>) -> isize {
let skips: isize = (0..(nucl + 1))
.map(|n| vstrands[helix].skip[n as usize])
.sum();
nucl as isize + skips
}