diff --git a/src/obikmer/src/cmd/filter.rs b/src/obikmer/src/cmd/filter.rs index 993f80a..ae3f45f 100644 --- a/src/obikmer/src/cmd/filter.rs +++ b/src/obikmer/src/cmd/filter.rs @@ -2,14 +2,14 @@ use std::path::PathBuf; use clap::Args; use obikindex::{KmerIndex, MergeMode}; -use obikpartitionner::filter::{MaxTotalCount, MinTotalCount}; +use obikpartitionner::filter::{MaxTotalCount, MinComplexity, MinTotalCount}; use obisys::Reporter; use tracing::info; use super::predicate::FilterArgs as KmerFilterArgs; #[derive(Args)] -pub struct FilterArgs { +pub struct FilterCmdArgs { /// Source index directory pub source: PathBuf, @@ -28,6 +28,18 @@ pub struct FilterArgs { #[arg(long)] pub max_total_count: Option, + /// Minimum normalized entropy (complexity) to keep a k-mer — same metric + /// as `obikmer index`'s --theta, applied here to k-mers already committed + /// to the source index (reconstructed from unitigs.bin). K-mers scoring + /// below this are removed. + #[arg(long)] + pub min_complexity: Option, + + /// Maximum sub-word size for the complexity computation (see `obikmer + /// index`'s --level-max). Only used when --min-complexity is set. + #[arg(long, default_value_t = 6)] + pub complexity_level_max: usize, + /// Output as presence/absence instead of counts #[arg(long)] pub presence: bool, @@ -37,7 +49,7 @@ pub struct FilterArgs { pub force: bool, } -pub fn run(args: FilterArgs) { +pub fn run(args: FilterCmdArgs) { let src = KmerIndex::open(&args.source).unwrap_or_else(|e| { eprintln!("error opening source index: {e}"); std::process::exit(1); @@ -62,6 +74,9 @@ pub fn run(args: FilterArgs) { if let Some(v) = args.max_total_count { filters.push(Box::new(MaxTotalCount { total: v })); } + if let Some(theta) = args.min_complexity { + filters.push(Box::new(MinComplexity { level_max: args.complexity_level_max, theta })); + } let mut rep = Reporter::new(); KmerIndex::rebuild(&args.output, &src, &filters, mode, args.force, &mut rep) diff --git a/src/obikmer/src/main.rs b/src/obikmer/src/main.rs index a0b270b..cb701e6 100644 --- a/src/obikmer/src/main.rs +++ b/src/obikmer/src/main.rs @@ -21,7 +21,7 @@ enum Commands { /// Merge multiple built indexes into one Merge(cmd::merge::MergeArgs), /// Apply row-level selection (σ) to an index: retain only k-mers matching the predicates - Filter(cmd::filter::FilterArgs), + Filter(cmd::filter::FilterCmdArgs), /// Project and/or aggregate genome columns into a new or in-place index Select(cmd::select::SelectArgs), /// Query an index with sequences and annotate matches diff --git a/src/obikpartitionner/Cargo.toml b/src/obikpartitionner/Cargo.toml index 1c42cab..a0c810e 100644 --- a/src/obikpartitionner/Cargo.toml +++ b/src/obikpartitionner/Cargo.toml @@ -6,7 +6,6 @@ edition = "2024" [dev-dependencies] tempfile = "3" obikseq = { path = "../obikseq", features = ["test-utils"] } -obiskbuilder = { path = "../obiskbuilder" } obiread = { path = "../obiread" } obikrope = { path = "../obikrope" } @@ -14,6 +13,7 @@ obikrope = { path = "../obikrope" } niffler = "3.0.0" remove_dir_all = "0.8" obikseq = { path = "../obikseq" } +obiskbuilder = { path = "../obiskbuilder" } obiskio = { path = "../obiskio" } obidebruinj = { path = "../obidebruinj" } obilayeredmap = { path = "../obilayeredmap" } diff --git a/src/obikpartitionner/src/dump_layer.rs b/src/obikpartitionner/src/dump_layer.rs index e1c8226..00a8507 100644 --- a/src/obikpartitionner/src/dump_layer.rs +++ b/src/obikpartitionner/src/dump_layer.rs @@ -62,7 +62,7 @@ impl KmerPartition { for (kmer, _, _) in reader.iter_indexed_canonical_kmers() { if let Some(slot) = mphf.find(kmer) { let row = mat.row(slot); - if passes_all(filters, &row, n_genomes) { + if passes_all(filters, kmer, &row, n_genomes) { cont = cb(kmer, row); if !cont { break; } } @@ -75,7 +75,7 @@ impl KmerPartition { for (kmer, _, _) in reader.iter_indexed_canonical_kmers() { if let Some(slot) = mphf.find(kmer) { let row: Box<[u32]> = mat.row(slot).iter().map(|&b| b as u32).collect(); - if passes_all(filters, &row, n_genomes) { + if passes_all(filters, kmer, &row, n_genomes) { cont = cb(kmer, row); if !cont { break; } } @@ -83,16 +83,17 @@ impl KmerPartition { } cont } else { - // No data matrix: implicit presence — all values are 1. - // The filter result is identical for every kmer, so evaluate once. + // No data matrix: implicit presence — all values are 1. `row` + // is identical for every kmer, but a filter can still depend + // on the kmer's own sequence (e.g. MinComplexity), so this + // cannot be evaluated once for the whole layer — filters must + // still be tested per kmer. let all_present: Box<[u32]> = vec![1u32; n_genomes].into(); let mut cont = true; - if passes_all(filters, &all_present, n_genomes) { - for (kmer, _, _) in reader.iter_indexed_canonical_kmers() { - if mphf.find(kmer).is_some() { - cont = cb(kmer, all_present.clone()); - if !cont { break; } - } + for (kmer, _, _) in reader.iter_indexed_canonical_kmers() { + if mphf.find(kmer).is_some() && passes_all(filters, kmer, &all_present, n_genomes) { + cont = cb(kmer, all_present.clone()); + if !cont { break; } } } cont @@ -140,7 +141,7 @@ impl KmerPartition { for (kmer, _, _) in reader.iter_indexed_canonical_kmers() { if let Some(slot) = mphf.find(kmer) { let row = mat.row(slot); - if passes_all(filters, &row, n_genomes) { + if passes_all(filters, kmer, &row, n_genomes) { cont = cb(part, layer, kmer, row); if !cont { break; } } @@ -153,7 +154,7 @@ impl KmerPartition { for (kmer, _, _) in reader.iter_indexed_canonical_kmers() { if let Some(slot) = mphf.find(kmer) { let row: Box<[u32]> = mat.row(slot).iter().map(|&b| b as u32).collect(); - if passes_all(filters, &row, n_genomes) { + if passes_all(filters, kmer, &row, n_genomes) { cont = cb(part, layer, kmer, row); if !cont { break; } } @@ -161,14 +162,15 @@ impl KmerPartition { } cont } else { + // Same as iter_partition_kmers: row is constant but a filter + // may still depend on the kmer's own sequence, so this must + // be tested per kmer, not once for the whole layer. let all_present: Box<[u32]> = vec![1u32; n_genomes].into(); let mut cont = true; - if passes_all(filters, &all_present, n_genomes) { - for (kmer, _, _) in reader.iter_indexed_canonical_kmers() { - if mphf.find(kmer).is_some() { - cont = cb(part, layer, kmer, all_present.clone()); - if !cont { break; } - } + for (kmer, _, _) in reader.iter_indexed_canonical_kmers() { + if mphf.find(kmer).is_some() && passes_all(filters, kmer, &all_present, n_genomes) { + cont = cb(part, layer, kmer, all_present.clone()); + if !cont { break; } } } cont diff --git a/src/obikpartitionner/src/filter.rs b/src/obikpartitionner/src/filter.rs index 00f3b03..db69aba 100644 --- a/src/obikpartitionner/src/filter.rs +++ b/src/obikpartitionner/src/filter.rs @@ -1,17 +1,24 @@ use obicompactvec::FilterMask; +use obikseq::CanonicalKmer; -/// Trait for kmer row filters. +/// Trait for kmer filters. /// +/// `kmer` is the k-mer's own canonical sequence, reconstructed from the +/// source index's `unitigs.bin` (always present — see `rebuild_layer.rs`); /// `row` contains raw per-genome counts (or 0/1 for presence/absence data). -/// `n_genomes` equals `row.len()`. +/// `n_genomes` equals `row.len()`. Most filters only need `row` — `kmer` is +/// there for filters that reason about the k-mer's sequence itself (e.g. +/// [`MinComplexity`]). pub trait KmerFilter: Send + Sync { - fn passes(&self, row: &[u32], n_genomes: usize) -> bool; + fn passes(&self, kmer: CanonicalKmer, row: &[u32], n_genomes: usize) -> bool; /// Express this filter as a [`FilterMask`] column-operation expression. /// /// Returns `Some(expr)` if the filter can be evaluated solely from matrix /// column aggregates (no per-kmer row scan needed). Returns `None` if the - /// filter requires row-level inspection. + /// filter requires row-level inspection — always the case for a filter + /// that needs the k-mer's sequence, since a `FilterMask` only expresses + /// per-genome column aggregates, never per-slot sequence data. /// /// `threshold` semantics in the returned mask use `>= threshold`, matching /// [`obicompactvec::MatrixGroupOps`]. Implementations must add 1 to any @@ -23,8 +30,13 @@ pub trait KmerFilter: Send + Sync { /// True when `row` passes every filter in `filters`. /// Returns `true` if `filters` is empty. -pub fn passes_all(filters: &[Box], row: &[u32], n_genomes: usize) -> bool { - filters.iter().all(|f| f.passes(row, n_genomes)) +pub fn passes_all( + filters: &[Box], + kmer: CanonicalKmer, + row: &[u32], + n_genomes: usize, +) -> bool { + filters.iter().all(|f| f.passes(kmer, row, n_genomes)) } // ── Quorum filters ───────────────────────────────────────────────────────────── @@ -40,7 +52,7 @@ pub struct MinGenomeFraction { } impl KmerFilter for MinGenomeFraction { - fn passes(&self, row: &[u32], n_genomes: usize) -> bool { + fn passes(&self, _kmer: CanonicalKmer, row: &[u32], n_genomes: usize) -> bool { let p = present_count(row, self.threshold); p as f64 / n_genomes as f64 >= self.frac } @@ -63,7 +75,7 @@ pub struct MaxGenomeFraction { } impl KmerFilter for MaxGenomeFraction { - fn passes(&self, row: &[u32], n_genomes: usize) -> bool { + fn passes(&self, _kmer: CanonicalKmer, row: &[u32], n_genomes: usize) -> bool { let p = present_count(row, self.threshold); p as f64 / n_genomes as f64 <= self.frac } @@ -86,7 +98,7 @@ pub struct MinGenomeCount { } impl KmerFilter for MinGenomeCount { - fn passes(&self, row: &[u32], _n_genomes: usize) -> bool { + fn passes(&self, _kmer: CanonicalKmer, row: &[u32], _n_genomes: usize) -> bool { present_count(row, self.threshold) >= self.count } @@ -107,7 +119,7 @@ pub struct MaxGenomeCount { } impl KmerFilter for MaxGenomeCount { - fn passes(&self, row: &[u32], _n_genomes: usize) -> bool { + fn passes(&self, _kmer: CanonicalKmer, row: &[u32], _n_genomes: usize) -> bool { present_count(row, self.threshold) <= self.count } @@ -129,7 +141,7 @@ pub struct MinTotalCount { } impl KmerFilter for MinTotalCount { - fn passes(&self, row: &[u32], _n_genomes: usize) -> bool { + fn passes(&self, _kmer: CanonicalKmer, row: &[u32], _n_genomes: usize) -> bool { row.iter().sum::() >= self.total } @@ -147,7 +159,7 @@ pub struct MaxTotalCount { } impl KmerFilter for MaxTotalCount { - fn passes(&self, row: &[u32], _n_genomes: usize) -> bool { + fn passes(&self, _kmer: CanonicalKmer, row: &[u32], _n_genomes: usize) -> bool { row.iter().sum::() <= self.total } @@ -212,7 +224,7 @@ impl GroupQuorumFilter { } impl KmerFilter for GroupQuorumFilter { - fn passes(&self, row: &[u32], _n_genomes: usize) -> bool { + fn passes(&self, _kmer: CanonicalKmer, row: &[u32], _n_genomes: usize) -> bool { if !self.ingroup_idx.is_empty() { let n = self.ingroup_idx.iter() .filter(|&&i| row.get(i).copied().unwrap_or(0) > self.threshold) @@ -260,3 +272,27 @@ impl KmerFilter for GroupQuorumFilter { Some(FilterMask::And(parts)) } } + +// ── Complexity filter (post-hoc, sequence-based) ────────────────────────────── + +/// Reject k-mers with normalized entropy below `theta` — the same complexity +/// metric `obikmer index`'s `--theta`/`--level-max` apply *during* superkmer +/// construction (see [`obiskbuilder::KmerEntropy`]), applied here after the +/// fact, to k-mers already committed to a built index. +/// +/// Unlike every other filter in this module, this one needs the k-mer's own +/// sequence, not its per-genome row — `column_mask_expr` is never overridden +/// (stays `None`), so this filter always forces the row-level scan path in +/// `rebuild_layer.rs` (which reconstructs the sequence from `unitigs.bin` +/// regardless, so no extra I/O beyond what filtering already requires). +pub struct MinComplexity { + pub level_max: usize, + pub theta: f64, +} + +impl KmerFilter for MinComplexity { + fn passes(&self, kmer: CanonicalKmer, _row: &[u32], _n_genomes: usize) -> bool { + use obiskbuilder::KmerEntropy; + kmer.entropy(self.level_max) >= self.theta + } +} diff --git a/src/obikpartitionner/src/rebuild_layer.rs b/src/obikpartitionner/src/rebuild_layer.rs index b8893ef..e608c49 100644 --- a/src/obikpartitionner/src/rebuild_layer.rs +++ b/src/obikpartitionner/src/rebuild_layer.rs @@ -126,7 +126,7 @@ fn iter_src_kmers_masked( Some(m) => m.get(slot), None => { let row = src_data.fill_row_by_slot(slot, n_genomes); - filters.iter().all(|f| f.passes(&row, n_genomes)) + filters.iter().all(|f| f.passes(kmer, &row, n_genomes)) } }; if passes { cb(kmer); } @@ -165,7 +165,7 @@ fn iter_src_layers( cb(kmer, row.into_boxed_slice()); } else { let row = src_data.fill_row_by_slot(slot, n_genomes); - if filters.iter().all(|f| f.passes(&row, n_genomes)) { + if filters.iter().all(|f| f.passes(kmer, &row, n_genomes)) { cb(kmer, row.into_boxed_slice()); } } diff --git a/src/obiskbuilder/src/iter.rs b/src/obiskbuilder/src/iter.rs index 0839b97..8a33123 100644 --- a/src/obiskbuilder/src/iter.rs +++ b/src/obiskbuilder/src/iter.rs @@ -149,157 +149,5 @@ impl Iterator for SuperKmerIter<'_> { // ── tests ───────────────────────────────────────────────────────────────────── #[cfg(test)] -mod tests { - use super::*; - use obikrope::Rope; - - fn setup() { - obikseq::params::set_k(K); - obikseq::params::set_m(5); - } - - fn make_rope(data: &[u8]) -> Rope { - let mut r = Rope::new(None); - r.push(data.to_vec()); - r - } - - fn run_nofilter(data: &[u8], k: usize) -> Vec> { - let rope = make_rope(data); - SuperKmerIter::new(&rope, k, 1, 0.0) - .map(|rsk| rsk.superkmer().to_ascii()) - .collect() - } - - // k=11, m=5 — valeurs minimales du projet (k ∈ [11,31]) - const K: usize = 11; - - /// Collect the set of canonical k-mers from a raw ASCII sequence (no NUL). - fn direct_canonical_kmers(seq: &[u8]) -> std::collections::HashSet> { - (0..seq.len().saturating_sub(K - 1)) - .map(|i| obikseq::SuperKmer::from_ascii(&seq[i..i + K]).to_ascii()) - .collect() - } - - /// Collect the set of canonical k-mers emitted by SuperKmerIter over a rope. - fn iter_canonical_kmers(rope: &Rope) -> std::collections::HashSet> { - SuperKmerIter::new(rope, K, 1, 0.0) - .flat_map(|rsk| { - rsk.superkmer() - .iter_canonical_kmers() - .map(|km| km.to_ascii()) - .collect::>() - }) - .collect() - } - - #[test] - fn coverage_single_segment() { - setup(); - let seq = b"ACGTACGTACGTACGTACGT"; - let rope = make_rope(&[seq.as_ref(), b"\x00"].concat()); - let direct = direct_canonical_kmers(seq); - let from_iter = iter_canonical_kmers(&rope); - let missing: Vec<_> = direct.difference(&from_iter).collect(); - assert!( - missing.is_empty(), - "k-mers perdus dans segment unique : {missing:?}" - ); - } - - #[test] - fn coverage_two_segments() { - setup(); - let seg1 = b"ACGTACGTACGTACGTACGT"; - let seg2 = b"TGCATGCATGCATGCATGCA"; - let rope = make_rope(&[seg1.as_ref(), b"\x00", seg2.as_ref(), b"\x00"].concat()); - let mut direct = direct_canonical_kmers(seg1); - direct.extend(direct_canonical_kmers(seg2)); - let from_iter = iter_canonical_kmers(&rope); - let missing: Vec<_> = direct.difference(&from_iter).collect(); - assert!( - missing.is_empty(), - "k-mers perdus dans deux segments : {missing:?}" - ); - } - - #[test] - fn coverage_minimizer_boundary() { - setup(); - // sequence assez longue pour forcer plusieurs changements de minimiseur - let seq: Vec = (0..80).map(|i| b"ACGT"[i % 4]).collect(); - let rope = make_rope(&[seq.as_slice(), b"\x00"].concat()); - let direct = direct_canonical_kmers(&seq); - let from_iter = iter_canonical_kmers(&rope); - let missing: Vec<_> = direct.difference(&from_iter).collect(); - assert!( - missing.is_empty(), - "k-mers perdus à la frontière de minimiseur : {missing:?}" - ); - } - - #[test] - fn single_segment_one_superkmer() { - setup(); - let out = run_nofilter(b"ACGTACGTACGTACGTACGT\x00", K); - assert!(!out.is_empty()); - let total: Vec = out.into_iter().flatten().collect(); - assert!(total.len() >= K); - } - - #[test] - fn segment_shorter_than_k_emits_nothing() { - setup(); - let out = run_nofilter(b"ACGTACGT\x00", K); - assert_eq!(out, Vec::>::new()); - } - - #[test] - fn empty_input_emits_nothing() { - setup(); - let out = run_nofilter(b"", K); - assert_eq!(out, Vec::>::new()); - } - - #[test] - fn two_segments_both_emitted() { - setup(); - let out = run_nofilter(b"ACGTACGTACGTACGT\x00TGCATGCATGCATGCA\x00", K); - assert!(!out.is_empty()); - } - - #[test] - fn low_complexity_kmer_is_rejected() { - setup(); - let out_pass = run_nofilter(b"AAAAAAAAAAAACGTACGTACGT\x00", K); - assert!(!out_pass.is_empty()); - - let rope = make_rope(b"AAAAAAAAAAAAAAAAAAAA\x00"); - let out_reject: Vec> = SuperKmerIter::new(&rope, K, 6, 0.9) - .map(|rsk| rsk.superkmer().to_ascii()) - .collect(); - assert!(out_reject.is_empty()); - } - - #[test] - fn multi_slice_rope() { - setup(); - let data = b"ACGTACGTACGTACGTACGT\x00"; - let mid = data.len() / 2; - let mut rope = Rope::new(None); - rope.push(data[..mid].to_vec()); - rope.push(data[mid..].to_vec()); - let out: Vec> = SuperKmerIter::new(&rope, K, 1, 0.0) - .map(|rsk| rsk.superkmer().to_ascii()) - .collect(); - assert!(!out.is_empty()); - } - - #[test] - fn yields_minimizer_value() { - setup(); - let rope = make_rope(b"ACGTACGTACGTACGTACGT\x00"); - let results: Vec = SuperKmerIter::new(&rope, K, 1, 0.0).collect(); - assert!(!results.is_empty()); - } -} +#[path = "tests/iter.rs"] +mod tests; diff --git a/src/obiskbuilder/src/kmer_entropy.rs b/src/obiskbuilder/src/kmer_entropy.rs new file mode 100644 index 0000000..f8fcd7f --- /dev/null +++ b/src/obiskbuilder/src/kmer_entropy.rs @@ -0,0 +1,49 @@ +//! Normalized entropy of an isolated, already-built k-mer. +//! +//! [`SuperKmerIter`](crate::SuperKmerIter) uses [`RollingStat`] to reject +//! low-complexity k-mers *during* superkmer construction, incrementally, over +//! a streaming window. [`KmerEntropy`] exposes the same metric for a single +//! k-mer taken in isolation (e.g. one already reconstructed from an index's +//! `unitigs.bin`, with no surrounding sequence) — built on the identical +//! `RollingStat` code path, not a re-derived formula, so a `theta` threshold +//! chosen for index-build-time filtering means the same thing when applied +//! after the fact (e.g. `obikmer filter`). + +use obikseq::CanonicalKmer; + +use crate::rolling_stat::RollingStat; + +/// Extension trait: compute the normalized entropy of a single canonical +/// k-mer, independent of any surrounding sequence. +pub trait KmerEntropy { + /// Normalized entropy across sub-word orders `1..=level_max` (the + /// minimum is taken across orders, same as + /// [`RollingStat::normalized_entropy`]). Lower means less complex; + /// `theta` in `index`/`filter` rejects k-mers with a score `< theta`. + /// + /// # Panics + /// + /// Requires both `obikseq::params::k()` *and* `params::m()` to already be + /// set (via `set_k`/`set_m`), even though the minimizer length plays no + /// conceptual role in an entropy score: `RollingStat` is a shared, + /// general-purpose struct (also used for minimizer selection during + /// superkmer decomposition) and unconditionally sizes an `m`-dependent + /// mask in its constructor. Callers must call `set_m` with a valid value + /// (e.g. the index's own stored `minimizer_size`) even when only the + /// entropy score is needed. + fn entropy(&self, level_max: usize) -> f64; +} + +impl KmerEntropy for CanonicalKmer { + fn entropy(&self, level_max: usize) -> f64 { + let mut stat = RollingStat::new(level_max); + for &base in self.to_ascii().iter() { + stat.push(base); + } + stat.normalized_entropy().unwrap_or(1.0) + } +} + +#[cfg(test)] +#[path = "tests/kmer_entropy.rs"] +mod tests; diff --git a/src/obiskbuilder/src/lib.rs b/src/obiskbuilder/src/lib.rs index 6dc2de1..61d8cf1 100644 --- a/src/obiskbuilder/src/lib.rs +++ b/src/obiskbuilder/src/lib.rs @@ -6,6 +6,7 @@ #![deny(missing_docs)] pub mod iter; +pub mod kmer_entropy; pub mod stream_iter; mod scratch; @@ -14,6 +15,7 @@ pub(crate) mod entropy_table; pub(crate) mod rolling_stat; pub use iter::SuperKmerIter; +pub use kmer_entropy::KmerEntropy; pub use scratch::SuperKmerScratch; pub use stream_iter::SuperKmerStreamIter; diff --git a/src/obiskbuilder/src/tests/iter.rs b/src/obiskbuilder/src/tests/iter.rs new file mode 100644 index 0000000..4e03036 --- /dev/null +++ b/src/obiskbuilder/src/tests/iter.rs @@ -0,0 +1,152 @@ +use super::*; +use obikrope::Rope; + +fn setup() { + obikseq::params::set_k(K); + obikseq::params::set_m(5); +} + +fn make_rope(data: &[u8]) -> Rope { + let mut r = Rope::new(None); + r.push(data.to_vec()); + r +} + +fn run_nofilter(data: &[u8], k: usize) -> Vec> { + let rope = make_rope(data); + SuperKmerIter::new(&rope, k, 1, 0.0) + .map(|rsk| rsk.superkmer().to_ascii()) + .collect() +} + +// k=11, m=5 — valeurs minimales du projet (k ∈ [11,31]) +const K: usize = 11; + +/// Collect the set of canonical k-mers from a raw ASCII sequence (no NUL). +fn direct_canonical_kmers(seq: &[u8]) -> std::collections::HashSet> { + (0..seq.len().saturating_sub(K - 1)) + .map(|i| obikseq::SuperKmer::from_ascii(&seq[i..i + K]).to_ascii()) + .collect() +} + +/// Collect the set of canonical k-mers emitted by SuperKmerIter over a rope. +fn iter_canonical_kmers(rope: &Rope) -> std::collections::HashSet> { + SuperKmerIter::new(rope, K, 1, 0.0) + .flat_map(|rsk| { + rsk.superkmer() + .iter_canonical_kmers() + .map(|km| km.to_ascii()) + .collect::>() + }) + .collect() +} + +#[test] +fn coverage_single_segment() { + setup(); + let seq = b"ACGTACGTACGTACGTACGT"; + let rope = make_rope(&[seq.as_ref(), b"\x00"].concat()); + let direct = direct_canonical_kmers(seq); + let from_iter = iter_canonical_kmers(&rope); + let missing: Vec<_> = direct.difference(&from_iter).collect(); + assert!( + missing.is_empty(), + "k-mers perdus dans segment unique : {missing:?}" + ); +} + +#[test] +fn coverage_two_segments() { + setup(); + let seg1 = b"ACGTACGTACGTACGTACGT"; + let seg2 = b"TGCATGCATGCATGCATGCA"; + let rope = make_rope(&[seg1.as_ref(), b"\x00", seg2.as_ref(), b"\x00"].concat()); + let mut direct = direct_canonical_kmers(seg1); + direct.extend(direct_canonical_kmers(seg2)); + let from_iter = iter_canonical_kmers(&rope); + let missing: Vec<_> = direct.difference(&from_iter).collect(); + assert!( + missing.is_empty(), + "k-mers perdus dans deux segments : {missing:?}" + ); +} + +#[test] +fn coverage_minimizer_boundary() { + setup(); + // sequence assez longue pour forcer plusieurs changements de minimiseur + let seq: Vec = (0..80).map(|i| b"ACGT"[i % 4]).collect(); + let rope = make_rope(&[seq.as_slice(), b"\x00"].concat()); + let direct = direct_canonical_kmers(&seq); + let from_iter = iter_canonical_kmers(&rope); + let missing: Vec<_> = direct.difference(&from_iter).collect(); + assert!( + missing.is_empty(), + "k-mers perdus à la frontière de minimiseur : {missing:?}" + ); +} + +#[test] +fn single_segment_one_superkmer() { + setup(); + let out = run_nofilter(b"ACGTACGTACGTACGTACGT\x00", K); + assert!(!out.is_empty()); + let total: Vec = out.into_iter().flatten().collect(); + assert!(total.len() >= K); +} + +#[test] +fn segment_shorter_than_k_emits_nothing() { + setup(); + let out = run_nofilter(b"ACGTACGT\x00", K); + assert_eq!(out, Vec::>::new()); +} + +#[test] +fn empty_input_emits_nothing() { + setup(); + let out = run_nofilter(b"", K); + assert_eq!(out, Vec::>::new()); +} + +#[test] +fn two_segments_both_emitted() { + setup(); + let out = run_nofilter(b"ACGTACGTACGTACGT\x00TGCATGCATGCATGCA\x00", K); + assert!(!out.is_empty()); +} + +#[test] +fn low_complexity_kmer_is_rejected() { + setup(); + let out_pass = run_nofilter(b"AAAAAAAAAAAACGTACGTACGT\x00", K); + assert!(!out_pass.is_empty()); + + let rope = make_rope(b"AAAAAAAAAAAAAAAAAAAA\x00"); + let out_reject: Vec> = SuperKmerIter::new(&rope, K, 6, 0.9) + .map(|rsk| rsk.superkmer().to_ascii()) + .collect(); + assert!(out_reject.is_empty()); +} + +#[test] +fn multi_slice_rope() { + setup(); + let data = b"ACGTACGTACGTACGTACGT\x00"; + let mid = data.len() / 2; + let mut rope = Rope::new(None); + rope.push(data[..mid].to_vec()); + rope.push(data[mid..].to_vec()); + let out: Vec> = SuperKmerIter::new(&rope, K, 1, 0.0) + .map(|rsk| rsk.superkmer().to_ascii()) + .collect(); + assert!(!out.is_empty()); +} + +#[test] +fn yields_minimizer_value() { + setup(); + let rope = make_rope(b"ACGTACGTACGTACGTACGT\x00"); + let results: Vec = SuperKmerIter::new(&rope, K, 1, 0.0).collect(); + assert!(!results.is_empty()); +} diff --git a/src/obiskbuilder/src/tests/kmer_entropy.rs b/src/obiskbuilder/src/tests/kmer_entropy.rs new file mode 100644 index 0000000..1456c3a --- /dev/null +++ b/src/obiskbuilder/src/tests/kmer_entropy.rs @@ -0,0 +1,54 @@ +use super::*; +use obikseq::Sequence; +use obikseq::kmer::Kmer; + +const K: usize = 21; +const M: usize = 9; // RollingStat also tracks the minimizer window internally +const LEVEL_MAX: usize = 6; + +fn kmer_from_ascii(seq: &[u8]) -> CanonicalKmer { + obikseq::set_k(K); + obikseq::set_m(M); + Kmer::from_ascii(seq).expect("valid k-mer sequence").canonical() +} + +#[test] +fn homopolymer_scores_lower_than_diverse_sequence() { + let homopolymer = kmer_from_ascii(b"AAAAAAAAAAAAAAAAAAAAA"); // 21 bases + let diverse = kmer_from_ascii(b"CATTAGCGTACCTGATCAGGT"); // 21 bases, same as used elsewhere in this workspace's tests + + let e_homopolymer = homopolymer.entropy(LEVEL_MAX); + let e_diverse = diverse.entropy(LEVEL_MAX); + + assert!( + e_homopolymer < e_diverse, + "homopolymer ({e_homopolymer}) should score lower than a diverse sequence ({e_diverse})" + ); + // A pure homopolymer is the most degenerate case representable — its + // score should sit near the bottom of the range, not just "somewhat lower". + assert!(e_homopolymer < 0.3, "homopolymer entropy unexpectedly high: {e_homopolymer}"); +} + +#[test] +fn entropy_is_deterministic_for_the_same_kmer() { + let a = kmer_from_ascii(b"CATTAGCGTACCTGATCAGGT"); + let b = kmer_from_ascii(b"CATTAGCGTACCTGATCAGGT"); + assert_eq!(a.entropy(LEVEL_MAX), b.entropy(LEVEL_MAX)); +} + +#[test] +fn entropy_is_within_zero_one_range() { + let mut repeat = "AT".repeat(K / 2 + 1); + repeat.truncate(K); + + for seq in [ + "AAAAAAAAAAAAAAAAAAAAA".to_string(), + repeat, + "CATTAGCGTACCTGATCAGGT".to_string(), + ] { + assert_eq!(seq.len(), K, "test sequence must be exactly K bases: {seq:?}"); + let kmer = kmer_from_ascii(seq.as_bytes()); + let e = kmer.entropy(LEVEL_MAX); + assert!((0.0..=1.0).contains(&e), "entropy {e} out of [0,1] for {seq:?}"); + } +}