#![deny(missing_docs)]
#[macro_use]
extern crate serde_derive;
extern crate serde;
pub use anyhow::Error;
use std::borrow::Cow;
use std::collections::{HashMap, HashSet, VecDeque};
use std::f64::consts::PI;
use std::fmt;
pub mod geometry;
use crate::geometry::*;
pub mod sequences;
pub mod cadnano;
pub mod mrdna;
pub mod autoroll;
pub use autoroll::{RollSystem, SpringSystem};
pub mod patch;
use patch::Patch;
pub mod guess3d;
#[derive(Serialize, Deserialize, Clone)]
pub struct Design<StrandLabel, DomainLabel> {
pub version: String,
pub helices: Vec<Helix>,
pub strands: Vec<Strand<StrandLabel, DomainLabel>>,
#[serde(skip_serializing_if = "Option::is_none", default)]
pub parameters: Option<Parameters>,
}
type Tnucl = (isize, isize, bool);
fn is_false(x: &bool) -> bool {
!*x
}
fn none<Label>() -> Option<Label> {
None
}
#[derive(Debug, Default, Clone, Serialize, Deserialize)]
pub struct Strand<Label, DomainLabel> {
pub domains: Vec<Domain<DomainLabel>>,
#[serde(skip_serializing_if = "Option::is_none", default)]
pub sequence: Option<Cow<'static, str>>,
#[serde(skip_serializing_if = "is_false", default)]
pub cyclic: bool,
#[serde(skip_serializing_if = "Option::is_none", default)]
pub color: Option<Color>,
#[serde(skip_serializing_if = "Option::is_none", default = "none")]
pub label: Option<Label>,
}
#[derive(Debug, Clone, Serialize, Deserialize)]
#[serde(untagged)]
pub enum Color {
Int(u32),
Hex(String),
Rgb {
r: u8,
g: u8,
b: u8
},
}
impl Color {
pub fn as_int(&self) -> u32 {
match *self {
Color::Int(n) => n,
Color::Hex(ref s) => {
let s = s.trim_start_matches("0x");
let s = s.trim_start_matches("#");
u32::from_str_radix(s, 16).unwrap()
}
Color::Rgb { r, g, b } => {
((r as u32) << 16) | ((g as u32) << 8) | (b as u32)
}
}
}
}
impl<Label, DomainLabel> Strand<Label, DomainLabel> {
pub fn find_nucl(&self, nucl: isize) -> (isize, isize, bool) {
let mut current = 0;
let mut i = 0;
while current + &self.domains[i].length() <= nucl {
current += &self.domains[i].length();
i += 1;
}
let position = if self.domains[i].forward {
self.domains[i].start + nucl - current
} else {
self.domains[i].end - 1 - nucl + current
};
(self.domains[i].helix, position, self.domains[i].forward)
}
}
#[derive(Debug, Default, Clone, Serialize, Deserialize)]
pub struct Domain<Label> {
pub helix: isize,
pub start: isize,
pub end: isize,
pub forward: bool,
#[serde(skip_serializing_if = "Option::is_none", default = "none")]
pub label: Option<Label>,
pub sequence: Option<Cow<'static, str>>,
}
impl<Label> Domain<Label> {
pub fn iter(&self) -> DomainIter {
DomainIter {
start: self.start,
end: self.end,
forward: self.forward,
}
}
pub fn translate(self, dx: isize, dy: isize) -> Self {
use std::convert::TryFrom;
Domain {
start: self.start + dx,
end: self.end + dx,
helix: usize::try_from(self.helix as isize + dy).unwrap() as isize,
..self
}
}
pub fn shift_x(self, dx: isize) -> Self {
Domain {
start: self.start + dx,
end: self.end + dx,
..self
}
}
pub fn shift_y(self, dy: isize) -> Self {
use std::convert::TryFrom;
Domain {
helix: usize::try_from(self.helix as isize + dy).unwrap() as isize,
..self
}
}
pub fn length(&self) -> isize {
self.end - self.start
}
pub fn pseudo_copy(&self) -> Self {
Domain {
helix: self.helix,
start: self.start,
end: self.end,
forward: self.forward,
label: None,
sequence: None,
}
}
pub fn contains(&self, h: isize, x: isize, b: bool) -> bool {
self.helix == h && self.forward == b && self.start <= x && self.end > x
}
pub fn first_nucl(&self) -> Option<(isize, isize, bool)> {
if self.start >= self.end {
None
} else {
if self.forward {
Some((self.helix, self.start, self.forward))
} else {
Some((self.helix, self.end - 1, self.forward))
}
}
}
pub fn last_nucl(&self) -> Option<(isize, isize, bool)> {
if self.start >= self.end {
None
} else {
if self.forward {
Some((self.helix, self.end - 1, self.forward))
} else {
Some((self.helix, self.start, self.forward))
}
}
}
}
pub struct DomainIter {
start: isize,
end: isize,
forward: bool,
}
impl Iterator for DomainIter {
type Item = isize;
fn next(&mut self) -> Option<Self::Item> {
if self.start >= self.end {
None
} else {
if self.forward {
let s = self.start;
self.start += 1;
Some(s)
} else {
let s = self.end;
self.end -= 1;
Some(s - 1)
}
}
}
}
pub const M13_7249: &'static str = include_str!("m13");
pub const M13_594: &'static str = include_str!("m13_594");
#[derive(Copy, Clone, Debug, Serialize, Deserialize)]
pub struct Parameters {
pub z_step: f64,
pub helix_radius: f64,
pub bases_per_turn: f64,
pub groove_angle: f64,
pub inter_helix_gap: f64,
}
impl Parameters {
pub const DEFAULT: Parameters = Parameters {
z_step: 0.332,
helix_radius: 1.,
bases_per_turn: 10.44,
groove_angle: -24. * PI / 34.,
inter_helix_gap: 0.65,
};
}
#[derive(Serialize, Deserialize, Clone)]
pub struct Helix {
#[serde(default="zero_point")]
pub position: Point<f64>,
#[serde(default="zero_f64")]
pub roll: f64,
#[serde(default="zero_f64")]
pub yaw: f64,
#[serde(default="zero_f64")]
pub pitch: f64,
pub max_offset: Option<isize>,
pub major_ticks: Option<Vec<isize>>,
}
fn zero_point() -> Point<f64> {
Point { x: 0., y: 0. ,z: 0. }
}
fn zero_f64() -> f64 {
0.
}
impl fmt::Debug for Helix {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
f.debug_tuple("").field(&self.position).finish()
}
}
impl Helix {
pub fn theta(&self, n: isize, forward: bool, cst: &Parameters) -> f64 {
let shift = if forward { cst.groove_angle } else { 0. };
n as f64 * 2. * PI / cst.bases_per_turn + shift + self.roll
}
pub fn space_pos(&self, p: &Parameters, n: isize, forward: bool) -> [f64; 3] {
let theta = self.theta(n, forward, p);
let mut ret = [
n as f64 * p.z_step,
theta.cos() * p.helix_radius,
theta.sin() * p.helix_radius,
];
ret = self.rotate_point(ret);
ret[0] += self.position.x;
ret[1] += self.position.y;
ret[2] += self.position.z;
ret
}
pub(crate) fn rotate_point(&self, ret: [f64; 3]) -> [f64; 3] {
let forward = [
self.yaw.cos() * self.pitch.cos(),
self.pitch.sin(),
-self.yaw.sin() * self.pitch.cos(),
];
let right = [self.yaw.sin(), 0., self.yaw.cos()];
let up = [
right[1] * forward[2] - right[2] * forward[1],
right[2] * forward[0] - right[0] * forward[2],
right[0] * forward[1] - right[1] * forward[0],
];
[
ret[0] * forward[0] + ret[1] * up[0] + ret[2] * right[0],
ret[0] * forward[1] + ret[1] * up[1] + ret[2] * right[1],
ret[0] * forward[2] + ret[1] * up[2] + ret[2] * right[2],
]
}
pub fn basis(&self) -> [[f64; 3]; 3] {
[
self.rotate_point([1., 0., 0.]),
self.rotate_point([0., self.roll.cos(), self.roll.sin()]),
self.rotate_point([0., -self.roll.sin(), self.roll.cos()]),
]
}
pub fn axis_pos(&self, p: &Parameters, n: isize) -> [f64; 3] {
let mut ret = [n as f64 * p.z_step, 0., 0.];
ret = self.rotate_point(ret);
ret[0] += self.position.x;
ret[1] += self.position.y;
ret[2] += self.position.z;
ret
}
pub fn overlap(&self, other: &Helix, p: &Parameters) -> bool {
let dir_vec = Point::from_coord(self.axis_pos(p, 1)) - self.position;
let vec1 = Point::from_coord(other.axis_pos(p, 30)) - self.position;
if vec1.cross_product(&dir_vec).norm() / dir_vec.norm() > p.helix_radius {
false
} else {
let vec2 = Point::from_coord(other.axis_pos(p, -30)) - self.position;
vec2.cross_product(&dir_vec).norm() / dir_vec.norm() < p.helix_radius
}
}
pub fn clone_up(&self, p: &Parameters) -> Self {
let mut new_position = [0., p.helix_radius * 2. + p.inter_helix_gap, 0.];
new_position = self.rotate_point(new_position);
new_position[0] += self.position.x;
new_position[1] += self.position.y;
new_position[2] += self.position.z;
Helix {
position: Point::from_coord(new_position),
.. self.clone()
}
}
pub fn clone_down(&self, p: &Parameters) -> Self {
let mut new_position = [0., -p.helix_radius * 2. - p.inter_helix_gap, 0.];
new_position = self.rotate_point(new_position);
new_position[0] += self.position.x;
new_position[1] += self.position.y;
new_position[2] += self.position.z;
Helix {
position: Point::from_coord(new_position),
.. self.clone()
}
}
pub fn clone_left(&self, p: &Parameters) -> Self {
let mut new_position = [0., 0., -p.helix_radius * 2. - p.inter_helix_gap];
new_position = self.rotate_point(new_position);
new_position[0] += self.position.x;
new_position[1] += self.position.y;
new_position[2] += self.position.z;
Helix {
position: Point::from_coord(new_position),
.. self.clone()
}
}
pub fn clone_forward(&self, p: &Parameters) -> Self {
let mut new_position = [0., 0., p.helix_radius * 2. + p.inter_helix_gap];
new_position = self.rotate_point(new_position);
new_position[0] += self.position.x;
new_position[1] += self.position.y;
new_position[2] += self.position.z;
Helix {
position: Point::from_coord(new_position),
.. self.clone()
}
}
pub fn closest_nucl<F: Into<f64> + Copy>(&self, point: [F; 3], p: &Parameters) -> isize {
let point = Point::from_coord([point[0].into(), point[1].into(), point[2].into()]);
let mut up = 10000;
let mut low = -10000;
while up - low > 1 {
let point_low = Point::from_coord(self.axis_pos(p, low));
let point_up = Point::from_coord(self.axis_pos(p, up));
let dist_low = (point_low - point).norm();
let dist_up = (point_up - point).norm();
if dist_low > dist_up {
low = (up + low) / 2;
} else {
up = (up + low) / 2;
}
}
let point_low = Point::from_coord(self.axis_pos(p, low));
let point_up = Point::from_coord(self.axis_pos(p, up));
let dist_low = (point_low - point).norm();
let dist_up = (point_up - point).norm();
if dist_low > dist_up {
up
} else {
low
}
}
}
#[derive(Debug, Copy, Clone, PartialEq, Eq, PartialOrd, Ord, Hash)]
pub struct HelixId(pub isize);
impl From<isize> for HelixId {
fn from(e: isize) -> Self {
HelixId(e)
}
}
impl From<i32> for HelixId {
fn from(e: i32) -> Self {
HelixId(e as isize)
}
}
impl From<usize> for HelixId {
fn from(e: usize) -> Self {
HelixId(e as isize)
}
}
pub struct StrandId(pub usize);
#[doc(hidden)]
#[derive(Serialize, Deserialize)]
pub struct Res {
pub out: String,
pub err: String,
}
pub struct StrandRef<'a, StrandLabel, DomainLabel> {
builder: &'a mut Design<StrandLabel, DomainLabel>,
strand_id: usize,
x: Option<isize>,
x0: Option<isize>,
}
impl<StrandLabel: serde::Serialize, DomainLabel: serde::Serialize>
Design<StrandLabel, DomainLabel>
{
pub fn new() -> Self {
Design {
version: env!("CARGO_PKG_VERSION").to_string(),
helices: Vec::new(),
strands: Vec::new(),
parameters: Some(Parameters::DEFAULT),
}
}
pub fn square_grid_helix(&self, h: isize, v: isize) -> Helix {
Helix {
position: Point {
x: 0.,
y: h as f64
* (self.parameters.unwrap().helix_radius * 2.
+ self.parameters.unwrap().inter_helix_gap),
z: v as f64
* (self.parameters.unwrap().helix_radius * 2.
+ self.parameters.unwrap().inter_helix_gap),
},
roll: 0.,
pitch: 0.,
yaw: 0.,
major_ticks: None,
max_offset: None,
}
}
pub fn honneycomb_grid_helix(&self, h: isize, v: isize) -> Helix {
let param = &self.parameters.expect("parameters is none");
let r = param.helix_radius + param.inter_helix_gap / 2.;
let upper = -3. * h as f64 * r;
let lower = upper - r;
Helix {
position: Point {
x: 0.,
y: if h % 2 == v % 2 { upper } else { lower },
z: v as f64 * r * 3f64.sqrt(),
},
roll: 0.,
pitch: 0.,
yaw: 0.,
major_ticks: None,
max_offset: None,
}
}
pub fn adjusted_grid_helix(&self, h: isize, v: isize) -> Helix {
let roll = if (h + v) % 2 == 0 {
self.parameters.unwrap().groove_angle - PI
} else {
0.
};
Helix {
position: Point {
x: 0.,
y: -h as f64
* (self.parameters.unwrap().helix_radius * 2.
+ self.parameters.unwrap().inter_helix_gap),
z: v as f64
* (self.parameters.unwrap().helix_radius * 2.
+ self.parameters.unwrap().inter_helix_gap),
},
roll,
pitch: 0.,
yaw: 0.,
major_ticks: None,
max_offset: None,
}
}
pub fn with_parameters(parameters: Parameters) -> Self {
let mut b = Design::new();
b.parameters = Some(parameters);
b
}
pub fn add_grid_helix(&mut self, h: isize, v: isize) -> HelixId {
self.helices.push(self.square_grid_helix(h, v));
HelixId((self.helices.len() - 1) as isize)
}
pub fn add_honneycomb_helix(&mut self, h: isize, v: isize) -> HelixId {
self.helices.push(self.honneycomb_grid_helix(h, v));
HelixId((self.helices.len() - 1) as isize)
}
pub fn add_adjusted_grid_helix(&mut self, h: isize, v: isize) -> HelixId {
self.helices.push(self.adjusted_grid_helix(h, v));
HelixId((self.helices.len() - 1) as isize)
}
pub fn add_helix(
&mut self,
x: f64,
y: f64,
z: f64,
roll: f64,
pitch: f64,
yaw: f64,
) -> HelixId {
self.helices.push(Helix {
position: Point { x, y, z },
roll,
pitch,
yaw,
major_ticks: None,
max_offset: None,
});
HelixId((self.helices.len() - 1) as isize)
}
pub fn strand<H: Into<HelixId>>(
&mut self,
helix: H,
start: isize,
) -> StrandRef<StrandLabel, DomainLabel> {
let strand_id = self.strands.len();
self.strands.push(Strand {
domains: vec![Domain {
helix: helix.into().0,
start,
end: start,
forward: true,
label: None,
sequence: None,
}],
label: None,
sequence: None,
cyclic: false,
color: None,
});
StrandRef {
strand_id,
builder: self,
x: Some(start),
x0: Some(start),
}
}
pub fn empty_strand(&mut self) -> StrandRef<StrandLabel, DomainLabel> {
let strand_id = self.strands.len();
self.strands.push(Strand {
domains: Vec::new(),
cyclic: false,
color: None,
sequence: None,
label: None,
});
StrandRef {
strand_id,
builder: self,
x: None,
x0: None,
}
}
pub fn get_strand_ref(&mut self, strand_id: usize) -> StrandRef<StrandLabel, DomainLabel> {
let strand = &self.strands[strand_id];
if strand.cyclic {
panic!("Attempt to keep building on a cyclic strand");
}
let domain = &strand.domains[strand.domains.len() - 1];
let x0 = if domain.forward {
domain.start
} else {
domain.end - 1
};
let x = if domain.forward {
domain.end - 1
} else {
domain.start
};
StrandRef {
strand_id,
builder: self,
x: Some(x),
x0: Some(x0),
}
}
pub fn fix_colors(&mut self) {
for s in self.strands.iter_mut() {
if s.color.is_none() {
s.color = Some(s.default_color());
}
}
}
pub fn write(&self) -> Result<(), Error> {
serde_json::to_writer(std::io::stdout(), self)?;
Ok(())
}
pub fn write_to<P: AsRef<std::path::Path>>(&self, p: P) -> Result<(), Error> {
if std::env::var("REMOTE").is_ok() {
serde_json::to_writer_pretty(std::io::stdout(), &self)?
} else {
serde_json::to_writer_pretty(std::fs::File::create(p)?, &self)?
}
Ok(())
}
pub fn json_string(&self) -> String {
serde_json::to_string_pretty(&self).unwrap()
}
pub fn unique_positions(&self) -> Option<(isize, isize, bool, usize, usize)> {
let mut positions = HashMap::new();
for (i, s) in self.strands.iter().enumerate() {
for dom in s.domains.iter() {
for pos in dom.iter() {
if let Some(prev) = positions.insert((dom.helix, pos, dom.forward), i) {
if s.cyclic {
let j = positions.get(&(
s.domains[0].helix,
s.domains[0].iter().next().unwrap(),
s.domains[0].forward,
));
if let Some(x) = j {
if *x == prev {
continue;
}
}
}
return Some((dom.helix, pos, dom.forward, i, prev));
}
}
}
}
None
}
pub fn check_and_write_to<P: AsRef<std::path::Path>>(&self, p: P) -> Result<(), Error> {
if let Some((helix, pos, forward, i, prev)) = self.unique_positions() {
let sens = if forward { "Sens" } else { "AntiSens" };
println!("WARNING: SOME STRANDS ARE CONFLICTING");
println!(
"First conflict found: On helix {}, at position {}, on the {} strand",
helix, pos, sens
);
println!("Strands number {} and {} are in conflict", i, prev);
}
self.write_to(p)
}
pub fn get_domain_builder(&mut self, h: isize, x: isize, b: bool) -> &mut Domain<DomainLabel> {
let mut must_add = true;
let mut i_found = 0;
let mut j_found = 0;
for s in 0..self.strands.len() {
for j in 0..self.strands[s].domains.len() {
let d = &self.strands[s].domains[j];
if d.forward == b && d.helix == h && x >= d.start && x < d.end {
must_add = false;
i_found = s;
j_found = j;
break;
}
}
}
if must_add {
self.strands.push(Strand {
domains: vec![Domain {
helix: h,
start: x,
end: x + 1,
forward: b,
label: None,
sequence: None,
}],
label: None,
sequence: None,
cyclic: false,
color: None,
});
self.strands.last_mut().unwrap().domains.last_mut().unwrap()
} else {
&mut self.strands[i_found].domains[j_found]
}
}
pub fn get_domain_bound(&self, h: isize, x: isize, b: bool) -> (isize, isize, usize, usize) {
for s in 0..self.strands.len() {
for j in 0..self.strands[s].domains.len() {
let d = &self.strands[s].domains[j];
if d.forward == b && d.helix == h && x >= d.start && x < d.end {
return (d.start, d.end, s, j);
}
}
}
(x, x + 1, self.strands.len(), 0)
}
pub fn get_strand_nucl(&self, h: isize, x: isize, b: bool) -> Option<usize> {
for s in 0..self.strands.len() {
for d in &self.strands[s].domains {
if d.contains(h, x, b) {
return Some(s);
}
}
}
None
}
pub fn cross_nucl(&mut self, h1: isize, x1: isize, b1: bool, h2: isize, x2: isize, b2: bool) {
let patch = self.patch_cross_nucl(h1, x1, b1, h2, x2, b2);
if let Some(patch) = patch {
patch.apply(self)
}
}
pub fn patch_cross_nucl(
&self,
h1: isize,
x1: isize,
b1: bool,
h2: isize,
x2: isize,
b2: bool,
) -> Option<Patch> {
let mut position = None;
let mut dest = None;
for i in 0..self.strands.len() {
if self.strands[i].len() > 0 {
let (h, x, b) = self.strands[i].first_nucl().unwrap();
if h == h1 && b == b1 && x == x1 && dest.is_none() {
dest = Some((h, x, b));
continue;
}
if h == h2 && b == b2 && x == x2 && dest.is_none() {
dest = Some((h, x, b));
continue;
}
let (h, x, b) = self.strands[i].last_nucl().unwrap();
if h == h1 && b == b1 && x == x1 && position.is_none() {
position = Some((h, x, b));
continue;
}
if h == h2 && b == b2 && x == x2 && position.is_none() {
position = Some((h, x, b));
continue;
}
}
}
if !dest.is_none() && !position.is_none() {
let position = position.unwrap();
let dest = dest.unwrap();
Some(Patch::MakeBound(
position.0 as usize,
position.2,
position.1,
dest.0 as usize,
dest.2,
dest.1,
))
} else {
None
}
}
pub fn overlap(&self, h: &Helix) -> bool {
for h1 in &self.helices {
if h1.overlap(h, &self.parameters.unwrap()) {
return true;
}
}
false
}
pub fn add_link_helix(
&mut self,
helix_1: usize,
pos_1: isize,
forward_1: bool,
helix_2: usize,
pos_2: isize,
forward_2: bool,
) -> usize {
let param = &self
.parameters
.expect("Failed to unwrap parameters in add_link_helix");
let axis_1 = Point::from_coord(self.helices[helix_1].axis_pos(param, 1))
- Point::from_coord(self.helices[helix_1].axis_pos(param, 0));
let axis_2 = Point::from_coord(self.helices[helix_2].axis_pos(param, 1))
- Point::from_coord(self.helices[helix_2].axis_pos(param, 0));
let axis_1 = axis_1 / axis_1.norm();
let axis_2 = axis_2 / axis_2.norm();
let shift_1 = if forward_1 {
param.helix_radius * axis_1
} else {
-param.helix_radius * axis_1
};
let shift_2 = if forward_2 {
param.helix_radius * axis_2
} else {
-param.helix_radius * axis_2
};
let point_1 = Point::from_coord(self.helices[helix_1].axis_pos(param, pos_1)) + shift_1;
let point_2 = Point::from_coord(self.helices[helix_2].axis_pos(param, pos_2)) + shift_2;
let axis = point_2 - point_1;
let pitch = (axis.y / axis.norm()).asin();
let yaw = (-axis.z).atan2(axis.x);
self.add_helix(point_1.x, point_1.y, point_1.z, 0., pitch, yaw);
self.helices.len() - 1
}
pub fn add_strand_link(
&mut self,
support_helix: usize,
from_helix: usize,
from_nucl: isize,
to_helix: usize,
to_nucl: isize,
) {
let param = &self
.parameters
.expect("parameters is None in add_strand_link");
let point_from = self.helices[from_helix].axis_pos(param, from_nucl);
let point_to = self.helices[to_helix].axis_pos(param, to_nucl);
let start = self.helices[support_helix].closest_nucl(point_from, param);
let end = self.helices[support_helix].closest_nucl(point_to, param);
self.strand(support_helix, start).to(end);
self.strand(support_helix, end).to(start);
}
}
impl<S, T> Design<S, T> {
pub fn auto_roll(&mut self) {
let mut system = RollSystem::from_design(self);
system.solve(self, 0.05);
}
pub fn clean_up(&mut self) {
let mut strand_remove: Vec<usize> = Vec::new();
for s_id in 0..self.strands.len() {
let s = &mut self.strands[s_id];
let mut domain_remove: Vec<usize> = Vec::new();
if s.len() == 0 {
strand_remove.push(s_id);
continue;
}
for d_id in 0..s.domains.len() {
if s.domains[d_id].length() <= 0 {
domain_remove.push(d_id);
}
}
for d_id in domain_remove.iter().rev() {
s.domains.remove(*d_id);
}
}
for s_id in strand_remove.iter().rev() {
self.strands.remove(*s_id);
}
}
pub(crate) fn get_cross_overs(&self) -> Vec<(isize, isize, bool, isize, isize, bool)> {
let mut ret = Vec::new();
for s in &self.strands {
for i in 0..(s.domains.len() - 1) {
if s.domains[i].helix != s.domains[i + 1].helix {
let position = if s.domains[i].forward {
s.domains[i].end - 1
} else {
s.domains[i].start
};
let destination = if s.domains[i + 1].forward {
s.domains[i + 1].start
} else {
s.domains[i + 1].end - 1
};
ret.push((
s.domains[i].helix,
position,
s.domains[i].forward,
s.domains[i + 1].helix,
destination,
s.domains[i + 1].forward,
));
}
}
}
ret
}
pub fn optimize_cross_over(&mut self) {
let mut xovers_map: HashMap<(isize, isize, bool), (isize, isize, bool)> = HashMap::new();
let mut nucl_map: HashMap<(isize, isize, bool), (usize, usize)> = HashMap::new();
let mut optimized = HashSet::new();
let mut nucl_set = HashSet::new();
let mut cyclic_strands = Vec::new();
for s_id in 0..self.strands.len() {
let s = &mut self.strands[s_id];
for i in 0..(s.domains.len() - 1) {
let dom = &s.domains[i];
let dom_ = &s.domains[i + 1];
if dom.helix != dom_.helix {
let position = if dom.forward {
(dom.helix, dom.end - 1, true)
} else {
(dom.helix, dom.start, false)
};
let destination = if dom_.forward {
(dom_.helix, dom_.start, true)
} else {
(dom_.helix, dom_.end - 1, false)
};
xovers_map.insert(position, destination);
xovers_map.insert(destination, position);
nucl_map.insert(position, (s_id, i));
nucl_map.insert(destination, (s_id, i + 1));
}
}
if s.cyclic {
cyclic_strands.push(s_id);
s.cyclic = false;
}
for i in 0..(s.domains.len() - 1) {
let dom = &s.domains[i];
let dom_ = &s.domains[i + 1];
if dom.helix != dom_.helix {
let position = if dom.forward {
(dom.helix, dom.end - 1, true)
} else {
(dom.helix, dom.start, false)
};
let destination = if dom_.forward {
(dom_.helix, dom_.start, true)
} else {
(dom_.helix, dom_.end - 1, false)
};
xovers_map.insert(position, destination);
xovers_map.insert(destination, position);
nucl_map.insert(position, (s_id, i));
nucl_map.insert(destination, (s_id, i + 1));
}
for j in dom.start..dom.end {
nucl_set.insert((dom.helix, j, dom.forward));
}
for j in dom_.start..dom_.end {
nucl_set.insert((dom_.helix, j, dom_.forward));
}
}
if !s.cyclic && s.len() > 0 {
let last = s.last_nucl().unwrap();
println!("last {:?}", last);
nucl_map.insert(last, (s_id, s.domains.len() - 1));
let first = s.first_nucl().unwrap();
nucl_map.insert(first, (s_id, 0));
println!("first {:?}", first);
}
}
for nucl in nucl_map.keys() {
println!("nucl = {:?}", nucl);
if !xovers_map.contains_key(nucl) {
println!("continue");
continue;
}
if optimized.insert(*nucl) {
let mut chain = VecDeque::new();
chain.push_back(*nucl);
let (h, x, b) = *nucl;
println!("nucl = {:?}", nucl);
let mut front;
if nucl_map.contains_key(&(h, x - 1, b)) {
front = (h, x - 1, b);
} else if nucl_map.contains_key(&(h, x + 1, b)) {
front = (h, x + 1, b);
} else {
front = (h, x, b);
}
println!("front = {:?}", front);
while optimized.insert(front) {
chain.push_front(front);
println!("pushed");
let (h, x, b) = front;
let left = (h, x - 1, b);
let forward = (h, x + 1, b);
if nucl_map.contains_key(&left) && !optimized.contains(&left) {
front = left;
println!("front = left {:?}", front);
} else if nucl_map.contains_key(&forward) && !optimized.contains(&forward) {
front = forward;
println!("front = forward {:?}", front);
} else if xovers_map.contains_key(&front) {
front = *xovers_map.get(&front).unwrap();
println!("front = xovers {:?}", front);
}
}
let mut back = *xovers_map.get(&nucl).unwrap();
println!("back = {:?}", back);
while optimized.insert(back) {
println!("pushed");
chain.push_back(back);
let (h, x, b) = back;
let left = (h, x - 1, b);
let forward = (h, x + 1, b);
if nucl_map.contains_key(&left) && !optimized.contains(&left) {
back = left;
} else if nucl_map.contains_key(&forward) && !optimized.contains(&forward) {
back = forward;
} else if xovers_map.contains_key(&back) {
back = *xovers_map.get(&back).unwrap();
}
println!("back = {:?}", back);
}
if chain.len() > 1 {
self.optimize_chain(&chain, &xovers_map, &nucl_map, &nucl_set);
}
}
}
for s_id in cyclic_strands {
let ret = &mut self.strands[s_id];
ret.domains.pop();
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 forward = dom0.forward;
ret.domains.push(Domain {
helix,
start,
end,
forward,
label: None,
sequence: None,
});
} else {
let len = ret.domains.len();
let forward = dom0.forward;
let start = dom0.start;
let end = dom0.end;
let last_dom = &mut ret.domains[len - 1];
if forward {
last_dom.end = start + 1;
} else {
last_dom.start = end - 1;
}
}
ret.cyclic = true;
}
}
fn optimize_chain(
&mut self,
chain: &VecDeque<Tnucl>,
xovers_map: &HashMap<Tnucl, Tnucl>,
nucl_map: &HashMap<Tnucl, (usize, usize)>,
nucl_set: &HashSet<Tnucl>,
) {
let mut helices_map: HashMap<isize, Vec<(isize, bool)>> = HashMap::new();
let mut helices = Vec::new();
for (h, x, b) in chain {
if let Some(v) = helices_map.get_mut(h) {
v.push((*x, *b));
} else {
helices_map.insert(*h, vec![(*x, *b)]);
helices.push(*h);
}
}
let delta: Vec<isize> = (-5..5).collect();
let delta_len = delta.len();
let default_delta = 5;
let mut opt = vec![vec![(0., default_delta); delta_len]; helices.len()];
for k in 1..helices.len() {
let h = helices[k];
let d = |alpha: isize, beta: isize| {
let mut ret: f64 = 0.;
for (x_0, b) in helices_map.get(&h).unwrap() {
if !nucl_set.contains(&(h, x_0 + alpha, !(*b))) {
return std::f64::INFINITY;
}
if xovers_map.get(&(h, *x_0, *b)).is_none() {
continue;
}
if xovers_map.get(&(h, *x_0, *b)).unwrap().0 == helices[k - 1] {
let (h_, y_0, b_) = xovers_map.get(&(h, *x_0, *b)).unwrap();
if !nucl_set.contains(&(*h_, y_0 + beta, !(*b_))) {
return std::f64::INFINITY;
}
let point_a = Point::from_coord(self.space_pos(h, x_0 + alpha, *b));
let point_b = Point::from_coord(self.space_pos(*h_, y_0 + beta, *b_));
let norm = (point_a - point_b).norm();
ret = ret.max(norm);
}
}
ret
};
for x in 0..delta.len() {
let mut arg_min = default_delta;
let mut val_min = std::f64::INFINITY;
for y in 0..delta.len() {
let val = d(delta[x], delta[y]).max(opt[(k - 1) as usize][y].0);
if val < val_min {
arg_min = y;
val_min = val;
}
}
opt[k][x] = (val_min, arg_min);
}
}
println!("{:?}", chain);
println!("helices {:?}", helices);
println!("{:?}", opt);
let k = helices.len() - 1;
let mut arg_min = default_delta;
let mut val_min = std::f64::INFINITY;
for x in 0..delta.len() {
if opt[k][x].0 < val_min {
arg_min = x;
val_min = opt[k][x].0;
}
}
let mut z = arg_min;
for k in (0..helices.len()).rev() {
let h = helices[k];
for (x, b) in helices_map.get(&h).unwrap() {
let nucl = (h, *x, *b);
let (s, d) = nucl_map.get(&nucl).unwrap();
let domain = &mut self.strands[*s].domains[*d];
println!("shift {}, {}, {}", h, x, delta[z]);
println!(
"before domain {}, {}, {}, {}",
domain.helix, domain.start, domain.end, domain.forward
);
if domain.start == *x {
domain.start += delta[z];
} else {
domain.end += delta[z];
}
println!(
"after domain {}, {}, {}, {}",
domain.helix, domain.start, domain.end, domain.forward
);
}
z = opt[k][z].1;
}
}
pub fn space_pos(&self, h: isize, x: isize, b: bool) -> [f64; 3] {
self.helices[h as usize].space_pos(&self.parameters.unwrap(), x, b)
}
pub fn get_cross_overs_helices(&self, h1: isize, h2: isize) -> Vec<(isize, bool, isize, bool)> {
let cross_overs = self.get_cross_overs();
let mut ret = Vec::new();
for (g1, x1, b1, g2, x2, b2) in cross_overs {
if g1 == h1 && g2 == h2 {
ret.push((x1, b1, x2, b2));
} else if g2 == h1 && g1 == h2 {
ret.push((x2, b2, x1, b1));
}
}
ret
}
fn get_total_dist_helices(&self, h1: isize, h2: isize, theta: f64, phi: f64) -> f64 {
let cross_overs = self.get_cross_overs_helices(h1, h2);
let helix1 = Helix {
roll: theta,
.. self.helices[h1 as usize].clone()
};
let helix2 = Helix {
roll: phi,
.. self.helices[h2 as usize].clone()
};
let mut ret = 0.;
for (x1, b1, x2, b2) in cross_overs {
let pos_1 = helix1.space_pos(&self.parameters.unwrap(), x1, b1);
let pos_2 = helix2.space_pos(&self.parameters.unwrap(), x2, b2);
ret += (pos_1[0] - pos_2[0]) * (pos_1[0] - pos_2[0]);
ret += (pos_1[1] - pos_2[1]) * (pos_1[1] - pos_2[1]);
ret += (pos_1[2] - pos_2[2]) * (pos_1[2] - pos_2[2]);
}
ret
}
pub fn adjust_helix(&mut self, h1: usize, h2: usize) {
let mut theta = self.helices[h1].roll;
let mut phi = self.helices[h2].roll;
loop {
let best = self.get_total_dist_helices(h1 as isize, h2 as isize, theta, phi);
if self.get_total_dist_helices(h1 as isize, h2 as isize, theta - 0.1, phi - 0.1) < best
{
theta -= 0.1;
phi -= 0.1;
} else if self.get_total_dist_helices(h1 as isize, h2 as isize, theta + 0.1, phi - 0.1)
< best
{
theta += 0.1;
phi -= 0.1;
} else if self.get_total_dist_helices(h1 as isize, h2 as isize, theta - 0.1, phi + 0.1)
< best
{
theta -= 0.1;
phi += 0.1;
} else if self.get_total_dist_helices(h1 as isize, h2 as isize, theta + 0.1, phi + 0.1)
< best
{
theta += 0.1;
phi += 0.1;
} else if self.get_total_dist_helices(h1 as isize, h2 as isize, theta, phi + 0.1) < best
{
phi += 0.1;
} else if self.get_total_dist_helices(h1 as isize, h2 as isize, theta, phi - 0.1) < best
{
phi -= 0.1;
} else if self.get_total_dist_helices(h1 as isize, h2 as isize, theta - 0.1, phi) < best
{
theta -= 0.1;
} else if self.get_total_dist_helices(h1 as isize, h2 as isize, theta + 0.1, phi) < best
{
theta += 0.1;
} else {
self.helices[h1].roll = theta;
self.helices[h2].roll = phi;
return;
}
}
}
pub fn get_overlapping_nucl(&self) -> HashMap<(isize, isize, bool), Vec<usize>> {
let nucl_map = self.get_nucl_map();
let mut ret: HashMap<(isize, isize, bool), Vec<usize>> = HashMap::new();
for (nucl, v) in nucl_map.iter() {
if v.len() > 1 {
ret.insert(*nucl, v.to_vec());
}
}
ret
}
pub fn get_nucl_map(&self) -> HashMap<(isize, isize, bool), Vec<usize>> {
let mut nucl_map: HashMap<(isize, isize, bool), Vec<usize>> = HashMap::new();
for i in 0..self.strands.len() {
let s = &self.strands[i];
for d in &s.domains {
for n in d.iter() {
let nucl = (d.helix, n, d.forward);
if let Some(v) = nucl_map.get_mut(&nucl) {
v.push(i);
} else {
nucl_map.insert(nucl, vec![i]);
}
}
}
}
nucl_map
}
}
const KELLY: [Color; 19] = [
Color::Int(0xF3C300),
Color::Int(0x875692),
Color::Int(0xA1CAF1),
Color::Int(0xBE0032),
Color::Int(0xC2B280),
Color::Int(0x848482),
Color::Int(0x008856),
Color::Int(0xE68FAC),
Color::Int(0x0067A5),
Color::Int(0xF99379),
Color::Int(0x604E97),
Color::Int(0xF6A600),
Color::Int(0xB3446C),
Color::Int(0xDCD300),
Color::Int(0x882D17),
Color::Int(0x8DB600),
Color::Int(0x654522),
Color::Int(0xE25822),
Color::Int(0x2B3D26),
];
impl<StrandLabel, DomainLabel> Strand<StrandLabel, DomainLabel> {
pub fn default_color(&self) -> Color {
for domain in self.domains.iter() {
let x1 = if domain.forward {
domain.end - 1
} else {
domain.start
};
let h = domain.helix as isize;
let x = x1 + (x1 % 11) + 5 * h;
let n = KELLY.len() as isize;
return KELLY[(((x % n) + n) % n) as usize].clone();
}
Color::Int(0)
}
pub fn set_kelly_color(&mut self, color: usize) {
self.color = Some(KELLY[color % KELLY.len()].clone())
}
pub fn is_empty(&self) -> bool {
for d in &self.domains {
if d.length() > 0 {
return false;
}
}
true
}
pub fn len(&self) -> usize {
self.domains
.iter()
.map(|d| ((d.end - d.start).max(0)) as usize)
.sum()
}
pub fn from_domain(domain: Domain<DomainLabel>) -> Self {
Strand {
domains: vec![domain],
label: None,
sequence: None,
cyclic: false,
color: None,
}
}
pub fn first_nucl(&self) -> Option<(isize, isize, bool)> {
if self.len() == 0 {
None
} else {
let mut i = 0;
while self.domains[i].length() == 0 {
i += 1;
}
let x = if self.domains[i].forward {
self.domains[i].start
} else {
self.domains[i].end - 1
};
Some((self.domains[i].helix, x, self.domains[i].forward))
}
}
pub fn last_nucl(&self) -> Option<(isize, isize, bool)> {
if self.len() == 0 {
None
} else {
let mut i = self.domains.len() - 1;
while self.domains[i].length() == 0 {
i -= 1;
}
let x = if self.domains[i].forward {
self.domains[i].end - 1
} else {
self.domains[i].start
};
Some((self.domains[i].helix, x, self.domains[i].forward))
}
}
pub fn nth_nucl(&self, n: usize) -> Option<(isize, isize, bool)> {
if !self.cyclic && self.len() <= n {
None
} else {
let mut nb_skipped = 0;
let mut i = 0;
while self.domains[i].length() + nb_skipped <= n as isize {
nb_skipped += self.domains[i].length();
i += 1;
i %= self.domains.len();
}
let dom = &self.domains[i];
if dom.forward {
Some((dom.helix, dom.start + n as isize - nb_skipped, true))
} else {
Some((dom.helix, dom.end - 1 - n as isize + nb_skipped, false))
}
}
}
pub fn set_sequence(&mut self, seq: String) {
self.sequence = Some(Cow::Owned(seq));
}
}
impl<'a: 'b, 'b, StrandLabel, DomainLabel> StrandRef<'a, StrandLabel, DomainLabel> {
pub fn id(&self) -> StrandId {
StrandId(self.strand_id)
}
pub fn design(&mut self) -> &mut Design<StrandLabel, DomainLabel> {
self.builder
}
pub fn len(&self) -> usize {
self.domains()
.iter()
.map(|d| (d.end - d.start) as usize)
.sum()
}
pub fn with_color(self, color: u32) -> Self {
self.builder.strands[self.strand_id].color = Some(Color::Int(color));
self
}
pub fn with_kelly_color(self, color: usize) -> Self {
self.builder.strands[self.strand_id].color = Some(KELLY[color].clone());
self
}
pub fn domains(&self) -> &[Domain<DomainLabel>] {
&self.builder.strands[self.strand_id].domains
}
pub fn domains_mut(&mut self) -> &mut Vec<Domain<DomainLabel>> {
&mut self.builder.strands[self.strand_id].domains
}
pub fn push_domain(&mut self, domain: Domain<DomainLabel>) {
self.x = Some(if domain.forward {
domain.end - 1
} else {
domain.start
});
self.builder.strands[self.strand_id].domains.push(domain);
}
pub fn extend<I: Iterator<Item = Domain<DomainLabel>>>(&mut self, i: I) {
for i in i {
self.push_domain(i)
}
}
pub fn with_sequence<I: Into<Cow<'static, str>>>(self, sequence: I) -> Self {
self.builder.strands[self.strand_id].sequence = Some(sequence.into());
self
}
pub fn with_domain_sequence<I: Into<Cow<'static, str>>>(self, sequence: I) -> Self {
self.builder.strands[self.strand_id]
.domains
.last_mut()
.unwrap()
.sequence = Some(sequence.into());
self
}
pub fn to(mut self, x: isize) -> Self {
let strand = self.builder.strands[self.strand_id]
.domains
.last_mut()
.unwrap();
self.x = Some(x);
let x0 = self.x0.unwrap();
if x > x0 {
strand.forward = true;
strand.start = x0;
strand.end = x + 1;
} else {
strand.forward = false;
strand.start = x;
strand.end = x0 + 1;
}
self
}
pub fn to_point(mut self, point: [f32; 3]) -> Self {
let strand = self.builder.strands[self.strand_id]
.domains
.last_mut()
.unwrap();
let x = self.builder.helices[strand.helix as usize].closest_nucl(
point,
&self
.builder
.parameters
.expect("parameters is None in to_point"),
);
self.x = Some(x);
let x0 = self.x0.unwrap();
if x > x0 {
strand.forward = true;
strand.start = x0;
strand.end = x + 1;
} else {
strand.forward = false;
strand.start = x;
strand.end = x0 + 1;
}
self
}
pub fn pop(mut self) -> Self {
self.builder.strands[self.strand_id].domains.pop();
if let Some(last) = self.builder.strands[self.strand_id].domains.last() {
self.x = Some(if last.forward { last.end - 1 } else { last.start })
} else {
self.x = None
}
self
}
pub fn rev_to(mut self, x: isize) -> Self {
if self.builder.strands[self.strand_id].domains.is_empty() {
return self;
}
let ref mut strand = self.builder.strands[self.strand_id].domains[0];
self.x = Some(x);
let x0 = self.x0.unwrap();
if x > x0 {
strand.forward = true;
strand.start = x0;
strand.end = x + 1;
} else {
strand.forward = false;
strand.start = x;
strand.end = x0 + 1;
}
self
}
pub fn with_label<L: Into<StrandLabel>>(self, label: L) -> Self {
self.builder.strands[self.strand_id].label = Some(label.into());
self
}
pub fn with_domain_label<L: Into<DomainLabel>>(self, label: L) -> Self {
self.builder.strands[self.strand_id]
.domains
.last_mut()
.unwrap()
.label = Some(label.into());
self
}
pub fn next_domain_to(mut self, to: isize) -> Self {
let helix = {
let last = self.builder.strands[self.strand_id].domains.last().unwrap();
if last.end == last.start {
return self.to(to);
}
last.helix
};
let (start, end, forward) = if to < self.x.unwrap() {
let n = (to, self.x.unwrap(), false);
self.x = Some(to);
n
} else {
let n = (self.x.unwrap() + 1, to + 1, true);
self.x = Some(to);
n
};
let domain = Domain {
helix,
start,
end,
forward,
label: None,
sequence: None,
};
self.builder.strands[self.strand_id].domains.push(domain);
self
}
pub fn cross<H: Into<HelixId>>(self, h: H) -> Self {
let x = self.x.unwrap();
self.cross_to(h, x)
}
pub fn cross_to<H: Into<HelixId>>(mut self, h: H, start: isize) -> Self {
let helix = h.into().0;
if helix >= self.builder.helices.len() as isize {
panic!("Crossing to an undefined helix")
}
self.builder.strands[self.strand_id].domains.push(Domain {
helix,
start,
end: start,
forward: true,
label: None,
sequence: None,
});
self.x0 = Some(start);
self.x = Some(start);
self
}
pub fn hflip(self, x: isize) -> Self {
for dom in self.builder.strands[self.strand_id].domains.iter_mut() {
let start = dom.start;
dom.start = x - dom.end;
dom.end = x - start;
}
self.builder.strands[self.strand_id].domains.reverse();
self
}
pub fn vflip<H: Into<HelixId>>(self, y: H) -> Self {
let y = y.into();
for dom in self.builder.strands[self.strand_id].domains.iter_mut() {
assert!(y.0 >= dom.helix);
dom.helix = 2 * y.0 - dom.helix
}
self
}
pub fn translate(self, x: isize, y: isize) -> Self {
use std::convert::TryFrom;
for dom in self.builder.strands[self.strand_id].domains.iter_mut() {
dom.helix = usize::try_from(dom.helix as isize + y).unwrap() as isize;
dom.start += x;
dom.end += x;
}
self
}
pub fn cycle(self) -> () {
let id = self.strand_id;
let helix = self.builder.strands[id].domains[0].helix;
let first_orientation;
let start = if self.builder.strands[id].domains[0].forward {
first_orientation = true;
self.builder.strands[id].domains[0].start
} else {
first_orientation = false;
self.builder.strands[id].domains[0].end - 1
};
self.builder.strands[id].cyclic = true;
let curent_dom = self.builder.strands[id].domains.last().unwrap();
let builder = if curent_dom.helix != helix {
if curent_dom.forward && curent_dom.end - 1 >= start
|| !curent_dom.forward && curent_dom.end - 1 <= start
{
self.cross(helix).to(start)
} else if curent_dom.forward != first_orientation {
self.to(start).cross(helix)
} else {
self.cross_to(helix, start)
}
} else {
self.to(start)
};
let last_dom = builder.builder.strands[id].domains.last_mut().unwrap();
last_dom.forward = first_orientation;
}
pub fn reverse(self) -> Self {
self.builder.strands[self.strand_id].domains.reverse();
for dom in self.builder.strands[self.strand_id].domains.iter_mut() {
dom.forward = !dom.forward
}
self
}
}
fn sequence<StrandLabel, DomainLabel>(
positions: &HashMap<(usize, isize, bool), char>,
strand: &Strand<StrandLabel, DomainLabel>,
) -> String {
let mut seq = String::new();
for dom in strand.domains.iter() {
if let Some(ref s) = dom.sequence {
seq.push_str(s);
continue;
}
assert!(dom.helix >= 0);
for pos in dom.iter() {
seq.push(sequences::complement(
*positions
.get(&(dom.helix as usize, pos, !dom.forward))
.unwrap(),
))
}
}
seq
}
pub mod transforms {
use super::Domain;
pub fn rotate<Label>(domains: &mut [Domain<Label>]) {
let max_helix = domains.iter().map(|x| x.helix).max().unwrap();
for dom in domains.iter_mut() {
let start = dom.start;
dom.start = -dom.end;
dom.end = -start;
dom.helix = max_helix - dom.helix;
dom.forward = !dom.forward;
}
}
pub fn v_flip<Label>(domains: &mut [Domain<Label>], max_helix: isize) {
for dom in domains.iter_mut() {
dom.helix = max_helix - dom.helix;
}
}
pub fn reverse<Label>(domains: &mut [Domain<Label>]) {
for dom in domains.iter_mut() {
dom.forward = !dom.forward;
}
}
pub fn h_flip<Label>(domains: &mut [Domain<Label>]) {
for dom in domains.iter_mut() {
let start = dom.start;
dom.start = -dom.end;
dom.end = -start;
}
}
pub fn translate<Label>(domains: &mut [Domain<Label>], dh: isize, dx: isize) {
use std::convert::TryFrom;
for dom in domains.iter_mut() {
dom.helix = usize::try_from(dom.helix as isize + dh).unwrap() as isize;
dom.start += dx;
dom.end += dx;
}
}
}
impl<StrandLabel, DomainLabel> Design<StrandLabel, DomainLabel> {
fn positions(&self) -> HashMap<(usize, isize, bool), char> {
let mut positions = HashMap::new();
for s in self.strands.iter() {
if let Some(ref seq) = s.sequence {
let mut seq = seq.chars();
for dom in s.domains.iter() {
if dom.helix < 0 {
continue;
}
for pos in dom.iter() {
let c = seq.next().unwrap();
positions.insert((dom.helix as usize, pos, dom.forward), c);
}
}
}
}
positions
}
pub fn sequences(&self) -> Vec<String> {
let positions = self.positions();
let mut sequences = Vec::new();
for s in self.strands.iter() {
if s.sequence.is_none() {
sequences.push(sequence(&positions, s))
}
}
sequences
}
}
#[cfg(feature = "excel")]
impl<StrandLabel: Ord + std::fmt::Display, DomainLabel> Design<StrandLabel, DomainLabel> {
pub fn make_plates_96<
H: FnMut(&Strand<StrandLabel, DomainLabel>) -> usize,
F: FnMut(&Strand<StrandLabel, DomainLabel>) -> Option<bool>,
G: FnMut(usize) -> String,
>(
&mut self,
file: &str,
mut seq_class: H,
mut filter: F,
mut format: G,
) {
self.strands
.sort_by(|a, b| seq_class(&a).cmp(&seq_class(b)));
let mut current_plate = 0;
let mut current_class = 0;
let mut plate = Plate96::new();
let mut plates = vec![Vec::new()];
let positions = self.positions();
for s in self.strands.iter() {
let class = seq_class(s);
if class != current_class && plate.row > 0 {
plate.incr_col();
}
if plate.n != current_plate {
plates.push(Vec::new());
current_plate = plate.n;
}
if let Some(x) = filter(s) {
if x && (plate.row != 0 || plate.column != 0) {
plate.incr_plate();
plates.push(Vec::new());
current_plate = plate.n;
}
plates.last_mut().unwrap().push((
plate.column,
plate.row,
&s.label,
sequence(&positions, s),
));
plate.incr()
}
current_class = class;
}
use simple_excel_writer::*;
let mut wb = Workbook::create(file);
for (i, seqs) in plates.iter().enumerate() {
let mut sheet = wb.create_sheet(&format(i));
sheet.add_column(Column { width: 30.0 });
sheet.add_column(Column { width: 30.0 });
sheet.add_column(Column { width: 30.0 });
wb.write_sheet(&mut sheet, |sw| {
sw.append_row(row!["Well Position", "Name", "Sequence"])?;
for &(col, row, label, ref seq) in seqs.iter() {
let row = (b'A' + row as u8) as char;
let label = if let Some(ref label) = label {
format!("{}", label)
} else {
format!("{}{}", row, col + 1)
};
sw.append_row(row![
format!("{}{}", row, col + 1).as_str(),
label.as_str(),
seq.as_str()
])?
}
Ok(())
})
.unwrap()
}
wb.close().unwrap();
}
}
#[cfg(feature = "excel")]
#[derive(Debug, Copy, Clone)]
pub struct Plate96 {
pub n: usize,
pub column: usize,
pub row: usize,
}
#[cfg(feature = "excel")]
impl Plate96 {
pub fn new() -> Self {
Plate96 {
n: 0,
column: 0,
row: 0,
}
}
pub fn incr(&mut self) {
self.row += 1;
if self.row >= 8 {
self.incr_col()
}
}
pub fn incr_col(&mut self) {
self.column += 1;
self.row = 0;
if self.column >= 12 {
self.incr_plate()
}
}
pub fn incr_plate(&mut self) {
self.column = 0;
self.row = 0;
self.n += 1;
}
}