Push qowsvpqmoukq #61
@@ -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<u32>,
|
||||
|
||||
/// 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<f64>,
|
||||
|
||||
/// 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)
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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" }
|
||||
|
||||
@@ -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,18 +83,19 @@ 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() {
|
||||
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,16 +162,17 @@ 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() {
|
||||
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
|
||||
};
|
||||
|
||||
|
||||
@@ -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<dyn KmerFilter>], row: &[u32], n_genomes: usize) -> bool {
|
||||
filters.iter().all(|f| f.passes(row, n_genomes))
|
||||
pub fn passes_all(
|
||||
filters: &[Box<dyn KmerFilter>],
|
||||
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::<u32>() >= 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::<u32>() <= 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
|
||||
}
|
||||
}
|
||||
|
||||
@@ -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());
|
||||
}
|
||||
}
|
||||
|
||||
@@ -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<Vec<u8>> {
|
||||
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<Vec<u8>> {
|
||||
(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<Vec<u8>> {
|
||||
SuperKmerIter::new(rope, K, 1, 0.0)
|
||||
.flat_map(|rsk| {
|
||||
rsk.superkmer()
|
||||
.iter_canonical_kmers()
|
||||
.map(|km| km.to_ascii())
|
||||
.collect::<Vec<_>>()
|
||||
})
|
||||
.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<u8> = (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<u8> = 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::<Vec<u8>>::new());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn empty_input_emits_nothing() {
|
||||
setup();
|
||||
let out = run_nofilter(b"", K);
|
||||
assert_eq!(out, Vec::<Vec<u8>>::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<Vec<u8>> = 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<Vec<u8>> = 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<RoutableSuperKmer> = SuperKmerIter::new(&rope, K, 1, 0.0).collect();
|
||||
assert!(!results.is_empty());
|
||||
}
|
||||
}
|
||||
#[path = "tests/iter.rs"]
|
||||
mod tests;
|
||||
|
||||
@@ -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;
|
||||
@@ -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;
|
||||
|
||||
|
||||
@@ -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<Vec<u8>> {
|
||||
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<Vec<u8>> {
|
||||
(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<Vec<u8>> {
|
||||
SuperKmerIter::new(rope, K, 1, 0.0)
|
||||
.flat_map(|rsk| {
|
||||
rsk.superkmer()
|
||||
.iter_canonical_kmers()
|
||||
.map(|km| km.to_ascii())
|
||||
.collect::<Vec<_>>()
|
||||
})
|
||||
.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<u8> = (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<u8> = 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::<Vec<u8>>::new());
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn empty_input_emits_nothing() {
|
||||
setup();
|
||||
let out = run_nofilter(b"", K);
|
||||
assert_eq!(out, Vec::<Vec<u8>>::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<Vec<u8>> = 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<Vec<u8>> = 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<RoutableSuperKmer> = SuperKmerIter::new(&rope, K, 1, 0.0).collect();
|
||||
assert!(!results.is_empty());
|
||||
}
|
||||
@@ -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:?}");
|
||||
}
|
||||
}
|
||||
Reference in New Issue
Block a user