Feat/introduce revcomp - #201
Conversation
There was a problem hiding this comment.
Code Review
This pull request introduces a new revcomp command to the bqtools CLI, enabling users to reverse complement sequences in BINSEQ files. It also adds a safety check to prevent writing binary BINSEQ data to stdout without the --pipe flag. The review feedback highlights three key areas for improvement: first, a bug in on_thread_complete where thread-local record counts are not added to the global count, leading to underreported totals; second, a performance optimization for the reverse_complement function to process buffers in a single pass; and third, a bug where single-end files are left untouched instead of being reverse complemented when an invalid mate flag is provided.
| fn reverse_complement(buf: &mut [u8]) { | ||
| buf.reverse(); | ||
| for base in buf.iter_mut() { | ||
| *base = match *base { | ||
| b'A' => b'T', | ||
| b'C' => b'G', | ||
| b'G' => b'C', | ||
| b'T' => b'A', | ||
| b'a' => b't', | ||
| b'c' => b'g', | ||
| b'g' => b'c', | ||
| b't' => b'a', | ||
| other => other, | ||
| }; | ||
| } | ||
| } |
There was a problem hiding this comment.
The current implementation of reverse_complement performs two passes over the buffer: one to reverse the bytes (buf.reverse()) and another to complement them. Since bqtools is designed for high-performance processing, we can optimize this to a single pass by swapping and complementing from both ends of the buffer simultaneously. This reduces memory bandwidth usage and cache misses.
fn reverse_complement(buf: &mut [u8]) {
let mut i = 0;
let mut j = buf.len().saturating_sub(1);
while i < j {
let left = complement(buf[i]);
let right = complement(buf[j]);
buf[i] = right;
buf[j] = left;
i += 1;
j -= 1;
}
if i == j && !buf.is_empty() {
buf[i] = complement(buf[i]);
}
}
#[inline]
fn complement(base: u8) -> u8 {
match base {
b'A' => b'T',
b'C' => b'G',
b'G' => b'C',
b'T' => b'A',
b'a' => b't',
b'c' => b'g',
b'g' => b'c',
b't' => b'a',
other => other,
}
}| pub fn run(args: &RevcompCommand) -> Result<()> { | ||
| let reader = BinseqReader::new(args.input.path())?; | ||
| if !reader.is_paired() && args.mate != Mate::Both { | ||
| warn!("Ignoring `--mate/-M` flag as only single channel found in file"); | ||
| } | ||
|
|
||
| let builder = get_builder(args)?; | ||
| let ohandle = args.output.as_writer()?; | ||
| let writer = builder.build(ohandle)?; | ||
| let mut processor = RevCompProcessor::new(writer, args.mate)?; |
There was a problem hiding this comment.
When a single-end file is processed and the user specifies a specific mate (e.g., -M 2), a warning is logged stating that the --mate/-M flag is ignored. However, the code still passes args.mate (which is Mate::Two) to RevCompProcessor. In process_record, rc_primary is computed as matches!(self.mate, Mate::One | Mate::Both). Since self.mate is Mate::Two, rc_primary evaluates to false, and the single channel is left untouched (not reverse complemented). This contradicts the warning and results in silent incorrect behavior.
To fix this, we should override the mate variable to Mate::Both when the file is single-end so that the single channel is actually reverse complemented as warned.
pub fn run(args: &RevcompCommand) -> Result<()> {
let reader = BinseqReader::new(args.input.path())?;
let mut mate = args.mate;
if !reader.is_paired() && mate != Mate::Both {
warn!("Ignoring `--mate/-M` flag as only single channel found in file");
mate = Mate::Both;
}
let builder = get_builder(args)?;
let ohandle = args.output.as_writer()?;
let writer = builder.build(ohandle)?;
let mut processor = RevCompProcessor::new(writer, mate)?;
No description provided.