proper offset implementation, non offsets encoded as (0,0)

Author: WardDeb

pydeeptools/deeptools/bamCoverage2.py

@@ -141,8 +141,11 @@ def main(args=None):
         args.effectiveGenomeSize = 0
     if not args.normalizeUsing:
         args.normalizeUsing = 'None'
+    print(args.Offset)
     if not args.Offset:
-        args.Offset = [1, -1]
+        args.Offset = [0, 0]
+    elif len(args.Offset) == 1:
+        args.Offset = [args.Offset[0], 0]
     if not args.extendReads:
         args.extendReads = 0
     if not args.filterRNAstrand:

src/bamcoverage.rs

@@ -2,6 +2,7 @@ use pyo3::prelude::*;
 use pyo3::types::PyList;
 use rayon::prelude::*;
 use rayon::ThreadPoolBuilder;
+use core::panic;
 use std::io::prelude::*;
 use std::io::BufReader;
 use std::fs::File;
@@ -25,7 +26,7 @@ pub fn r_bamcoverage(
     scalefactor: f32, // default 1.0
     // processing options
     mnase: bool,
-    _offset: Py<PyList>, // list of max 2 [offset 5', offset 3'], if no offset is required we have [1, -1]
+    _offset: Py<PyList>, // list of 2 [offset 5', offset 3'], if no offset is required we have [0, 0]
     extendreads: u32, // if 0, no extension
     centerreads: bool,
     filterrnastrand: &str, // forward, reverse or 'None'
@@ -45,14 +46,14 @@ pub fn r_bamcoverage(
     verbose: bool,
     collapse: bool,
 ) -> PyResult<()> {
-    let mut offset: (i32, i32) = (1, -1);
+    let mut offset: (i32, i32) = (0, 0);
     let mut ignorechr: Vec<String> = Vec::new();
     Python::with_gil(|py| {
         let offsetv: Vec<i32> = _offset.extract(py).expect("Failed to retrieve offset.");
-        if offsetv.len() == 1 {
-            offset = (offsetv[0], -1);
-        } else if offsetv.len() == 2 {
+        if offsetv.len() == 2 {
             offset = (offsetv[0], offsetv[1]);
+        } else {
+            panic!("Error: Offset should be a list of 2. Received: {:?}", offsetv);
         }
         ignorechr = _ignorechr.extract(py).expect("Failed to retrieve ignorechr.");
     });
@@ -71,10 +72,7 @@ pub fn r_bamcoverage(
     if mnase && offset != (1, -1) {
         println!("Warning: Both MNase and offset are set. The offset will be ignored !");
     }
-    if offset != (1, -1) {
-        if offset.0 == 0 {
-            panic!("Offsets cannot be 0 as they are 1-based positions inside each alignment.");
-        }
+    if offset != (0, 0) {
         if offset.1 > 0 && offset.1 < offset.0 {
             panic!("Right side offset cannot be smaller than the left side offset.");
         }

src/covcalc.rs

@@ -281,15 +281,16 @@ pub fn bam_pileup<'a>(
                     }
                 }
                 if filters.manipulate {
-                    let manipulated_blocks = filters.manipulate_record(&mut record);
-                    if manipulated_blocks.is_none() {
+                    let manipulated_blockpos = filters.manipulate_record(&mut record);
+                    if manipulated_blockpos.is_none() {
                         continue;
                     }
-                    let block = manipulated_blocks.unwrap();
-                    let ixstart = ((block[0] - region.1) / binsize) as usize;
-                    let ixend = ((block[1] - region.1) / binsize) as usize;
+                    let indices: HashSet<usize> = manipulated_blockpos
+                        .unwrap()
+                        .into_iter()
+                        .map(|x| ((x - region.1) / binsize) as usize)
+                        .collect();
 
-                    let indices: HashSet<usize> = (ixstart..ixend).collect();
                     indices.into_iter()
                         .for_each(|ix| {
                             if ix < counts.len() {

src/filtering.rs

@@ -40,13 +40,13 @@ impl Alignmentfilters {
         let _mifl = minfraglen.unwrap_or(0);
         let _mafl = maxfraglen.unwrap_or(0);
         let _mnase = mnase.unwrap_or(false);
-        let _offset =  offset.unwrap_or((1, -1));
+        let _offset =  offset.unwrap_or((0, 0));
         let _frs = filterrnastrand.unwrap_or(String::from("None"));
         let _extend = extendreads.unwrap_or(0);
         let _center = centerreads.unwrap_or(false);
 
         // Set the manipulate bool for a quick escape in case manipulation is not needed.
-        if _offset != (1, -1) || _mnase || _extend > 0 || _center {
+        if _offset != (0, 0) || _mnase || _extend > 0 || _center {
             manipulate = true;
         }
 
@@ -144,7 +144,7 @@ impl Alignmentfilters {
         }
     }
 
-    pub fn manipulate_record(&self, rec: &Record) -> Option<[u32; 2]> {
+    pub fn manipulate_record(&self, rec: &Record) -> Option<Vec<u32>> {
         // Not just simple yes - no filtering, but actually manipulating the record.
         // In general, this is the case for MNase mode, offset, extendreads and centerreads.
         if self.mnase {
@@ -154,26 +154,65 @@ impl Alignmentfilters {
             let rinsertsize = rec.insert_size().abs() as u32;
             if rec.is_proper_pair() && !rec.is_reverse() && rinsertsize > 1 {
                 // Not sure why the insert size test is here, but note we always filter prior to manipulating.
-                // insertsize even
                 let recpos: u32 = rec.pos() as u32;
                 let frag_start = recpos - 1 + rinsertsize / 2;
 
                 if rinsertsize % 2 == 0 {
-                    // 
                     return Some(
-                        [frag_start, frag_start + 2]
+                        (frag_start..frag_start + 2).collect()
                     );
                 } else {
-                    // Uneven middle, frag start + 4 to keep output consistent with deeptools 3.
                     return Some(
-                        [frag_start, frag_start + 4]
+                        (frag_start..frag_start+4).collect()
                     );
                 }
             }
             return None;
         }
-        if self.offset != (1, -1) {
-            println!("Offset implementation");
+        if self.offset != (0, 0) {
+            // Collect blocks and flatten them out.
+            let mut blockvec: Vec<u32> = rec
+                .aligned_blocks()
+                .flat_map(|x| x[0] as u32..x[1] as u32)
+                .collect();
+            let blocklen = blockvec.len() as i32;
+
+            // Convert potential negative indices to positive indices
+            // It could be that for the offset only one value is given, in which case we only use that site
+            if self.offset.1 == 0 {
+                // 
+                let pos = if self.offset.0 < 0 {blocklen + self.offset.0 } else {self.offset.0 - 1};
+                if pos < 0 || pos >= blocklen {
+                    return None;
+                }
+
+                if rec.is_reverse() {
+                    blockvec.reverse();
+                    blockvec = blockvec[pos as usize..pos as usize + 1].to_vec();
+                    blockvec.reverse();
+                    return Some(blockvec);
+                } else {
+                    blockvec = blockvec[pos as usize..pos as usize + 1].to_vec();
+                    return Some(blockvec);
+                }
+            } else {
+                let start = if self.offset.0 < 0 { blocklen + self.offset.0 } else { self.offset.0 -1};
+                let end = if self.offset.1 < 0 { blocklen + self.offset.1 + 1 } else { self.offset.1 };
+
+                // if the range falls outside the vec, return none (retain deeptools 3 behavior)
+                if start < 0 || end < 0 || start >= blocklen || end >= blocklen || start >= end {
+                    return None;
+                }
+                if rec.is_reverse() {
+                    blockvec.reverse();
+                    blockvec = blockvec[start as usize..end as usize].to_vec();
+                    blockvec.reverse();
+                    return Some(blockvec);
+                } else {
+                    blockvec = blockvec[start as usize..end as usize].to_vec();
+                    return Some(blockvec);
+                }
+            }
         }
         if self.extendreads > 0 {
             println!("extendreads implementation");