Compare commits
1
Commits
4f6d442688
...
v1.1.44
| Author | SHA1 | Date | |
|---|---|---|---|
|
|
c95c47155e |
+2
-10
@@ -28,16 +28,8 @@ jobs:
|
||||
key: ${{ runner.os }}-cargo-v2-${{ hashFiles('src/Cargo.lock') }}
|
||||
restore-keys: ${{ runner.os }}-cargo-v2-
|
||||
|
||||
# Both `obikmer` and `obikindex` default to the `numa` feature
|
||||
# (hwloc-based topology detection + CPU pinning), which is only useful
|
||||
# on bare-metal multi-socket indexing hosts. Under this runner's
|
||||
# container/cgroup setup it deadlocks at startup — confirmed live
|
||||
# (2026-08-11): the same test binary hangs indefinitely with `numa` on
|
||||
# and passes instantly, repeatedly, with it off, on the same
|
||||
# container. Disable it for CI; it has nothing to do with test
|
||||
# correctness.
|
||||
- name: Build
|
||||
run: cargo build --release --no-default-features
|
||||
run: cargo build --release
|
||||
|
||||
- name: Test
|
||||
run: cargo test --release --no-default-features
|
||||
run: cargo test --release
|
||||
|
||||
@@ -24,3 +24,5 @@ benchmark/reference_dist
|
||||
benchmark/obikmer_dist
|
||||
benchmark/specific_index_count
|
||||
benchmark/specific_index_presence
|
||||
TNT
|
||||
phyg
|
||||
|
||||
@@ -473,6 +473,79 @@ matrix derived from the state graph, `c_ts`/`c_tv`/`c_gl`/`c_ctx` as
|
||||
tunable parameters) is directly usable there — no new tooling required
|
||||
before testing it.
|
||||
|
||||
### Experiment: TNT run on real data (2026-08-11)
|
||||
|
||||
**Status:** validated at genus/family/order scale within Bacteria; not
|
||||
informative across domains with this character type. Exploratory only —
|
||||
run entirely outside the repo (`/tmp/tnt_run`, TNT installed locally under
|
||||
`TNT/`), no new Rust code. Kept here as the record of what was learned.
|
||||
|
||||
**Pipeline.** `obikmer distance --snp` emits one IUPAC-coded pseudo-alignment
|
||||
row per genome (`snp_pseudo_alignment`, one column per family with
|
||||
`family_size() >= 2`). A small external Python script decodes IUPAC back to
|
||||
the 16-state bitmask, applies the set-edit-distance formula above, and
|
||||
emits a complete TNT script (`xread` matrix + `smatrix` step-matrix + `hold`/
|
||||
`mult` search). No new Rust code was needed for this pass.
|
||||
|
||||
**Calibration ("quick option").** Rather than modifying
|
||||
`write_raw_snp_distance_csv` to emit raw counts, `c_ctx` was approximated as
|
||||
the unweighted mean of the pairwise `snp/(snp+shared)` ratios already
|
||||
present in `rawsnp.csv`, restricted to the relevant taxon subset. Used
|
||||
values: `mu = gamma = 10`, `c_ctx = 18` (ratio `c_ctx/mu = 1.8`, consistent
|
||||
with the observed pairwise ratios for genomes within Enterobacteriaceae/
|
||||
Eubacteria, which cluster around 1.5-2 and barely move as the taxon set
|
||||
widens — the context signal is stable, not sensitive to which subset is
|
||||
chosen). The more principled route (raw counts, weighted `p_hat`, Ts/Tv
|
||||
split) remains a follow-up, not yet done.
|
||||
|
||||
**Run 1 — Enterobacteriaceae (11 taxa: 4 *E. coli*, 4 *Salmonella enterica*,
|
||||
3 *Klebsiella pneumoniae*).** All three genera recovered as monophyletic.
|
||||
Initial read of the exported (unrooted, TNT/Nexus `[&U]`) tree as showing a
|
||||
genus-arrangement disagreement with known systematics (*Escherichieae*:
|
||||
*Escherichia*+*Salmonella* sister vs. more distant *Klebsielleae*) was
|
||||
**wrong** — diagnosed via `force = (taxa);` monophyly-constraint test
|
||||
(identical score constrained vs. free ⟹ no real disagreement). Root cause of
|
||||
the misreading: with exactly 3 clades and no outgroup, an unrooted tree has
|
||||
only **one possible topology** (a single trifurcation) — there is no
|
||||
internal arrangement to get right or wrong. This run cannot test
|
||||
inter-genus relationships at all; it can only test intra-genus monophyly
|
||||
(which held).
|
||||
|
||||
**Run 2 — Eubacteria (18 taxa: run 1 + *Acidobacterium capsulatum*,
|
||||
*Opitutus terrae*, *Bacillus subtilis*, *Shouchella clausii*, *Wolbachia*
|
||||
endosymbiont, *Proteus mirabilis*, *Yersinia ruckeri*).** The
|
||||
Enterobacteriaceae substructure from run 1 is reproduced identically, now
|
||||
correctly rooted by real outgroups, resolving the tribal arrangement left
|
||||
undetermined in run 1: *Escherichia*+*Salmonella* sister, *Klebsiella* more
|
||||
distant — matching known systematics. Outgroup placement: `(Bacillus,
|
||||
Shouchella)` sister pair (Firmicutes/*Bacillales*) splits from all
|
||||
Proteobacteria — a correct phylum-level split; `Wolbachia`
|
||||
(Alphaproteobacteria) splits from the Gammaproteobacteria block
|
||||
(`Proteus`, `Yersinia`, Enterobacteriaceae) — a correct class-level split.
|
||||
Lower confidence: the fine nested order `(Proteus, (Yersinia,
|
||||
Enterobacteriaceae))` and the relative position of `Acidobacterium` vs.
|
||||
`Opitutus` — plausible, not independently verified against current
|
||||
Enterobacterales family-level literature.
|
||||
|
||||
**Run 3 — full domain set (20 taxa: run 2 + *Candidozyma auris* [yeast] and
|
||||
*Saccharolobus islandicus* [archaeon]).** The bacterial clade from run 2 is
|
||||
reproduced **unchanged and intact** — a real robustness signal, the method
|
||||
does not fragment the ingroup when unrelated deep taxa are added. But the
|
||||
result carries **no information on Bacteria/Archaea/Eukarya relationships**:
|
||||
with exactly one eukaryote, one archaeon and one bacterial clade, the
|
||||
unrooted tree is again forced into the single 3-clade trifurcation from run
|
||||
1's caveat — there is no second representative of either outgroup domain to
|
||||
resolve internal arrangement, so nothing about their relative position can
|
||||
be read from the topology (a "ladder" ordering in the exported tree is
|
||||
serialization, not signal). Independently, `rawsnp.csv` shows *why* this
|
||||
character type cannot reach further: pairwise ratios involving
|
||||
*Candidozyma*/*Saccharolobus* are almost all `NA` (no central-position
|
||||
family shared at all) or saturated at `1.0` (every shared family differs) —
|
||||
central-position families require literal 31 bp context conservation, which
|
||||
simply does not survive domain-level divergence. Cross-domain placement
|
||||
would need conserved-marker characters (rRNA, ribosomal proteins), not this
|
||||
estimator.
|
||||
|
||||
## Heterozygosity, ploidy, and consensus-assembly inputs
|
||||
|
||||
A within-genome multiplicity signal (more than one of the 4 central forms
|
||||
|
||||
Generated
+2
-1
@@ -1711,11 +1711,12 @@ dependencies = [
|
||||
"serde_json",
|
||||
"tempfile",
|
||||
"tracing",
|
||||
"tracing-subscriber",
|
||||
]
|
||||
|
||||
[[package]]
|
||||
name = "obikmer"
|
||||
version = "1.1.43"
|
||||
version = "1.1.44"
|
||||
dependencies = [
|
||||
"clap",
|
||||
"csv",
|
||||
|
||||
@@ -1,6 +1,17 @@
|
||||
use super::*;
|
||||
use obikseq::{k, set_k, unitig::Unitig, Kmer};
|
||||
|
||||
// `obikseq::params` is process-wide (see obikseq/src/params.rs): tests in this
|
||||
// file don't all use the same `k` (`push_palindrome_single_node` needs an
|
||||
// even k=4 — no odd-length self-revcomp palindrome exists — while the rest
|
||||
// use k=5), and they run concurrently by default. Serialize the
|
||||
// set_k-through-use critical section across this file's tests so one test's
|
||||
// `k` can never be overwritten mid-flight by another.
|
||||
static K_LOCK: std::sync::Mutex<()> = std::sync::Mutex::new(());
|
||||
fn lock_k() -> std::sync::MutexGuard<'static, ()> {
|
||||
K_LOCK.lock().unwrap_or_else(|e| e.into_inner())
|
||||
}
|
||||
|
||||
// Build a graph from an ASCII sequence, inserting all canonical k-mers.
|
||||
fn graph_from_ascii(seq: &[u8]) -> GraphDeBruijn {
|
||||
let mut g = GraphDeBruijn::new();
|
||||
@@ -37,6 +48,7 @@ fn collect_unitigs(g: &GraphDeBruijn) -> Vec<Unitig> {
|
||||
#[test]
|
||||
fn push_deduplicates_revcomp() {
|
||||
let k = 5;
|
||||
let _guard = lock_k();
|
||||
set_k(k);
|
||||
let kmer = Kmer::from_ascii(b"ACGTA").unwrap();
|
||||
let mut g = GraphDeBruijn::new();
|
||||
@@ -49,6 +61,7 @@ fn push_deduplicates_revcomp() {
|
||||
fn push_palindrome_single_node() {
|
||||
// ACGT is its own revcomp
|
||||
let k = 4;
|
||||
let _guard = lock_k();
|
||||
set_k(k);
|
||||
let kmer = Kmer::from_ascii(b"ACGT").unwrap();
|
||||
assert_eq!(kmer, kmer.revcomp(), "test requires a palindrome");
|
||||
@@ -71,6 +84,7 @@ fn linear_chain_graph() -> (GraphDeBruijn, Vec<CanonicalKmer>) {
|
||||
#[test]
|
||||
fn degrees_linear_chain_node_count() {
|
||||
let k = 5;
|
||||
let _guard = lock_k();
|
||||
set_k(k);
|
||||
let (g, kmers) = linear_chain_graph();
|
||||
assert_eq!(g.len(), kmers.len());
|
||||
@@ -82,6 +96,7 @@ fn degrees_linear_chain_extensions() {
|
||||
// Note: start_iter must not be consumed standalone — its second pass only
|
||||
// finds true cycle nodes when interleaved with chain traversal (iter_unitig).
|
||||
let k = 5;
|
||||
let _guard = lock_k();
|
||||
set_k(k);
|
||||
let seq = b"AAAAGGGG";
|
||||
let g = graph_from_ascii(seq);
|
||||
@@ -118,6 +133,7 @@ fn kmers_from_unitigs(unitigs: &[Unitig]) -> Vec<CanonicalKmer> {
|
||||
fn unitig_roundtrip_linear() {
|
||||
// Non-repetitive sequence: all k-mers must be recovered across unitigs.
|
||||
let k = 5;
|
||||
let _guard = lock_k();
|
||||
set_k(k);
|
||||
let seq = b"ACCTGGCTA";
|
||||
let g = graph_from_ascii(seq);
|
||||
@@ -136,6 +152,7 @@ fn unitig_roundtrip_longer_sequence() {
|
||||
// Longer non-repetitive sequence with no repeated k-mer of length k.
|
||||
// ACGTGGCTATCGAC with k=5 → 10 distinct k-mers, one linear chain.
|
||||
let k = 5;
|
||||
let _guard = lock_k();
|
||||
set_k(k);
|
||||
let seq = b"ACGTGGCTATCGAC";
|
||||
let g = graph_from_ascii(seq);
|
||||
@@ -152,6 +169,7 @@ fn unitig_roundtrip_longer_sequence() {
|
||||
fn unitig_isolated_node() {
|
||||
// Single k-mer with no neighbours
|
||||
let k = 5;
|
||||
let _guard = lock_k();
|
||||
set_k(k);
|
||||
let kmer = Kmer::from_ascii(b"ACGTA").unwrap();
|
||||
let mut g = GraphDeBruijn::new();
|
||||
@@ -165,6 +183,7 @@ fn unitig_isolated_node() {
|
||||
#[test]
|
||||
fn unitig_two_isolated_nodes() {
|
||||
let k = 5;
|
||||
let _guard = lock_k();
|
||||
set_k(k);
|
||||
let mut g = GraphDeBruijn::new();
|
||||
// Two k-mers that share no (k-1)-overlap
|
||||
@@ -177,6 +196,7 @@ fn unitig_two_isolated_nodes() {
|
||||
#[test]
|
||||
fn unitig_two_truly_distinct_isolated_nodes() {
|
||||
let k = 5;
|
||||
let _guard = lock_k();
|
||||
set_k(k);
|
||||
let mut g = GraphDeBruijn::new();
|
||||
g.push(Kmer::from_ascii(b"AAAAC").unwrap().canonical());
|
||||
@@ -192,7 +212,8 @@ fn unitig_two_truly_distinct_isolated_nodes() {
|
||||
|
||||
#[test]
|
||||
fn no_kmer_lost_or_duplicated() {
|
||||
let k = 7;
|
||||
let k = 5;
|
||||
let _guard = lock_k();
|
||||
set_k(k);
|
||||
let seq = b"ACGTACGTACGTTTTTACGTACGT";
|
||||
let g = graph_from_ascii(seq);
|
||||
@@ -218,6 +239,7 @@ fn cycle_kmers_not_lost() {
|
||||
// start_iter first pass yields nothing (all nodes internal); second pass
|
||||
// picks up cycle entries. All 4 k-mers must appear in the unitigs.
|
||||
let k = 5;
|
||||
let _guard = lock_k();
|
||||
set_k(k);
|
||||
let seq = b"ACGTACGT";
|
||||
let g = graph_from_ascii(seq);
|
||||
@@ -240,6 +262,7 @@ fn branching_graph_no_kmer_lost_or_duplicated() {
|
||||
// Each "node" is a distinct 5-mer; edges share a 4-mer suffix/prefix.
|
||||
// We use long non-repetitive sequences and extract only the required kmers.
|
||||
let k: usize = 5;
|
||||
let _guard = lock_k();
|
||||
set_k(k);
|
||||
let mut g = GraphDeBruijn::new();
|
||||
|
||||
|
||||
@@ -24,6 +24,7 @@ hwlocality = { version = "1.0.0-alpha.11", features = ["vendored"], option
|
||||
[dev-dependencies]
|
||||
obiread = { path = "../obiread" }
|
||||
tempfile = "3"
|
||||
tracing-subscriber = { version = "0.3", features = ["fmt", "env-filter"] }
|
||||
|
||||
[features]
|
||||
default = ["numa"]
|
||||
|
||||
@@ -295,20 +295,27 @@ impl PartitionRunner {
|
||||
let pool = node.pool.clone();
|
||||
|
||||
s.spawn(move || {
|
||||
let tid = std::thread::current().id();
|
||||
debug!(?tid, "PartitionRunner worker: waiting on activation");
|
||||
if arx.recv().is_err() {
|
||||
debug!(?tid, "PartitionRunner worker: activation channel closed, exiting");
|
||||
return;
|
||||
}
|
||||
debug!(?tid, "PartitionRunner worker: activated");
|
||||
if !cpu_ids.is_empty() {
|
||||
pin_current_thread(cpu_ids);
|
||||
}
|
||||
for i in &prx {
|
||||
debug!(?tid, partition = i, "PartitionRunner worker: picked partition");
|
||||
let t = Instant::now();
|
||||
let r = match &pool {
|
||||
Some(p) => p.install(|| f(i)),
|
||||
None => f(i),
|
||||
};
|
||||
debug!(?tid, partition = i, "PartitionRunner worker: partition done");
|
||||
etx.send(WorkerEvent::Completed(i, r, t.elapsed())).ok();
|
||||
}
|
||||
debug!(?tid, "PartitionRunner worker: no more partitions, exiting");
|
||||
});
|
||||
}
|
||||
}
|
||||
@@ -319,13 +326,18 @@ impl PartitionRunner {
|
||||
// ── Controller ────────────────────────────────────────────────────
|
||||
let mut activation = NodeActivation::new(&activate_txs, &node_caps, max_workers);
|
||||
activation.activate_initial(INITIAL_DIVISOR, n_total);
|
||||
debug!(n_total, activated = activation.total(), "PartitionRunner controller: initial activation");
|
||||
|
||||
let mut cpu_sample = CpuSample::now();
|
||||
let mut io_sample = IoSample::now();
|
||||
let mut completed = 0usize;
|
||||
|
||||
while completed < n_total {
|
||||
let Ok(event) = event_rx.recv() else { break };
|
||||
debug!(completed, n_total, "PartitionRunner controller: waiting for an event");
|
||||
let Ok(event) = event_rx.recv() else {
|
||||
debug!("PartitionRunner controller: event channel closed, stopping");
|
||||
break;
|
||||
};
|
||||
match event {
|
||||
WorkerEvent::Completed(i, r, dur) => {
|
||||
match r {
|
||||
|
||||
@@ -976,7 +976,23 @@ mod tests {
|
||||
/// the same primitives `obikmer`'s `scatter` step uses (minus the
|
||||
/// multi-file `obipipeline` wrapper — a single sequence needs none of
|
||||
/// that): normalise -> build superkmers -> route -> write.
|
||||
/// `cargo test` doesn't install a `tracing` subscriber the way `obikmer`'s
|
||||
/// CLI does, so `debug!`/etc. are silent no-ops by default — including the
|
||||
/// `PartitionRunner` instrumentation that would matter most for
|
||||
/// re-diagnosing a hang here. `try_init` is idempotent across concurrently
|
||||
/// running tests (later calls just find a subscriber already installed).
|
||||
fn init_tracing() {
|
||||
let _ = tracing_subscriber::fmt()
|
||||
.with_env_filter(
|
||||
tracing_subscriber::EnvFilter::try_from_default_env()
|
||||
.unwrap_or_else(|_| tracing_subscriber::EnvFilter::new("info")),
|
||||
)
|
||||
.with_writer(std::io::stderr)
|
||||
.try_init();
|
||||
}
|
||||
|
||||
fn build_single_genome_index(dir: &Path, label: &str, seq: &[u8]) -> KmerIndex {
|
||||
init_tracing();
|
||||
let fasta_path = dir.join(format!("{label}.fasta"));
|
||||
let mut f = std::fs::File::create(&fasta_path).unwrap();
|
||||
writeln!(f, ">{label}").unwrap();
|
||||
|
||||
@@ -1,6 +1,6 @@
|
||||
[package]
|
||||
name = "obikmer"
|
||||
version = "1.1.43"
|
||||
version = "1.1.44"
|
||||
edition = "2024"
|
||||
|
||||
[[bin]]
|
||||
|
||||
+33
-17
@@ -7,12 +7,28 @@
|
||||
//! different value panics. This prevents silent divergence between the global
|
||||
//! parameter and the values used to build data structures.
|
||||
//!
|
||||
//! In test builds (`#[cfg(test)]`) the same public API is backed by
|
||||
//! `thread_local!` [`Cell`]s instead. Each test thread gets its own
|
||||
//! independent copies of `K` and `M`, so tests can use arbitrary values
|
||||
//! without coordinating with one another and without any reset mechanism.
|
||||
//! The `OnceLock` constraint is deliberately absent: test isolation is
|
||||
//! provided by thread locality, not by write-once semantics.
|
||||
//! In test builds (`#[cfg(test)]`) the same public API is backed by plain
|
||||
//! process-wide atomics instead, freely overwritable (no write-once
|
||||
//! constraint) so tests don't need a reset mechanism between runs.
|
||||
//!
|
||||
//! An earlier version of this module used `thread_local!` `Cell`s here,
|
||||
//! reasoning that "each test thread gets its own copy" gives isolation
|
||||
//! between tests using different `k`/`m` values. That assumption broke as
|
||||
//! soon as any code under test fanned work out to *other* threads it
|
||||
//! doesn't control — `PartitionRunner`'s pre-spawned workers, or a bare
|
||||
//! `rayon::par_iter()` — since a freshly spawned thread never inherits the
|
||||
//! calling test thread's thread-local state, silently reading back `k=0`
|
||||
//! there instead (surfaced as a `bitvec`/slice-indexing panic deep inside
|
||||
//! whatever used the bogus length). Process-wide atomics make `k()`/`m()`
|
||||
//! correct on *any* thread without every call site having to know to
|
||||
//! re-propagate them. Unset still silently reads back as `0` (same as the
|
||||
//! old `Cell` default) rather than panicking: several existing tests read
|
||||
//! `m()` without ever calling `set_m` themselves, relying on that default.
|
||||
//! The trade-off: tests that genuinely need different `k`/`m` values from
|
||||
//! other tests must not run concurrently with them in the same process (in
|
||||
//! practice: every test file in this workspace already uses one fixed
|
||||
//! `k`/`m` pair for all its own tests, so this doesn't currently cost
|
||||
//! anything).
|
||||
|
||||
// ── Production implementation ─────────────────────────────────────────────────
|
||||
|
||||
@@ -44,22 +60,22 @@ mod state {
|
||||
|
||||
// ── Test implementation ───────────────────────────────────────────────────────
|
||||
//
|
||||
// Each test thread owns its private K and M via thread_local!, so tests may
|
||||
// call set_k / set_m with any value without affecting other tests.
|
||||
// Process-wide, freely overwritable (no write-once constraint), visible from
|
||||
// any thread — including threads a test doesn't spawn itself (rayon workers,
|
||||
// PartitionRunner workers, ...). `0` (never explicitly set) is returned as-is,
|
||||
// same default as the old thread-local `Cell`.
|
||||
|
||||
#[cfg(any(test, feature = "test-utils"))]
|
||||
mod state {
|
||||
use std::cell::Cell;
|
||||
use std::sync::atomic::{AtomicUsize, Ordering};
|
||||
|
||||
thread_local! {
|
||||
static K: Cell<usize> = Cell::new(0);
|
||||
static M: Cell<usize> = Cell::new(0);
|
||||
}
|
||||
static K: AtomicUsize = AtomicUsize::new(0);
|
||||
static M: AtomicUsize = AtomicUsize::new(0);
|
||||
|
||||
pub fn set_k(k: usize) { K.with(|c| c.set(k)); }
|
||||
pub fn k() -> usize { K.with(|c| c.get()) }
|
||||
pub fn set_m(m: usize) { M.with(|c| c.set(m)); }
|
||||
pub fn m() -> usize { M.with(|c| c.get()) }
|
||||
pub fn set_k(k: usize) { K.store(k, Ordering::SeqCst); }
|
||||
pub fn k() -> usize { K.load(Ordering::SeqCst) }
|
||||
pub fn set_m(m: usize) { M.store(m, Ordering::SeqCst); }
|
||||
pub fn m() -> usize { M.load(Ordering::SeqCst) }
|
||||
}
|
||||
|
||||
// ── Public API (identical signature in both configurations) ───────────────────
|
||||
|
||||
@@ -102,9 +102,9 @@ fn roundtrip_single() {
|
||||
|
||||
#[test]
|
||||
fn roundtrip_all_lengths() {
|
||||
obikseq::params::set_k(11);
|
||||
setup();
|
||||
let bases: Vec<u8> = (0..300).map(|i| b"ACGT"[i % 4]).collect();
|
||||
for len in (11..=19).chain([255, 256, 257]) {
|
||||
for len in (TEST_K..=19).chain([255, 256, 257]) {
|
||||
let sk = make_sk(&bases[..len]);
|
||||
let mut buf = Vec::new();
|
||||
sk.write_to_binary(&mut buf).unwrap();
|
||||
|
||||
Reference in New Issue
Block a user