Push zunrplorkwkt #70
@@ -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<Item = obikseq::Unitig> + '_ {
|
||||
self.iter().flat_map(|l| l.iter_unitigs())
|
||||
}
|
||||
|
||||
#[inline]
|
||||
pub fn hash(&self, partition: usize, layer: usize, kmer: CanonicalKmer) -> Option<usize> {
|
||||
Some(self.get_layer(partition, layer)?.hash(kmer))
|
||||
|
||||
@@ -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"),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@@ -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};
|
||||
|
||||
@@ -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<UnitigFileReader>` clone, `Send + 'static`,
|
||||
/// streamed from disk.
|
||||
pub struct UnitigIter {
|
||||
inner: Box<dyn Iterator<Item = (usize, obikseq::Unitig)> + 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::Item> {
|
||||
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)
|
||||
|
||||
@@ -223,6 +223,12 @@ impl<D: LayerData> TypedLayer<D> {
|
||||
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<UnitigFileWriter> {
|
||||
MphfLayer::unitig_writer(out_dir)
|
||||
}
|
||||
|
||||
@@ -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<u8>,
|
||||
|
||||
/// Fingerprint bits per slot (b). FP per z-window = 1/2^(b·z).
|
||||
#[arg(long, default_value = None)]
|
||||
pub evidence_bits: Option<u8>,
|
||||
|
||||
/// Target false-positive rate per z-window (e.g. 0.01).
|
||||
#[arg(long, default_value = None)]
|
||||
pub fp: Option<f64>,
|
||||
}
|
||||
|
||||
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);
|
||||
}
|
||||
@@ -1,3 +1,4 @@
|
||||
pub mod estimate;
|
||||
pub mod index;
|
||||
pub mod merge;
|
||||
pub mod superkmer;
|
||||
|
||||
@@ -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),
|
||||
}
|
||||
}
|
||||
|
||||
@@ -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<Self>,
|
||||
) -> impl Iterator<Item = (usize, Unitig)> + 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<Item = Kmer> + '_ {
|
||||
self.iter_chunks_sequential()
|
||||
.flat_map(|(_, u)| u.into_kmers())
|
||||
|
||||
Reference in New Issue
Block a user