Feat/introduce revcomp

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.

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),
 }

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;

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<String>,
 
     #[clap(flatten)]
@@ -159,6 +159,11 @@ pub struct OutputBinseq {
 }
 impl OutputBinseq {
     pub fn as_writer(&self) -> Result<Box<dyn Write + Send>> {
+        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<OutputBinseqOptions> 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());
+    }
+}

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,
+}

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;

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<BinseqWriterBuilder> {
+    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<Vec<String>> {
+        let content = std::fs::read_to_string(path)?;
+        let mut seqs: Vec<String> = 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(())
+    }
+}