Merge pull request #85 from ArcInstitute/binseq-0.9.0

Binseq 0.9.0

Author: noamteyssier

.github/workflows/ci.yml

@@ -26,41 +26,44 @@ jobs:
       - name: Linting
         run: cargo clippy --verbose
 
-  example_read_write:
-    runs-on: ubuntu-latest
-    steps:
-      - uses: actions/checkout@v3
-      - name: run example
-        run: cargo run --release --example read_write
-
-  example_parallel:
+  example_grep:
     runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        ext: [bq, vbq, cbq]
     steps:
       - uses: actions/checkout@v3
-      - name: run example
-        run: cargo run --release --example parallel_processing
+      - name: run example ${{ matrix.ext }}
+        run: cargo run --release --example grep -- ./data/subset.${{ matrix.ext }} "ACGTACGT"
 
-  example_example:
+  example_range:
     runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        ext: [bq, vbq, cbq]
     steps:
       - uses: actions/checkout@v3
-      - name: run example
-        run: cargo run --release --example example
+      - name: run example ${{ matrix.ext }}
+        run: cargo run --release --example parallel_range -- ./data/subset.${{ matrix.ext }} 4 30 200
 
-  example_grep:
+  example_write:
     runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        ext: [bq, vbq, cbq]
     steps:
       - uses: actions/checkout@v3
-      - name: run example bq
-        run: cargo run --release --example grep ./data/subset.bq
-      - name: run example vbq
-        run: cargo run --release --example grep ./data/subset.vbq
+      - name: run example (single) ${{ matrix.ext }}
+        run: cargo run --release --example write -- ./data/subset_R1.fastq.gz -o ./output.${{ matrix.ext }}
+      - name: run example (paired) ${{ matrix.ext }}
+        run: cargo run --release --example write -- ./data/subset_R1.fastq.gz ./data/subset_R2.fastq.gz -o ./output.${{ matrix.ext }}
 
-  example_range:
+  example_read:
     runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        ext: [bq, vbq, cbq]
     steps:
       - uses: actions/checkout@v3
-      - name: run example (bq)
-        run: cargo run --release --example parallel_range -- ./data/subset.bq 4 30 200
-      - name: run example (vbq)
-        run: cargo run --release --example parallel_range -- ./data/subset.vbq 4 30 200
+      - name: run example ${{ matrix.ext }}
+        run: cargo run --release --example read -- ./data/subset.${{ matrix.ext }}

Cargo.toml

@@ -1,7 +1,7 @@
 [package]
 name = "binseq"
 version = "0.8.3"
-edition = "2021"
+edition = "2024"
 description = "A high efficiency binary format for sequencing data"
 license = "MIT"
 authors = ["Noam Teyssier <noam.teyssier@arcinstitute.org>"]
@@ -11,25 +11,32 @@ categories = ["science::bioinformatics", "encoding", "data-structures"]
 keywords = ["bioinformatics", "nucleotide", "sequencing", "genomics", "fastq"]
 
 [dependencies]
-anyhow = "1.0.100"
+anyhow = {version = "1.0.100", optional = true}
 auto_impl = "1.3.0"
-bitnuc = "0.3.2"
-bytemuck = "1.24.0"
+bitnuc = "0.4.0"
+bytemuck = { version = "1.24.0", features = ["derive", "extern_crate_alloc"] }
 byteorder = "1.5.0"
-itoa = "1.0.15"
+itoa = "1.0.17"
+memchr = "2.7.6"
 memmap2 = "0.9.9"
 num_cpus = "1.17.0"
+paraseq = { version = "0.4.8", optional = true }
+parking_lot = {version = "0.12.5", optional = true }
 rand = { version = "0.9.2", features = ["small_rng"] }
+sucds = "0.8.3"
 thiserror = "2.0.17"
 zstd = { version = "0.13.3", features = ["zstdmt"] }
 
 [dev-dependencies]
-nucgen = "0.2.0"
-niffler = "3.0.0"
-seq_io = "0.3.4"
+anyhow = "1.0.100"
 parking_lot = "0.12.5"
-itoa = "1.0.15"
-memchr = "2.7.6"
+clap = { version = "4.5.54", features = ["derive"] }
+paraseq = "0.4.8"
+
+[features]
+default = ["paraseq", "anyhow"]
+anyhow = ["dep:anyhow"]
+paraseq = ["dep:paraseq", "dep:parking_lot"]
 
 [lints.clippy]
 pedantic = { level = "warn", priority = -1 }

README.md

@@ -10,12 +10,18 @@
 BINSEQ is a binary file format family designed for efficient storage and processing of DNA sequences.
 They make use of two-bit encoding for nucleotides and are optimized for high-performance parallel processing.
 
-BINSEQ currently has two flavors:
+BINSEQ has three variants:
 
 1. **BQ**: (`*.bq`) files are for _fixed-length_ records **without** quality scores.
 2. **VBQ**: (`*.vbq`) files are for _variable-length_ records **with optional** quality scores and headers.
