Merge pull request #164 from ArcInstitute/dev-0.5.1

Dev 0.5.1

Author: noamteyssier

CLAUDE.md

@@ -0,0 +1,74 @@
+# CLAUDE.md
+
+This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository.
+
+## Project Overview
+
+bqtools is a Rust CLI for working with BINSEQ files — a binary format family for high-performance DNA sequence processing. It encodes, decodes, greps, concatenates, samples, and pipes BINSEQ files (`.bq`, `.vbq`, `.cbq`). CBQ is the recommended format for most applications.
+
+## Build & Test Commands
+
+```bash
+cargo build                    # Debug build
+cargo build --release          # Optimized build (uses LTO, slow)
+cargo install --path .         # Install binary locally
+
+cargo test --verbose           # Run all tests
+cargo test --verbose -F fuzzy  # Run tests including fuzzy feature
+cargo test <test_name>         # Run a single test by name
+
+cargo fmt --check              # Check formatting
+cargo clippy --verbose         # Lint (pedantic clippy enabled)
+```
+
+Logging is controlled via `BQTOOLS_LOG` environment variable (uses `env_logger`).
+
+## Feature Flags
+
+- `htslib` (default): SAM/BAM/CRAM support via rust-htslib
+- `gcs` (default): Google Cloud Storage file reading
+- `fuzzy` (optional): Fuzzy matching via `sassy` — requires `RUSTFLAGS="-C target-cpu=native"`
+
+Build without defaults: `cargo build --no-default-features -F fuzzy,gcs`
+
+## Architecture
+
+### Module Layout
+
+- **`src/cli/`** — Clap derive-based argument definitions. `cli.rs` has the top-level `Commands` enum. `input.rs` and `output.rs` handle complex input/output argument parsing (file formats, compression, paired-end, spans).
+- **`src/commands/`** — Command implementations, each in its own subdirectory. `utils.rs` has shared compression helpers.
+- **`src/types.rs`** — Type aliases (`BoxedReader`, `BoxedWriter`).
+- **`src/main.rs`** — CLI dispatch and SIGPIPE handling.
+
+### Key Patterns
+
+**Parallel processing**: Commands use the `paraseq` crate's `ParallelProcessor` trait for embarrassingly parallel batch processing. Each command has a `processor.rs` implementing this trait with thread-local buffers and `Arc<Mutex<T>>` for shared global state.
+
+**Grep backends**: The grep command uses a `PatternMatcher` enum dispatching to three backends — `regex`, `aho-corasick` (fixed-string, multi-pattern), and `sassy` (fuzzy, feature-gated). The same pattern applies to `PatternCounter` for the `-P` pattern-count mode.
+
+**Encode modes**: Encoding dispatches across atomic (single/paired files), recursive (directory walk via `walkdir`), manifest (file list), and batch (multi-file thread distribution) modes.
+
+**Writer abstraction**: `SplitWriter` supports interleaved (single file) and split (separate R1/R2) output modes with polymorphic writers (file, stdout, compressed, chunked).
+
+### Core Dependencies
+
+| Crate | Role |
+|-------|------|
+| `binseq` | BINSEQ format read/write |
+| `bitnuc` | 2-bit/4-bit nucleotide encoding |
+| `paraseq` | Parallel FASTX/BINSEQ processing |
+| `clap` | CLI argument parsing (derive) |
+| `anyhow` | Error handling throughout |
+
+### Testing
+
+Integration tests live in `tests/`. `tests/common.rs` provides a builder (`write_fastx()`) for generating random FASTQ/FASTA test data with configurable compression (none, gzip, zstd). Tests use cartesian products over format/compression/mode combinations. Dev dependencies: `bon` (builder macro), `nucgen` (random sequences), `tempfile`, `itertools`.
+
+## Contribution Guide
+
+When making changes, keep the following documentation in sync:
+
+1. **CLAUDE.md** — Update this file when adding new commands, changing architecture, or modifying build/test workflows.
+2. **README.md** — Update usage examples and feature descriptions when adding or changing user-facing functionality (new commands, flags, behavior changes).
+3. **Clap doc comments** — All CLI arguments, flags, and subcommands use clap derive macros with `/// doc comments` and `#[clap(long_about)]` attributes. When adding or modifying flags, write clear help text directly on the struct fields in `src/cli/`. These doc comments are the `--help` output users see.
+4. **New feature flags** — If adding a Cargo feature flag, document it in both `CLAUDE.md` (Feature Flags section) and `README.md` (Feature Flags / Installation section).

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "bqtools"
-version = "0.5.0"
+version = "0.5.1"
 edition = "2021"
 license = "MIT"
 authors = ["Noam Teyssier <noam.teyssier@arcinstitute.org>"]

