From f1f2ba14a47cb2f528c4a3992abe7e39ef36e88e Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 1 Jul 2026 09:41:24 +0200 Subject: [PATCH 1/3] fix(cbq): restore N-position backfill in serial decode (#94) CBQ stores N bases as a 2-bit placeholder A plus an Elias-Fano index of the N-positions that is backfilled on decode. In the serial Reader path, decompress_columns pre-filled ef_bytes with `resize(len_nef, 0)` and then called copy_decode, whose io::Write for Vec appends. The serialized EliasFano bytes therefore landed after len_nef zero bytes, so deserialize_from read the leading zeros and produced an empty index (num_ones == 0). backfill_npos then panicked in debug (sucds debug_assert) or silently no-oped in release, leaving every N decoded as A. Clear ef_bytes before copy_decode so the buffer holds exactly the decompressed bytes, mirroring the headers/qual columns. Also reset len_nef in ColumnarBlock::clear() for consistent block reuse. The mmap process_parallel path (dctx.decompress writes at offset 0) was already correct, and the on-disk bytes are unchanged, so no format bump is needed. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/cbq/core/block.rs | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/cbq/core/block.rs b/src/cbq/core/block.rs index 96b75b9..69d91ff 100644 --- a/src/cbq/core/block.rs +++ b/src/cbq/core/block.rs @@ -112,6 +112,7 @@ impl ColumnarBlock { self.num_records = 0; self.current_size = 0; self.num_npos = 0; + self.len_nef = 0; } // clear spans @@ -450,7 +451,12 @@ impl ColumnarBlock { // decompress npos if !self.z_npos.is_empty() { - self.ef_bytes.resize(self.len_nef, 0); + // Clear (do NOT `resize(len_nef, 0)`): `copy_decode`'s destination is a + // `&mut Vec`, whose `io::Write` impl *appends*. Pre-filling with zeros would + // leave `ef_bytes == [0; len_nef] ++ [real EF bytes]`, so `deserialize_from` would + // read the leading zeros and reconstruct an empty EliasFano (`num_ones() == 0`), + // silently dropping every N on decode. See ArcInstitute/binseq#94. + self.ef_bytes.clear(); copy_decode(self.z_npos.as_slice(), &mut self.ef_bytes)?; let ef = EliasFano::deserialize_from(self.ef_bytes.as_slice())?; From 3d2d4d153989f0598a232678251beb9520aae75d Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 1 Jul 2026 09:41:29 +0200 Subject: [PATCH 2/3] test(cbq): add N-nucleotide round-trip regression tests (#94) Existing CBQ tests were ACGT-only, which is why the N-decode bug shipped. Add coverage for N at leading, interior, trailing, and all-N positions, across single and multiple blocks, with and without quality plus headers, over both the serial Reader path and the parallel MmapReader / process_parallel path. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/cbq/read.rs | 75 +++++++++++++++++++++++++ src/cbq/write.rs | 141 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 216 insertions(+) diff --git a/src/cbq/read.rs b/src/cbq/read.rs index ab07913..80438b1 100644 --- a/src/cbq/read.rs +++ b/src/cbq/read.rs @@ -634,4 +634,79 @@ mod tests { let final_count = *count.lock().unwrap(); assert_eq!(final_count, num_records); } + + // ==================== N-nucleotide round-trip (issue #94) ==================== + + /// Collects each record's `(global index, decoded sequence)` during parallel + /// processing so the result can be re-ordered and compared to the input. + #[derive(Clone)] + struct SeqCollector { + seqs: Arc)>>>, + } + impl ParallelProcessor for SeqCollector { + fn process_record(&mut self, record: R) -> Result<()> { + let mut buf = Vec::new(); + record.decode_s(&mut buf)?; + self.seqs.lock().unwrap().push((record.index(), buf)); + Ok(()) + } + } + + /// The mmap / `process_parallel` decode path must restore every `N`. + /// Complements the serial-path tests in `write.rs` for ArcInstitute/binseq#94. + #[test] + fn test_n_roundtrip_parallel_mmap() -> Result<()> { + use crate::{BinseqWriterBuilder, SequencingRecordBuilder, write::Format}; + + // N-containing sequences; a small block size forces multiple blocks so + // several threads each decode N-bearing blocks. + let seqs: Vec> = (0..128u32) + .map(|i| { + if i % 7 == 0 { + vec![b'N'; 8] + } else { + let mut s = vec![b'A', b'C', b'G', b'T', b'A', b'C', b'G', b'T']; + let pos = i as usize % s.len(); + s[pos] = b'N'; + s + } + }) + .collect(); + + let path = + std::env::temp_dir().join(format!("binseq_n94_parallel_{}.cbq", std::process::id())); + + // Write to a real file (MmapReader requires a path). + { + let mut writer = BinseqWriterBuilder::new(Format::Cbq) + .block_size(64) + .build(std::fs::File::create(&path)?)?; + for seq in &seqs { + writer.push(SequencingRecordBuilder::default().s_seq(seq).build()?)?; + } + writer.finish()?; + } + + // Decode back in parallel via mmap. + let collected = Arc::new(std::sync::Mutex::new(Vec::<(u64, Vec)>::new())); + { + let reader = MmapReader::new(&path)?; + assert!(reader.num_blocks() > 1, "expected multiple blocks"); + let processor = SeqCollector { + seqs: collected.clone(), + }; + reader.process_parallel(processor, 4)?; + } + + let _ = std::fs::remove_file(&path); + + let mut got = Arc::try_unwrap(collected).unwrap().into_inner().unwrap(); + got.sort_by_key(|(idx, _)| *idx); + let got_seqs: Vec> = got.into_iter().map(|(_, s)| s).collect(); + assert_eq!( + got_seqs, seqs, + "parallel mmap decode must round-trip N positions" + ); + Ok(()) + } } diff --git a/src/cbq/write.rs b/src/cbq/write.rs index 6301f7d..9de6ad8 100644 --- a/src/cbq/write.rs +++ b/src/cbq/write.rs @@ -451,4 +451,145 @@ mod tests { Ok(()) } + + // ==================== N-nucleotide round-trip (issue #94) ==================== + // + // CBQ stores `N` as a 2-bit placeholder `A` plus an Elias-Fano index of the + // N-positions that is used to backfill `N` on decode. These tests guard the + // serial `Reader`/`decompress_columns` path, which previously dropped every N + // (decoding it back as `A`). See ArcInstitute/binseq#94. + + /// Build a `FileHeader` with quality scores and headers enabled. + fn header_with_quality_headers(block_size: usize) -> FileHeader { + FileHeaderBuilder::default() + .is_paired(false) + .with_headers(true) + .with_qualities(true) + .with_flags(false) + .with_block_size(block_size) + .build() + } + + /// Read back every record's `(sequence, quality, header)` from a finished CBQ buffer. + fn read_all_records(bytes: Vec) -> Vec<(Vec, Vec, Vec)> { + let mut reader = Reader::new(Cursor::new(bytes)).expect("failed to open reader"); + let mut out = Vec::new(); + let mut cumulative = 0u64; + while let Some(block_header) = reader.read_block().expect("failed to read block") { + cumulative += block_header.num_records; + reader + .block + .decompress_columns() + .expect("failed to decompress block"); + let range = BlockRange::new(0, cumulative); + for rec in reader.block.iter_records(range) { + out.push(( + rec.sseq().to_vec(), + rec.squal().to_vec(), + rec.sheader().to_vec(), + )); + } + } + out + } + + /// Sequences with `N` at every notable position: leading, interior, trailing, + /// all-N, and a no-N control. Direct reproduction of ArcInstitute/binseq#94. + #[test] + fn test_n_roundtrip_serial_edge_positions() -> Result<()> { + let seqs: Vec<&[u8]> = vec![ + b"NACGTACGT", // leading + b"ACGTNACGT", // interior + b"ACGTACGTN", // trailing + b"NNNNNNNN", // all-N + b"ACGTACGT", // no-N control + ]; + + // Large block size: all records land in a single block. + let mut writer = ColumnarBlockWriter::new(Vec::new(), header(1 << 20))?; + for &seq in &seqs { + writer.push(record(seq))?; + } + writer.finish()?; + + let read_back = read_all_sequences(writer.inner); + let expected: Vec> = seqs.iter().map(|s| s.to_vec()).collect(); + assert_eq!(read_back, expected, "N positions must round-trip verbatim"); + Ok(()) + } + + /// Many N-containing sequences with a small block size, forcing multiple + /// blocks. Guards per-block Elias-Fano correctness and `ef_bytes`/`len_nef` + /// reuse across `clear()`. + #[test] + fn test_n_roundtrip_serial_multi_block() -> Result<()> { + let seqs: Vec> = (0..64u32) + .map(|i| { + if i % 5 == 0 { + // occasionally all-N + vec![b'N'; 8] + } else { + // vary N placement so per-block position sets differ + let mut s = vec![b'A', b'C', b'G', b'T', b'A', b'C', b'G', b'T']; + let pos = i as usize % s.len(); + s[pos] = b'N'; + s + } + }) + .collect(); + + let mut writer = ColumnarBlockWriter::new(Vec::new(), header(64))?; + for seq in &seqs { + writer.push(record(seq))?; + } + writer.finish()?; + + // Sanity: the small block size must actually span more than one block. + assert!( + writer.headers.len() > 1, + "expected multiple blocks, got {}", + writer.headers.len() + ); + + let read_back = read_all_sequences(writer.inner); + assert_eq!( + read_back, seqs, + "N round-trip must hold across multiple blocks" + ); + Ok(()) + } + + /// Mirrors the exact reproduction in ArcInstitute/binseq#94: quality + headers + /// enabled, with N at leading / all / interior positions. + #[test] + fn test_n_roundtrip_with_quality_and_headers() -> Result<()> { + let records: [(&[u8], &[u8], &[u8]); 3] = [ + (b"r1", b"NACGTACGT", b"IIIIIIIII"), + (b"r2", b"NNNNNNNN", b"IIIIIIII"), + (b"r3", b"ACGTNACGT", b"IIIIIIIII"), + ]; + + let mut writer = + ColumnarBlockWriter::new(Vec::new(), header_with_quality_headers(1 << 20))?; + for (h, s, q) in records { + let rec = SequencingRecordBuilder::default() + .s_header(h) + .s_seq(s) + .s_qual(q) + .build()?; + writer.push(rec)?; + } + writer.finish()?; + + let read_back = read_all_records(writer.inner); + let expected: Vec<(Vec, Vec, Vec)> = records + .iter() + .map(|(h, s, q)| (s.to_vec(), q.to_vec(), h.to_vec())) + .collect(); + assert_eq!( + read_back, expected, + "seq/qual/header must round-trip with N present" + ); + Ok(()) + } } From 03e108b88a17ee6186e6e6815db901e55ad50920 Mon Sep 17 00:00:00 2001 From: Benjamin Demaille Date: Wed, 1 Jul 2026 10:29:00 +0200 Subject: [PATCH 3/3] perf(cbq): reserve ef_bytes capacity before npos decode (#94) Address PR review: pre-size ef_bytes to the known decompressed length (len_nef) after clearing it, so copy_decode's appends avoid incremental reallocations. Capacity-only change; decode semantics are unchanged. Co-Authored-By: Claude Opus 4.8 (1M context) --- src/cbq/core/block.rs | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/cbq/core/block.rs b/src/cbq/core/block.rs index 69d91ff..3737d45 100644 --- a/src/cbq/core/block.rs +++ b/src/cbq/core/block.rs @@ -457,6 +457,9 @@ impl ColumnarBlock { // read the leading zeros and reconstruct an empty EliasFano (`num_ones() == 0`), // silently dropping every N on decode. See ArcInstitute/binseq#94. self.ef_bytes.clear(); + // Pre-size (capacity only, not length): the decompressed EF is exactly + // `len_nef` bytes, so this avoids reallocations while `copy_decode` appends. + self.ef_bytes.reserve(self.len_nef); copy_decode(self.z_npos.as_slice(), &mut self.ef_bytes)?; let ef = EliasFano::deserialize_from(self.ef_bytes.as_slice())?;