Commit a5fcc567 authored by nlevy's avatar nlevy
Browse files

Fix checking of sequence coverage

parent 5ec54913
......@@ -7,7 +7,7 @@ const INSERTION_MARKER: &'static str = "**INSERTION**";
const TODO_MARKER: &'static str = "**TODO**";
use colored::Colorize;
use std::collections::BTreeSet;
use std::collections::{BTreeMap, BTreeSet};
use std::env::Args;
fn parse_command_line_arguments(mut args: Args) -> OrigamiInput {
......@@ -18,7 +18,7 @@ fn parse_command_line_arguments(mut args: Args) -> OrigamiInput {
let stapples_file_name = args.next().expect(USAGE);
let scaffold_sequence = match std::fs::read_to_string(scaffold_file_name) {
Ok(content) => content,
Ok(content) => content.trim_end().to_string(),
Err(e) => {
panic!("Could not read scaffold sequence {e}");
}
......@@ -66,12 +66,13 @@ fn parse_staple_file(content: String) -> Vec<Staple> {
fn check_coverage(input: &OrigamiInput) {
let mut nucls: BTreeSet<usize> = (0..input.scaffold_sequence.len()).collect();
let mut used_by: BTreeMap<usize, (usize, Staple)> = Default::default();
let mut success = true;
print!("Checking coverage...");
for staple in input.staples.iter() {
for (staple_id, staple) in input.staples.iter().enumerate() {
for (a, b) in staple.domains.iter().cloned() {
let range: Box<dyn Iterator<Item = isize>> = if b < a {
let range: Box<dyn Iterator<Item = isize>> = if b <= a {
Box::new((b as isize)..=(a as isize))
} else {
Box::new((0..=(a as isize)).chain(b..(input.scaffold_sequence.len() as isize)))
......@@ -80,6 +81,17 @@ fn check_coverage(input: &OrigamiInput) {
if !nucls.remove(&(n as usize)) {
success = false;
let msg = format!("The nucleotide {n} is used several times");
println!("{}", msg.red());
let dbg_msg = format!("it appears in domain [{a}, {b}]");
println!("{}", dbg_msg.red())
}
if let Some(previous_staple) =
used_by.insert(n as usize, (staple_id, staple.clone()))
{
let msg = format!(
"Used by {staple_id} \\ Previously used by staple {:?}",
previous_staple
);
println!("{}", msg.red())
}
}
......@@ -123,7 +135,7 @@ fn check_sequences(input: &OrigamiInput) {
let stapple_chars: Vec<char> = trimed_seq.chars().collect();
let mut stap_idx = 0;
for (a, b) in staple.domains.iter().cloned() {
let range: Box<dyn Iterator<Item = isize>> = if b < a {
let range: Box<dyn Iterator<Item = isize>> = if b <= a {
Box::new(((b as isize)..=(a as isize)).rev())
} else {
Box::new(
......@@ -139,6 +151,7 @@ fn check_sequences(input: &OrigamiInput) {
if stap_char != expected {
success = false;
println!("{:?}", staple);
let msg = format!("At pos {scaff_idx}, expected {expected}, got {stap_char}");
println!("{}", msg.red());
}
......@@ -162,6 +175,7 @@ struct OrigamiInput {
staples: Vec<Staple>,
}
#[derive(Clone, Debug)]
struct Staple {
domains: Vec<(isize, isize)>,
sequence: String,
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment