extern crate palette;
extern crate serde;
extern crate serde_json;
#[macro_use]
extern crate serde_derive;
use palette::{Hsl, Srgb};
use std::collections::HashMap;
use std::collections::HashSet;
mod cadnano;
pub use cadnano::*;
#[derive(Debug)]
pub enum Error {
IO(std::io::Error),
Json(serde_json::Error),
}
impl std::convert::From<std::io::Error> for Error {
fn from(e: std::io::Error) -> Self {
Error::IO(e)
}
}
impl std::convert::From<serde_json::Error> for Error {
fn from(e: serde_json::Error) -> Self {
Error::Json(e)
}
}
pub fn color_rgb(y: usize, x: usize) -> Srgb {
Hsl::new(((y * 31 + x * 3) % 360) as f32, 1., 0.4).into()
}
pub fn color(y: usize, x: usize) -> isize {
let rgb = color_rgb(y, x);
let (r, g, b) = rgb.into_components();
let r = (r * 255.).round() as isize;
let g = (g * 255.).round() as isize;
let b = (b * 255.).round() as isize;
(((r << 8) | g) << 8) | b
}
#[derive(Debug)]
pub struct Point {
pub helix: isize,
pub x: isize,
}
#[derive(Debug)]
pub struct Strand {
pub path: Vec<Point>,
}
#[derive(Debug)]
pub struct Strands {
pub scaffolds: Vec<Strand>,
pub staples: Vec<Strand>,
pub edits: Vec<(Point, isize)>,
}
pub enum Lattice {
Square,
Hexagonal,
}
impl Cadnano {
fn follow_scaffold(
&self,
points: &mut HashSet<(isize, isize)>,
nums: &HashMap<isize, usize>,
helix0: isize,
x0: isize,
) -> Strand {
let mut path = Vec::new();
let (mut helix, mut x) = (helix0, x0);
while helix != -1 {
let (a, b, _, _) = self.vstrands[*nums.get(&helix).unwrap() as usize].scaf[x as usize];
let n = path.len();
if n < 2 {
path.push(Point { helix, x });
} else {
if path[n-1].helix == path[n-2].helix && path[n-1].helix == helix {
path.last_mut().unwrap().x = x
} else {
path.push(Point { helix, x });
}
}
points.remove(&(helix, x));
helix = a;
x = b;
}
path.reverse();
path.pop();
println!("scaffold reversed: {:?}", path);
helix = helix0;
x = x0;
while helix != -1 {
let (_, _, c, d) = self.vstrands[*nums.get(&helix).unwrap() as usize].scaf[x as usize];
let n = path.len();
if n < 2 {
path.push(Point { helix, x });
} else {
if path[n-1].helix == path[n-2].helix && path[n-1].helix == helix {
path.last_mut().unwrap().x = x
} else {
path.push(Point { helix, x });
}
}
points.remove(&(helix, x));
helix = c;
x = d;
}
Strand {
path
}
}
fn follow_staple(
&self,
points: &mut HashSet<(isize, isize)>,
nums: &HashMap<isize, usize>,
helix0: isize,
x0: isize,
) -> Strand {
let mut path = Vec::new();
let (mut helix, mut x) = (helix0, x0);
while helix != -1 {
let (a, b, _, _) = self.vstrands[*nums.get(&helix).unwrap() as usize].stap[x as usize];
let n = path.len();
if n < 2 {
path.push(Point { helix, x });
} else {
if path[n-1].helix == path[n-2].helix && path[n-1].helix == helix {
path.last_mut().unwrap().x = x
} else {
path.push(Point { helix, x });
}
}
points.remove(&(helix, x));
helix = a;
x = b;
}
println!("reversed staple: {:?}", path);
path.reverse();
path.pop();
helix = helix0;
x = x0;
while helix != -1 {
let (_, _, c, d) = self.vstrands[*nums.get(&helix).unwrap() as usize].stap[x as usize];
let n = path.len();
if n < 2 {
path.push(Point { helix, x });
} else {
if path[n-1].helix == path[n-2].helix && path[n-1].helix == helix {
path.last_mut().unwrap().x = x
} else {
path.push(Point { helix, x });
}
}
points.remove(&(helix, x));
helix = c;
x = d;
}
Strand {
path
}
}
pub fn to_strands(&self) -> Strands {
let mut scaf_points = HashSet::new();
let mut stap_points = HashSet::new();
let mut nums = HashMap::new();
for (vs, vstrand) in self.vstrands.iter().enumerate() {
nums.insert(vstrand.num, vs);
for (x, st) in vstrand.scaf.iter().enumerate() {
if *st != (-1, -1, -1, -1) {
scaf_points.insert((vstrand.num, x as isize));
}
}
for (x, st) in vstrand.stap.iter().enumerate() {
if *st != (-1, -1, -1, -1) {
stap_points.insert((vstrand.num, x as isize));
}
}
}
let mut scaffolds = Vec::new();
loop {
let (num, x) = {
let mut it = scaf_points.iter();
if let Some(&(num, x)) = it.next() {
(num, x)
} else {
break
}
};
scaffolds.push(self.follow_scaffold(&mut scaf_points, &nums, num, x))
}
let mut staples = Vec::new();
loop {
let (num, x) = {
let mut it = stap_points.iter();
if let Some(&(num, x)) = it.next() {
(num, x)
} else {
break
}
};
staples.push(self.follow_staple(&mut stap_points, &nums, num, x))
}
let mut edits = Vec::new();
for (h, vs) in self.vstrands.iter().enumerate() {
for &x in vs.skip.iter().chain(vs.loop_.iter()) {
if x != 0 {
edits.push((
Point { helix: h as isize, x: x as isize },
x
))
}
}
}
Strands {
scaffolds,
staples,
edits,
}
}
pub fn from_strands(
name: String,
lattice: Lattice,
scaffolds: &[Strand],
staples: &[Strand],
edits: &[(Point, isize)],
) -> Cadnano {
let ref p0 = scaffolds[0].path[0];
let (min_x, max_x, min_y, max_y) = scaffolds
.iter().flat_map(|s| s.path.iter())
.chain(staples.iter().flat_map(|s| s.path.iter()))
.fold(
(p0.x, p0.x, p0.helix, p0.helix),
|(min_x, max_x, min_y, max_y), p| {
(
std::cmp::min(min_x, p.x),
std::cmp::max(max_x, p.x),
std::cmp::min(min_y, p.helix),
std::cmp::max(max_y, p.helix),
)
},
);
let min_x = min_x;
let columns = (max_x - min_x + 1) as usize;
let mut vstrands = Vec::new();
let columns = match lattice {
Lattice::Square => {
let mut cols = (columns / 32) * 32;
while cols % 21 == 0 || cols < columns {
cols += 32
}
cols
}
Lattice::Hexagonal => {
let mut cols = (columns / 21) * 21;
while cols % 32 == 0 || cols < columns {
cols += 21
}
cols
}
};
for y in min_y..=max_y {
let num = y - min_y;
let stap: Vec<(isize, isize, isize, isize)> = vec![(-1, -1, -1, -1); columns];
let scaf = stap.clone();
vstrands.push(VStrand {
stap_colors: vec![],
num,
scaf_loop: vec![],
stap,
skip: vec![0; columns as usize],
scaf,
stap_loop: vec![],
col: 1,
loop_: vec![0; columns as usize],
row: num as isize + 1,
})
}
for (is_scaffold, path) in
scaffolds.iter().map(|x| (true, x))
.chain(staples.iter().map(|x| (false, x)))
{
let (mut x0, mut y0): (isize, isize) = (-1, -1);
let it0 = path.path.iter().map(|p| Point {
helix: p.helix - min_y,
x: p.x - min_x,
});
let mut it1 = path.path.iter().map(|p| Point {
helix: p.helix - min_y,
x: p.x - min_x,
});
it1.next();
for (p0, p1) in it0.zip(it1) {
if p0.helix == p1.helix {
if p0.x < p1.x {
for x in p0.x..p1.x {
if is_scaffold {
vstrands[p0.helix as usize].scaf[x as usize] =
(y0, x0, p0.helix, x + 1);
} else {
vstrands[p0.helix as usize].stap[x as usize] =
(y0, x0, p0.helix, x + 1);
}
y0 = p0.helix;
x0 = x;
}
if is_scaffold {
vstrands[p1.helix as usize].scaf[p1.x as usize] = (y0, x0, -1, -1);
} else {
vstrands[p1.helix as usize].stap[p1.x as usize] = (y0, x0, -1, -1);
}
} else {
for x in ((p1.x + 1)..=p0.x).rev() {
if is_scaffold {
vstrands[p0.helix as usize].scaf[x as usize] =
(y0, x0, p1.helix, x - 1);
} else {
vstrands[p0.helix as usize].stap[x as usize] =
(y0, x0, p1.helix, x - 1);
}
y0 = p0.helix;
x0 = x;
}
if is_scaffold {
vstrands[p1.helix as usize].scaf[p1.x as usize] = (y0, x0, -1, -1);
} else {
vstrands[p1.helix as usize].stap[p1.x as usize] = (y0, x0, -1, -1);
}
}
} else {
if is_scaffold {
vstrands[p0.helix as usize].scaf[p0.x as usize] = (y0, x0, p1.helix, p1.x);
} else {
vstrands[p0.helix as usize].stap[p0.x as usize] = (y0, x0, p1.helix, p1.x);
}
x0 = p0.x;
y0 = p0.helix;
}
}
}
for &(ref p, ref n) in edits.iter() {
println!("{:?}", p);
if *n > 0 {
vstrands[p.helix as usize].loop_[p.x as usize] = *n
} else if *n < 0 {
vstrands[p.helix as usize].skip[p.x as usize] = *n
}
}
Cadnano { name, vstrands }
}
}
#[test]
fn test_cad() {
let scaf = vec![Strand {
path: vec![
Point { x: 0, helix: 0 },
Point { x: 20, helix: 0 },
Point { x: 20, helix: 1 },
Point { x: 0, helix: 1 },
],
}];
let staples = vec![Strand {
path: vec![
Point { x: 20, helix: 0 },
Point { x: 0, helix: 0 },
Point { x: 0, helix: 1 },
Point { x: 20, helix: 1 },
],
}];
let cad = Cadnano::from_strands("blabla".to_string(), Lattice::Square, &scaf, &staples, &[(Point { x: 4, helix: 0 }, 1), (Point { x: 4, helix: 1 }, -1)]);
cad.to_file("/tmp/test.json").unwrap();
println!("{:?}", cad.to_strands());
}