Add new parameter --max-read-depth to the vcf-report (#76)

* Add max-read-depth to CLI

* Change rows to HashSet

* fmt

* Update src/bcf/report/table_report/static_reader.rs

Co-authored-by: Till Hartmann <till.hartmann@udo.edu>

* Update src/bcf/report/table_report/static_reader.rs

Co-authored-by: Till Hartmann <till.hartmann@udo.edu>

* Addition to suggestion from @tedil

Co-authored-by: Till Hartmann <till.hartmann@udo.edu>

Author: fxwiegand

Cargo.toml

@@ -21,7 +21,8 @@ rustc-serialize = "0.3"
 quick-error = "1.2"
 newtype_derive = "0.1"
 custom_derive = "0.1"
-rand = "0.4"
+rand = "0.7.3"
+rand_core = "0.5.1"
 cogset = "0.2"
 num-bigint = "0.2"
 serde = "1.0"

src/bcf/report/table_report/create_report_table.rs

@@ -37,6 +37,7 @@ pub struct Report {
     vis: HashMap<String, String>,
 }
 
+#[allow(clippy::too_many_arguments)]
 pub(crate) fn make_table_report(
     vcf_path: &Path,
     fasta_path: &Path,
@@ -45,6 +46,7 @@ pub(crate) fn make_table_report(
     formats: Option<Vec<&str>>,
     sample: String,
     output_path: &str,
+    max_read_depth: u32,
 ) -> Result<(), Box<dyn Error>> {
     // HashMap<gene: String, Vec<Report>>, Vec<ann_field_identifiers: String>
     let mut reports = HashMap::new();
@@ -292,6 +294,7 @@ pub(crate) fn make_table_report(
                             chrom.clone(),
                             0,
                             end_position as u64 + 75,
+                            max_read_depth,
                         );
                         visualization =
                             manipulate_json(content, 0, end_position as u64 + 75, max_rows);
@@ -303,6 +306,7 @@ pub(crate) fn make_table_report(
                             chrom.clone(),
                             pos as u64 - 75,
                             fasta_length - 1,
+                            max_read_depth,
                         );
                         visualization =
                             manipulate_json(content, pos as u64 - 75, fasta_length - 1, max_rows);
@@ -314,6 +318,7 @@ pub(crate) fn make_table_report(
                             chrom.clone(),
                             pos as u64 - 75,
                             end_position as u64 + 75,
+                            max_read_depth,
                         );
                         visualization = manipulate_json(
                             content,
@@ -408,6 +413,7 @@ fn create_report_data(
     chrom: String,
     from: u64,
     to: u64,
+    max_read_depth: u32,
 ) -> (Json, usize) {
     let mut data = Vec::new();
 
@@ -416,7 +422,8 @@ fn create_report_data(
         data.push(nucleobase);
     }
 
-    let (bases, matches, max_rows) = get_static_reads(bam_path, fasta_path, chrom, from, to);
+    let (bases, matches, max_rows) =
+        get_static_reads(bam_path, fasta_path, chrom, from, to, max_read_depth);
 
     for b in bases {
         let base = json!(b);

src/bcf/report/table_report/mod.rs

@@ -10,6 +10,7 @@ use std::error::Error;
 use std::fs;
 use std::path::Path;
 
+#[allow(clippy::too_many_arguments)]
 pub fn table_report(
     vcf: &str,
     fasta: &str,
@@ -18,6 +19,7 @@ pub fn table_report(
     sample: &str,
     info: Option<Values>,
     format: Option<Values>,
+    max_read_depth: u32,
 ) -> Result<(), Box<dyn Error>> {
     let detail_path = output_path.to_owned() + "/details/" + sample;
     fs::create_dir(Path::new(&detail_path)).unwrap_or_else(|_| {
@@ -49,5 +51,6 @@ pub fn table_report(
         format_strings,
         sample.to_owned(),
         output_path,
+        max_read_depth,
     )?)
 }

src/bcf/report/table_report/static_reader.rs

@@ -2,8 +2,11 @@ use crate::bcf::report::table_report::alignment_reader::{
     make_nucleobases, read_indexed_bam, AlignmentMatch, AlignmentNucleobase,
 };
 use crate::bcf::report::table_report::create_report_table::VariantType;
+use rand::rngs::StdRng;
+use rand::seq::IteratorRandom;
+use rand_core::SeedableRng;
 use serde::Serialize;
-use std::collections::BTreeMap;
+use std::collections::{BTreeMap, HashSet};
 use std::path::Path;
 
 #[derive(Serialize, Clone, Debug)]
@@ -34,12 +37,13 @@ pub struct Variant {
 fn calc_rows(
     reads: Vec<AlignmentNucleobase>,
     matches: Vec<AlignmentMatch>,
+    max_read_depth: u32,
 ) -> (
     Vec<StaticAlignmentNucleobase>,
     Vec<StaticAlignmentMatch>,
     usize,
 ) {
-    let mut row_ends = vec![0; 1000];
+    let mut row_ends = vec![0; 10000];
 
     let mut read_names: BTreeMap<String, u16> = BTreeMap::new();
 
@@ -54,7 +58,7 @@ fn calc_rows(
         if read_names.contains_key(&r.name) {
             row = *read_names.get(&r.name).unwrap();
         } else {
-            for (i, _) in row_ends.iter().enumerate().take(1000).skip(1) {
+            for (i, _) in row_ends.iter().enumerate().take(10000).skip(1) {
                 if r.read_start > row_ends[i] {
                     if i > max_row {
                         max_row = i;
@@ -102,6 +106,23 @@ fn calc_rows(
         reads_wr.push(base);
     }
 
+    if max_row > max_read_depth as usize {
+        let mut rng = StdRng::seed_from_u64(42);
+        let random_rows: HashSet<_> = (0..max_row as u32)
+            .choose_multiple(&mut rng, max_read_depth as usize)
+            .into_iter()
+            .collect();
+        reads_wr = reads_wr
+            .into_iter()
+            .filter(|b| random_rows.contains(&&(b.row as u32)))
+            .collect();
+        matches_wr = matches_wr
+            .into_iter()
+            .filter(|b| random_rows.contains(&&(b.row as u32)))
+            .collect();
+        max_row = max_read_depth as usize;
+    }
+
     (reads_wr, matches_wr, max_row)
 }
 
@@ -111,12 +132,13 @@ pub fn get_static_reads(
     chrom: String,
     from: u64,
     to: u64,
+    max_read_depth: u32,
 ) -> (
     Vec<StaticAlignmentNucleobase>,
     Vec<StaticAlignmentMatch>,
     usize,
 ) {
     let alignments = read_indexed_bam(path, chrom.clone(), from, to);
     let (msm, m) = make_nucleobases(fasta_path, chrom, alignments, from, to);
-    calc_rows(msm, m)
+    calc_rows(msm, m, max_read_depth)
 }

src/cli.yaml

@@ -239,6 +239,11 @@ subcommands:
             short: c
             default_value: "1000"
             help: Set the maximum number of cells in the oncoprint per page. Lowering max-cells should improve the performance of the plots in the browser. Default value is 1000.
+        - max-read-depth:
+            required: false
+            short: d
+            default_value: "500"
+            help: Set the maximum read depth for plots of the third stage. If the maximum depth is exceeded, rows are randomly (with seed for reproducibility) chosen. Default value is 500.
         - info:
             multiple: true
             required: false

src/main.rs

@@ -86,6 +86,8 @@ fn main() -> Result<(), Box<dyn Error>> {
                 });
             }
             let max_cells = u32::from_str(matches.value_of("max-cells").unwrap()).unwrap();
+            let max_read_depth =
+                u32::from_str(matches.value_of("max-read-depth").unwrap()).unwrap();
             let custom_js = matches.value_of("custom-js");
             let tsv_data = matches.value_of("tsv");
             bcf::report::embed_js(output_path, custom_js)?;
@@ -123,6 +125,7 @@ fn main() -> Result<(), Box<dyn Error>> {
                     sample,
                     infos.clone(),
                     formats.clone(),
+                    max_read_depth,
                 )?;
             }