From 643f56db9b7fabdf48a2c9b6093b50f444789872 Mon Sep 17 00:00:00 2001 From: Jean Pierre ELBERS Date: Wed, 27 Oct 2021 16:14:09 +0200 Subject: [PATCH 1/2] added initial bedcov_rs1_rustbio.rs but coverage values do not match other implementations --- Cargo.toml | 6 +- bedcov/Makefile | 12 ++- bedcov/bedcov_rs1_rustbio.rs | 164 +++++++++++++++++++++++++++++++++++ 3 files changed, 180 insertions(+), 2 deletions(-) create mode 100644 bedcov/bedcov_rs1_rustbio.rs diff --git a/Cargo.toml b/Cargo.toml index e5ee09c..160dffe 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -6,7 +6,7 @@ edition = "2018" [dependencies] flate2 = "1.0.14" -bio = "0.31.0" +bio = "0.39.0" needletail = "0.4" [[bin]] @@ -16,3 +16,7 @@ path = "fqcnt/fqcnt_rs1_rustbio.rs" [[bin]] name = "fqcnt_needletail" path = "fqcnt/fqcnt_rs2_needletail.rs" + +[[bin]] +name = "bedcov_rustbio" +path = "bedcov/bedcov_rs1_rustbio.rs" diff --git a/bedcov/Makefile b/bedcov/Makefile index 458f41e..dfb672b 100644 --- a/bedcov/Makefile +++ b/bedcov/Makefile @@ -1,10 +1,12 @@ CFLAGS=-g -Wall -O2 LIBS=-lz -PROG=bedcov_c1_cgr bedcov_nim1_klib bedcov_cr1_klib bedcov_go1 bedcov_d_dgr bedcov_fsharp1 +PROG=bedcov_c1_cgr bedcov_nim1_klib bedcov_cr1_klib bedcov_go1 bedcov_rs bedcov_d_dgr bedcov_fsharp1 NIM=nim CRYSTAL=crystal GO=go +CARGO=cargo LDC=ldc2 +PROJECT_ROOT := $(shell realpath ../) .PHONY:all clean @@ -22,6 +24,14 @@ bedcov_cr1_klib:bedcov_cr1_klib.cr ../lib/klib.cr bedcov_go1:bedcov_go1.go $(GO) build $< +.PHONY: bedcov_rs +bedcov_rs: + $(CARGO) install --path $(PROJECT_ROOT) --force --bin bedcov_rustbio --root . + +.PHONY: rust +rust: fqcnt_rs + @echo "Rust binaries are in ./bin/" + bedcov_d_dgr:bedcov_d_dgr.d $(LDC) -O3 -release -flto=full $< diff --git a/bedcov/bedcov_rs1_rustbio.rs b/bedcov/bedcov_rs1_rustbio.rs new file mode 100644 index 0000000..5691a99 --- /dev/null +++ b/bedcov/bedcov_rs1_rustbio.rs @@ -0,0 +1,164 @@ +use std::convert::TryInto; +use std::path::PathBuf; +use bio::data_structures::interval_tree::IntervalTree; +use std::collections::HashMap; +use bio::io::bed; +use std::env; +use std::fs::File; + +fn main() { + // Get the first bed file's path + // if missing, then provide error + let path1: PathBuf = match env::args().nth(1) { + Some(p1) => PathBuf::from(p1), + _ => { + eprintln!("Usage: {} ", file!()); + std::process::exit(1) + } + }; + + // Get the second bed file's path + // if missing, then provide error + let path2: PathBuf = match env::args().nth(2) { + Some(p2) => PathBuf::from(p2), + _ => { + eprintln!("Usage: {} ", file!()); + std::process::exit(1) + } + }; + + // Open the first bed file's path + let file1 = match File::open(&path1) { + Ok(fh1) => fh1, + Err(err) => { + eprintln!("Failed to open file: {}", err); + std::process::exit(1) + } + }; + + // Open the second bed file's path + let file2 = match File::open(&path2) { + Ok(fh2) => fh2, + Err(err) => { + eprintln!("Failed to open file: {}", err); + std::process::exit(1) + } + }; + + let mut reader1 = bed::Reader::new(file1); + + // make the HashMap + let mut trees = HashMap::new(); + + // Make a vector of contigs + + let mut contigs = Vec::new(); + + // Get the unique contigs + for record in reader1.records() { + let rec = record.expect("Error reading record."); +// println!("rec.start() = {}, rec.end() = {}, rec.chrom() = {}", rec.start(), rec.end(), rec.chrom().to_string()); + contigs.push(rec.chrom().to_string()); + } + + contigs.sort_by(|a, b| a.partial_cmp(b).expect("NaN in vector")); + contigs.dedup(); + + // Open the first bed file's path + let file1 = match File::open(&path1) { + Ok(fh1) => fh1, + Err(err) => { + eprintln!("Failed to open file: {}", err); + std::process::exit(1) + } + }; + + let mut reader1 = bed::Reader::new(file1); + + for contig in contigs { +// println!("contig = {}", contig); + let mut currentcontig = Vec::new(); + for record1 in reader1.records() { + let rec1 = record1.expect("Error reading record."); +// println!("rec1.start() = {}, rec1.end() = {}, contig = {}", rec1.start(), rec1.end(), contig); + if rec1.chrom().to_string() == contig { + currentcontig.push(rec1.start()-1); + currentcontig.push(rec1.end()+1); +// println!("rec1.start() = {}, rec1.end() = {}, contig = {}", rec1.start(), rec1.end(), contig); + } + } + let mut tree = IntervalTree::new(); + let mut i=0; + for s in currentcontig.windows(2) { + i+=1; + if i % 2 == 0 { + // n is even + } + else { + let [a, b]: [u64; 2] = s.try_into().unwrap(); +// println!("{}\t{}", a, b); + tree.insert(a..b, contig.to_string()); + } + } + trees.insert(contig.to_string(),tree); + } + + let mut reader2 = bed::Reader::new(file2); + + + // initialize coverage vectors + let mut currentcontig2 = Vec::new(); + let mut currentcontig2start = Vec::new(); + let mut currentcontig2end = Vec::new(); + let mut currentcontig2length = Vec::new(); + + // Read through second bed file entry by entry and then + // and try to find overlaps + + let mut count=0; + for record2 in reader2.records() { + let rec2 = record2.expect("Error reading record."); + currentcontig2.push(rec2.chrom().to_string()); + currentcontig2start.push(rec2.start()); + currentcontig2end.push(rec2.end()); + currentcontig2length.push(rec2.end()-rec2.start()); + + if let Some(tree) = trees.get(&rec2.chrom().to_string()) { + tree.find(rec2.start()..rec2.end()); +/* currentcontig2.push(rec2.chrom().to_string()); + currentcontig2start.push(rec2.start()); + currentcontig2end.push(rec2.end()); + currentcontig2length.push(rec2.end()-rec2.start()); + println!("{:?}", currentcontig2); + println!("{:?}", currentcontig2start); + println!("{:?}", currentcontig2end); + println!("{:?}", currentcontig2length); + println!("count={}; index={}", count, currentcontig2.len());*/ + if currentcontig2.len() == 1 { + if currentcontig2[currentcontig2.len()-1] == rec2.chrom().to_string() { + if currentcontig2start[currentcontig2.len()-1] == rec2.start() { + if currentcontig2end[currentcontig2.len()-1] == rec2.end() { + if currentcontig2length[currentcontig2.len()-1] == rec2.end() - rec2.start() { + count+=1; + } + } + } + } + } + else { + if currentcontig2[currentcontig2.len()-2] == rec2.chrom().to_string() { + if currentcontig2start[currentcontig2.len()-2] == rec2.start() { + if currentcontig2end[currentcontig2.len()-2] == rec2.end() { + if currentcontig2length[currentcontig2.len()-2] == rec2.end() - rec2.start() { + count+=1; + } + } + } + } + } + println!("{}\t{}\t{}\t{}\t{}", currentcontig2[currentcontig2.len()-1], currentcontig2start[currentcontig2.len()-1], currentcontig2end[currentcontig2.len()-1], count, currentcontig2length[currentcontig2.len()-1]); + count=1; + } + + } +} From d3fd450f48208ddbe054c60e5787ed43903df7de Mon Sep 17 00:00:00 2001 From: Jean Elbers Date: Thu, 9 Dec 2021 15:10:43 +0100 Subject: [PATCH 2/2] fixed counting to match bedcov_c1_cgr --- bedcov/bedcov_rs1_rustbio.rs | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/bedcov/bedcov_rs1_rustbio.rs b/bedcov/bedcov_rs1_rustbio.rs index 5691a99..67e1903 100644 --- a/bedcov/bedcov_rs1_rustbio.rs +++ b/bedcov/bedcov_rs1_rustbio.rs @@ -115,8 +115,9 @@ fn main() { // Read through second bed file entry by entry and then // and try to find overlaps - let mut count=0; +// let mut count=0; for record2 in reader2.records() { + let mut count=0; let rec2 = record2.expect("Error reading record."); currentcontig2.push(rec2.chrom().to_string()); currentcontig2start.push(rec2.start()); @@ -124,12 +125,12 @@ fn main() { currentcontig2length.push(rec2.end()-rec2.start()); if let Some(tree) = trees.get(&rec2.chrom().to_string()) { - tree.find(rec2.start()..rec2.end()); -/* currentcontig2.push(rec2.chrom().to_string()); + for r in tree.find(rec2.start()..rec2.end()) { + currentcontig2.push(rec2.chrom().to_string()); currentcontig2start.push(rec2.start()); currentcontig2end.push(rec2.end()); currentcontig2length.push(rec2.end()-rec2.start()); - println!("{:?}", currentcontig2); +/* println!("{:?}", currentcontig2); println!("{:?}", currentcontig2start); println!("{:?}", currentcontig2end); println!("{:?}", currentcontig2length); @@ -156,8 +157,8 @@ fn main() { } } } + } println!("{}\t{}\t{}\t{}\t{}", currentcontig2[currentcontig2.len()-1], currentcontig2start[currentcontig2.len()-1], currentcontig2end[currentcontig2.len()-1], count, currentcontig2length[currentcontig2.len()-1]); - count=1; } }