Skip to content

Commit

Permalink
logging in pair overlap
Browse files Browse the repository at this point in the history
  • Loading branch information
brentp committed Aug 24, 2023
1 parent 6201ea1 commit 7a2e646
Show file tree
Hide file tree
Showing 3 changed files with 27 additions and 4 deletions.
10 changes: 9 additions & 1 deletion src/bin/commands/trimmer.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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,

Expand Down Expand Up @@ -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,
Expand All @@ -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);
}
Expand Down
7 changes: 4 additions & 3 deletions src/lib/base_quality.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
14 changes: 14 additions & 0 deletions src/lib/pair_overlap.rs
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down

0 comments on commit 7a2e646

Please sign in to comment.