|
|
@@ -303,48 +303,12 @@ pub fn par_whole_scan(dict_file: &str, bam_path: &str, out_dir: &str) -> anyhow:
|
|
|
Ok(())
|
|
|
}
|
|
|
|
|
|
-// pub fn mrd_count_ratio(mrd_bam: &str, count_file: &str) -> anyhow::Result<()> {
|
|
|
-// let count_file = File::open(count_file)?;
|
|
|
-// let reader = BufReader::new(count_file);
|
|
|
-// let mut tumoral_counts = Vec::new();
|
|
|
-// for line in reader.lines() {
|
|
|
-// let line_split: Vec<&str> = line?.split("\t").collect();
|
|
|
-// let pos = *line_split.get(0).context("Can't parse the tsv file")?;
|
|
|
-// let (contig, pos) = pos.split_once(':').context("Can't parse the tsv file")?;
|
|
|
-// let (start, end) = pos.split_once('-').context("Can't parse the tsv file")?;
|
|
|
-// let start: u32 = start.parse()?;
|
|
|
-// let end: u32 = end.parse()?;
|
|
|
-// let n: u32 = line_split
|
|
|
-// .get(1)
|
|
|
-// .context("Can't parse the tsv file")?
|
|
|
-// .parse::<u32>()?;
|
|
|
-// tumoral_counts.push((contig.to_string(), start, end, n));
|
|
|
-// }
|
|
|
-//
|
|
|
-// let ratios: Vec<(usize, String, u32, u32, u32, u32, Option<f64>)> = tumoral_counts
|
|
|
-// .par_iter()
|
|
|
-// .enumerate()
|
|
|
-// .map(|(i, (contig, start, end, n))| {
|
|
|
-// let bin = Bin::new(mrd_bam, contig, *start, *end - *start + 1).unwrap();
|
|
|
-// let n_mrd = bin.n_reads() as u32;
|
|
|
-// let r = if n_mrd == 0 {
|
|
|
-// None
|
|
|
-// } else {
|
|
|
-// Some((*n as f64 / n_mrd as f64).log10())
|
|
|
-// };
|
|
|
-// (i, contig.to_string(), *start, *end, *n, n_mrd, r)
|
|
|
-// })
|
|
|
-// .collect();
|
|
|
-//
|
|
|
-// // filter_outliers_modified_z_score_with_indices(ratios.par_iter().map(|(_, _, _, _, _, _, r)))
|
|
|
-// Ok(())
|
|
|
-// }
|
|
|
-
|
|
|
#[cfg(test)]
|
|
|
mod tests {
|
|
|
|
|
|
- use counts::{ranges_over, ranges_under, Counts};
|
|
|
- use plotly::{color::NamedColor, common::Marker, Bar, Plot, Scatter};
|
|
|
+ use counts::{ranges_between, ranges_over, ranges_under, Counts};
|
|
|
+ use pandora_lib_graph::cytoband::{svg_chromosome, AdditionalRect, RectPosition};
|
|
|
+ use plotly::{common::Marker, Bar, Plot};
|
|
|
use rust_htslib::bam::Reader;
|
|
|
|
|
|
use super::*;
|
|
|
@@ -377,9 +341,6 @@ mod tests {
|
|
|
let bam_path = format!("/data/longreads_basic_pipe/{id}/diag/{id}_diag_hs1.bam");
|
|
|
let out_dir = format!("/data/longreads_basic_pipe/{id}/diag/scan");
|
|
|
par_whole_scan("/data/ref/hs1/chm13v2.0.dict", &bam_path, &out_dir).unwrap();
|
|
|
- // let bam_path = format!("/data/longreads_basic_pipe/{id}/mrd/{id}_mrd_hs1.bam");
|
|
|
- // let out_dir = format!("/data/longreads_basic_pipe/{id}/mrd/scan");
|
|
|
- // par_whole_scan("/data/ref/hs1/chm13v2.0.dict", &bam_path, &out_dir).unwrap();
|
|
|
}
|
|
|
|
|
|
#[test]
|
|
|
@@ -448,10 +409,11 @@ mod tests {
|
|
|
fn load() -> anyhow::Result<()> {
|
|
|
init();
|
|
|
info!("loading");
|
|
|
- let count_file = "/data/longreads_basic_pipe/ROBIN/diag/scan/chr1_counts.tsv";
|
|
|
+ let contig = "chr22";
|
|
|
+ let count_file = &format!("/data/longreads_basic_pipe/ROBIN/diag/scan/{contig}_counts.tsv");
|
|
|
let counts = Counts::from_files(vec![count_file]);
|
|
|
|
|
|
- let chr1_nd_reads = counts.nd_reads("chr1")?;
|
|
|
+ let chr1_nd_reads = counts.nd_reads(contig)?;
|
|
|
println!(
|
|
|
"Percentiles: 1% {}, 50% {}, 99% {}",
|
|
|
chr1_nd_reads.percentile(1.0).unwrap(),
|
|
|
@@ -460,12 +422,61 @@ mod tests {
|
|
|
);
|
|
|
println!("< 6x: {:.2}%", chr1_nd_reads.proportion_under(6) * 100.0);
|
|
|
println!("> 15x: {:.2}%", chr1_nd_reads.proportion_above(15) * 100.0);
|
|
|
- let d = counts.get("chr1")?;
|
|
|
-
|
|
|
- let all_ranges = ranges_over(&d, 6, 10);
|
|
|
- println!("Ranges over 1 : {:?}", all_ranges.len());
|
|
|
- let all_ranges = ranges_under(&d, 6, 10);
|
|
|
- println!("Ranges under 1 : {:?}", all_ranges.len());
|
|
|
+ let d = counts.get(contig)?;
|
|
|
+
|
|
|
+ let tol = 20;
|
|
|
+
|
|
|
+ let under_6_rects: Vec<AdditionalRect> = ranges_under(&d, 6, tol)
|
|
|
+ .iter()
|
|
|
+ .filter(|(s, e)| e > s)
|
|
|
+ // .filter(|(s, e)| e - s > tol)
|
|
|
+ .map(|(start, end)| AdditionalRect {
|
|
|
+ start: *start as u32 * 1000,
|
|
|
+ end: *end as u32 * 1000,
|
|
|
+ color: String::from("red"),
|
|
|
+ position: RectPosition::Below(0),
|
|
|
+ })
|
|
|
+ .collect();
|
|
|
+
|
|
|
+ let over_6_rects: Vec<AdditionalRect> = ranges_between(&d, 6, 15, tol)
|
|
|
+ .iter()
|
|
|
+ .filter(|(s, e)| e > s)
|
|
|
+ // .filter(|(s, e)| e - s > tol)
|
|
|
+ .map(|(start, end)| AdditionalRect {
|
|
|
+ start: *start as u32 * 1000,
|
|
|
+ end: *end as u32 * 1000,
|
|
|
+ color: String::from("green"),
|
|
|
+ position: RectPosition::Below(1),
|
|
|
+ })
|
|
|
+ .collect();
|
|
|
+
|
|
|
+ let over15: Vec<AdditionalRect> = ranges_over(&d, 15, tol)
|
|
|
+ .iter()
|
|
|
+ .filter(|(s, e)| e > s)
|
|
|
+ // .filter(|(s, e)| e - s > tol)
|
|
|
+ .map(|(start, end)| AdditionalRect {
|
|
|
+ start: *start as u32 * 1000,
|
|
|
+ end: *end as u32 * 1000,
|
|
|
+ color: String::from("blue"),
|
|
|
+ position: RectPosition::Below(2),
|
|
|
+ })
|
|
|
+ .collect();
|
|
|
+
|
|
|
+ let mut all = Vec::new();
|
|
|
+ all.extend(under_6_rects);
|
|
|
+ all.extend(over_6_rects);
|
|
|
+ all.extend(over15);
|
|
|
+
|
|
|
+ svg_chromosome(
|
|
|
+ contig,
|
|
|
+ 1000,
|
|
|
+ 50,
|
|
|
+ "/data/ref/hs1/cytoBandMapped.bed",
|
|
|
+ "/data/chr1.svg",
|
|
|
+ &all,
|
|
|
+ &Vec::new(),
|
|
|
+ )
|
|
|
+ .unwrap();
|
|
|
|
|
|
let mut plot = Plot::new();
|
|
|
|
|
|
@@ -497,50 +508,9 @@ mod tests {
|
|
|
|
|
|
plot.add_trace(bars);
|
|
|
|
|
|
- // let trace = Scatter::new(
|
|
|
- // data.iter().map(|(x, _)| *x).collect::<Vec<f64>>(),
|
|
|
- // data.iter().map(|(_, y)| *y).collect::<Vec<f64>>(),
|
|
|
- // );
|
|
|
- // plot.add_trace(trace);
|
|
|
- //
|
|
|
plot.use_local_plotly();
|
|
|
-
|
|
|
plot.write_image("/data/test2.svg", plotly::ImageFormat::SVG, 800, 600, 1.0);
|
|
|
|
|
|
- // plot.write_html("out.html");
|
|
|
Ok(())
|
|
|
}
|
|
|
-
|
|
|
- // #[test]
|
|
|
- // fn phasing() -> anyhow::Result<()> {
|
|
|
- // let id = "SALICETTO";
|
|
|
- // let min_records = 2;
|
|
|
- //
|
|
|
- // let logger =
|
|
|
- // env_logger::Builder::from_env(env_logger::Env::default().default_filter_or("info"))
|
|
|
- // .build();
|
|
|
- // let multi = MultiProgress::new();
|
|
|
- // LogWrapper::new(multi.clone(), logger).try_init().unwrap();
|
|
|
- //
|
|
|
- // let config = PhaserConfig::new(id, "/data/longreads_basic_pipe", min_records, 0.33);
|
|
|
- // phase(config, multi)
|
|
|
- // }
|
|
|
- //
|
|
|
- // #[test]
|
|
|
- // fn load_phase() -> anyhow::Result<()> {
|
|
|
- // init();
|
|
|
- // let id = "SALICETTO";
|
|
|
- // let contig = "chr7";
|
|
|
- // let phases_dir = format!("/data/longreads_basic_pipe/{id}/diag/phases");
|
|
|
- // let phase_path = format!("{phases_dir}/{id}_{contig}_phases.postcard.gz");
|
|
|
- // let p = load_phases(&phase_path)?;
|
|
|
- // info!("{} phases", p.len());
|
|
|
- //
|
|
|
- // for phase in p {
|
|
|
- // if let Some(phase_id) = &phase.id {
|
|
|
- // info!("{}\t{}", phase_id, phase.mean_vaf());
|
|
|
- // }
|
|
|
- // }
|
|
|
- // Ok(())
|
|
|
- // }
|
|
|
}
|
Your Name