g

src/lib.rs

@@ -155,7 +155,7 @@ pub fn get_start_end_qual(
     stop: i32,
     mapq: u8,
 ) -> Result<Vec<i32>> {
-    bam.fetch((chr, start, stop))?;
+    bam.fetch((chr, start - 1, stop))?;
     let length = stop - start;
     let mut results = Vec::new();
     results.resize(length as usize, 0);
@@ -163,9 +163,11 @@ pub fn get_start_end_qual(
     for read in bam.records() {
         let record = read.context(format!("Error while parsing record"))?;
         let rstart = record.pos() as i32;
+        // let rstart = record.reference_start() as i32;
         let rend = record.reference_end() as i32;
 
         if rstart >= start && rstart < stop {
+            println!("start {rstart}");
             if record.mapq() >= mapq {
                 let index = rstart - start;
                 *results.get_mut(index as usize).unwrap() += 1;
@@ -173,6 +175,7 @@ pub fn get_start_end_qual(
         }
 
         if rend >= start && rend < stop {
+            println!("{rend}");
             if record.mapq() >= mapq {
                 let index = rend - start;
                 *results.get_mut(index as usize).unwrap() += 1;
@@ -190,7 +193,7 @@ pub fn get_start_end_qual_rec(
     stop: i32,
     mapq: u8,
 ) -> Result<(Vec<Vec<Record>>, Vec<Vec<Record>>)> {
-    bam.fetch((chr, start, stop))?;
+    bam.fetch((chr, start - 1, stop))?;
     let length = stop - start;
     let mut results_start: Vec<Vec<Record>> = Vec::new();
     results_start.resize(length as usize, vec![]);
@@ -205,9 +208,7 @@ pub fn get_start_end_qual_rec(
         if rstart >= start && rstart < stop {
             if record.mapq() >= mapq {
                 let index = rstart - start;
-                let u = results_start
-                    .get_mut(index as usize)
-                    .unwrap();
+                let u = results_start.get_mut(index as usize).unwrap();
                 u.push(record.clone());
             }
         }
@@ -215,9 +216,7 @@ pub fn get_start_end_qual_rec(
         if rend >= start && rend < stop {
             if record.mapq() >= mapq {
                 let index = rend - start;
-                let u = results_end
-                    .get_mut(index as usize)
-                    .unwrap();
+                let u = results_end.get_mut(index as usize).unwrap();
                 u.push(record.clone());
             }
         }
@@ -314,3 +313,27 @@ pub fn range_len_qual(
     }
     Ok(depths)
 }
+
+#[cfg(test)]
+mod tests {
+    // Note this useful idiom: importing names from outer (for mod tests) scope.
+    use super::*;
+
+    #[test]
+    fn test() {
+        let bam_path = "/data/longreads_basic_pipe/CAMARA/diag/CAMARA_diag_hs1.bam";
+        let chr = "chr9";
+        let start = 21839796;
+        let mapq = 50;
+
+
+        let mut bam = rust_htslib::bam::IndexedReader::from_path(bam_path)
+            .unwrap();
+
+        let a = get_start_end_qual(&mut bam, chr, start, start + 1, mapq).unwrap();
+        println!("{a:?}");
+
+        let res = get_start_end_qual_rec(&mut bam, chr, start, start + 1, mapq).unwrap();
+        println!("{res:?}");
+    }
+}