From fb31a35c761d93d12c0dfce79778d1b538ad31ba Mon Sep 17 00:00:00 2001 From: Eric Coissac Date: Sat, 22 Aug 2026 12:48:49 +0200 Subject: [PATCH] Add unitig streaming iterators and estimate CLI subcommand Introduce `iter_unitigs` methods across the index cache, content layer, MphfLayer, and typed layer to stream whole reconstructed sequences directly from underlying storage without decomposing into k-mers. Add a corresponding low-level streaming iterator in obiskio for thread-safe, lazy reads from memory-mapped files. Include a new CLI subcommand to compute and display approximate false-positive rates based on provided indexing parameters. --- src/obikidxcache/src/index_cache.rs | 10 +++++++ src/obikindex/src/layer/content_layer.rs | 12 ++++++++ src/obikindex/src/layer/mod.rs | 2 +- src/obikindex/src/layer/mphf_layer.rs | 28 +++++++++++++++++ src/obikindex/src/layer/typed_layer.rs | 6 ++++ src/obikmer2/src/cmd/estimate/mod.rs | 38 ++++++++++++++++++++++++ src/obikmer2/src/cmd/mod.rs | 1 + src/obikmer2/src/main.rs | 3 ++ src/obiskio/src/unitig_index/reader.rs | 25 ++++++++++++++++ 9 files changed, 124 insertions(+), 1 deletion(-) create mode 100644 src/obikmer2/src/cmd/estimate/mod.rs diff --git a/src/obikidxcache/src/index_cache.rs b/src/obikidxcache/src/index_cache.rs index bd6a517a..f4ed54d0 100644 --- a/src/obikidxcache/src/index_cache.rs +++ b/src/obikidxcache/src/index_cache.rs @@ -119,6 +119,16 @@ impl<'a> IndexCache<'a> { }) } + /// Iterate over the unitigs of every cached layer, chained — the + /// partition/index level of `KmerLayer::iter_unitigs`'s layer/partition/ + /// index chain, covering whichever scope this cache was opened with + /// (one partition, several, or the whole index — see [`new`](Self::new)). + /// No extra opening: every layer here is already open, this only + /// streams what [`iter`](Self::iter) already gives out. + pub fn iter_unitigs(&self) -> impl Iterator + '_ { + self.iter().flat_map(|l| l.iter_unitigs()) + } + #[inline] pub fn hash(&self, partition: usize, layer: usize, kmer: CanonicalKmer) -> Option { Some(self.get_layer(partition, layer)?.hash(kmer)) diff --git a/src/obikindex/src/layer/content_layer.rs b/src/obikindex/src/layer/content_layer.rs index f38784d3..10e7dbc1 100644 --- a/src/obikindex/src/layer/content_layer.rs +++ b/src/obikindex/src/layer/content_layer.rs @@ -311,4 +311,16 @@ impl KmerLayer { } } } + + /// Iterate over the layer's unitigs (whole sequences, not decomposed + /// into k-mers) — the "layer level" step of the layer/partition/index + /// chain of unitig iterators; see `KmerPartition::iter_unitigs`/ + /// `KmerIndex::iter_unitigs`. + pub fn iter_unitigs(&self) -> crate::layer::mphf_layer::UnitigIter { + match self { + KmerLayer::Count { layer, .. } => layer.iter_unitigs(), + KmerLayer::Presence { layer, .. } => layer.iter_unitigs(), + KmerLayer::Empty { .. } => panic!("Layer::iter_unitigs() called on an Empty layer"), + } + } } diff --git a/src/obikindex/src/layer/mod.rs b/src/obikindex/src/layer/mod.rs index fb0b1334..c06950d1 100644 --- a/src/obikindex/src/layer/mod.rs +++ b/src/obikindex/src/layer/mod.rs @@ -9,5 +9,5 @@ pub(crate) mod utils; pub use crate::index::error::{OKIError, OKIResult}; pub use content_layer::KmerLayer; pub use meta::IndexMode; -pub use mphf_layer::{EvidenceKind, KmerBatchIter, KmerIter, MphfLayer, MphfOnly}; +pub use mphf_layer::{EvidenceKind, KmerBatchIter, KmerIter, MphfLayer, MphfOnly, UnitigIter}; pub use typed_layer::{HasLayerContent, HasStorageKind, Hit, LayerContent, LayerData, TypedLayer}; diff --git a/src/obikindex/src/layer/mphf_layer.rs b/src/obikindex/src/layer/mphf_layer.rs index 7721ebac..0ce61819 100644 --- a/src/obikindex/src/layer/mphf_layer.rs +++ b/src/obikindex/src/layer/mphf_layer.rs @@ -341,6 +341,16 @@ impl MphfLayer { (base, batch) }) } + + /// Iterate over the layer's unitigs (whole reconstructed sequences, not + /// decomposed into k-mers), in the physical order of `unitigs.bin`. Same + /// `Send + 'static` shape as [`iter_kmers`](Self::iter_kmers), via + /// `UnitigFileReader::iter_unitigs_owned`. + pub fn iter_unitigs(&self) -> UnitigIter { + UnitigIter { + inner: Box::new(self.unitigs.iter_unitigs_owned()), + } + } } // ── Iterator types ──────────────────────────────────────────────────────────── @@ -365,6 +375,24 @@ impl Iterator for KmerIter { } } +/// Iterator over the unitigs stored in a layer, whole (not decomposed into +/// k-mers). Produced by [`MphfLayer::iter_unitigs`]. Same ownership shape as +/// [`KmerIter`] — owns an `Arc` clone, `Send + 'static`, +/// streamed from disk. +pub struct UnitigIter { + inner: Box + Send>, +} + +impl Iterator for UnitigIter { + type Item = obikseq::Unitig; + + /// Return the next unitig in iteration order (its `chunk_id` dropped — + /// use [`Iterator::enumerate`] if the index within the layer is needed). + fn next(&mut self) -> Option { + self.inner.next().map(|(_, unitig)| unitig) + } +} + /// Iterator over batches of canonical kmers stored in a layer. /// /// Produced by [`MphfLayer::iter_kmers_batch`]. Each call to [`next`](Self::next) diff --git a/src/obikindex/src/layer/typed_layer.rs b/src/obikindex/src/layer/typed_layer.rs index 99cf5a81..5df93beb 100644 --- a/src/obikindex/src/layer/typed_layer.rs +++ b/src/obikindex/src/layer/typed_layer.rs @@ -223,6 +223,12 @@ impl TypedLayer { self.mphf.enumerate_kmers_batch(n) } + /// Iterate over the layer's unitigs (whole sequences, not decomposed + /// into k-mers). + pub fn iter_unitigs(&self) -> crate::layer::mphf_layer::UnitigIter { + self.mphf.iter_unitigs() + } + pub fn unitig_writer(out_dir: &Path) -> OKIResult { MphfLayer::unitig_writer(out_dir) } diff --git a/src/obikmer2/src/cmd/estimate/mod.rs b/src/obikmer2/src/cmd/estimate/mod.rs new file mode 100644 index 00000000..0d8fe982 --- /dev/null +++ b/src/obikmer2/src/cmd/estimate/mod.rs @@ -0,0 +1,38 @@ +use clap::Args; + +use super::index::resolve_approx_params; + +#[derive(Args)] +pub struct EstimateArgs { + /// k-mer size used for querying (same as --kmer-size in index) + #[arg(short = 'k', long, default_value_t = 31)] + pub kmer_size: usize, + + /// Findere z parameter: number of consecutive k-mers that must all match. + /// Effective indexed k-mer size is kmer_size - z + 1. + #[arg(short = 'z', long, default_value = None)] + pub findere_z: Option, + + /// Fingerprint bits per slot (b). FP per z-window = 1/2^(b·z). + #[arg(long, default_value = None)] + pub evidence_bits: Option, + + /// Target false-positive rate per z-window (e.g. 0.01). + #[arg(long, default_value = None)] + pub fp: Option, +} + +pub fn run(args: EstimateArgs) { + let (z, b, fp_window) = resolve_approx_params(args.findere_z, args.evidence_bits, args.fp); + + let k_query = args.kmer_size; + let k_index = k_query.saturating_sub(z as usize - 1); + let fp_kmer = 1.0_f64 / 2_f64.powi(b as i32); + + println!("{:<22} {}", "k (query):", k_query); + println!("{:<22} {}", "k (indexed):", k_index); + println!("{:<22} {}", "z:", z); + println!("{:<22} {}", "evidence bits (b):", b); + println!("{:<22} {:.3e} (1/2^{})", "FP per k-mer:", fp_kmer, b); + println!("{:<22} {:.3e} (1/2^{})", "FP per z-window:", fp_window, b as u32 * z as u32); +} diff --git a/src/obikmer2/src/cmd/mod.rs b/src/obikmer2/src/cmd/mod.rs index 307cf14d..8419a2d8 100644 --- a/src/obikmer2/src/cmd/mod.rs +++ b/src/obikmer2/src/cmd/mod.rs @@ -1,3 +1,4 @@ +pub mod estimate; pub mod index; pub mod merge; pub mod superkmer; diff --git a/src/obikmer2/src/main.rs b/src/obikmer2/src/main.rs index 501f9901..7cc80788 100644 --- a/src/obikmer2/src/main.rs +++ b/src/obikmer2/src/main.rs @@ -19,6 +19,8 @@ enum Commands { Superkmer(cmd::superkmer::SuperkmerArgs), /// Merge multiple genome indexes into one Merge(cmd::merge::MergeArgs), + /// Estimate approximate-evidence false-positive rates for given parameters + Estimate(cmd::estimate::EstimateArgs), } fn main() { @@ -34,5 +36,6 @@ fn main() { Commands::Index(args) => cmd::index::run(args), Commands::Superkmer(args) => cmd::superkmer::run(args), Commands::Merge(args) => cmd::merge::run(args), + Commands::Estimate(args) => cmd::estimate::run(args), } } diff --git a/src/obiskio/src/unitig_index/reader.rs b/src/obiskio/src/unitig_index/reader.rs index 54226d53..d8dd461d 100644 --- a/src/obiskio/src/unitig_index/reader.rs +++ b/src/obiskio/src/unitig_index/reader.rs @@ -208,6 +208,31 @@ impl UnitigFileReader { self.iter_chunks_sequential() } + /// Same streamed sequence as [`iter_unitigs`](Self::iter_unitigs), but + /// owning a clone of `self` instead of borrowing it — `Send + 'static`, + /// so it can be handed to a threaded consumer without first collecting + /// the layer's unitigs into memory. Same shape as + /// [`iter_indexed_canonical_kmers_owned`](Self::iter_indexed_canonical_kmers_owned), + /// stopping short of decomposing into k-mers. + pub fn iter_unitigs_owned( + self: &Arc, + ) -> impl Iterator + Send + 'static { + let this = Arc::clone(self); + let k = this.k; + let n = this.n_unitigs; + let mut offset = 0usize; + (0..n).map(move |chunk_id| { + let mmap = &*this.mmap; + let seql = mmap[offset] as usize + k; + let byte_len = (seql + 3) / 4; + let bytes = mmap[offset + 1..offset + 1 + byte_len] + .to_vec() + .into_boxed_slice(); + offset += 1 + byte_len; + (chunk_id, Unitig::new((seql % 4) as u8, bytes)) + }) + } + pub fn iter_kmers(&self) -> impl Iterator + '_ { self.iter_chunks_sequential() .flat_map(|(_, u)| u.into_kmers())