README.md

@@ -342,6 +342,18 @@ bqtools grep input.bq --xfile patterns.txt
 bqtools grep input.bq --file patterns.txt -x
 ```
 
+You can count the number of matching records with `-C` or get the fraction of matching records with `--frac`:
+
+```bash
+# Count the number of matching records
+bqtools grep input.bq "ACGTACGT" -C
+
+# Count matching records and show fraction of total
+bqtools grep input.bq "ACGTACGT" -F
+```
+
+The output of `--frac` is a TSV with three columns: [Count, Total, Fraction]
+
 `bqtools` also introduces a new feature for the counting the occurrences of individual patterns.
 This is useful for seeing how many times each pattern occurs across a sequencing dataset without having to iterate over the dataset multiple times using traditional methods.
 

src/cli/grep.rs

@@ -52,6 +52,13 @@ pub struct GrepArgs {
     #[clap(short = 'C', long, conflicts_with = "pattern_count")]
     pub count: bool,
 
+    /// Show match count as a fraction of total records
+    ///
+    /// Implies --count (-C). Displays the number of matches,
+    /// total records, and the fraction of records matching.
+    #[clap(short = 'F', long, conflicts_with = "pattern_count")]
+    pub frac: bool,
+
     /// Only match patterns that are within this range.
     ///
     /// Will not match if the pattern is outside the range or if

src/cli/input.rs

@@ -1,8 +1,8 @@
-use std::path::PathBuf;
+use std::{path::PathBuf, str::FromStr};
 
 use anyhow::{bail, Result};
 use clap::Parser;
-use log::debug;
+use log::{debug, error, warn};
 use paraseq::fastx;
 
 #[cfg(not(feature = "gcs"))]
@@ -34,7 +34,8 @@ pub struct InputFile {
     #[clap(short, long)]
     pub batch_size: Option<usize>,
 
-    #[clap(short = 'I', long, help = "Interleaved input file format")]
+    /// Input is paired-interleaved
+    #[clap(short = 'I', long, conflicts_with = "paired")]
     pub interleaved: bool,
 
     /// Apply encoding to all fasta/fastq files in the provided directory input.
@@ -242,6 +243,10 @@ pub struct BatchEncodingOptions {
 pub struct InputBinseq {
     #[clap(help = "Input binseq file")]
     pub input: String,
+
+    /// Span of records to process. If not specified, all records will be processed.
+    #[clap(long)]
+    pub span: Option<Span>,
 }
 impl InputBinseq {
     pub fn path(&self) -> &str {
@@ -256,3 +261,67 @@ pub struct MultiInputBinseq {
     #[clap(num_args = 1..)]
     pub input: Vec<String>,
 }
+
+#[derive(Debug, Clone, Copy)]
+pub struct Span {
+    start: Option<usize>,
+    end: Option<usize>,
+}
+impl Span {
+    fn validate(&mut self, max_records: usize) -> Result<()> {
+        if let Some(start) = self.start {
+            if start > max_records {
+                error!(
+                    "Provided start ({}) exceeds maximum number of records ({})",
+                    start, max_records
+                );
+                bail!("Maximum number of records exceeded")
+            }
+        }
+        if let Some(end) = self.end {
+            if end > max_records {
+                warn!(
+                    "Clipping provided endpoint ({}) to maximum number of records ({})",
+                    end, max_records
+                );
+            }
+            self.end = Some(end.min(max_records));
+        }
+        Ok(())
+    }
+    pub fn get_range(&mut self, max_records: usize) -> Result<std::ops::Range<usize>> {
+        self.validate(max_records)?;
+        match (self.start, self.end) {
+            (Some(start), Some(end)) => Ok(start..end),
+            (Some(start), None) => Ok(start..max_records),
+            (None, Some(end)) => Ok(0..end),
+            (None, None) => Ok(0..max_records),
+        }
+    }
+}
+
+impl FromStr for Span {
+    type Err = String;
+
+    fn from_str(s: &str) -> Result<Self, Self::Err> {
+        let (start_str, end_str) = s
+            .split_once("..")
+            .ok_or_else(|| format!("expected range like '10..20', got '{s}'"))?;
+
+        let parse_bound = |bound_str: &str, name: &str| {
+            if bound_str.is_empty() {
+                Ok(None)
+            } else {
+                bound_str
+                    .parse()
+                    .map(Some)
+                    .map_err(|_| format!("invalid {name}: '{bound_str}'"))
+            }
+        };
+
+        let start = parse_bound(start_str, "start")?;
+        let end = parse_bound(end_str, "end")?;
+
+        Ok(Self { start, end })
+    }
+}

src/commands/decode/mod.rs

@@ -59,7 +59,16 @@ pub fn run(args: &DecodeCommand) -> Result<()> {
         None
     };
     let proc = Decoder::new(writer, format, mate);
-    reader.process_parallel(proc.clone(), args.output.threads())?;
+    if let Some(mut span) = args.input.span {
+        let num_records = reader.num_records()?;
+        reader.process_parallel_range(
+            proc.clone(),
+            args.output.threads(),
+            span.get_range(num_records)?,
+        )?
+    } else {
+        reader.process_parallel(proc.clone(), args.output.threads())?;
+    }
     let num_records = proc.num_records();
     info!("Processed {num_records} records...");
     Ok(())

src/commands/grep/filter/processor.rs

@@ -25,12 +25,18 @@ pub struct FilterProcessor<Pm: PatternMatch> {
     /// Only count the number of matches
     count: bool,
 
+    /// Show count as fraction of total records
+    frac: bool,
+
     /// Match within range
     range: Option<SimpleRange>,
 
     /// Local count
     local_count: usize,
 
+    /// Local total records processed
+    local_total: usize,
+
     /// Local primary/extended sequence match indices
     smatches: MatchRanges,
     xmatches: MatchRanges,
@@ -54,6 +60,7 @@ pub struct FilterProcessor<Pm: PatternMatch> {
     /// Global values
     global_writer: Arc<Mutex<SplitWriter>>,
     global_count: Arc<Mutex<usize>>,
+    global_total: Arc<Mutex<usize>>,
 }
 impl<Pm: PatternMatch> FilterProcessor<Pm> {
     #[allow(clippy::fn_params_excessive_bools)]
@@ -63,6 +70,7 @@ impl<Pm: PatternMatch> FilterProcessor<Pm> {
         and_logic: bool,
         invert: bool,
         count: bool,
+        frac: bool,
         range: Option<SimpleRange>,
         writer: SplitWriter,
         format: FileFormat,
@@ -82,14 +90,17 @@ impl<Pm: PatternMatch> FilterProcessor<Pm> {
             and_logic,
             invert,
             count,
+            frac,
             range,
             format,
             mate,
             color,
             is_split: writer.is_split(),
             global_writer: Arc::new(Mutex::new(writer)),
             local_count: 0,
+            local_total: 0,
             global_count: Arc::new(Mutex::new(0)),
+            global_total: Arc::new(Mutex::new(0)),
         }
     }
     pub fn clear_matches(&mut self) {
@@ -132,13 +143,26 @@ impl<Pm: PatternMatch> FilterProcessor<Pm> {
         }
     }
     pub fn pprint_counts(&self) {
-        println!("{}", self.global_count.lock());
+        let count = *self.global_count.lock();
+        if self.frac {
+            let total = *self.global_total.lock();
+            let frac = if total > 0 {
+                count as f64 / total as f64
+            } else {
+                0.0
+            };
+            println!("count\ttotal\tfrac");
+            println!("{count}\t{total}\t{frac:.4}");
+        } else {
+            println!("{count}");
+        }
     }
 }
 
 impl<Pm: PatternMatch> ParallelProcessor for FilterProcessor<Pm> {
     fn process_record<B: BinseqRecord>(&mut self, record: B) -> binseq::Result<()> {
         self.clear_matches();
+        self.local_total += 1;
 
         let sbuf = record.sseq();
         let xbuf = record.xseq();
@@ -232,6 +256,10 @@ impl<Pm: PatternMatch> ParallelProcessor for FilterProcessor<Pm> {
         *self.global_count.lock() += self.local_count;
         self.local_count = 0;
 
+        // Increment the global total and reset local
+        *self.global_total.lock() += self.local_total;
+        self.local_total = 0;
+
         Ok(())
     }
 }

src/commands/grep/mod.rs

@@ -61,7 +61,16 @@ fn build_counter(args: &GrepCommand) -> Result<PatternCounter> {
 fn run_pattern_count(args: &GrepCommand, reader: BinseqReader) -> Result<()> {
     let counter = build_counter(args)?;
     let proc = PatternCountProcessor::new(counter, args.grep.range);
-    reader.process_parallel(proc.clone(), args.output.threads())?;
+    if let Some(mut span) = args.input.span {
+        let num_records = reader.num_records()?;
+        reader.process_parallel_range(
+            proc.clone(),
+            args.output.threads(),
+            span.get_range(num_records)?,
+        )?
+    } else {
+        reader.process_parallel(proc.clone(), args.output.threads())?;
+    }
     proc.pprint_pattern_counts()?;
     Ok(())
 }
@@ -110,20 +119,32 @@ fn run_grep(
     format: FileFormat,
     mate: Option<Mate>,
 ) -> Result<()> {
+    let count = args.grep.count || args.grep.frac;
     let matcher = build_matcher(args)?;
     let proc = FilterProcessor::new(
         matcher,
         args.grep.and_logic(),
         args.grep.invert,
-        args.grep.count,
+        count,
+        args.grep.frac,
         args.grep.range,
         writer,
         format,
         mate,
         args.should_color(),
     );
-    reader.process_parallel(proc.clone(), args.output.threads())?;
-    if args.grep.count {
+
+    if let Some(mut span) = args.input.span {
+        let num_records = reader.num_records()?;
+        reader.process_parallel_range(
+            proc.clone(),
+            args.output.threads(),
+            span.get_range(num_records)?,
+        )?;
+    } else {
+        reader.process_parallel(proc.clone(), args.output.threads())?;
+    }
+    if count {
         proc.pprint_counts();
     }
 

src/commands/pipe/mod.rs

@@ -6,7 +6,7 @@ use std::thread;
 
 use anyhow::Result;
 use binseq::BinseqReader;
-use log::info;
+use log::{info, warn};
 
 use crate::cli::{FileFormat, PipeCommand};
 use processor::PipeProcessor;
@@ -23,6 +23,10 @@ pub enum RecordPair {
 }
 
 pub fn run(args: &PipeCommand) -> Result<()> {
+    if args.input.span.is_some() {
+        warn!("Span is ignored when using pipe subcommand");
+    }
+
     let format = args.format()?;
     let reader = BinseqReader::new(args.input.path())?;
     let num_records = reader.num_records()?;