Remove stderr gag if BAM reading fails so the error message is visible

src/hidive/src/fetch.rs

@@ -24,31 +24,20 @@ use skydive;
 ///
 /// Panics if any locus in `loci_list` cannot be parsed.
 pub fn start(output: &PathBuf, loci_list: &Vec<String>, bam_paths: &Vec<PathBuf>, require_spanning_reads: bool) {
-    // Initialize a HashSet to store unique loci after parsing
-    let mut loci = HashSet::new();
-
-    // Iterate over each locus in the provided list
-    for locus in loci_list {
-        // Attempt to parse the locus using a function from the skydive module
-        match skydive::utils::parse_locus(locus.to_owned()) {
-            Ok(l_fmt) => {
-                // If parsing is successful, insert the formatted locus into the HashSet
-                loci.insert(l_fmt);
-            }
-            Err(_) => {
-                // If parsing fails, panic and terminate the program, providing an error message
-                panic!("Could not parse locus '{}'.", locus);
-            }
-        }
-    }
+    let loci = skydive::utils::parse_loci(loci_list);
 
     // Convert the list of BAM file paths into a HashSet of URLs
     let reads_urls: HashSet<Url> = bam_paths
         .iter()
         // Use filter_map to attempt to parse each path as a URL, and collect the successful ones
-        .filter_map(
-            |path| Url::parse(&path.to_string_lossy()).ok()
-        )
+        .filter_map(|path| {
+            let path_str = path.to_string_lossy();
+            if path_str.starts_with("gs://") {
+                Url::parse(&path_str).ok()
+            } else {
+                Url::from_file_path(path.absolutize().unwrap()).ok()
+            }
+        })
         .collect();
 
     // Get the system's temporary directory path
@@ -61,7 +50,7 @@ pub fn start(output: &PathBuf, loci_list: &Vec<String>, bam_paths: &Vec<PathBuf>
     let r = skydive::stage::stage_data(&output_path, &loci, &reads_urls, &cache_path, require_spanning_reads);
 
     match r {
-        Ok(_) => { eprintln!("Done!") },
+        Ok(_) => {},
         Err(_) => { panic!("Failed to write multi-sample locus BAM.") },
     }
 }

src/hidive/src/trim.rs

@@ -45,6 +45,47 @@ pub fn start(output: &PathBuf, loci_list: &Vec<String>, bam_path: &PathBuf) {
         .iter()
         .for_each(|(chr, start, stop)| {
             let _ = bam.fetch(((*chr).as_bytes(), *start, *stop));
+
+            let mut bmap = HashMap::new();
+
+            for p in bam.pileup() {
+                let pileup = p.unwrap();
+
+                if *start <= (pileup.pos() as u64) && (pileup.pos() as u64) < *stop {
+                    for alignment in pileup.alignments() {
+                        let qname = String::from_utf8_lossy(alignment.record().qname()).into_owned();
+                        if !bmap.contains_key(&qname) {
+                            bmap.insert(qname.to_owned(), String::new());
+                        }
+
+                        if !alignment.is_del() && !alignment.is_refskip() {
+                            let a = alignment.record().seq()[alignment.qpos().unwrap()];
+
+                            bmap.get_mut(&qname).unwrap().push(a as char);
+                        }
+
+                        match alignment.indel() {
+                            bam::pileup::Indel::Ins(len) =>  {
+                                let pos1 = alignment.qpos().unwrap() as usize;
+                                let pos2 = pos1 + (len as usize);
+                                for pos in pos1..pos2 {
+                                    let a = alignment.record().seq()[pos];
+
+                                    bmap.get_mut(&qname).unwrap().push(a as char);
+                                }
+                            },
+                            _ => {}
+                        }
+                    }
+                }
+            }
+
+            for kv in bmap {
+                let output_str = format!(">{}\n{}", kv.0, kv.1);
+                writeln!(output_file, "{}", output_str).expect("Unable to write to file");
+            }
+
+            /*
             for (_, r) in bam.records().enumerate() {
                 let record = r.unwrap();
 
@@ -127,6 +168,7 @@ pub fn start(output: &PathBuf, loci_list: &Vec<String>, bam_path: &PathBuf) {
                     writeln!(output_file, "{}", output_str).expect("Unable to write to file");
                 }
             }
+            */
         });
 }
 

src/skydive/src/stage.rs

@@ -159,11 +159,18 @@ pub fn stage_data(
     // that the Jupyter notebook user doesn't get confused by intermediate
     // error messages that are nothing to worry about. The gag will end
     // automatically when it goes out of scope at the end of the function.
-    let _stderr_gag = Gag::stderr().unwrap();
+    let mut _stderr_gag = Gag::stderr().unwrap();
 
     // Get the header from the first BAM file.
     let first_url = reads_urls.iter().next().unwrap();
-    let first_header = get_bam_header(first_url, cache_path)?;
+    let first_header = match get_bam_header(first_url, cache_path) {
+        Ok(first_header) => first_header,
+        Err(e) => {
+            drop(_stderr_gag);
+
+            panic!("Error: {}", e);
+        }
+    };
 
     // Create a new header based on the first header.
     let mut header = bam::Header::from_template(&first_header);
@@ -175,13 +182,21 @@ pub fn stage_data(
         .collect();
 
     // Stage data from all BAM files.
-    let all_data= stage_data_from_all_files(reads_urls, loci, cache_path)?;
+    let all_data = match stage_data_from_all_files(reads_urls, loci, cache_path) {
+        Ok(all_data) => all_data,
+        Err(e) => {
+            drop(_stderr_gag);
+
+            panic!("Error: {}", e);
+        }
+    };
+
+    // Populate the output header with read group information, avoiding the read group that's already in the header.
     all_data
         .iter()
         .for_each(|(h, _)| {
             let mut hr = HeaderRecord::new("RG".as_bytes());
 
-            // Populate the header with read group information, avoiding the read group that's already in the header.
             h
                 .iter()
                 .filter(|a| {

src/skydive/src/utils.rs

@@ -1,5 +1,29 @@
+use std::collections::HashSet;
+
 use anyhow::Result;
 
+pub fn parse_loci(loci_list: &Vec<String>) -> HashSet<(String, u64, u64)> {
+    // Initialize a HashSet to store unique loci after parsing
+    let mut loci = HashSet::new();
+
+    // Iterate over each locus in the provided list
+    for locus in loci_list {
+        // Attempt to parse the locus using a function from the skydive module
+        match parse_locus(locus.to_owned()) {
+            Ok(l_fmt) => {
+                // If parsing is successful, insert the formatted locus into the HashSet
+                loci.insert(l_fmt);
+            }
+            Err(_) => {
+                // If parsing fails, panic and terminate the program, providing an error message
+                panic!("Could not parse locus '{}'.", locus);
+            }
+        }
+    }
+
+    loci
+}
+
 pub fn parse_locus(locus: String) -> Result<(String, u64, u64)> {
     let l_fmt = locus.replace(",", "");
     let parts: Vec<&str> = l_fmt.split(|c| (c == ':' || c == '-')).collect();