CBQ mis-decodes N-containing records (silent N→A in release, panic in debug)

[BenjaminDEMAILLE]

Summary

A .cbq file containing N bases does not round-trip: on decode, every N is emitted as A. In a debug build this panics; in a release build it is silent (the placeholder base is returned as real signal). Tested with binseq 0.9.2 + sucds 0.8.3 (both current-latest).

Reproduction

Write a CBQ with the canonical writer and N-containing sequences, then read it back:

use binseq::{BinseqRecord, prelude::*};
use binseq::cbq::{Reader, BlockRange};

// write
let mut w = BinseqWriterBuilder::new(Format::Cbq)
    .paired(false).quality(true).headers(true)
    .build(std::fs::File::create("n.cbq").unwrap()).unwrap();
for (h, s, q) in [("r1", b"NACGTACGT".as_slice(), b"IIIIIIIII".as_slice()),
                  ("r2", b"NNNNNNNN".as_slice(),  b"IIIIIIII".as_slice()),
                  ("r3", b"ACGTNACGT".as_slice(), b"IIIIIIIII".as_slice())] {
    let rec = SequencingRecordBuilder::default()
        .s_header(h.as_bytes()).s_seq(s).s_qual(q).build().unwrap();
    w.push(rec).unwrap();
}
w.finish().unwrap();

// read (serial streaming path)
let mut r = Reader::new(std::fs::File::open("n.cbq").unwrap()).unwrap();
let mut cumulative = 0u64;
while let Some(bh) = r.read_block().unwrap() {
    cumulative += bh.num_records;
    r.block.decompress_columns().unwrap();
    for rec in r.block.iter_records(BlockRange::new(0, cumulative)) {
        let mut seq = Vec::new();
        rec.decode_s(&mut seq).unwrap();
        println!("{}", String::from_utf8_lossy(&seq));
    }
}

Observed

• Debug build: panic at sucds-0.8.3/src/mii_sequences/elias_fano/iter.rs:23 — debug_assert_ne!(high_bits.num_ones(), 0).
• Release build: silent corruption:

stored	decoded
NACGTACGT	AACGTACGT
NNNNNNNN	AAAAAAAA
ACGTNACGT	ACGTAACGT

The parallel process_parallel path and bqtools decode on the same file corrupt identically (consistent with the defect being in decode, not the caller).

Expected

N bases round-trip as N.

Mechanism (best guess)

CBQ stores bases 2-bit (N as a placeholder A) plus an Elias-Fano index of the N-positions; ColumnarBlock::decompress_columns() calls backfill_npos() to restore the Ns. The deserialized npos Elias-Fano is degenerate — EliasFano::len() >= 1 but high_bits.num_ones() == 0 — so ef.iter(0) asserts (debug) or yields nothing (release), and the placeholder As are never overwritten. Root cause looks like the sucds EliasFano serialize/deserialize round-trip for the N-position set (or how binseq encodes/serializes that EF).

Impact

For a bisulfite/methylation consumer this is the worst failure class: a decoded A where the read had a no-call N is treated as real signal, silently, in release. It blocks adopting CBQ as an N-supporting input format until the round-trip is fixed. (Reported from the Bismark Rust port; happy to share a self-contained repro crate.)

[noamteyssier]

@BenjaminDEMAILLE

Thank you for noticing this and submitting a PR. I patched this with the latest 0.9.3 release of binseq.

To be clear - this bug only affected the streaming reader path for CBQ which is not the standard path and was built originally for development or single-thread usage and is generally not mature.

If you want to really make use of BINSEQ's parallel features I suggest you use the BinseqReader or cbq::MmapReader structs with the process_parallel API as is shown in the examples