fq trimming

src/bin/commands/trimmer.rs

@@ -109,6 +109,11 @@ fn create_writer<P: AsRef<Path>>(name: P) -> Result<BufWriter<File>, Error> {
     Ok(BufWriter::new(File::create(name)?))
 }
 
+fn check_extension(p: &Path) -> bool {
+    let ext = p.extension().map_or("", |v| v.to_str().unwrap_or(""));
+    ["bgz", "gz"].contains(&ext)
+}
+
 impl TrimmerOpts {
     fn prepare(&self) -> Result<(Pool, Vec<PooledWriter>, VecOfFqReaders), Error> {
         let fgio = Io::new(5, BUFFER_SIZE);
@@ -125,7 +130,7 @@ impl TrimmerOpts {
             .output
             .iter()
             .map(|p| {
-                if !(p.ends_with(".gz") || p.ends_with(".bgz")) {
+                if !check_extension(p) {
                     log::warn!(
                         "Output file {} does not end with .gz or .bgz. Writing bgzipped output to {} instead.",
                         p.display(),
@@ -260,3 +265,29 @@ impl Command for TrimmerOpts {
         Ok(())
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use std::path::PathBuf;
+
+    #[test]
+    fn test_check_extension_gz() {
+        let path = PathBuf::from("test_file.gz");
+        assert_eq!(
+            check_extension(&path),
+            true,
+            "The function should return true for .gz extension"
+        );
+    }
+
+    #[test]
+    fn test_check_extension_txt() {
+        let path = PathBuf::from("test_file.txt");
+        assert_eq!(
+            check_extension(&path),
+            false,
+            "The function should return false for .txt extension"
+        );
+    }
+}

src/lib/pair_overlap.rs

@@ -51,6 +51,41 @@ impl PairOverlap {
         let mut r2_quals = r2.qual.clone();
         r2_quals.reverse();
 
+        if log::log_enabled!(log::Level::Debug) {
+            let r1_seq = String::from_utf8(r1_seq.clone()).unwrap();
+            let r2_seq = String::from_utf8(r2_seq.clone()).unwrap();
+            let r1_space = if self.adapter { " ".repeat(self.shift) } else { "".to_string() };
+            let r2_space = if self.adapter { "".to_string() } else { " ".repeat(self.shift) };
+
+            let r1_string = format!("{}{}", r1_space, r1_seq);
+            let r2_string = format!("{}{}", r2_space, r2_seq);
+
+            let chunk_size = 80;
+            // now wrap strings at chunk_size chars and print each pair of lines.
+            let r1_chunks = r1_string
+                .as_bytes()
+                .chunks(chunk_size)
+                .map(|x| String::from_utf8(x.to_vec()).unwrap())
+                .collect::<Vec<String>>();
+
+            let r2_chunks = r2_string
+                .as_bytes()
+                .chunks(chunk_size)
+                .map(|x| String::from_utf8(x.to_vec()).unwrap())
+                .collect::<Vec<String>>();
+
+            for i in 0..r2_chunks.len() {
+                let diffs = r1_chunks[i]
+                    .chars()
+                    .zip(r2_chunks[i].chars())
+                    .map(|(a, b)| if a != b && a != ' ' && b != ' ' { '*' } else { ' ' })
+                    .collect::<Vec<char>>();
+                eprintln!("r1: {}", r1_chunks[i]);
+                eprintln!("r2: {}", r2_chunks[i]);
+                eprintln!("dx: {}", diffs.iter().collect::<String>());
+            }
+        }
+
         if self.adapter {
             std::mem::swap(&mut r1_seq, &mut r2_seq);
             std::mem::swap(&mut r1_quals, &mut r2_quals);
@@ -97,19 +132,30 @@ impl PairOverlap {
             std::mem::swap(&mut r1_seq, &mut r2_seq);
             std::mem::swap(&mut r1_quals, &mut r2_quals);
             correction_counts = (correction_counts.1, correction_counts.0);
+            // clip right end of r1
+            // but we must check for scenario like:
+            // r1:      ACGAAATA
+            // r2: GGGGGACGAAATA
+            // shift: 5
+            // where we don't need to clip r1.
+            let r1_overhang =
+                if self.overlap >= r1_seq.len() { 0 } else { r1_seq.len() - self.overlap };
             if let Some(maskq) = mask_quality {
                 // mask qualities
                 let l = r1_quals.len();
-                for q in &mut r1_quals[l - self.shift..] {
-                    *q = maskq;
+                if r1_overhang > 0 {
+                    for q in &mut r1_quals[l - r1_overhang as usize..] {
+                        *q = maskq;
+                    }
                 }
                 for q in &mut r2_quals[..self.shift] {
                     *q = maskq;
                 }
             } else {
-                // clip right end of r1
-                r1_seq = r1_seq[..r1_seq.len() - self.shift].to_vec();
-                r1_quals = r1_quals[..r1_quals.len() - self.shift].to_vec();
+                if r1_overhang > 0 {
+                    r1_seq = r1_seq[..r1_seq.len() - r1_overhang as usize].to_vec();
+                    r1_quals = r1_quals[..r1_quals.len() - r1_overhang as usize].to_vec();
+                }
                 // we're still reversed here, so we also clip of left end of r2
                 r2_seq = r2_seq[self.shift..].to_vec();
                 r2_quals = r2_quals[self.shift..].to_vec();
@@ -138,11 +184,11 @@ pub fn find_overlap(
     min_overlap: usize,
     max_overlap_error_rate: f64,
 ) -> Option<PairOverlap> {
-    // note that it's possible to do this without allocating, but then we must (re)complement many times.
     if s1.len() < min_overlap || s2.len() < min_overlap {
         return None;
     }
 
+    // note that it's possible to do this without allocating, but then we must (re)complement many times.
     let s2 = s2.iter().rev().map(|x| complement(*x)).collect::<Vec<u8>>();
 
     'shift_loop: for shift in 0..s1.len().max(min_overlap) - min_overlap {
@@ -186,6 +232,42 @@ pub fn find_overlap(
 mod tests {
     use super::*;
 
+    #[test]
+    fn test_with_adapter() {
+        let r1_seq = /* 1234567890 */ b"GCGCGTGCGTTCTGC";
+        let r2_seq = b"AAAAAAAAAAAACACGTGCGCGTGCGTTCTGC"
+            .iter()
+            .rev()
+            .map(|x| complement(*x))
+            .collect::<Vec<u8>>();
+
+        let overlap = find_overlap(r1_seq, &r2_seq, 4, 0.1);
+        assert!(overlap.is_some());
+        let overlap = overlap.unwrap();
+        assert!(overlap.shift == 17);
+        assert!(overlap.adapter);
+        assert!(overlap.overlap == r1_seq.len());
+        eprintln!("overlap: {:?}", overlap);
+
+        let mut r1 = OwnedRecord {
+            head: b"r1".to_vec(),
+            seq: r1_seq.to_vec(),
+            qual: vec![b'I'; r1_seq.len()],
+        };
+        eprintln!("r1-quals: {:?}", r1.qual);
+        let mut r2 = OwnedRecord {
+            head: b"r2".to_vec(),
+            seq: r2_seq.to_vec(),
+            qual: vec![b'I'; r2_seq.len()],
+        };
+
+        let _corrections = overlap.correct(&mut r1, &mut r2, 1, None);
+        eprintln!("r1-seq-after: {}", String::from_utf8(r1.seq.clone()).unwrap());
+        eprintln!("r2-seq-after: {}", String::from_utf8(r2.seq.clone()).unwrap());
+        assert!(r1.seq == r1_seq);
+        assert!(r2.seq == r2_seq[..r2_seq.len() - overlap.shift].to_vec());
+    }
+
     #[test]
     fn test_simple_overlap() {
         // shift: 0
@@ -253,8 +335,8 @@ mod tests {
     #[test]
     fn test_overlap_correction_with_adapter() {
         let overlap = PairOverlap { shift: 3, overlap: 5, diffs: 0, adapter: true };
-        //       IIIQIIII
-        // r1:   AAAAAGGC
+        //        IIIQIIII
+        // r1:    AAAAAGGC
         // r2: ACGAAATA
 
         let mut r1 = OwnedRecord {