diff --git a/README.md b/README.md index b4a4e80..8c4b98c 100644 --- a/README.md +++ b/README.md @@ -45,6 +45,7 @@ It is the _fastest_ variant but _is lossy_ by design. - **Grep**: Search for fixed-string, regex, or fuzzy matches in BINSEQ files. - **Split**: Split a BINSEQ file into multiple files based on matching patterns. - **Pipe**: Create named-pipes for efficient data processing with legacy tools that don't support BINSEQ, optionally spawning and supervising the consumer commands directly (`-x`/`-X`). +- **Revcomp**: Reverse complement the sequences in a BINSEQ file. ## Installation @@ -118,6 +119,7 @@ bqtools grep --help bqtools split --help bqtools pipe --help bqtools qc --help +bqtools revcomp --help ``` ### Encoding @@ -257,6 +259,28 @@ Combine multiple BINSEQ files: bqtools cat file1.bq file2.bq file3.bq -o combined.bq ``` +> Note: `cat`, `revcomp`, and other commands that write BINSEQ output require either `-o/--output` +> or an explicit `--pipe` flag; binary BINSEQ data is never written to stdout implicitly. + +### Reverse Complementing + +Reverse complement the sequences in a BINSEQ file, preserving its format and configuration: + +```bash +bqtools revcomp input.cbq -o output.cbq +``` + +For paired files, both mates are reverse complemented by default. Use `-M/--mate` to +reverse complement only one of the two mates (the other is left untouched): + +```bash +# Only reverse complement mate 1 +bqtools revcomp input.cbq -o output.cbq -M 1 + +# Only reverse complement mate 2 +bqtools revcomp input.cbq -o output.cbq -M 2 +``` + ### Information and Statistics Show information and statistics about a BINSEQ file. diff --git a/src/cli/cli.rs b/src/cli/cli.rs index ebd2154..fda0d5e 100644 --- a/src/cli/cli.rs +++ b/src/cli/cli.rs @@ -8,7 +8,7 @@ use clap::{ use super::{ CatCommand, DecodeCommand, EncodeCommand, GrepCommand, InfoCommand, PipeCommand, QcCommand, - SampleCommand, SplitCommand, + RevcompCommand, SampleCommand, SplitCommand, }; // Configures Clap v3-style help menu colors @@ -45,4 +45,6 @@ pub enum Commands { Pipe(PipeCommand), Qc(QcCommand), + + Revcomp(RevcompCommand), } diff --git a/src/cli/mod.rs b/src/cli/mod.rs index c241bfa..ebd3476 100644 --- a/src/cli/mod.rs +++ b/src/cli/mod.rs @@ -9,6 +9,7 @@ mod input; mod output; mod pipe; mod qc; +mod revcomp; mod sample; mod split; @@ -25,5 +26,6 @@ pub use input::{InputBinseq, InputFile, MultiInputBinseq}; pub use output::{BinseqConfig, BinseqMode, Mate, OutputBinseq, OutputFile}; pub use pipe::PipeCommand; pub use qc::{QcCommand, QcOptions}; +pub use revcomp::RevcompCommand; pub use sample::SampleCommand; pub use split::SplitCommand; diff --git a/src/cli/output.rs b/src/cli/output.rs index dd5e98c..c967ffd 100644 --- a/src/cli/output.rs +++ b/src/cli/output.rs @@ -147,7 +147,7 @@ pub struct OutputBinseq { #[clap(short = 'o', long)] /// Output binseq file /// - /// To output to stdout, use the `-P/--pipe` flag. + /// To output to stdout, use the `--pipe` flag. pub output: Option, #[clap(flatten)] @@ -159,6 +159,11 @@ pub struct OutputBinseq { } impl OutputBinseq { pub fn as_writer(&self) -> Result> { + if self.output.is_none() && !self.pipe { + bail!( + "Refusing to write binary BINSEQ data to stdout. Provide an output path with `-o/--output`, or pass `--pipe` to write to stdout explicitly." + ); + } let writer = match_output(self.output.as_deref())?; Ok(writer) } @@ -440,3 +445,32 @@ impl From for BinseqConfig { } } } + +#[cfg(test)] +mod tests { + use clap::Parser; + + use super::OutputBinseq; + + /// Without `-o` or `--pipe`, writing binary BINSEQ data to stdout must be + /// refused rather than silently dumping binary into the terminal. + #[test] + fn test_as_writer_rejects_bare_stdout() { + let args = OutputBinseq::try_parse_from(["output"]).unwrap(); + assert!(args.as_writer().is_err()); + } + + #[test] + fn test_as_writer_allows_explicit_pipe() { + let args = OutputBinseq::try_parse_from(["output", "--pipe"]).unwrap(); + assert!(args.as_writer().is_ok()); + } + + #[test] + fn test_as_writer_allows_output_path() { + let tmp = tempfile::NamedTempFile::new().unwrap(); + let args = + OutputBinseq::try_parse_from(["output", "-o", tmp.path().to_str().unwrap()]).unwrap(); + assert!(args.as_writer().is_ok()); + } +} diff --git a/src/cli/revcomp.rs b/src/cli/revcomp.rs new file mode 100644 index 0000000..61196fb --- /dev/null +++ b/src/cli/revcomp.rs @@ -0,0 +1,24 @@ +use clap::Parser; + +use super::{InputBinseq, Mate, OutputBinseq}; + +/// Reverse complement the sequences in a BINSEQ file. +#[derive(Parser, Debug)] +pub struct RevcompCommand { + #[clap(flatten)] + pub input: InputBinseq, + + #[clap(flatten)] + pub output: OutputBinseq, + + /// Which mate(s) to reverse complement + /// + /// Only relevant for paired BINSEQ files. Defaults to reverse + /// complementing both mates; ignored (with a warning) on single-end + /// files. + /// + /// Note: `-m` is already used by `--mode` (BINSEQ output format), so + /// this flag uses `-M` instead. + #[clap(short = 'M', long, default_value = "both")] + pub mate: Mate, +} diff --git a/src/commands/mod.rs b/src/commands/mod.rs index cd91d02..6572a06 100644 --- a/src/commands/mod.rs +++ b/src/commands/mod.rs @@ -5,6 +5,7 @@ pub mod grep; pub mod info; pub mod pipe; pub mod qc; +pub mod revcomp; pub mod sample; pub mod split; mod utils; diff --git a/src/commands/revcomp/mod.rs b/src/commands/revcomp/mod.rs new file mode 100644 index 0000000..7af05bc --- /dev/null +++ b/src/commands/revcomp/mod.rs @@ -0,0 +1,227 @@ +mod processor; + +use anyhow::Result; +use binseq::{bq, cbq, vbq, BinseqReader, BinseqWriterBuilder, ParallelReader}; +use log::{info, warn}; + +use crate::cli::{BinseqMode, Mate, RevcompCommand}; +use processor::RevCompProcessor; + +/// Builds a writer that mirrors the input file's own header/configuration, +/// since reverse complementing changes sequence content but not schema. +fn get_builder(args: &RevcompCommand) -> Result { + let builder = match args.input.mode()? { + BinseqMode::Bq => { + let reader = bq::MmapReader::new(args.input.path())?; + BinseqWriterBuilder::from_bq_header(reader.header()) + } + BinseqMode::Vbq => { + let reader = vbq::MmapReader::new(args.input.path())?; + BinseqWriterBuilder::from_vbq_header(reader.header()) + } + BinseqMode::Cbq => { + let reader = cbq::MmapReader::new(args.input.path())?; + BinseqWriterBuilder::from_cbq_header(reader.header()) + } + }; + Ok(builder) +} + +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)?; + + if let Some(mut span) = args.input.span { + let num_records = reader.num_records()?; + reader.process_parallel_range( + processor.clone(), + args.output.threads(), + span.get_range(num_records)?, + )?; + } else { + reader.process_parallel(processor.clone(), args.output.threads())?; + } + processor.finish()?; + + info!( + "Wrote {} reverse complemented records", + processor.get_global_record_count() + ); + + Ok(()) +} + +#[cfg(test)] +mod tests { + use anyhow::Result; + use clap::Parser; + use itertools::iproduct; + use tempfile::NamedTempFile; + + use crate::cli::BinseqMode; + use crate::testutils::{count_binseq, write_fastx, DEFAULT_NUM_RECORDS}; + + fn encode(in_path: &std::path::Path, out_path: &std::path::Path) -> Result<()> { + let cmd = crate::cli::EncodeCommand::try_parse_from([ + "encode", + in_path.to_str().unwrap(), + "-o", + out_path.to_str().unwrap(), + ])?; + crate::commands::encode::run(&cmd) + } + + fn decode_to_fasta(bq_path: &std::path::Path, out_path: &std::path::Path) -> Result<()> { + let cmd = crate::cli::DecodeCommand::try_parse_from([ + "decode", + bq_path.to_str().unwrap(), + "-o", + out_path.to_str().unwrap(), + "-f", + "a", + ])?; + crate::commands::decode::run(&cmd) + } + + fn revcomp( + in_path: &std::path::Path, + out_path: &std::path::Path, + extra: &[&str], + ) -> Result<()> { + let mut cmd_args = vec![ + "revcomp".to_string(), + in_path.to_str().unwrap().to_string(), + "-o".to_string(), + out_path.to_str().unwrap().to_string(), + ]; + cmd_args.extend(extra.iter().map(std::string::ToString::to_string)); + let cmd = crate::cli::RevcompCommand::try_parse_from(cmd_args)?; + super::run(&cmd) + } + + fn reverse_complement_str(seq: &str) -> String { + seq.chars() + .rev() + .map(|c| match c { + 'A' => 'T', + 'C' => 'G', + 'G' => 'C', + 'T' => 'A', + other => other, + }) + .collect() + } + + /// Extracts just the sequence lines from a FASTA file, sorted, so + /// comparisons are insensitive to reordering from parallel processing. + fn sorted_sequences(path: &std::path::Path) -> Result> { + let content = std::fs::read_to_string(path)?; + let mut seqs: Vec = content + .lines() + .filter(|l| !l.starts_with('>')) + .map(std::string::ToString::to_string) + .collect(); + seqs.sort_unstable(); + Ok(seqs) + } + + /// Round-tripping revcomp twice must recover the original sequences. + #[test] + fn test_revcomp_double_application_is_identity() -> Result<()> { + for mode in BinseqMode::enum_iter() { + let in_tmp = write_fastx().call()?; + let bq_tmp = NamedTempFile::with_suffix(mode.extension())?; + encode(in_tmp.path(), bq_tmp.path())?; + + let rc_once = NamedTempFile::with_suffix(mode.extension())?; + revcomp(bq_tmp.path(), rc_once.path(), &[])?; + + let rc_twice = NamedTempFile::with_suffix(mode.extension())?; + revcomp(rc_once.path(), rc_twice.path(), &[])?; + + let original_fa = NamedTempFile::with_suffix(".fasta")?; + decode_to_fasta(bq_tmp.path(), original_fa.path())?; + let roundtrip_fa = NamedTempFile::with_suffix(".fasta")?; + decode_to_fasta(rc_twice.path(), roundtrip_fa.path())?; + + assert_eq!( + sorted_sequences(original_fa.path())?, + sorted_sequences(roundtrip_fa.path())?, + "double revcomp should be identity for {mode:?}" + ); + + assert_eq!( + count_binseq(rc_once.path())?, + DEFAULT_NUM_RECORDS, + "revcomp record count wrong for {mode:?}" + ); + } + Ok(()) + } + + /// Reverse complementing a known sequence should produce the expected result. + #[test] + fn test_revcomp_known_sequence() -> Result<()> { + let seq = "ACGTACGTGATTACAACGTACGT"; + let in_tmp = NamedTempFile::with_suffix(".fastq")?; + { + use std::io::Write as _; + let mut f = std::fs::File::create(in_tmp.path())?; + writeln!(f, "@read1")?; + writeln!(f, "{seq}")?; + writeln!(f, "+")?; + writeln!(f, "{}", "I".repeat(seq.len()))?; + } + let bq_tmp = NamedTempFile::with_suffix(".cbq")?; + encode(in_tmp.path(), bq_tmp.path())?; + + let rc_tmp = NamedTempFile::with_suffix(".cbq")?; + revcomp(bq_tmp.path(), rc_tmp.path(), &[])?; + + let out_fa = NamedTempFile::with_suffix(".fasta")?; + decode_to_fasta(rc_tmp.path(), out_fa.path())?; + + let content = std::fs::read_to_string(out_fa.path())?; + assert!( + content.contains(&reverse_complement_str(seq)), + "expected reverse complement of {seq} in output: {content}" + ); + + Ok(()) + } + + /// With `-M 1`/`-M 2`, only the targeted mate should be reverse complemented. + #[test] + fn test_revcomp_paired_single_mate() -> Result<()> { + for (mode, mate_flag) in iproduct!(BinseqMode::enum_iter(), ["1", "2"]) { + let r1 = write_fastx().call()?; + let r2 = write_fastx().call()?; + let bq_tmp = NamedTempFile::with_suffix(mode.extension())?; + let cmd = crate::cli::EncodeCommand::try_parse_from([ + "encode", + r1.path().to_str().unwrap(), + r2.path().to_str().unwrap(), + "-o", + bq_tmp.path().to_str().unwrap(), + ])?; + crate::commands::encode::run(&cmd)?; + + let rc_tmp = NamedTempFile::with_suffix(mode.extension())?; + revcomp(bq_tmp.path(), rc_tmp.path(), &["-M", mate_flag])?; + + assert_eq!( + count_binseq(rc_tmp.path())?, + DEFAULT_NUM_RECORDS, + "revcomp single-mate record count wrong for {mode:?} mate={mate_flag}" + ); + } + Ok(()) + } +} diff --git a/src/commands/revcomp/processor.rs b/src/commands/revcomp/processor.rs new file mode 100644 index 0000000..54b9c7c --- /dev/null +++ b/src/commands/revcomp/processor.rs @@ -0,0 +1,206 @@ +use std::{io::Write, sync::Arc}; + +use binseq::{BinseqRecord, BinseqWriter, ParallelProcessor, SequencingRecordBuilder}; +use parking_lot::Mutex; + +use crate::cli::Mate; + +/// Reverse complements a nucleotide sequence buffer in place. +/// +/// Any byte outside `ACGTacgt` (e.g. `N`) is left untouched, matching the +/// behavior of 4-bit decoding, which collapses all ambiguity codes to `N`. +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, + }; + } +} + +pub struct RevCompProcessor { + /// Which mate(s) to reverse complement + mate: Mate, + + /// Thread-local writer for the processor + t_writer: BinseqWriter>, + /// Thread-local record count + t_count: usize, + + /// Thread-local scratch buffers for the transformed primary sequence/quality + sseq: Vec, + squal: Vec, + /// Thread-local scratch buffers for the transformed extended sequence/quality + xseq: Vec, + xqual: Vec, + + /// Global writer for the processor + writer: Arc>>, + /// Global record count + count: Arc>, +} +impl Clone for RevCompProcessor { + fn clone(&self) -> Self { + Self { + mate: self.mate, + t_writer: self.t_writer.clone(), + t_count: 0, + sseq: Vec::new(), + squal: Vec::new(), + xseq: Vec::new(), + xqual: Vec::new(), + writer: self.writer.clone(), + count: self.count.clone(), + } + } +} +impl RevCompProcessor { + pub fn new(writer: BinseqWriter, mate: Mate) -> binseq::Result { + let t_writer = writer.new_headless_buffer()?; + Ok(Self { + mate, + t_writer, + t_count: 0, + sseq: Vec::new(), + squal: Vec::new(), + xseq: Vec::new(), + xqual: Vec::new(), + writer: Arc::new(Mutex::new(writer)), + count: Arc::new(Mutex::new(0)), + }) + } + + fn write_batch(&mut self) -> binseq::Result<()> { + self.writer.lock().ingest_completed(&mut self.t_writer) + } + + fn write_final(&mut self) -> binseq::Result<()> { + self.writer.lock().ingest(&mut self.t_writer) + } + + pub fn finish(&mut self) -> binseq::Result<()> { + self.writer.lock().finish() + } + + pub fn get_global_record_count(&self) -> usize { + *self.count.lock() + } +} + +impl ParallelProcessor for RevCompProcessor { + fn process_record(&mut self, record: B) -> binseq::Result<()> { + let is_paired = record.is_paired(); + let has_quality = record.has_quality(); + let rc_primary = matches!(self.mate, Mate::One | Mate::Both); + let rc_extended = is_paired && matches!(self.mate, Mate::Two | Mate::Both); + + if rc_primary { + self.sseq.clear(); + self.sseq.extend_from_slice(record.sseq()); + reverse_complement(&mut self.sseq); + if has_quality { + self.squal.clear(); + self.squal.extend_from_slice(record.squal()); + self.squal.reverse(); + } + } + if rc_extended { + self.xseq.clear(); + self.xseq.extend_from_slice(record.xseq()); + reverse_complement(&mut self.xseq); + if has_quality { + self.xqual.clear(); + self.xqual.extend_from_slice(record.xqual()); + self.xqual.reverse(); + } + } + + let s_seq: &[u8] = if rc_primary { + &self.sseq + } else { + record.sseq() + }; + let s_qual: Option<&[u8]> = if !has_quality { + None + } else if rc_primary { + Some(&self.squal) + } else { + Some(record.squal()) + }; + + let rec = if is_paired { + let x_seq: &[u8] = if rc_extended { + &self.xseq + } else { + record.xseq() + }; + let x_qual: Option<&[u8]> = if !has_quality { + None + } else if rc_extended { + Some(&self.xqual) + } else { + Some(record.xqual()) + }; + SequencingRecordBuilder::default() + .s_seq(s_seq) + .opt_s_qual(s_qual) + .s_header(record.sheader()) + .x_seq(x_seq) + .opt_x_qual(x_qual) + .x_header(record.xheader()) + .build()? + } else { + SequencingRecordBuilder::default() + .s_seq(s_seq) + .opt_s_qual(s_qual) + .s_header(record.sheader()) + .build()? + }; + + if self.t_writer.push(rec)? { + self.t_count += 1; + } + Ok(()) + } + + fn on_batch_complete(&mut self) -> binseq::Result<()> { + *self.count.lock() += self.t_count; + self.t_count = 0; + self.write_batch() + } + + fn on_thread_complete(&mut self) -> binseq::Result<()> { + self.write_final() + } +} + +#[cfg(test)] +mod tests { + use super::reverse_complement; + + #[test] + fn test_reverse_complement_basic() { + let mut seq = b"ACGTACGT".to_vec(); + reverse_complement(&mut seq); + assert_eq!(seq, b"ACGTACGT"); + + let mut seq = b"GATTACA".to_vec(); + reverse_complement(&mut seq); + assert_eq!(seq, b"TGTAATC"); + } + + #[test] + fn test_reverse_complement_preserves_n() { + let mut seq = b"ACGTN".to_vec(); + reverse_complement(&mut seq); + assert_eq!(seq, b"NACGT"); + } +} diff --git a/src/main.rs b/src/main.rs index cde407f..6a77efb 100644 --- a/src/main.rs +++ b/src/main.rs @@ -52,6 +52,7 @@ fn main() -> Result<()> { Commands::Split(ref split) => commands::split::run(split), Commands::Pipe(ref pipe) => commands::pipe::run(pipe), Commands::Qc(ref qc) => commands::qc::run(qc), + Commands::Revcomp(ref revcomp) => commands::revcomp::run(revcomp), }?; trace!("done"); Ok(())