+3. **CBQ**: (`*.cbq`) files are for _columnar variable-length_ records **with optional** quality scores and headers.
 
-Both flavors support both single and paired sequences.
+All variants support both single and paired sequences.
+
+**Note:** For most use cases, the newest variant _CBQ_ is recommended due to its flexibility, storage efficiency, and decoding speed.
+It supersedes _VBQ_ in terms of performance and storage efficiency, at a small cost in encoding speed.
+VBQ will still be supported but newer projects should consider using _CBQ_ instead.
+For information on the structure of _CBQ_ files, see the [documentation](https://docs.rs/binseq/latest/binseq/cbq/).
 
 ## Getting Started
 
@@ -24,4 +30,4 @@ This is a **library** for reading and writing BINSEQ files, for a **command-line
 To get started please refer to our [documentation](https://docs.rs/binseq/latest/binseq/).
 For example programs which make use of the library check out our [examples directory](https://github.com/arcinstitute/binseq/tree/main/examples).
 
-For more information about the BINSEQ file family, please refer to our [preprint](https://www.biorxiv.org/content/10.1101/2025.04.08.647863v1).
+For more information about the BINSEQ file family, please refer to our [preprint](https://www.biorxiv.org/content/10.1101/2025.04.08.647863v2).

data/subset.cbq

@@ -0,0 +1,122 @@
+use std::{fs::File, io::BufWriter};
+
+use anyhow::Result;
+use binseq::{BinseqWriterBuilder, write::Format};
+use bitnuc::BitSize;
+use clap::Parser;
+
+type BoxedWriter = Box<dyn std::io::Write + Send>;
+
+#[derive(Parser)]
+struct Args {
+    /// Input FASTX to encode into BINSEQ format
+    #[clap(required = true)]
+    input: String,
+
+    /// Input FASTX to encode into BINSEQ format (R2)
+    #[clap(required = false)]
+    input2: Option<String>,
+
+    /// Output file path for BINSEQ format
+    #[clap(short = 'o', long)]
+    output: Option<String>,
+
+    /// Default prefix for writing BINSEQ: `<prefix>.<ext>`
+    #[clap(short = 'p', long, default_value = "output")]
+    prefix: String,
+
+    /// Format of the output BINSEQ file
+    ///
+    /// [bq: bq|BQ|b, vbq: vbq|VBQ|v, cbq: cbq|CBQ|c]
+    #[clap(short = 'f', long)]
+    format: Option<Format>,
+
+    /// Exclude quality information in BINSEQ output
+    ///
+    /// (bq ignores quality always)
+    #[clap(short = 'Q', long)]
+    exclude_quality: bool,
+
+    /// Exclude sequence headers in BINSEQ output
+    ///
+    /// (bq ignores headers always)
+    #[clap(short = 'H', long)]
+    exclude_headers: bool,
+
+    /// Compression level for BINSEQ output (0: auto)
+    #[clap(long, default_value_t = 0)]
+    compression_level: i32,
+
+    /// Default BITSIZE for BINSEQ output (2: 2bit, 4: 4bit)
+    #[clap(long, default_value_t = 2)]
+    bitsize: u8,
+
+    /// Default BLOCKSIZE in KB for BINSEQ output (vbq,cbq)
+    #[clap(long, default_value_t = 128)]
+    blocksize: usize,
+
+    /// Number of threads to use for parallel processing, 0: all available
+    #[clap(short = 'T', long, default_value = "0")]
+    threads: usize,
+}
+impl Args {
+    /// Determines the output format based on the file extension or the provided format
+    fn format(&self) -> Format {
+        if let Some(format) = self.format {
+            format
+        } else {
+            if let Some(output) = &self.output {
+                match output.split(".").last() {
+                    Some("bq") => Format::Bq,
+                    Some("vbq") => Format::Vbq,
+                    Some("cbq") => Format::Cbq,
+                    _ => Format::default(),
+                }
+            } else {
+                Format::default()
+            }
+        }
+    }
+    fn bitsize(&self) -> BitSize {
+        match self.bitsize {
+            4 => BitSize::Four,
+            _ => BitSize::Two,
+        }
+    }
+
+    /// Creates an output file handle
+    fn ohandle(&self) -> Result<BoxedWriter> {
+        let path = if let Some(output) = &self.output {
+            output.to_string()
+        } else {
+            format!("{}{}", &self.prefix, self.format().extension())
+        };
+        let ofile = File::create(path).map(BufWriter::new)?;
+        Ok(Box::new(ofile))
+    }
+
+    fn is_paired(&self) -> bool {
+        self.input2.is_some()
+    }
+}
+
+fn main() -> Result<()> {
+    let args = Args::parse();
+    let handle = args.ohandle()?;
+    let builder = BinseqWriterBuilder::new(args.format())
+        .bitsize(args.bitsize())
+        .block_size(args.blocksize * 1024)
+        .headers(!args.exclude_headers)
+        .quality(!args.exclude_quality)
+        .compression_level(args.compression_level)
+        .encode_fastx(handle);
+    if args.is_paired() {
+        builder.input_paired(&args.input, args.input2.as_ref().unwrap())
+    } else {
+        builder.input(&args.input)
+    }
+    .threads(args.threads)
+    .run()?;
+
+    Ok(())
+}

data/subset.vbq

@@ -1,14 +1,14 @@
 use std::sync::Arc;
 
 use anyhow::Result;
-use binseq::{context::SeqCtx, prelude::*};
+use binseq::prelude::*;
+use clap::Parser;
 use memchr::memmem::Finder;
 use parking_lot::Mutex;
 
 #[derive(Clone)]
 pub struct GrepCounter {
     // (thread) local variables
-    ctx: SeqCtx,
     local_count: usize,
 
     // search pattern (using memchr::memmem::Finder for fast searching)
@@ -21,7 +21,6 @@ impl GrepCounter {
     #[must_use]
     pub fn new(pattern: &[u8]) -> Self {
         Self {
-            ctx: SeqCtx::default(),
             pattern: Finder::new(pattern).into_owned(),
             local_count: 0,
             count: Arc::new(Mutex::new(0)),
@@ -38,9 +37,7 @@ impl GrepCounter {
 }
 impl ParallelProcessor for GrepCounter {
     fn process_record<R: binseq::BinseqRecord>(&mut self, record: R) -> binseq::Result<()> {
-        self.ctx.fill(&record)?;
-
-        if self.match_sequence(&self.ctx.sbuf()) || self.match_sequence(&self.ctx.xbuf()) {
+        if self.match_sequence(&record.sseq()) || self.match_sequence(&record.xseq()) {
             self.local_count += 1;
         }
 
@@ -54,21 +51,26 @@ impl ParallelProcessor for GrepCounter {
     }
 }
 
-fn main() -> Result<()> {
-    let path = std::env::args()
-        .nth(1)
-        .unwrap_or("./data/subset.bq".to_string());
-    let pattern = std::env::args()
-        .nth(2)
-        .unwrap_or("ACGT".to_string())
-        .as_bytes()
-        .to_vec();
-    let n_threads = std::env::args().nth(3).unwrap_or("1".to_string()).parse()?;
+#[derive(Parser)]
+struct Args {
+    /// Input BINSEQ path to grep
+    #[clap(required = true)]
+    input: String,
 
-    let reader = BinseqReader::new(&path)?;
-    let counter = GrepCounter::new(&pattern);
-    reader.process_parallel(counter.clone(), n_threads)?;
-    counter.pprint();
+    /// Pattern to search for (either sseq or xseq)
+    #[clap(required = true)]
+    pattern: String,
 
+    /// Threads to use [0: auto]
+    #[clap(short = 'T', long, default_value_t = 0)]
+    threads: usize,
+}
+
+fn main() -> Result<()> {
+    let args = Args::parse();
+    let reader = BinseqReader::new(&args.input)?;
+    let counter = GrepCounter::new(args.pattern.as_bytes());
+    reader.process_parallel(counter.clone(), args.threads)?;
+    counter.pprint();
     Ok(())
 }

examples/auto-write.rs

@@ -2,10 +2,10 @@ use std::io::{BufReader, BufWriter};
 use std::net::{TcpListener, TcpStream};
 use std::thread;
 
-use binseq::bq::{BinseqHeader, BinseqHeaderBuilder, StreamReader, StreamWriterBuilder};
+use binseq::bq::{FileHeader, FileHeaderBuilder, StreamReader, StreamWriterBuilder};
 use binseq::{BinseqRecord, Policy, Result};
 
-fn server(header: BinseqHeader, sequence: &[u8]) -> Result<()> {
+fn server(header: FileHeader, sequence: &[u8]) -> Result<()> {
     // Create a listener on localhost:3000
     let listener = TcpListener::bind("127.0.0.1:3000").expect("Failed to bind to address");
     println!("Server listening on 127.0.0.1:3000");
@@ -25,6 +25,7 @@ fn server(header: BinseqHeader, sequence: &[u8]) -> Result<()> {
 
     // Write sequences in a loop
     for i in 0..10 {
+        #[allow(deprecated)]
         writer.write_record(Some(i), sequence)?;
         println!("Server: Sent record {i}");
 
@@ -79,7 +80,7 @@ fn client() -> Result<()> {
 
 fn main() -> Result<()> {
     // Create a header for sequences of length 100
-    let header = BinseqHeaderBuilder::new().slen(100).build()?;
+    let header = FileHeaderBuilder::new().slen(100).build()?;
 
     // Create some example sequence data
     let sequence = b"ACGT".repeat(25); // 100 nucleotides