g

src/lib.rs

@@ -1,5 +1,8 @@
 use anyhow::{Context, Ok, Result};
-use rust_htslib::{bam, bam::Read};
+use rust_htslib::{
+    bam,
+    bam::{ext::BamRecordExtensions, record, Read},
+};
 
 use average::Mean;
 pub fn get_hts_nt_pileup(
@@ -130,13 +133,16 @@ pub fn range_depths(
 
     let mut depths = vec![0];
     depths.resize((stop - start) as usize, 0);
-
     for p in bam.pileup() {
         let pileup = p.context(format!("eRR"))?;
         let rstart = pileup.pos() as i32;
+
         if rstart >= start && rstart < stop {
             // depths.push(pileup.depth());
-            let v = depths.get_mut((rstart - start) as usize).context(format!("Errrr {}", rstart - start))?;
+
+            let v = depths
+                .get_mut((rstart - start) as usize)
+                .context(format!("Errrr {}", rstart - start))?;
             *v = pileup.depth();
         } else if rstart >= stop {
             break;
@@ -145,3 +151,33 @@ pub fn range_depths(
 
     Ok(depths)
 }
+
+pub fn range_depths_qual(
+    bam: &mut rust_htslib::bam::IndexedReader,
+    chr: &str,
+    start: i32,
+    stop: i32,
+    mapq: u8,
+) -> Result<Vec<usize>> {
+    bam.fetch((chr, start, stop))?;
+
+    let mut depths = vec![0];
+    depths.resize((stop - start) as usize, 0);
+    for p in bam.pileup() {
+        let pileup = p.context(format!("eRR"))?;
+        let rstart = pileup.pos() as i32;
+
+        if rstart >= start && rstart < stop {
+            let v = depths
+                .get_mut((rstart - start) as usize)
+                .context(format!("Errrr {}", rstart - start))?;
+            *v = pileup
+                .alignments()
+                .filter(|e| e.record().mapq() >= mapq)
+                .count();
+        } else if rstart >= stop {
+            break;
+        }
+    }
+    Ok(depths)
+}