diff --git a/src/bin/commands/trimmer.rs b/src/bin/commands/trimmer.rs index 79c50e2..af56e38 100644 --- a/src/bin/commands/trimmer.rs +++ b/src/bin/commands/trimmer.rs @@ -56,7 +56,7 @@ pub(crate) struct TrimmerOpts { #[clap(long, default_value = "15")] osc_window: usize, - /// Required of oscillations in a window to trigger trimming/masking. + /// Required number of oscillations in a window to trigger trimming/masking. #[clap(long, default_value = "4")] osc_max_oscillations: usize, @@ -173,6 +173,13 @@ impl Command for TrimmerOpts { self.overlap_min_length, self.overlap_max_error_rate, ) { + log::debug!( + "found overlap in pair: {} shift: {}, overlap: {}, adapter: {}", + r1.id().unwrap_or("read"), + overlap.shift, + overlap.overlap, + overlap.adapter + ); stats.overlap_stats.update(overlap); let corrections = overlap.correct( &mut r1, @@ -184,6 +191,7 @@ impl Command for TrimmerOpts { Some(self.mask_quality) }, ); + log::debug!("corrections: {:?}", corrections); stats.overlap_stats.update_corrections(corrections.0, stats::ReadI::R1); stats.overlap_stats.update_corrections(corrections.1, stats::ReadI::R2); } diff --git a/src/lib/base_quality.rs b/src/lib/base_quality.rs index 56fa7cd..6a8c596 100644 --- a/src/lib/base_quality.rs +++ b/src/lib/base_quality.rs @@ -31,11 +31,12 @@ pub fn identify_trim_point( // here we have indices of oscillations. e.g.: [50, 52, 53, 54, 55, 57, 58, 59, 60, 62, 63, 64, 65, 67, 68, 69, 70, 72, 73, 74] // use a window of `max_oscilations` and check if there are at least `max_oscillations` oscillations in the window. - osc.windows(max_oscillations).find(|w| + osc.windows(max_oscillations) + .find(|w| // given a e.g. [50, 52, 55] means we found 3 oscillations in 5 bases (55 - 50) // and we had window_size to find that many. so if the last - first < window_size, we found an osc window. - w[w.len() - 1] - w[0] < window_size - ).map(|w| w[0]) + w[w.len() - 1] - w[0] < window_size) + .map(|w| w[0]) } // Indicates which tail(s) to clip. diff --git a/src/lib/pair_overlap.rs b/src/lib/pair_overlap.rs index 8b28315..b3af56b 100644 --- a/src/lib/pair_overlap.rs +++ b/src/lib/pair_overlap.rs @@ -68,10 +68,24 @@ impl PairOverlap { } else { // when bases disagree, set the lower quality base to the higher quality base. if r1_quals[i + self.shift] >= r2_quals[i] + bq_delta { + log::debug!( + "correcting read2 at {} with read1 at {}, base: {} to {}", + i, + i + self.shift, + r2_seq[i] as char, + r1_seq[i + self.shift] as char + ); r2_seq[i] = r1_seq[i + self.shift]; r2_quals[i] = r1_quals[i + self.shift]; correction_counts.1 += 1; } else if r2_quals[i] >= r1_quals[i + self.shift] + bq_delta { + log::debug!( + "correcting read1 at {} with read2 at {}, base: {} to {}", + i + self.shift, + i, + r1_seq[i + self.shift] as char, + r2_seq[i] as char + ); r1_seq[i + self.shift] = r2_seq[i]; r1_quals[i + self.shift] = r2_quals[i]; correction_counts.0 += 1;