filtering options finalized

WardDeb · WardDeb · commit b6b58462e631 · 2025-04-03T17:08:47.000+02:00
diff --git a/docs/content/changelog.rst b/docs/content/changelog.rst
@@ -24,34 +24,6 @@ Changes in deepTools4.0
 - options label, smartLabels, genomeChunkLength are removed.
 - ignoreDuplicates is removed, and (if wanted) should be set by the SamFlagExclude setting.
 
-
-
-
-
-
-
-
-
-
-
-
-# Good to know
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
 # Testing
  
 ## computeMatrix
diff --git a/src/alignmentsieve.rs b/src/alignmentsieve.rs
@@ -75,6 +75,7 @@ pub fn r_alignmentsieve(
         Some(filterrnastrand.to_string()),
         None,
         None,
+        None,
     );
     let pool = ThreadPoolBuilder::new().num_threads(1).build().unwrap();
     let (sieve, filtersieve, totalreads, filteredreads) = pool.install(|| {
diff --git a/src/bamcompare.rs b/src/bamcompare.rs
@@ -28,7 +28,8 @@ pub fn r_bamcompare(
     operation: &str,
     pseudocount: f32,
     // filtering options
-    extendreads: u32, // if 0, no extension
+    extendreads: bool, // if 0, no extension
+    extendreadslen: u32, // length of extension (0 if PE or if not extending)
     centerreads: bool,
     blacklist: &str, // path to blacklist filename, or 'None'
     minmappingquality: u8, // 
@@ -74,8 +75,9 @@ pub fn r_bamcompare(
             _ => panic!("Error: Cannot determine filetype of blacklist file.")
         }
     }
+    // 
     // Set alignment filters
-    let filters = Alignmentfilters::new(
+    let mut filter1 = Alignmentfilters::new(
         blacklistregions,
         Some(minmappingquality),
         Some(samflaginclude),
@@ -86,18 +88,41 @@ pub fn r_bamcompare(
         None, // No offset
         None, // No strand filtering.
         Some(extendreads),
+        Some(extendreadslen),
         Some(centerreads),
     );
+    let mut filter2 = filter1.clone();
+    if filter1.extendreads && filter1.extendreadslen == 0 && !ispe1 {
+        panic!("Error: Extendreads is set, but not to a specific length (and library is single end). Specify the length with the --extendReads parameter.");
+    }
+    if filter2.extendreads && filter2.extendreadslen == 0 && !ispe2 {
+        panic!("Error: Extendreads is set, but not to a specific length (and library is single end). Specify the length with the --extendReads parameter.");
+    }
+    if filter1.extendreads && filter1.extendreadslen == 0 && ispe1 {
+        // We need a pass over the bamfile already to get the mean fragment length.
+        filter1.set_extendreadslen(bamifile1, nproc, &regions);
+        if verbose {
+            println!("fragment length for read extension set as: {} for bamfile 1", filter1.extendreadslen);
+        }
+    }
+    if filter2.extendreads && filter2.extendreadslen == 0 && ispe2 {
+        // We need a pass over the bamfile already to get the mean fragment length.
+        filter2.set_extendreadslen(bamifile2, nproc, &regions);
+        if verbose {
+            println!("fragment length for read extension set as: {} for bamfile 2", filter2.extendreadslen);
+        }
+    }
+
     let pool = ThreadPoolBuilder::new().num_threads(nproc).build().unwrap();
     
     // Set up the bam files in a Vec.
-    let bamfiles = vec![(bamifile1, ispe1), (bamifile2, ispe2)];
+    let bamfiles: Vec<(&str, bool, &Alignmentfilters)> = vec![(bamifile1, ispe1, &filter1), (bamifile2, ispe2, &filter2)];
 
     let mut covcalcs: Vec<ParsedBamFile> = pool.install(|| {
         bamfiles.par_iter()
-            .map(|(bamfile, ispe)| {
+            .map(|(bamfile, ispe, alfilter)| {
                 let (bg, mapped, unmapped, readlen, fraglen) = regionblocks.par_iter()
-                    .map(|i| bam_pileup(bamfile, &i, &binsize, &ispe, &ignorechr, &filters, false, false, true))
+                    .map(|i| bam_pileup(bamfile, &i, &binsize, &ispe, &ignorechr, alfilter , false, false, true))
                     .reduce(
                         || (vec![], 0, 0, vec![], vec![]),
                         |(mut _bg, mut _mapped, mut _unmapped, mut _readlen, mut _fraglen), (bg, mapped, unmapped, readlen, fraglen)| {
diff --git a/src/bamcoverage.rs b/src/bamcoverage.rs
@@ -27,8 +27,9 @@ pub fn r_bamcoverage(
     // processing options
     mnase: bool,
     _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,
+    extendreads: bool, // true for extension of reads
+    extendreadslen: u32, // if extendreads is set, and SE, this length is used for extension.
+    centerreads: bool, // to center the reads or not.
     filterrnastrand: &str, // forward, reverse or 'None'
     blacklist: &str, // path to blacklist filename, or 'None'
     _ignorechr: Py<PyList>, // list of chromosomes to ignore. Is empty if none.
@@ -69,14 +70,21 @@ pub fn r_bamcoverage(
     if norm != "None" && scalefactor != 1.0 {
         println!("Warning: You have set a normalization option ({}), but also a scale factor. Only the scale factor will be used", norm);
     }
-    if mnase && offset != (1, -1) {
+    if mnase && offset != (0, 0) {
         println!("Warning: Both MNase and offset are set. The offset will be ignored !");
     }
     if offset != (0, 0) {
         if offset.1 > 0 && offset.1 < offset.0 {
             panic!("Right side offset cannot be smaller than the left side offset.");
         }
     }
+    if extendreads && extendreadslen == 0 && !ispe {
+        panic!("Error: Extendreads is set, but not to a specific length (and library is single end). Specify the length with the --extendReads parameter.");
+    }
+    if centerreads && !extendreads {
+        println!("Warning: Centerreads is set, but extendreads is not. Centering will do nothing in this case.");
+    }
+
     if verbose {
         println!("Chromosomes to ignore for normalization: {:?}", ignorechr);
     }
@@ -120,7 +128,7 @@ pub fn r_bamcoverage(
         }
     }
     // Set alignment filters
-    let filters = Alignmentfilters::new(
+    let mut filters = Alignmentfilters::new(
         backlistregions,
         Some(minmappingquality),
         Some(samflaginclude),
@@ -131,8 +139,18 @@ pub fn r_bamcoverage(
         Some(offset),
         Some(filterrnastrand.to_string()),
         Some(extendreads),
+        Some(extendreadslen),
         Some(centerreads),
     );
+    // If extendreads, extendreadslen & ispe, we need a pass over the bamfile already now to get the mean fraglen.
+    if filters.extendreads && filters.extendreadslen == 0 && ispe {
+        // We need a pass over the bamfile already to get the mean fragment length.
+        filters.set_extendreadslen(bamifile, nproc, &regions);
+        if verbose {
+            println!("fragment length for read extension set as: {}", filters.extendreadslen);
+        }
+    }
+
     let pool = ThreadPoolBuilder::new().num_threads(nproc).build().unwrap();
     let (bg, mapped, _unmapped, readlen, fraglen) = pool.install(|| {
         regionblocks.par_iter()
@@ -145,10 +163,12 @@ pub fn r_bamcoverage(
                     _fraglen.extend(fraglen);
                     _mapped += mapped;
                     _unmapped += unmapped;
+
                     (_bg, _mapped, _unmapped, _readlen, _fraglen)
                 }
             )
     });
+
     let readlen = median(readlen);
     let fraglen = median(fraglen);
     if verbose {
@@ -166,7 +186,7 @@ pub fn r_bamcoverage(
         scalefactor,
         &verbose
     );
-    
+
     // Create output stream
     let lines = bg.into_iter().flat_map(
         |bg| {
@@ -187,30 +207,7 @@ pub fn r_bamcoverage(
     if verbose {
         println!("Writing output to: {}", ofile);
     }
-    
-    write_covfile(lines, ofile, ofiletype, chromsizes);
-    
-    // // Create output stream
-    // let lines: Vec<(String, Value)> = bg.into_par_iter().flat_map(
-    //     |bg| {
-    //         let reader = BufReader::new(File::open(bg).unwrap());
-    //         reader.lines().map(
-    //             |l| {
-    //                 let l = l.unwrap();
-    //                 let fields: Vec<&str> = l.split('\t').collect();
-    //                 let chrom: String = fields[0].to_string();
-    //                 let start: u32 = fields[1].parse().unwrap();
-    //                 let end: u32 = fields[2].parse().unwrap();
-    //                 let cov: f32 = fields[3].parse().unwrap();
-    //                 (chrom, Value {start: start, end: end, value: cov * sf})
-    //             }
-    //         ).collect::<Vec<_>>()
-    //     }
-    // ).collect();
-    // if verbose {
-    //     println!("Writing output to: {}", ofile);
-    // }
-    // write_covfile(lines.into_iter(), ofile, ofiletype, chromsizes);
 
+    write_covfile(lines, ofile, ofiletype, chromsizes);
     Ok(())
 }
diff --git a/src/covcalc.rs b/src/covcalc.rs
@@ -151,8 +151,6 @@ pub fn bam_pileup<'a>(
     let mut binsize = &bs;
     // constant to check if read is first in pair (relevant later)
     const FREAD: u16 = 0x40;
-    
-    
 
     // init variables for mapping statistics and lengths
     let mut mapped_reads: u32 = 0;
@@ -201,7 +199,7 @@ pub fn bam_pileup<'a>(
                     for record in bam.records() {
                         let mut record = record.expect("Error parsing record.");
                         if filters.filter {
-                            if filters.filter_record(&record) {
+                            if filters.filter_record(&record, &region.0.as_str()) {
                                 continue;
                             }
                         }
@@ -240,7 +238,7 @@ pub fn bam_pileup<'a>(
                         for record in bam.records() {
                             let mut record = record.expect("Error parsing record.");
                             if filters.filter {
-                                if filters.filter_record(&record) {
+                                if filters.filter_record(&record, region.0.as_str()) {
                                     continue;
                                 }
                             }
@@ -272,15 +270,17 @@ pub fn bam_pileup<'a>(
             counts = vec![0.0; (region.2 - region.1).div_ceil(*binsize) as usize];
             // let mut binstart = region.1;
             let mut binix: u32 = 0;
+
             for record in bam.records() {
                 let mut record = record.expect("Error parsing record.");
 
                 if filters.filter {
-                    if filters.filter_record(&record) {
+                    if filters.filter_record(&record, region.0.as_str()) {
                         continue;
                     }
                 }
                 if filters.manipulate {
+
                     let manipulated_blockpos = filters.manipulate_record(&mut record);
                     if manipulated_blockpos.is_none() {
                         continue;
@@ -380,19 +380,10 @@ pub fn bam_pileup<'a>(
     }
     let bgpath = bg.into_temp_path();
     let tmpvec   = vec![bgpath];
+    
     return (tmpvec, mapped_reads, unmapped_reads, readlens, fraglens);
 }
 
-fn pos_in_blacklist(pos: i64, chrom: &str, blacklist: &Vec<Region>) -> bool {
-    for region in blacklist.iter() {
-        // Note that get_startu / getendu is used as they are guaranteed to be u32's to start with.
-        if region.get_startu() <= pos as u32 && pos as u32 <= region.get_endu() && region.chrom == chrom {
-            return true;
-        }
-    }
-    return false;
-}
-
 #[derive(Clone, Debug)]
 pub struct Region {
     pub chrom: String,
@@ -1330,33 +1321,9 @@ pub fn region_divider(regs: &Vec<Region>) -> Vec<Vec<Region>> {
                 tempregionvec = Vec::new();
                 bplen = 0
             }
-            // // our regions are rather large, so we can split these up (in case both start/end are Revalue:U)
-            // match (&reg.start, &reg.end) {
-            //     (Revalue::U(start), Revalue::U(end)) => {
-            //         let mut start: u32 = *start;
-            //         let mut end: u32 = *end;
-            //         while start < end {
-            //             let newend = std::cmp::min(start + 10000000, end);
-            //             let mut entryname = format!("{}:{}-{}", reg.chrom, start, newend);
-            //             tempregionvec.push( Region {
-            //                 chrom: reg.chrom.clone(),
-            //                 start: Revalue::U(start),
-            //                 end: Revalue::U(newend),
-            //                 score: reg.score.clone(),
-            //                 strand: reg.strand.clone(),
-            //                 name: entryname,
-            //                 regionlength: newend-start
-            //             } );
-            //             start = newend;
-            //         }
-            //         blocks.push(tempregionvec);
-            //         tempregionvec = Vec::new();
-            //     },
-            //     _ => {
-            //         blocks.push(vec![reg.clone()]);
-            //     }
+
             blocks.push(vec![reg.clone()]);
-            // }
+
         } else {
             tempregionvec.push(reg.clone());
             bplen += reg.regionlength;
diff --git a/src/filtering.rs b/src/filtering.rs
diff --git a/src/multibamsummary.rs b/src/multibamsummary.rs
