diff --git a/docmd/implementation/filtering.md b/docmd/implementation/filtering.md index ea6d4a2..b45938e 100644 --- a/docmd/implementation/filtering.md +++ b/docmd/implementation/filtering.md @@ -92,18 +92,48 @@ For each genome: | Flag | Applies to | Meaning | |------|-----------|---------| -| `--min-count N` | ingroup | k-mer present in at least N ingroup genomes | -| `--max-count N` | ingroup | k-mer present in at most N ingroup genomes | +| `--min-count N` | ingroup | k-mer present in at least N ingroup genomes (N may be negative, see below) | +| `--max-count N` | ingroup | k-mer present in at most N ingroup genomes (N may be negative, see below) | | `--min-frac F` | ingroup | k-mer present in at least fraction F of ingroup genomes | | `--max-frac F` | ingroup | k-mer present in at most fraction F of ingroup genomes | -| `--min-outgroup-count N` | outgroup | k-mer present in at least N outgroup genomes | -| `--max-outgroup-count N` | outgroup | k-mer present in at most N outgroup genomes | +| `--min-outgroup-count N` | outgroup | k-mer present in at least N outgroup genomes (N may be negative, see below) | +| `--max-outgroup-count N` | outgroup | k-mer present in at most N outgroup genomes (N may be negative, see below) | | `--min-outgroup-frac F` | outgroup | k-mer present in at least fraction F of outgroup genomes | | `--max-outgroup-frac F` | outgroup | k-mer present in at most fraction F of outgroup genomes | | `--min-total-count N` | all genomes | sum of per-genome counts ≥ N (`filter` only) | | `--max-total-count N` | all genomes | sum of per-genome counts ≤ N (`filter` only) | | `--presence-threshold N` | all | per-genome count > N to be considered "present" (default 0) | +### Negative counts — offset from group size + +The four integer count flags (`--min-count`, `--max-count`, `--min-outgroup-count`, +`--max-outgroup-count`) accept **negative** values, interpreted as an offset counted +down from the group size `n`, resolved at run time once `n` is known: + +| Value | Effective threshold | +|-------|---------------------| +| `N ≥ 0` | literal absolute count `N` | +| `-x` (x > 0) | `max(1, n − x)` — "all but x" | + +`-1` literally means *all but one*, `-2` *all but two*, and so on. This expresses +a quorum relative to the group size that a plain fraction cannot state exactly +(e.g. "present in every genome except at most one" is `n−1`, which is `0.9` for +`n = 10` but `0.857…` for `n = 7`). + +The threshold is **floored at 1**, never 0: the negative form always keeps +constraining the group. Without the floor, `--min-count -1` on a singleton +ingroup (`n = 1`) would resolve to `0` ("at least 0") and silently drop the +constraint; the floor makes it `1` ("present in that one genome") instead. + +To express a count of `0` (e.g. "absent from the ingroup"), use the literal `0`, +not a negative — `0` and `-0` are indistinguishable, so the offset form starts at +`-1`. + +> **Edge case** — on an *empty* group (`n = 0`, e.g. a predicate matching no +> genome), a negative count still resolves to `1`, an impossible constraint that +> rejects every k-mer. This is consistent with an empty group letting nothing +> through, but differs from the "no constraint" behaviour of the fraction flags. + **Conditional defaults** — the defaults for `--min-frac` and `--max-outgroup-count` depend on two conditions: whether the corresponding group was declared, **and** whether any quorum flag for that group was explicitly set. @@ -215,6 +245,17 @@ obikmer filter src --output dst \ --max-outgroup-count 0 ``` +Noise-tolerant core — keep k-mers present in *all but one* ingroup genome +(`-1` = `n−1`) and absent from *all but one* of the outgroup: + +```sh +obikmer filter src --output dst \ + --ingroup "genus=Betula" \ + --outgroup "*" \ + --min-count -1 \ + --max-outgroup-count -1 +``` + To dump only k-mers specific to *Betula nana*: ```sh diff --git a/docmd/theory/evolutionary_distances.md b/docmd/theory/evolutionary_distances.md index 3d6be69..6ff266b 100644 --- a/docmd/theory/evolutionary_distances.md +++ b/docmd/theory/evolutionary_distances.md @@ -7,12 +7,30 @@ for the implemented Jaccard/Mash distances). ## Motivation -Mash infers a mutation rate from a single scalar (Jaccard) via a model that -assumes independence across the k positions. For k odd, a substitution at the -**central base** of a k-mer, with the `2m` flanking bases (`m = (k-1)/2`) -otherwise conserved, is directly observable: it is a SNP. This gives access not -just to a rate but to the substitution's nature (transition/transversion), -enabling classical corrected distances. +**Primary intent: restrict the comparison to what is actually comparable.** +Mash's Jaccard is computed over the **union** of both genomes' k-mer content: +anything not identically shared is folded into a single undifferentiated +mass, whether the cause is a point substitution, a genuinely absent +homologous region (lineage-specific content, gene-family expansion, HGT, +genome-size asymmetry), or a diverged paralogous copy. The model then +back-infers a single mutation rate from that mass, silently attributing +non-homology to mutation. The central-SNP approach instead conditions every +comparison on local, positive evidence of homology: a locus only enters the +statistic if its `2m` flanking bases (`m = (k-1)/2`) are found intact in +*both* genomes — genuinely absent or non-homologous content is excluded from +the comparison entirely (neither numerator nor denominator), rather than +silently counted as divergence. This is a conditioning on comparability, not +just a richer summary statistic — see "Statistic and correspondence with +`shared`" below for how it plays out against genome-size asymmetry and +diverged gene families, and "Heterozygosity, ploidy, and consensus-assembly +inputs" for the corresponding paralogy/heterozygosity filter. + +**Secondary benefit: access to the substitution's nature.** Because the +central base of an odd-k window is directly observable once the flanks are +confirmed conserved, this also yields more than a rate — the +transition/transversion split — enabling classical corrected distances +(Jukes-Cantor, Kimura 2-parameter, LogDet) that a single Jaccard scalar +cannot support. ## Statistic and correspondence with `shared` @@ -21,24 +39,125 @@ the substitution to sit at the window's **center** makes exactly one window per SNP eligible — a 1:1 correspondence between SNP and center-neighbor k-mer pair, avoiding the ~k-fold overcount of an any-position neighbor search. -A locus with a fully conserved `k`-window (flanks **and** center) is precisely -an exact-shared k-mer — already produced by the existing `shared_kmers` matrix -(`--shared-kmers`, `BitPartials::partial_jaccard` / `CountPartials::partial_threshold_jaccard`). -A locus with conserved flanks but a substituted center is a "central SNP". Both -count each locus exactly once, in matching units: +A locus with a fully conserved `k`-window (flanks **and** center) is an +exact-shared k-mer at that locus; a locus with conserved flanks but a +substituted center is a "central SNP". Both count each locus exactly once, in +matching units: ``` p_hat[i,j] = SNP[i,j] / (SNP[i,j] + shared[i,j]) ``` -`p_hat` is `P(center substituted | 2m flanks conserved)`. Only the numerator -`SNP[i,j]` needs computing; the denominator term is already available. +`p_hat` is `P(center substituted | 2m flanks conserved)`. `shared[i,j]` here +is **not** the general-purpose `shared_kmers` matrix used by Jaccard/Mash +(`--shared-kmers`, `BitPartials::partial_jaccard` / +`CountPartials::partial_threshold_jaccard`) — that matrix counts raw k-mer +identity with no per-genome copy-number constraint, whereas `p_hat`'s +denominator applies the eligibility rule defined below (raw or +paralogy-filtered). Both `SNP` and `shared` are accumulated by the same +sweep, from the same per-locus candidate set (source k-mer + 3 variants), +under the same eligibility rule — see "Locus eligibility" below, and +"Heterozygosity, ploidy, and consensus-assembly inputs" for why the +copy-number constraint matters and what it costs. **Canonical invariance**: for odd k, the central position maps to itself under reverse-complement (`m -> k-1-m = m`, base complemented). A transition maps to a transition, a transversion to a transversion — the transition/transversion split is well-defined in canonical space. +## Locus eligibility: raw definition vs. paralogy filter + +For each k-mer `x` observed in genome A (source, one MPHF slot; the 3 +central-position variants generated as in the sweep below): check whether +A's locus (flanks fixed) is resolvable in genome B under one of the 4 +central forms. + +**Raw / no model.** The locus counts in the denominator iff at least one of +the 4 forms is present in B; it counts in the numerator iff the form found in +B differs from A's own. No constraint on A's or B's own copy number at this +locus. Open question, not resolved: what if **more than one** of the 4 forms +is present in B simultaneously (ambiguous target — count once arbitrarily, +count all, or drop)? The stringent filter below sidesteps the question by +construction rather than answering it. + +**Stringent / paralogy-aware.** The locus counts only if exactly one of the +4 forms is present in A **and** exactly one is present in B (`count == 1` at +that slot too, when a count index is available, to also exclude same-allele +duplicates that presence alone cannot see). This drops the raw definition's +ambiguous-B case automatically, at the cost of also dropping heterozygous +sites indiscriminately alongside true duplications (see "Heterozygosity, +ploidy, and consensus-assembly inputs" below). + +**Rejected: parsimony-based multiset pairing for multiplicity > 1.** Rather +than dropping ambiguous loci, pair identical alleles between A and B first +(0-mutation explanation preferred), then take `min(unmatched_A, unmatched_B)` +as inferred SNP pairs. Rejected on two grounds: (1) circularity — selecting +pairs by minimal apparent divergence, then measuring divergence on those same +pairs, deflates the estimate by construction, not a neutral heuristic; (2) +the discriminating signal is a single base among 4 possible values, and the +flanks are *already* guaranteed identical for every candidate by +construction (that is how the locus was selected) — no information remains +in a k-mer window to tell which copy in B truly corresponds to which copy in +A once multiplicity > 1 on either side. Any pairing rule invents a +correspondence the data cannot support. Multiplicity > 1 is treated as +non-identifiable, not as a puzzle to solve with a heuristic. + +## Heterozygosity, ploidy, and consensus-assembly inputs + +A within-genome multiplicity signal (more than one of the 4 central forms +present at a locus) is produced identically by two distinct causes: +paralogous duplication and diploid/polyploid heterozygosity. K-mer data alone +cannot distinguish them. The one real discriminator is sequencing depth +(heterozygous site: total depth of the present forms ~= the genome's +single-copy average; duplication: ~2x or more) — but that signal only exists +if genome "counts" are raw-read depth (FASTQ input), not occurrence counts in +an assembled FASTA, where per-locus depth is not preserved. + +**Magnitude is taxon- and mating-system-dependent, not universal.** +Heterozygosity density: mammals ~1 site / 1-1.5 kb (~0.1%); highly +outcrossing plants (maize, poplar) reported an order of magnitude higher +(~1%); self-fertilising plants (*Arabidopsis thaliana*) near zero — but with +a documented failure mode where segmental duplication masquerades as +"pseudo-heterozygosity"; fungi split between haploid vegetative stages +(non-issue) and dikaryotic Basidiomycetes, where two long-diverged haploid +nuclei coexist without fusing. The estimator's target use case (closely +related genomes, k=31) is exactly where the stringent filter above costs the +least for low-heterozygosity taxa and the most for outcrossing/dikaryotic +ones — no universal threshold; this is a scope caveat to document, not a +problem to solve generically. + +**Why assembled-consensus inputs don't make measured distances wrong.** +Phylogenetic inputs are near-universally assemblies, not raw reads, and +assemblers collapse heterozygous sites to one consensus allele per +position — effectively an arbitrary, largely uncorrelated-between-assemblies +choice at each het site. This does not inject unbounded noise: standard +population genetics gives `d_xy = d_a + (pi_A + pi_B)/2` — the expected +pairwise difference between a random allele of population A and a random +allele of population B equals the net (fixed) divergence `d_a` plus the +average of the two populations' own within-population diversity `pi`. +Consensus flattening realises exactly this random-allele draw, so the +measured genome-to-genome distance is a `d_xy`-like quantity, not `d_a` — +inflated by heterozygosity by a well-characterised additive term, not +distorted unpredictably. The term is negligible when `pi << d_xy` (the common +case for cross-species comparisons), and becomes material precisely in the +two cases already flagged above: very closely related genomes (this +estimator's explicit target) and highly heterozygous outcrossing organisms, +where `pi` and `d_xy` are the same order of magnitude. + +Caveat: this assumes the flattening is uncorrelated with the phylogenetic +signal — plausible for de novo assembly, not guaranteed for reference-guided +assembly biased toward one allele (e.g. the reference's) at each het site, +which would turn the noise term into a systematic bias toward the reference +lineage. Not evaluated here. + +**Forward-looking implication, not part of the current design.** The +multiplicity > 1 signal discarded by the stringent filter is a crude +per-genome proxy for `pi` (under low background paralogy). If a `pi_hat` per +genome were tallied alongside `SnpTally`, a `d_a` correction +(`p_hat - mean(pi_hat_i, pi_hat_j)/2`, roughly) could recover an estimate +closer to net divergence instead of `d_xy` — a possible extension, not +scoped here. + ## Sufficient statistic: 4x4 base-pair tally Tabulating the joint distribution of `(center_i, center_j)` over conserved-flank @@ -51,8 +170,9 @@ loci, per genome pair, is sufficient for every downstream correction: | Kimura 2-parameter | transition rate P, transversion rate Q | `d = 1/2 ln(1/(1-2P-Q)) + 1/4 ln(1/(1-2Q))` | | LogDet/paralinear | full 4x4 + base-composition margins | `d ~= -1/4 ln det(F)`, robust to non-stationary base composition | -JC/K2P need only the total and the transition/transversion split. LogDet needs -the full 4x4 (adds a diagonal tally of shared k-mers' central base — cheap). +JC/K2P need only the total and the transition/transversion split (the +diagonal collapses to a single "shared" total). LogDet needs the full 4x4, +already populated at no extra cost (see Step 1/2 below). Memory for the 4x4 tally: `n^2 * 16` counters. Trivial for the project's genome-scale use case (tens to hundreds of genomes); ~13 GB at n=10^4 — outside @@ -67,6 +187,16 @@ scope but worth flagging if n grows. (universal in practice), conserved flanks correlate with slow centers, so `p_hat` underestimates the genome-wide average rate — it specifically estimates the substitution rate of **conserved regions**. + Two distinct factors are at play here, not one: `P(centre of a given + window is a SNP) = p` exactly, **independent of k** — a direct restatement + of the raw per-site rate via the bijective window<->centre-position + correspondence (Statistic section above), not a k-dependent quantity. + `(1-p)^{2m}` is the *separate*, genuinely k-dependent ascertainment factor + (are the flanks also intact). The two multiply: + `P(usable window showing a central SNP) = p * (1-p)^{2m}` — e.g. at + p=1/31 (~3.2%), k=31: `p * (1-p)^30 ~= 0.0323 * 0.374 ~= 1.2%`, i.e. about + 1 window in 83, not 1 in 31 (which is only the centre-mutated fraction, + before requiring intact flanks). 2. **Bias toward isolated SNPs.** Two SNPs within k of each other disqualify each other's flanks. Hypervariable regions are invisible by construction. 3. **Indels are invisible.** A frameshift destroys k-mer matches in a block; @@ -179,10 +309,13 @@ k-mer set. Sequential per-partition sweep (Route D): reuse the already-mmap'd per-partition MPHF/evidence/presence structures for O(1) variant lookups, dedup via the fixed sweep-order rule (`q >= p`, plus an in-partition -tie-break), no scratch files, no graph materialisation. Denominator reused -from the existing `shared_kmers` matrix. Distances (p, JC, K2P, LogDet) as -finalisations of the resulting 4x4 tally, mirroring the `partial_* -> *_dist_matrix` -pattern used for Jaccard/Mash/Bray-Curtis/etc. +tie-break), no scratch files, no graph materialisation. Both the SNP +(off-diagonal) and shared (diagonal) counts are accumulated by this same +sweep, under the locus-eligibility rule chosen (raw or paralogy-filtered) — +not reused from the general-purpose `shared_kmers` matrix, whose raw-identity +definition does not apply the same copy-number constraint. Distances (p, JC, +K2P, LogDet) as finalisations of the resulting 4x4 tally, mirroring the +`partial_* -> *_dist_matrix` pattern used for Jaccard/Mash/Bray-Curtis/etc. ## Detailed implementation plan @@ -218,9 +351,13 @@ A `SnpTally` holding, per genome pair, the 4x4 joint count of central bases: `n * n * 4 * 4` `u64` (or a packed lower-triangular form since it is symmetric). Provide `merge(&mut self, other: &SnpTally)` for the thread-local reduce, and accessors yielding, per pair `(i,j)`: total off-diagonal (SNP), -transition count `P`, transversion count `Q`. The diagonal need not be stored -for JC/K2P (the denominator comes from `shared_kmers`); store it only if -LogDet is wanted, by also tallying shared k-mers' central base. +diagonal (shared, i.e. `p_hat`'s denominator minus SNP), transition count +`P`, transversion count `Q`. The diagonal is always populated — it is not an +optional LogDet-only extra, since `p_hat`'s denominator is no longer sourced +from the external `shared_kmers` matrix (see "Locus eligibility" and +"Statistic and correspondence with `shared`" above): the source k-mer's own +presence/count vector, already in hand when it is enumerated, supplies the +diagonal entry directly, at no extra lookup cost. ### Step 2 — the sweep (`obikindex`, new `snp.rs`) @@ -233,6 +370,11 @@ for p in 0..n_partitions: # OUTER — sequential enumerate distinct canonical k-mers of p (one per MPHF slot) with their presence/count vectors # column-major, as query stage 2 par_iter over these source k-mers: # INNER — rayon, thread-local tally + apply eligibility rule to the source's own vector (raw: none; # diagonal + stringent: exactly one of the 4 forms present in each genome) # gate + for i in genomes eligible with source base a: + for j in genomes eligible with source base a: + thread_tally[i,j][a,a] += 1 # diagonal — no extra lookup for each of the 3 central variants: q = partition_of(variant) if q < p: continue # dedup: forward targets only @@ -240,8 +382,9 @@ for p in 0..n_partitions: # OUTER — sequential slot = layers[q].find_slot(variant) # MphfLayer::find, mmap'd if hit: vb = variant presence/count vector - for i in genomes with source base a: - for j in genomes with variant base b: + apply eligibility rule to vb (as above) + for i in eligible genomes with source base a: + for j in eligible genomes with variant base b: thread_tally[i,j][a,b] += 1 merge thread-local tallies into global SnpTally ``` @@ -257,11 +400,231 @@ column-major scan already implemented in `dump_layer.rs` / iff its count at that slot is `>= presence_threshold` (trivially `>= 1` for presence indexes). +### Open problem (unresolved, session end — not yet fully convinced) + +The `q >= p` / tie-break dedup rule in Step 2's pseudocode above is **flawed** +for the stringent (paralogy-filtered) eligibility rule: it only ever brings +two family members into view at once (the source and one looked-up variant), +never all four simultaneously, and which subset gets compared depends on +partition sweep order. "Exactly one of the 4 forms present in genome A" is a +whole-family property and cannot be decided correctly from a sequence of +pairwise, order-dependent glimpses — the pseudocode above needs revision, not +just the eligibility gate bolted onto it as written. + +Direction discussed, **not yet settled**: + +1. **Every distinct source k-mer looks up all 3 variants unconditionally** + (drop the `q < p` skip entirely) so that every observed family member + independently gathers all 4 vectors (its own + whichever of the 3 + variants exist) at once — a whole-family, order-independent view, computed + redundantly once per observed member. Same total lookup order of + magnitude as already budgeted (`3 * N_distinct`), just organised + differently (no lookup actually skipped, versus the original rule which + skipped roughly half). +2. **Tie-break after gathering, not before**: only the member whose own + canonical encoding is the smallest *among the members actually observed* + (now known, since all were just looked up) writes to `SnpTally`; the + others silently discard their redundant computation. Deterministic, + order-independent — as a side effect this also removes the "outer loop + must stay sequential" constraint from the cost/parallelism discussion + above, since no step depends on partition processing order any more. +3. **Proposed optimisation**: precompute, once at index build time, a + compact global (not per-genome) annex per MPHF slot — the count of + *other* family members observed anywhere in the dataset (0-3). Slots with + count 0 (majority under low divergence and few genomes, but see the + scaling caveat below) need no cross-lookup at all: eligibility reduces to + a local `count == 1` check at that single slot, and only slots with count + >= 1 enter the 3-lookup sweep machinery above. Revised (see Step 2b + below): minorant status *is* stored alongside the count after all, on 3 + bits rather than 2 — it comes for free from the same lookups needed to + count siblings, and storing it lets the sweep discard non-minorant slots + without re-fetching anything. + +**Minorant/sibling-count relationship, worked out precisely.** "Minorant" is +a one-way implication from sibling count, not an equivalence: `0 siblings +=> minorant` (trivially — with no other observed member, the k-mer is by +definition the smallest of the observed set, itself alone), and its +contrapositive `not minorant => >= 1 sibling`. The converse does not hold: +being the minorant says nothing about sibling count — a minorant can have 0, +1, 2 or 3 siblings, all with larger encodings than itself. Consequence: this +confirms, as a logical necessity rather than a heuristic, that a 0-sibling +slot can always write its diagonal contribution with zero ambiguity and no +lookup (it is unconditionally its own minorant) — but it gives no shortcut +for the >= 1-sibling case, where minorant status still requires the actual +comparison of gathered encodings; sibling count alone never determines it. + +**When to compute the annex, and cache invalidation.** Sibling count is a +property of the whole set of columns (genomes/groups) currently in the +index, not of any single genome — it cannot be computed correctly at +mono-genome build time (a family may gain siblings, or its minorant may +change, once more genomes are merged in later). Computing it eagerly at +every `merge` would also waste work on intermediate merged states nobody +ever queries. Instead: compute it lazily, on first `distance` call against a +given index, and persist the result alongside that index for subsequent +calls — the same lazy-derived-cache pattern `PersistentBitMatrix` already +uses for `Columnar` -> `Packed`. This requires no explicit invalidation for +`merge` or `filter` (`obikindex/src/merge.rs`, `obikmer/src/cmd/filter.rs`): +both only ever write to a fresh `--output` directory, never mutate an input +index in place, so a re-merged/re-filtered index is simply a new state with +no annex yet. `select --in-place` (`select_layer.rs:139-235`) is the +exception: it aggregates genome columns into groups (Any/All/None/Sum/Min/ +Max) by mutating the existing index's files without changing its location. +It does not remove k-mer rows, but it can still change eligibility and +sibling counts derived from those rows (e.g. a `Sum` over several +single-copy genomes can read as multi-copy at the group level). Because it +mutates in place, **`select --in-place` must explicitly invalidate (delete +or mark stale) any cached sibling-count annex for that index** — the one +operation in the current pipeline where this doesn't happen for free. + +Not yet convinced this is the right shape, and Step 2's pseudocode above has +not been rewritten to match — flagged for the next pass rather than resolved +here. + +### Step 2b — sibling-count / minorant annex (consolidated plan) + +Scope: only the precursor annex (sibling count 0-3 per slot, minorant +decided on demand) — not the SNP tally itself, whose Step 2 sweep remains +unresolved above. This piece is simpler than the sweep, because it writes to +an independent per-slot value, not a shared cross-k-mer accumulator, so it +needs no dedup/ownership logic at all at this stage. + +1. **Primitive.** Reuse `central_canonical_neighbors()` from Step 0 + unchanged — the 3 canonicalised central-substitution variants of a k-mer. +2. **New annex type** (`obicompactvec`, alongside `bitmatrix.rs`): a 3-bit- + per-slot packed array, one per partition — same on-disk shape family as + `PersistentBitMatrix`'s `Packed` variant, but simpler (no per-genome + columns, a single derived read-only value per slot). 3 bits, not 2: + revised to also store minorant status alongside sibling count, since it + comes for free from the same lookups (point 3 below) — 5 real states + (not-minorant; minorant with 0/1/2/3 siblings) fit in 3 bits (8 states, + 3 unused). This lets the future SNP sweep discard a non-minorant slot + **instantly**, with no lookup at all, instead of having to regenerate and + look up its siblings just to rediscover it isn't the designated writer — + moving that cost into this one-time, cached pass instead of repeating it + on every future sweep. The otherwise-unreachable combination + "not-minorant + 0 siblings" (impossible: 0 siblings always implies + minorant, see below) doubles as a free **"not yet computed" sentinel** — + annex files for all partitions/layers can be pre-initialised to this + value before the computation pass runs, distinguishing genuinely-computed + 0-sibling slots from not-yet-processed ones with no extra storage. +3. **Computation pass** (`obikindex`, new `siblings.rs`): **one + `obipipeline` run per layer, iterated sequentially over the index's + layers** — settled after two false starts, worth recording both. + - *False start 1*: "fully parallel over every partition/slot at once, + no ordering at all". Correctness is fine with this (sibling count and + minorant are order-independent, unlike the old `q >= p` dedup they + replace), but it reintroduces, at a larger scale, exactly the + memory-blowup the original Step 2 sweep's sequential-outer-loop + constraint existed to prevent: scattering every source partition at + once multiplies the in-flight outgoing-query volume by the number of + partitions. + - *False start 2*: push the layer loop itself into the pipeline (source + = the index's layers, a first `Flat` stage expands each layer into + its k-mers). `obipipeline`'s scheduler already bounds memory on its + own — it dispatches every item through a **shared** worker pool at + each stage boundary (`scheduler.rs:217-372`, `dispatch()` into a + common `worker_tx` queue, any free worker picks up any pending item; + not "one worker owns a chunk end to end"), with a biased `Select` + that prioritises draining items already advanced in the chain over + admitting new source items (`scheduler.rs:271-282`: stage results + outrank the source, "vider le pipeline en priorité" / "dernier + recours" for new data) — so bounded channel `capacity` plus this + drain-first bias already caps in-flight work without any external + sequential discipline. Correct, but it means k-mers from several + layers can be completing concurrently, so the sink would need to + track several open per-layer annex-file writers at once — real, + avoidable complexity. + - **Settled design**: keep the layer loop external and sequential — + not for memory (the pipeline's own `capacity`/priority mechanism + already provides that, for free, regardless), but so each pipeline + run's sink targets exactly one layer's annex file, no concurrent + multi-writer bookkeeping. Per layer: source = that layer's distinct + k-mers; a `Flat` (1->N) stage generates the 3 central variants of a + k-mer, each tagged with its origin (local slot); a transform stage + routes each variant to its target partition (unchanged per-k-mer + minimiser); a transform stage performs the lookup (existence-only — + `find_slot` hit/miss, cheaper than the SNP sweep's full column + fetch); a final stage/sink folds each answer into its origin's + running state (below) and, once a layer's k-mers are all resolved, + flushes the completed array to that layer's annex file. Many small, + single-purpose stages on purpose, to let the scheduler interleave + them finely across many in-flight items — this deliberately does + **not** mirror how `obipipeline` is used elsewhere today: `query.rs`'s + `process_chunk` lumps parse+route+query+serialise into one closure + (`query.rs:325,743-758`), and `scatter.rs` only pipelines file- + reading/superkmer construction, routing partitions afterwards in a + plain sequential loop (`KmerPartition::write_batch`, + `partition.rs:140`) — both under-use the fine-grained scheduling the + mechanism offers, so they are not precedents to copy, only existing + (and arguably improvable, out of scope here) usages. Cross-partition + lookups (querying another layer's MPHF for a variant) remain + necessary as before — only the *output* side is kept single-layer. + - **Reconciliation**: processed at the granularity of one *answer batch + per destination partition*, not one source k-mer at a time — this is a + proper shuffle, not a per-k-mer wait. Each source partition `p` holds a + small array of running states `(minorant = true, siblings = 0)`, one + per local slot, initialised at scatter time and **persisting across + however many destination-partition batches answer it** (up to 3, one + per variant, not necessarily all from the same `q`). Every scattered + query carries an origin tag (source partition + local slot) so its + answer can be routed back. When target partition `q` returns its batch + (all answers for every query that named `q`, regardless of which source + k-mer or which source partition they came from), that batch is walked + once, locally, and each answer updates — via its origin tag — the + matching entry in *its* source partition's array: a miss changes + nothing; a hit does `siblings += 1`, and if the found sibling's own + encoding is smaller than the source's, `minorant = false`. A given + source k-mer's state is final only once every destination batch + concerning it has been folded in; its partition's array is flushed to + the persistent annex once complete. Commutative per entry, so the order + in which destination batches arrive and get folded in doesn't matter. + + **Open optimisation, not adopted yet — real tradeoff, not a free win.** + Since looking up sibling `y` from `x`'s visit already yields everything + needed to fill `y`'s own annex entry too, one visit per *family* could in + principle replace one visit per *observed family member* — cutting this + pass's cost roughly by the average family size instead of paying + `3 * N_distinct` regardless. But it means threads processing different + source k-mers can end up writing the *same* sibling's slot concurrently — + the fully independent, ownership-free parallelism of the plan above is + deliberately traded away for this gain. It stays safe only because the + computed value for a given slot is deterministic regardless of who + computes it, so redundant concurrent writes converge to the same + value — correct as long as each write is atomic, no locking needed — but + it is a real design complexity increase over "every member redoes its + own 3 lookups independently," not a strict improvement to adopt by + default. +4. **Trigger and caching** (`obikindex::KmerIndex`/`distance.rs`): compute + lazily on first `distance` call for an SNP-family metric against a given + index; check for an existing annex file first (mirrors + `PersistentBitMatrix::open()`'s auto-detect-and-fall-back, + `bitmatrix.rs:264-287`); if absent, run step 3 and persist; if present, + mmap and reuse. +5. **Invalidation.** `merge` and `filter` always write to a fresh `--output` + directory (`obikindex/src/merge.rs`, `obikmer/src/cmd/filter.rs`) so a + re-merged/re-filtered index simply has no annex yet — nothing to + invalidate. `select --in-place` (`select_layer.rs:139-235`) mutates + columns of an existing index without changing its location, which can + change sibling counts without removing rows — it must explicitly delete + any cached annex for that index as part of its in-place rewrite. +6. **Testing**: hand-built tiny indexes with known sibling counts (0-3); + order-independence (recompute twice on a static index, identical + result, given the fully-parallel no-ownership design); invalidation + (annex absent/correctly recomputed after `select --in-place`); once + Step 2's sweep is fixed, a regression check that sibling_count == 0 + slots are never looked up cross-partition during the sweep. + +Cost: `3 * N_distinct` existence-only lookups, computed once per index +state and amortised over every subsequent `distance` call that reuses the +cached annex — cheaper per-lookup than the sweep itself (hit/miss only, no +column fetch). + ### Step 3 — finalisation (`obikindex`) -From the global `SnpTally` + the existing `shared_kmers` matrix, derive n x n -distance matrices, each a pure function of the accumulated counts (same shape -as `jaccard_to_mash`): +From the global `SnpTally` alone (diagonal and off-diagonal both populated by +the sweep, see Step 1/2 — no dependency on the external `shared_kmers` +matrix), derive n x n distance matrices, each a pure function of the +accumulated counts (same shape as `jaccard_to_mash`): - `p_hat[i,j] = SNP / (SNP + shared)` - Jukes-Cantor, Kimura-2P (from `P`, `Q`), optionally LogDet (needs the diff --git a/src/obikmer/src/cmd/predicate.rs b/src/obikmer/src/cmd/predicate.rs index 47baab9..b9ae40f 100644 --- a/src/obikmer/src/cmd/predicate.rs +++ b/src/obikmer/src/cmd/predicate.rs @@ -151,12 +151,14 @@ pub struct FilterArgs { pub outgroup: Vec, /// Minimum number of ingroup genomes containing the k-mer - #[arg(long)] - pub min_count: Option, + /// (negative: offset from group size, e.g. -1 = all but one) + #[arg(long, allow_hyphen_values = true)] + pub min_count: Option, /// Maximum number of ingroup genomes containing the k-mer - #[arg(long)] - pub max_count: Option, + /// (negative: offset from group size, e.g. -1 = all but one) + #[arg(long, allow_hyphen_values = true)] + pub max_count: Option, /// Minimum fraction of ingroup genomes containing the k-mer [0.0–1.0] /// (default 1.0 when --ingroup is set, 0.0 otherwise) @@ -168,13 +170,15 @@ pub struct FilterArgs { pub max_frac: Option, /// Minimum number of outgroup genomes containing the k-mer - #[arg(long)] - pub min_outgroup_count: Option, + /// (negative: offset from outgroup size, e.g. -1 = all but one) + #[arg(long, allow_hyphen_values = true)] + pub min_outgroup_count: Option, /// Maximum number of outgroup genomes containing the k-mer - /// (default 0 when --outgroup is set, no constraint otherwise) - #[arg(long)] - pub max_outgroup_count: Option, + /// (default 0 when --outgroup is set, no constraint otherwise; + /// negative: offset from outgroup size, e.g. -1 = all but one) + #[arg(long, allow_hyphen_values = true)] + pub max_outgroup_count: Option, /// Minimum fraction of outgroup genomes containing the k-mer [0.0–1.0] #[arg(long)] @@ -239,12 +243,12 @@ pub fn matching_genome_indices(pred_str: &str, genomes: &[GenomeInfo]) -> Result pub struct GroupFilterParams { pub threshold: u32, - pub min_count: Option, - pub max_count: Option, + pub min_count: Option, + pub max_count: Option, pub min_frac: Option, pub max_frac: Option, - pub min_outgroup_count: Option, - pub max_outgroup_count: Option, + pub min_outgroup_count: Option, + pub max_outgroup_count: Option, pub min_outgroup_frac: Option, pub max_outgroup_frac: Option, } @@ -279,12 +283,20 @@ pub fn build_group_filter( let default_min_frac = if !ingroup_preds.is_empty() && !ingroup_quorum_explicit { 1.0 } else { 0.0 }; let default_max_outgroup_count = if !outgroup_preds.is_empty() && !outgroup_quorum_explicit { 0 } else { out_size }; - let min_count = p.min_count.unwrap_or(0); - let max_count = p.max_count.unwrap_or(in_size); + // Resolve a signed count: negative means an offset from the group size + // (e.g. -1 = all but one), floored at 1 so the negative form always keeps + // constraining the group — even a singleton group, where n-1 would be 0 + // and would otherwise drop the constraint entirely. + let resolve = |v: isize, size: usize| -> usize { + if v < 0 { (size as isize + v).max(1) as usize } else { v as usize } + }; + + let min_count = p.min_count.map(|v| resolve(v, in_size)).unwrap_or(0); + let max_count = p.max_count.map(|v| resolve(v, in_size)).unwrap_or(in_size); let min_frac = p.min_frac.unwrap_or(default_min_frac); let max_frac = p.max_frac.unwrap_or(1.0); - let min_outgroup_count = p.min_outgroup_count.unwrap_or(0); - let max_outgroup_count = p.max_outgroup_count.unwrap_or(default_max_outgroup_count); + let min_outgroup_count = p.min_outgroup_count.map(|v| resolve(v, out_size)).unwrap_or(0); + let max_outgroup_count = p.max_outgroup_count.map(|v| resolve(v, out_size)).unwrap_or(default_max_outgroup_count); let min_outgroup_frac = p.min_outgroup_frac.unwrap_or(0.0); let max_outgroup_frac = p.max_outgroup_frac.unwrap_or(1.0);