OBIKmers/obikmer
1
0 Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 34 Wiki Activity

Implement advanced indexing architecture and query features 

Introduces a partitioned, NUMA-aware indexing system, supports kmer filtering, evidence conversion, index packing, and advanced query capabilities including estimation and phylogenetic analysis.
Eric Coissac committed 2026-09-12 08:51:20 +02:00
commit e381ef2fb1
24 files changed
+1203
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

architecture.md
+33
@@ -0,0 +1,33 @@
# Architecture notes for advanced use
This page describes execution-level behavior relevant to sizing and running `obikmer` on large datasets or multi-socket machines. It complements the [[formats-index_layout|index format]] and [[theory-indexing_architecture|theory]] pages.
## Sequence invariant
Every input sequence is treated purely as a compact representation of a set of overlapping kmers:
- Only the `A`/`C`/`G`/`T` alphabet (case-insensitive) is recognized; a sequence is cut at any other character (including IUPAC ambiguity codes), so runs containing them are not represented in the index.
- Sequences are internally processed in chunks of at most 256 nucleotides; a chunk shorter than k is dropped. This is invisible to the user beyond the ACGT-only, minimum-length-k constraints above.
- Kmers are always handled in canonical form (see [[theory-encoding|DNA encoding]]), so the tool is strand-agnostic throughout: a kmer and its reverse complement are always the same entry.
## Index dimensioning
An index directory is organized as `KmerIndex → partitions → layers`, with a canonical kmer belonging to exactly one (partition, layer) pair. This is what makes set operations (merge, filter, distance) parallel and coordination-free across partitions.
- **Partition count** (`-p`/`--partitions`, rounded up to a power of 2) is the main dimensioning knob: more partitions means more independent parallel units and a smaller working set per partition, at the cost of more open files during construction.
- **Layers** accumulate as an index grows through successive merges; per-partition query cost grows with the number of layers (worst case linear, expected constant since most kmer lookups resolve in the first layer they could plausibly be in).
- Genome columns (count or presence data) are kept at a consistent width across every layer and partition after a merge, which is what allows whole-index aggregate distances (Jaccard, Bray-Curtis, Euclidean, Hellinger, …) to be computed as a two-pass cascade (local partial sums per partition, then a global combination) with no double counting.
## Parallel execution and NUMA awareness
Partition-level work (index construction, `merge`, `filter`, `convert`, `select`, `phylo`'s sibling-annex/Sankoff computations) is dispatched by a partition runner that adapts to the machine's memory topology, detected automatically at startup via hwloc:
- On a multi-socket / multi-NUMA-node machine, one thread pool is pinned per NUMA node, and each partition is processed entirely by threads pinned to one node — keeping the memory a partition touches local to that node's DRAM. This matters because touching kmer data across NUMA nodes without pinning can degrade throughput by an order of magnitude or more on large multi-socket machines.
- On a single-socket machine, Apple Silicon, or if hwloc cannot report NUMA topology, all cores are treated as one node with no pinning and negligible overhead — this is the default behavior on macOS.
- Within a node, the number of active worker threads ramps up progressively rather than being fixed up front: it starts conservatively and grows in steps, but only as long as measured CPU efficiency or disk I/O throughput keeps improving. If neither improves after a growth step, the runner stops adding workers — avoiding oversubscription on stages that are memory-bandwidth-bound rather than CPU- or I/O-bound. Ramp speed scales with the number of cores per node, so a single-node machine ramps just as fast as a large multi-node one.
No CLI flag controls this directly; it is fully automatic at runtime. NUMA-aware pinning can be compiled out (Cargo feature `numa`, on by default), in which case a plain global thread pool is used instead.
## Kmer filtering (`filter`)
[[usage-filter|`filter`]] evaluates predicates against the genome metadata matrix directly whenever every active filter can be expressed as a column-level test (e.g. "any outgroup column non-zero"), producing a per-slot keep/drop decision without touching kmer sequence data at all. If any active filter cannot be expressed this way, evaluation falls back to a per-kmer, row-level check. Either way, the result is always written as a single, freshly compacted layer (`unitigs.bin` and the MPHF are rebuilt from the surviving kmers), never as an additional layer on top of the source index.
formats-index_layout.md
+56
@@ -0,0 +1,56 @@
# Index construction and on-disk layout
## Construction pipeline
Building an index ([[usage-index_command|`index`]]) proceeds through a fixed sequence of phases, each operating independently per partition (see [[theory-indexing_architecture|Partitioning and indexing architecture]]):
1. **Scatter.** A single streaming pass over the input. Each sequence fragment is cut at non-ACGT bases, passed through the low-complexity entropy filter (see [[theory-entropy_filter|Low-complexity kmer filter]]), and any resulting segment shorter than k is dropped. Surviving segments are decomposed into super-kmers, canonicalized, and routed by `hash(minimizer) mod n_partitions` into one file per partition.
2. **Dereplication.** Within each partition, identical super-kmer sequences are merged and their occurrence counts summed. This count is per super-kmer, not per kmer — a kmer's true abundance is the sum of the counts of every super-kmer containing it.
3. **Exact counting.** Every kmer in every dereplicated super-kmer is enumerated and its exact total count computed. A per-genome kmer frequency spectrum is produced at this stage.
4. **Quorum filtering.** Kmers outside the `--min-abundance`/`--max-abundance` range are dropped, and super-kmers are recompacted around the surviving kmer set.
5. **Local assembly.** The surviving kmers of each partition are assembled into unitigs — maximal non-branching runs of a local de Bruijn graph — such that every kmer appears exactly once, at one (unitig, offset) location.
6. **MPHF and evidence construction.** A minimal perfect hash function is built over the canonical kmers of each partition, together with the evidence structure needed to verify that a queried kmer was genuinely indexed (see below). Per-genome counts or presence bits are recorded alongside if requested.
Phases 1–5 are independent per partition and run in parallel; phase 6 finalizes each partition once its kmer set is fixed.
## Minimal perfect hash function (MPHF)
Each partition's surviving kmers are mapped to a dense range of integer slots by a minimal perfect hash function: no collisions, near-optimal space (a few bits per key), O(1) lookup. Because an MPHF maps *any* input to some slot — including kmers that were never indexed — a lookup alone cannot distinguish a genuinely indexed kmer from an arbitrary one; every lookup is followed by an evidence check.
## Evidence: exact vs. approximate
Two verification modes are available, selected at build time (`index --approx`) and convertible afterwards ([[usage-convert|`convert`]]):
- **Exact** (default): the hashed slot stores a pointer back into the partition's unitig data. At query time the kmer is reconstructed from that location and compared directly to the query. Zero false positives, at the cost of one extra random read per lookup.
- **Approximate** (`--approx`): the slot stores a short fingerprint (`--evidence-bits` bits) instead of a pointer; verification is a single fingerprint comparison. This trades a small, bounded false-positive rate ($1/2^b$ per kmer, reduced further to about $1/2^{b \cdot z}$ for a read requiring $z$ consecutive matching kmers via the `-z`/`--findere-z` parameter) for lower memory and disk usage, since no reconstruction index is needed. See [[usage-estimate|`estimate`]] to explore this trade-off before building.
## On-disk layout
```
<index_root>/
index.meta global configuration (k, minimizer size, partition count,
evidence mode, whether counts are stored) and genome list/metadata
scatter.done / count.done / index.done build-progress sentinels
spectrums/<label>.json per-genome kmer frequency histogram
partitions/
part_00000/ ... part_NNNNN/
index/
meta.json number of layers in this partition
layer_0/
unitigs.bin reconstructible kmer sequence data — always kept
unitigs.bin.idx random-access index into unitigs.bin (exact evidence only)
mphf.bin the minimal perfect hash function
evidence.bin exact evidence (exact mode only)
fingerprint.bin approximate evidence (approximate mode only)
counts/ per-genome kmer counts (if counts were requested)
presence/ per-genome presence/absence bits
layer_1/, layer_2/, ... added by later merges, same internal structure
```
`unitigs.bin` is the only file from which the indexed kmer content can be fully recovered; it is always retained. Every other file (MPHF, evidence, counts) is derived from it.
A **layer** corresponds to one increment of kmer content added to a partition — most commonly, one [[usage-merge|`merge`]] operation that introduces kmers not already present in the index. Genomes already present in the index simply gain new columns in the existing layers' count/presence data; only genuinely new kmer content is assembled into a new layer. Because of this, merging cost scales with the novel kmer content being added, not with the accumulated size of the index. A query against an index with several layers checks each layer's MPHF in turn.
Sources merged together must share the same kmer size, minimizer size, partition count, and evidence mode (including matching approximate-mode parameters); mismatches are rejected rather than silently reconciled — [[usage-convert|`convert`]] one of the sources first if needed.
`obikmer pack` consolidates a partition's per-column files (counts/presence) into a single file, reducing the number of file opens needed at query time.
index.md
+68
@@ -0,0 +1,68 @@
# obikmer
`obikmer` is a command-line tool for counting, indexing, querying and comparing DNA sequences represented as kmer sets. It targets individual genome datasets of tens of gigabases, with an emphasis on computational, memory, and disk efficiency.
All functionality is exposed through a single binary, `obikmer`, organized as subcommands.
## Core principles
- Kmers are of fixed, odd length $k$, chosen at index-construction time in the range $[11, 31]$ (see [[theory-kmers_and_superkmers|Kmers and super-kmers]]).
- Each kmer fits in a 64-bit word using a 2-bit-per-base encoding (see [[theory-encoding|DNA encoding]]).
- Kmers are handled in **canonical form** ($\text{canonical}(kmer) = \min(kmer, \text{revcomp}(kmer))$), making counting strand-independent.
- Sequences are decomposed into **super-kmers** before storage, anchored on a hash-selected **minimizer** (see [[theory-minimizer_selection|Minimizer selection]]), then routed to one of several **partitions** for parallel, memory-bounded processing (see [[theory-indexing_architecture|Partitioning and indexing architecture]]).
- Low-complexity kmers can be filtered out at index-construction time using an entropy-based score (see [[theory-entropy_filter|Low-complexity kmer filter]]).
## Commands
| Command | Purpose |
|---|---|
| [[usage-superkmer|`superkmer`]] | Extract super-kmers from a sequence file and write them to stdout |
| [[usage-index_command|`index`]] | Build a genome index |
| [[usage-merge|`merge`]] | Merge multiple indexes into one |
| [[usage-filter|`filter`]] | Retain only kmers matching ingroup/outgroup predicates |
| [[usage-select|`select`]] | Project and/or aggregate genome columns of an index |
| [[usage-query|`query`]] | Query an index with sequences and annotate matches |
| [[usage-dump|`dump`]] | Dump indexed kmers as CSV |
| [[usage-annotate|`annotate`]] | Add, update, or dump genome metadata |
| [[usage-phylo|`phylo`]] | Compute pairwise genome distances, trees, and phylogenetic exports |
| [[usage-unitig|`unitig`]] | Dump the unitigs of an index as FASTA |
| [[usage-estimate|`estimate`]] | Estimate approximate-index parameters before indexing |
| [[usage-convert|`convert`]] | Convert an index's evidence representation (exact/approximate/hybrid), in place |
| [[usage-utils|`utils`]] | Miscellaneous index maintenance and inspection utilities |
| [[usage-pack|`pack`]] | Pack per-column matrix files into a single-file format |
See [[usage-predicates|Genome predicates and taxonomy paths]] for the selection language shared by `filter`, `select`, `dump`, and `unitig`.
## Further reading
- [[formats-index_layout|Index construction and on-disk layout]]
- [[architecture|Architecture notes for advanced use]] — parallel execution, NUMA awareness, index dimensioning
## Input formats
- `superkmer` and `index`: FASTA (`.fa`, `.fasta`), FASTQ (`.fq`, `.fastq`), GenBank flat file (`.gb`, `.gbk`, `.gbff`), all optionally gzip-compressed; directories are expanded recursively; streaming stdin via `-` or when no input path is given.
- `query`: FASTA or FASTQ, optionally gzip-compressed; streaming stdin the same way.
## Parameter constraints
These constraints are checked at startup; an invalid value exits immediately with an error.
| Parameter | Constraint | Reason |
|---|---|---|
| $k$ (`--kmer-size`) | odd, $k \in [11, 31]$ | odd length guarantees the canonical form is always well defined; the range keeps a kmer within a 64-bit word while retaining specificity |
| $m$ (`--minimizer-size`) | odd, $3 \le m \le k-1$ | same palindrome argument as $k$; must be strictly shorter than the kmer |
| $z$ (`-z`, approximate evidence only) | $z \le k-1$ | the effective indexed kmer size is $k-z+1$ |
## Genome label constraints
Genome labels are arbitrary Unicode strings, with the following restrictions:
| Character | Forbidden | Reason |
|---|---|---|
| `/` | yes | filesystem path separator |
| `=` | yes | separator used by `--new-label` |
| `\0` | yes | null byte |
| `\n`, `\r`, `\t` | yes | would break CSV output |
| spaces | allowed | quote in the shell, e.g. `--new-label 'new label=old label'` |
Empty labels are rejected. A label derived automatically from the input file name (when `--label` is omitted) is not validated, since it is already filesystem-safe.
installation.md
+84
@@ -0,0 +1,84 @@
# Installation
## Prerequisites
### Rust toolchain
`obikmer` requires **Rust 1.85 or later** (edition 2024). Install or update via [rustup](https://rustup.rs):
```bash
curl --proto '=https' --tlsv1.2 -sSf https://sh.rustup.rs | sh
rustup update stable
```
### C build environment (required for hwloc)
`obikmer` embeds [hwloc](https://www.open-mpi.org/projects/hwloc/) (Hardware Locality) for NUMA-aware thread placement on multi-socket machines. hwloc is built from source at compile time, which requires a standard C build environment.
#### Linux (Debian/Ubuntu)
```bash
apt install build-essential automake libtool autoconf pkg-config
```
#### Linux (RHEL/Rocky/AlmaLinux)
```bash
dnf install gcc make automake libtool autoconf pkgconfig
```
#### HPC clusters
Most HPC clusters provide these tools via the module system:
```bash
module load gcc automake libtool autoconf
```
If in doubt, check that `autoreconf --version` and `libtool --version` return successfully.
#### macOS
```bash
brew install automake libtool autoconf pkg-config
```
## Building
```bash
git clone <repository-url>
cd obikmer/src
cargo build --release
```
The compiled binary is at `target/release/obikmer`.
### Building on HPC clusters (network filesystems)
HPC home directories are typically on a network filesystem (Lustre, NFS) optimized for large sequential reads, not for the many small file operations Cargo generates during compilation. Building directly on such a filesystem can be extremely slow.
Redirect the build directory to a local scratch disk:
```bash
CARGO_TARGET_DIR=/scratch/$USER/cargo-target cargo build --release
```
Adapt the path to the scratch space available on your cluster (`/var/tmp`, `/tmp`, `/scratch/local`, etc.). Once built, copy the binary to a permanent location:
```bash
cp /scratch/$USER/cargo-target/release/obikmer ~/bin/
```
## NUMA support
NUMA-aware thread placement is active automatically on multi-socket Linux machines, detected at runtime via hwloc. No build flag is required — it falls back gracefully to a single-pool strategy on:
- macOS (Apple Silicon, unified memory)
- single-socket Linux machines
- any system where hwloc reports only one NUMA node
## Verifying the installation
```bash
obikmer --help
```
theory-encoding.md
+28
@@ -0,0 +1,28 @@
# DNA encoding
## 2-bit nucleotide encoding
Every nucleotide is encoded on 2 bits, most-significant-bit first within each word:
| Base | Encoding |
|------|----------|
| A | `00` |
| C | `01` |
| G | `10` |
| T | `11` |
The Watson-Crick complement of a base is its bitwise NOT on 2 bits: $\text{complement}(base) = \lnot base \mathbin{\&} \texttt{0b11}$.
## Kmer encoding
A kmer of length $k$ ($k \le 31$) fits in a single 64-bit word. The first nucleotide occupies the two most significant bits, each following nucleotide occupies the next two bits, and unused low-order bits are zero. Extracting nucleotide i (0-indexed from the 5′ end) is a shift-and-mask operation.
Reverse complement is computed by bit manipulation directly on the packed word, without any lookup table: complement every base, reverse the byte order, then reverse the order of 2-bit groups within each byte in two more passes, and finally realign the result to the most-significant bits.
## Canonical form
The canonical form of a kmer is the lexicographic minimum of the kmer and its reverse complement:
$$\text{canonical}(kmer) = \min\big(kmer,\ \text{revcomp}(kmer)\big)$$
Using the canonical form halves the kmer space and makes counting strand-independent: a kmer and its reverse complement are always treated as the same entity, regardless of which DNA strand was sequenced.
theory-entropy_filter.md
+36
@@ -0,0 +1,36 @@
# Low-complexity kmer filter
Low-complexity kmers (homopolymer runs, tandem repeats) can dominate an index without carrying useful information. `obikmer` detects and excludes them during index construction using a normalized Shannon entropy score.
## Sub-word frequencies
For a kmer of length $k$ and a sub-word size $ws$ ($1 \le ws \le ws_{\max}$, default $ws_{\max} = 6$), the kmer is decomposed into its $k - ws + 1$ overlapping sub-words of length $ws$ by sliding a window across it. Each sub-word is tallied under its raw 2-bit-packed value, with no canonicalization.
## Corrected Shannon entropy
Let $f_j$ be the observed count of raw sub-word $j$, and $n_{\text{words}} = k - ws + 1$ the total number of sub-words. The entropy is:
$$H_{\text{corr}} = \log(n_{\text{words}}) - \frac{1}{n_{\text{words}}} \sum_j f_j \log f_j$$
## Small-sample correction
Because only $n_{\text{words}}$ sub-words are observed among up to $4^{ws}$ possible values, the achievable maximum entropy $H_{\max}$ is bounded below $\log(4^{ws})$ for small samples. $H_{\max}$ is computed from the most uniform integer distribution achievable with $n_{\text{words}}$ observations over $4^{ws}$ categories. The normalized entropy is:
$$\hat{H}(ws) = \frac{H_{\text{corr}}}{H_{\max}} \in [0, 1]$$
A value near 0 indicates low complexity (e.g. a homopolymer run); near 1 indicates high complexity, characteristic of a random sequence.
## Final score
The filter evaluates $\hat{H}(ws)$ for every word size from 1 to ws_max and keeps the minimum:
$$\text{entropy}(kmer) = \min_{ws=1}^{ws_{\max}} \hat{H}(ws)$$
Taking the minimum across word sizes ensures that repetition at any scale is detected: a homopolymer is caught at $ws=1$, a dinucleotide repeat at $ws=2$, and so on. A kmer is rejected if its entropy score falls below a threshold $\theta$ (default 0.7), a configurable collection parameter.
## Properties
The entropy score depends only on the kmer sequence itself, not on where or how many times it occurs:
- **Orientation invariance**: a kmer and its reverse complement always receive the same score.
- **Context independence**: a given kmer is always accepted or always rejected, regardless of which genome or read it appears in. The filter defines a fixed partition of the kmer space into low-complexity and valid kmers.
theory-indexing_architecture.md
+32
@@ -0,0 +1,32 @@
# Partitioning and indexing architecture
An index is split into a fixed number of **partitions**, each handling an independent, disjoint slice of the kmer space. Partitioning keeps the working set of each stage small enough to process efficiently and enables parallel construction and querying.
## Routing
The canonical minimizer of a super-kmer (see [[theory-minimizer_selection|Minimizer selection]]) is hashed to produce a $p$-bit routing value that selects the destination partition:
```
canonical minimizer → hash(minimizer) → p-bit value → partition index
```
The routing value is recomputed whenever it is needed (during construction and again at query time) rather than stored — it is not part of the on-disk super-kmer representation.
Within a partition, kmers are indexed as plain values via a minimal perfect hash function (see [[formats-index_layout|On-disk storage]]); the minimizer plays no further role once a super-kmer has reached its partition.
## Why hashing is necessary
A canonical minimizer is an m-mer ($m \in \{9, 11, 13, 15\}$), and its distribution over all possible m-mer values is not uniform — as the lexicographic minimum of a window, small values are systematically over-represented (@Zheng2020-ji; @Zheng2021-cc; @Pan2024-hb; @Kille2023-px; @Golan2025-xf). Routing directly on the raw minimizer value would therefore produce badly unbalanced partitions.
Hashing the minimizer before routing redistributes this skewed distribution uniformly across partitions. This works reliably because the number of partition-index bits $p$ is chosen well below the number of bits available in the minimizer ($2m$): even with strong bias in the minimizer distribution, the hash has enough entropy margin to absorb it, provided the number of distinct minimizers actually observed is much larger than the number of partitions.
## Parameter guidance
| Minimizer size $m$ | Minimizer bits ($2m$) | Typical partition-index bits $p$ | Partitions |
|----|-----------|-----------|------------|
| 9 | 18 | 6–8 | 64–256 |
| 11 | 22 | 8–10 | 256–1 024 |
| 13 | 26 | 10–12 | 1 024–4 096|
| 15 | 30 | 10–14 | 1 024–16 384|
The number of partitions must satisfy $p \le 2m$, and in practice $p$ is chosen well below that bound to leave a comfortable entropy margin. For $k=31$, $m=13$, $p=10$ (1024 partitions), partition load is well balanced on real genomic data.
theory-kmers_and_superkmers.md
+24
@@ -0,0 +1,24 @@
# Kmers and super-kmers
## Kmers
A **kmer** is a DNA subsequence of fixed length $k$. Two constraints apply to $k$, both enforced when a command starts (an invalid value exits immediately with an error):
- $k \in [11, 31]$: long enough to be specific, short enough to fit in a single 64-bit word at 2 bits/base ($k \le 32$ is the hard limit; $k < 11$ gives insufficient specificity).
- $k$ **is odd**: an odd-length sequence can never equal its own reverse complement, so the two orientations of any kmer are always distinct. This is required for the canonical form (see [[theory-encoding|DNA encoding]]) to be well defined.
## Super-kmers
A **super-kmer** is a maximal run of consecutive, overlapping kmers from a read that share the same canonical minimizer (see [[theory-minimizer_selection|Minimizer selection]]). Each kmer in the run overlaps the next by $k-1$ nucleotides. A super-kmer is capped at 256 nucleotides; a longer run is split at that boundary.
For a random minimizer of length $m$ over kmers of length $k$, the expected length of a super-kmer is approximately (@Zheng2020-ji; @Golan2025-xf):
$$L_{\text{nt}} \approx \frac{k-m+2}{2} + k - 1$$
For $k=31$, $m=13$ this is about 40 nucleotides; in practice super-kmers rarely exceed a few dozen nucleotides.
### Canonical super-kmers
A **canonical super-kmer** is the lexicographic minimum of a super-kmer and its reverse complement. When a read and its reverse complement are both encountered, they produce super-kmers that are reverse complements of each other; both reduce to the same canonical super-kmer, so a genomic region is represented once regardless of which strand was read.
Super-kmers are the unit of work used throughout construction and querying: sequences are decomposed into super-kmers first, and every downstream step (partition routing, deduplication, counting) operates on them rather than on individual kmers.
theory-minimizer_selection.md
+42
@@ -0,0 +1,42 @@
# Minimizer selection
## Definition
A **minimizer** of a kmer window is the m-mer ($m < k$) that is smallest, among all $k - m + 1$ overlapping m-mers in the window, under a chosen ordering. The minimizer is always taken in canonical form (lexicographic minimum of forward and reverse complement) so that selection is strand-independent.
The minimizer partitions a sequence into super-kmers: maximal runs of overlapping kmers that share the same minimizer (see [[theory-kmers_and_superkmers|Kmers and super-kmers]]).
## Hash-based ("random") minimizer
`obikmer` selects minimizers by hash order rather than plain lexicographic order. Ordering m-mers lexicographically on their 2-bit encoding systematically favors AT-rich m-mers (an all-A m-mer always encodes to 0), which causes low-complexity regions to dominate as minimizers and produces unbalanced partitions.
Instead, a well-distributed hash function $H$ is applied to the canonical (lexicographically minimal) form of each m-mer, and the m-mer with the smallest $H$ value wins. Because $H$ is a bijection with good avalanche properties, every distinct m-mer in a window has an equal chance of holding the minimum hash value, independent of its nucleotide composition.
The canonical form used as input to $H$ is still the lexicographic minimum of forward/reverse-complement — hashing is applied on top of it, not used to redefine it. Defining canonicity by hash value instead would bias the *distribution of hash values themselves* toward small values (the minimum of two independent hashes is not uniformly distributed), reintroducing a bias one layer down.
### Hash function
The hash function is a 64-bit mixing function (splitmix64-style finalizer) applied to the m-mer XORed with a fixed non-zero seed:
$$H(x) = \text{mix64}(x \oplus s), \quad s = \lfloor 2^{64}/\varphi \rfloor = \texttt{0x9e3779b97f4a7c15}$$
```
H(x):
x ← x ⊕ 0x9e3779b97f4a7c15
x ← x ⊕ (x >> 30)
x ← x × 0xbf58476d1ce4e5b9
x ← x ⊕ (x >> 27)
x ← x × 0x94d049bb133111eb
return x ⊕ (x >> 31)
```
The XOR seed avoids the finalizer's fixed point at 0 ($\text{mix64}(0) = 0$), which would otherwise make an all-A m-mer (canonical value 0) win every window comparison.
## Partition routing is independent of minimizer selection
The hash used to select a minimizer within a window (the minimum of several hash values) and the hash used to route a super-kmer to a storage partition are computed separately:
- **Selection** uses $H$ applied to every candidate m-mer in the window, keeping the minimum.
- **Partition routing** recomputes $H$ on the single selected minimizer only, once its position is fixed. This is a hash of one specific value, not the minimum of several, so it is uniformly distributed and safe to use directly for routing.
See [[theory-indexing_architecture|Partitioning and indexing architecture]] for how the routing value is turned into a partition index.
usage-annotate.md
+25
@@ -0,0 +1,25 @@
# annotate
Add or update genome metadata of an index from a CSV file, or dump the current metadata as CSV.
```bash
obikmer annotate INDEX --csv FILE [OPTIONS]
obikmer annotate INDEX --dump
```
## Arguments
| Argument | Description |
|---|---|
| `INDEX` | Index directory to annotate (modified in place) |
## Options
| Option | Default | Description |
|---|---|---|
| `--csv` | — | CSV file of metadata to apply (must contain an id column); required unless `--dump` is used |
| `--sep` | `,` | CSV field separator |
| `--id-col` | `id` | Name of the column containing genome labels |
| `--na-value` | `NA` | Value meaning "remove this field" (deletes the existing key if present) |
| `--no-overwrite` | off | Do not overwrite existing metadata keys |
| `--dump` | off | Print all genome metadata as CSV to stdout instead of applying a file |
usage-convert.md
+29
@@ -0,0 +1,29 @@
# convert
Convert an existing index's evidence representation in place, between exact, approximate, and hybrid.
```bash
obikmer convert INDEX (--exact-evidence | --approx-evidence BITS | --hybrid-evidence) [OPTIONS]
```
## Arguments
| Argument | Description |
|---|---|
| `INDEX` | Index directory to convert (modified in place) |
## Options
Exactly one of the first three is required:
| Option | Description |
|---|---|
| `--exact-evidence` | Convert to exact evidence (zero false positives) |
| `--approx-evidence BITS` | Convert to approximate (fingerprint-only) evidence; `BITS` = fingerprint bits per slot (b) |
| `--hybrid-evidence` | Convert to hybrid evidence (both exact and approximate bundles kept) |
| `--evidence-bits BITS` | Fingerprint bits per slot (b) — required with `--hybrid-evidence` when the source index is currently exact; rejected otherwise (the source already fixes `b`) |
| `-z, --findere-z Z` | Findere z parameter: number of consecutive stored kmers that must all match to confirm a hit. This does not shorten the indexed kmer length (fixed forever at `index` build time) — it extends the effective match window: on a k=31 index, `z=2` requires 32 consecutive matching bases, not 30 |
| `--fp FP` | Target false-positive rate per z-window (e.g. `0.01`); derives `b` or `z` when one of them isn't given directly |
| `--block-size N` | Block size for exact evidence's on-disk index (unitigs per block). Ignored when converting to pure approximate evidence. Default `1` |
See [[usage-index_command#exact-vs-approximate-evidence|`index`]] for the exact/approximate trade-off and the underlying false-positive model, and [[usage-estimate|`estimate`]] to explore parameters beforehand. The index directory is locked for exclusive access during conversion.
usage-dump.md
+25
@@ -0,0 +1,25 @@
# dump
Dump all kmers of an index as CSV, one row per kmer, with per-genome counts or presence.
```bash
obikmer dump INDEX [OPTIONS]
```
## Arguments
| Argument | Description |
|---|---|
| `INDEX` | Index directory to dump |
## Options
| Option | Default | Description |
|---|---|---|
| `--force-presence` | off | Output presence/absence (0/1) even if the index stores counts |
| `--debug` | off | Prefix each row with the partition and layer columns |
| `--head N` | none | Limit output to the first N kmers |
`dump` also accepts the shared [[usage-filter#predicate-options|predicate options]] (`--ingroup`, `--outgroup`, `--min-count`, etc.) to restrict which kmers are dumped.
Output is CSV on stdout.
usage-estimate.md
+18
@@ -0,0 +1,18 @@
# estimate
Estimate approximate-index parameters (z, evidence bits, false-positive rate) before building an index with `--approx`, without touching any files.
```bash
obikmer estimate [OPTIONS]
```
## Options
| Option | Default | Description |
|---|---|---|
| `-k, --kmer-size` | `31` | Kmer size used at query time (matches `index`'s `--kmer-size`) |
| `-z, --findere-z` | none | Findere z parameter |
| `--evidence-bits` | none | Fingerprint bits per slot (b) |
| `--fp` | none | Target false-positive rate per z-window |
Any two of `-z`, `--evidence-bits`, `--fp` may be given; the third is derived using the same model as `index --approx` and `convert --approx-evidence` ($FP = 1 / 2^{b \cdot z}$). The report printed to stdout includes: query $k$, effective indexed $k$ ($k-z+1$), $z$, evidence bits, per-kmer false-positive rate, and per-z-window false-positive rate.
usage-filter.md
+47
@@ -0,0 +1,47 @@
# filter
Apply row-level selection to an index: retain only kmers matching ingroup/outgroup predicates over genome membership, plus optional total-count and complexity thresholds. The output is a new, single-layer index.
```bash
obikmer filter SOURCE -o OUTPUT [OPTIONS]
```
## Arguments
| Argument | Description |
|---|---|
| `SOURCE` | Source index directory |
## Options
| Option | Default | Description |
|---|---|---|
| `-o, --output` | — (required) | Output index directory |
| `-f, --force` | off | Overwrite an existing output directory |
| `--presence` | off | Output presence/absence instead of counts |
| `--min-total-count` | none | Minimum total count across all genomes (count index only) |
| `--max-total-count` | none | Maximum total count across all genomes |
| `--min-complexity` | none | Minimum normalized entropy (same score as `--theta` at index build time), recomputed from the stored unitig sequences |
| `--complexity-level-max` | `6` | Maximum sub-word size for the complexity score (used only with `--min-complexity`) |
## Predicate options
| Option | Default | Description |
|---|---|---|
| `--ingroup` | none | Ingroup predicate (repeatable; each occurrence is ANDed) |
| `--outgroup` | none | Outgroup predicate (repeatable; each occurrence is ORed) |
| `--min-count` | 0, or group size + N if negative | Minimum number of ingroup genomes carrying the kmer |
| `--max-count` | ingroup group size | Maximum number of ingroup genomes carrying the kmer |
| `--min-frac` | `1.0` if `--ingroup` given without an explicit quorum, else `0.0` | Minimum fraction of ingroup genomes |
| `--max-frac` | `1.0` | Maximum fraction of ingroup genomes |
| `--min-outgroup-count` | `0` | Minimum number of outgroup genomes carrying the kmer |
| `--max-outgroup-count` | `0` if `--outgroup` given without an explicit quorum, else outgroup group size | Maximum number of outgroup genomes |
| `--min-outgroup-frac` | `0.0` | Minimum fraction of outgroup genomes |
| `--max-outgroup-frac` | `1.0` | Maximum fraction of outgroup genomes |
| `--presence-threshold` | `0` | Minimum count for a genome to be considered a carrier of a kmer |
See [[usage-predicates|Genome predicates and taxonomy paths]] for the predicate syntax used by `--ingroup`/`--outgroup`.
A negative `--min-count`/`--max-count` is interpreted as an offset from the group size — e.g. `--min-count=-1` means "all but one".
Declaring `--ingroup` with no explicit ingroup quorum flag implicitly sets `--min-frac 1.0` (present in every ingroup genome). Declaring `--outgroup` with no explicit outgroup quorum flag implicitly sets `--max-outgroup-count 0` (absent from every outgroup genome). Any explicit quorum flag for a group disables that group's implicit default.
usage-index_command.md
+50
@@ -0,0 +1,50 @@
# index
Build a genome index from one or more sequence files. Construction proceeds in phases (scatter → dereplicate → count → layered MPHF), described in [[formats-index_layout|On-disk storage]].
```bash
obikmer index -o OUTPUT [OPTIONS] [INPUTS...]
```
## Arguments
| Argument | Description |
|---|---|
| `INPUTS...` | Input sequence files or directories (FASTA/FASTQ/GenBank, gzip optional). If omitted, reads from stdin. |
## Options
| Option | Default | Description |
|---|---|---|
| `-o, --output` | — (required) | Output index directory |
| `--force` | off | Overwrite an existing output directory |
| `--label` | input file name without extension | Genome label stored in the index |
| `--meta KEY=VALUE` | none | Attach a categorical metadata field to the genome (repeatable) |
| `-k, --kmer-size` | `31` | Kmer size (odd, in [11, 31]) |
| `-m, --minimizer-size` | `11` | Minimizer size (odd, in $[3, k-1]$) |
| `--theta` | `0.7` | Entropy threshold for the low-complexity filter |
| `--level-max` | `6` | Maximum sub-word size for the entropy score |
| `-p, --partitions` | `256` | Number of partitions (rounded up to a power of 2) |
| `-T, --threads` | detected core count | Number of worker threads |
| `--max-open-files` | `threads / 4` (min 1) | Maximum number of input files open simultaneously |
| `--min-abundance` | `1` | Minimum abundance (inclusive) for a kmer to be retained |
| `--max-abundance` | none | Maximum abundance (inclusive) |
| `--with-counts` | off | Store per-kmer counts; otherwise only presence/absence is stored |
| `--keep-intermediate` | off | Keep intermediate build files instead of deleting them after construction |
| `--approx` | off | Use approximate evidence (Findere fingerprint) instead of exact evidence |
| `-z, --findere-z` | see below | Findere z parameter: number of consecutive kmers that must all match (approximate evidence only) |
| `--evidence-bits` | see below | Fingerprint bits per slot (b), approximate evidence only |
| `--fp` | see below | Target false-positive rate per z-window, approximate evidence only |
| `--block-size` | `1` | Block size, in unitigs, for the exact on-disk index (rounded up to a power of 2) |
## Exact vs. approximate evidence
By default, an index stores **exact** evidence: a kmer is either present or absent (or has an exact count with `--with-counts`), with no false positives.
With `--approx`, evidence is stored as a compact **fingerprint** instead, trading a small, tunable false-positive rate for reduced memory/disk usage. The false-positive model is:
$$FP = \frac{1}{2^{b \cdot z}}$$
where $b$ is `--evidence-bits` and $z$ is `--findere-z`. Any two of `-z`, `--evidence-bits`, `--fp` can be given and the third is derived; if none are given, defaults are $b=8$, $z=1$ ($FP \approx 1/256$). See [[usage-estimate|`estimate`]] to explore this trade-off before building an index, and [[usage-convert|`convert`]] to change an existing index's representation afterwards.
`z` must be strictly less than k: the effective indexed kmer length under approximate evidence is k−z+1.
usage-merge.md
+29
@@ -0,0 +1,29 @@
# merge
Merge multiple built indexes into a single index.
```bash
obikmer merge -o OUTPUT SOURCE... [OPTIONS]
```
## Arguments
| Argument | Description |
|---|---|
| `SOURCE...` | Index directories to merge (at least one required) |
## Options
| Option | Default | Description |
|---|---|---|
| `-o, --output` | — (required) | Output index directory |
| `--force` | off | Overwrite an existing output directory |
| `--force-presence` | off | Store the merged index as presence/absence even if all sources have counts |
| `--rename-duplicates` | off | Disambiguate duplicate genome labels (`.1`, `.2`, …) instead of failing |
| `--budget-fraction` | `0.5` | Fraction of available RAM reserved as the memory budget for parallel partition merging |
## Behaviour
The output mode is chosen automatically: if every source index stores counts, the merged index stores counts too; otherwise it is presence/absence. `--force-presence` forces presence/absence regardless of the sources.
By default, merging two indexes that share a genome label fails with an error; `--rename-duplicates` instead appends a numeric suffix to keep both copies.
usage-pack.md
+29
@@ -0,0 +1,29 @@
# pack
Pack an index's per-column matrix files into a single-file format to reduce query-time I/O (fewer file opens per query).
```bash
obikmer pack INDEX [--sparse]
```
## Arguments
| Argument | Description |
|---|---|
| `INDEX` | Index directory to pack (modified in place) |
## Options
| Option | Default | Description |
|---|---|---|
| `--sparse` | off | Pack presence/absence and count matrices into a sparse, deduplicated format instead of the dense one |
The index directory is locked for exclusive access while packing.
## `--sparse`
Matrix data (which genomes carry each kmer, or with what count) is often mostly empty — most kmers are present in only a handful of genomes out of the whole collection. The default (dense) packed format stores one entry per genome for every kmer regardless of how many genomes actually carry it; `--sparse` instead stores each kmer's genome list directly. For presence/absence matrices, identical genome lists shared by many kmers are also deduplicated (common in real data, since kmers from the same conserved region tend to be carried by the same genomes); for count matrices, the genome list is deduplicated the same way but each kmer's actual counts are kept per-kmer, since two kmers sharing the same genome list rarely carry the same counts.
On real genome collections this has measured at roughly 7x smaller on disk than the dense format for presence/absence, and single-kmer lookups (the shape `phylo`'s sibling-annex/entropy/Sankoff computations use) are typically faster too, since the smaller files mean less data to read from disk. The trade-off: reading a whole genome column at once (used by `--distance` matrix computations) is much slower on the sparse format than on the dense one, since there is no native column layout to read sequentially — prefer the dense format (the default, no `--sparse`) for indexes you mainly query with `phylo`'s `--distance` matrices.
`--sparse` applies to both presence/absence and count matrices — a count index (`--distance` matrix computations included) is packed sparse the same as a presence index.
usage-phylo.md
+344
@@ -0,0 +1,344 @@
# phylo
Compute pairwise distances between the genomes stored in an index, optionally build trees (NJ/UPGMA) from them, and optionally calibrate a 16-state parsimony model for a central-position SNP character with exports for external phylogenetic tools (TNT, PhyG, IQ-TREE).
```bash
obikmer phylo INDEX [OPTIONS]
```
## Arguments
| Argument | Description |
|---|---|
| `INDEX` | Index directory |
## Distance matrix (`--distance`)
| Option | Default | Description |
|---|---|---|
| `--distance` | `jaccard` | See the two tables below for the full list of accepted values |
| `--gamma-shape ALPHA\|auto` | none | Rate-heterogeneity correction, for `snp-*` values that support it (see below). Either a fixed $\alpha$ or `auto`/`estimate` to fit it from the data (see "Automatic $\alpha$ estimation" below). No effect on the other values; rejected if given together with a value that doesn't support it |
| `--presence-threshold` | `1` | Minimum count for a kmer to be considered present, for `jaccard`/`mash` on a count index |
| `--csv` | off | Write the matrix as plain CSV instead of the default relaxed-PHYLIP format |
| `--shared-kmers` | off | Also write the shared-kmer count matrix. Only valid with a whole-index metric, not a `snp-*` value |
| `--nj` | off | Compute and write a Neighbor-Joining tree (Newick) |
| `--upgma` | off | Compute and write a UPGMA tree (Newick) |
| `-o, --output` | none (stdout) | Output file prefix |
Every value routes to one of two independent computations:
### Whole-index metrics
| Value | Definition |
|---|---|
| `jaccard` | $D = 1 - \dfrac{\lvert A \cap B \rvert}{\lvert A \cup B \rvert}$ over the sets of kmers present in each genome |
| `mash` | derived from the Jaccard distance via $D = -\dfrac{1}{k} \ln\!\left(\dfrac{2J}{1+J}\right)$ where $J = 1 - D_{\text{jaccard}}$ and $k$ is the index's kmer size; clamped to 1.0 when $J \le 0$ |
| `hamming` | number of kmer positions where presence differs between the two genomes (presence index only, not normalized): $D = \sum_i \mathbb{1}[a_i \ne b_i]$ |
| `bray-curtis` | $D = 1 - \dfrac{2 \sum_i \min(c_i^A, c_i^B)}{\sum_i c_i^A + \sum_i c_i^B}$ on raw per-kmer counts |
| `relfreq-bray-curtis` | the same formula computed on per-genome relative frequencies $p_i = c_i / \sum_j c_j$ instead of raw counts |
| `euclidean` | $D = \sqrt{\sum_i (c_i^A - c_i^B)^2}$ on raw counts |
| `relfreq-euclidean` | the same formula on relative frequencies |
| `hellinger` | $D = \dfrac{1}{\sqrt{2}} \sqrt{\sum_i \left(\sqrt{p_i^A} - \sqrt{p_i^B}\right)^2}$ on relative frequencies, bounded in $[0, 1]$ |
| `hellinger-euclidean` | the unnormalized variant, $D = \sqrt{2} \times D_{\text{hellinger}}$ |
`hamming` requires a presence/absence index; the others work on either index type.
### `snp-*` corrections
Computed from the central-position SNP model (see "Central-position SNP model" below): a family is the set of up to 4 kmers sharing identical flanking sequence and differing only at the central base. These values require the sibling annex (`--sibling-annex`, below) and are, by default, computed exhaustively over every non-monomorphic family in the index; add `--subsample N` to bound the computation to approximately `N` families instead (see "Sampling at scale" below — the same flag `--pseudo-alignment`/`--sankoff` use, but optional here).
For a genome pair, let $L$ be its total number of eligible loci (both genomes single-copy at that family), $p$ the raw proportion of substitutions among those loci, $P$/$Q$ the transition/transversion proportions, $Q_1$/$Q_2$ Kimura's two transversion categories (A↔C & G↔T vs. A↔T & C↔G), $P_1$/$P_2$ the purine (A↔G) / pyrimidine (C↔T) transition proportions, and $\pi_A,\pi_C,\pi_G,\pi_T$ the pair's pooled base frequencies.
**`snp-raw`**
$$d = p$$
**`snp-jc`**
$$d = -\frac{3}{4}\ln\!\left(1-\frac{4p}{3}\right)$$
**`snp-k2p`**
$$
\begin{aligned}
a_1 &= 1-2P-Q \\
a_2 &= 1-2Q \\
d &= -\frac{1}{2}\ln a_1-\frac{1}{4}\ln a_2
\end{aligned}
$$
**`snp-k81`**
$$
\begin{aligned}
a_1 &= 1-2P-2Q_1 \\
a_2 &= 1-2P-2Q_2 \\
a_3 &= 1-2Q_1-2Q_2 \\
d &= -\frac{1}{4}\left(\ln a_1+\ln a_2+\ln a_3\right)
\end{aligned}
$$
**`snp-f81`**
$$
\begin{aligned}
E &= 1-\left(\pi_A^2+\pi_C^2+\pi_G^2+\pi_T^2\right) \\
d &= -E\ln\!\left(1-\frac{p}{E}\right)
\end{aligned}
$$
**`snp-t92`**
$$
\begin{aligned}
g &= \pi_C+\pi_G \\
w &= 2g(1-g) \\
a_1 &= 1-\frac{P}{w}-Q \\
a_2 &= 1-2Q \\
d &= -w\ln a_1-\frac{1}{2}(1-w)\ln a_2
\end{aligned}
$$
**`snp-tn93`**
$$
\begin{aligned}
g_R &= \pi_A+\pi_G \\
g_Y &= \pi_C+\pi_T \\
k_1 &= \frac{2\pi_A\pi_G}{g_R} \\
k_2 &= \frac{2\pi_C\pi_T}{g_Y} \\
k_3 &= 2\left(g_Rg_Y-\frac{\pi_A\pi_G\,g_Y}{g_R}-\frac{\pi_C\pi_T\,g_R}{g_Y}\right) \\
w_1 &= 1-\frac{P_1}{k_1}-\frac{Q}{2g_R} \\
w_2 &= 1-\frac{P_2}{k_2}-\frac{Q}{2g_Y} \\
w_3 &= 1-\frac{Q}{2g_Rg_Y} \\
d &= -k_1\ln w_1-k_2\ln w_2-k_3\ln w_3
\end{aligned}
$$
**`snp-tv`** — transversions only, deliberately uncorrected:
$$d = Q$$
`--gamma-shape ALPHA` applies to every value above except `snp-raw` and `snp-tv`: each $-\ln(x)$ term in the formulas above is replaced by $\alpha\left(x^{-1/\alpha}-1\right)$ (the same weight, same $x$).
### Automatic $\alpha$ estimation (`--gamma-shape auto`)
`--gamma-shape auto` (or the equivalent `--gamma-shape estimate`) fits $\alpha$ from the index itself instead of requiring a user-supplied value, using a method-of-moments estimator computed once, from the same sampling pass that builds the pairwise substitution tally — no extra scan of the index.
The estimator pools substitution counts by **partition** rather than by genome pair: for partition $i$, let $n_i$ be the total number of substitutions observed across every genome pair, and $L_i$ the total number of eligible loci across every genome pair, in that partition. Define the partition's observed substitution rate:
$$R_i = \frac{n_i}{L_i}$$
Under a single shared substitution rate with no among-site heterogeneity, each $R_i$ would vary only by Poisson sampling noise. Rate heterogeneity is modeled, as elsewhere in this correction, by a $\mathrm{Gamma}(\alpha,\alpha)$-distributed multiplicative rate (mean 1) shared by every locus in a partition — the classical Poisson–Gamma (negative-binomial) mixture. Under that model:
$$\mathbb{E}[R_i] = \mu \qquad \mathrm{Var}[R_i] = \frac{\mu}{L_i} + \frac{\mu^2}{\alpha}$$
where $\mu$ is the pooled substitution rate across every partition. Weighting each partition's squared deviation by its own $L_i$ removes the first (Poisson) term before attributing what's left to genuine rate heterogeneity:
$$\hat\mu = \frac{\sum_i n_i}{\sum_i L_i} \qquad V = \frac{\sum_i L_i\,(R_i-\hat\mu)^2}{\sum_i L_i} \qquad \bar L = \frac{\sum_i L_i}{\text{number of partitions}}$$
$$\hat\alpha = \frac{\hat\mu^2}{V - \hat\mu/\bar L}$$
If the measured variance $V$ doesn't exceed the Poisson floor $\hat\mu/\bar L$ (no detectable over-dispersion across partitions — the data are consistent with a single shared rate), $\alpha$ is left undefined: the correction is silently disabled for that run rather than applying a fabricated value, and a warning is logged. When an estimate is produced, it's logged at the `info` level before the distance matrix is computed.
Note: this is a method-of-moments estimator derived from the standard Poisson–Gamma relationship between substitution counts and gamma-distributed rate variation, applied per-partition — it is not part of Jin & Nei's (1990) original publication, which only defines the `+Γ` distance formula itself and, absent an estimate, recommends the fixed default $\alpha = 1$ (`--gamma-shape 1`) rather than proposing a way to estimate it from data. `alpha < 1` indicates strong among-site rate heterogeneity (many near-invariant loci, a few fast ones); `alpha` growing large makes the correction converge to the uncorrected formula.
### Output
Without `-o`, the matrix goes to stdout in relaxed-PHYLIP format (`n` on the first line, then one `label<TAB>value...` row per genome). With `--csv`, the format is instead a header row `genome,<label1>,<label2>,...` followed by one `<label>,<value1>,<value2>,...` row per genome, 6 decimals. Both formats are symmetric with a zero diagonal, except where noted below.
## `--exclude-genome`, `--min-shared-family`
| Option | Description |
|---|---|
| `--exclude-genome LABEL` | Exclude a genome (repeatable). Drops its row/column from the distance/shared-kmer matrix output, and removes it from the sampling used by `--pseudo-alignment`/`--sankoff`/a `snp-*` `--distance` value. Does not change the value computed for any remaining pair |
| `--min-shared-family N` | Auto-exclude, on top of `--exclude-genome`, any genome whose mean shared-family count against every other genome (see "Family Overlap" below) falls below `N`. Applies only to `--pseudo-alignment`/`--sankoff`/`snp-*` `--distance` — never to the whole-index metrics or their matrix/NJ/UPGMA output |
Neighbor-Joining and UPGMA trees (`--nj`/`--upgma`) are always built from every genome in the index, regardless of `--exclude-genome`/`--min-shared-family`.
## Central-position SNP model
Requires the sibling annex, built once per index:
| Option | Description |
|---|---|
| `--sibling-annex` | Build (or rebuild) the sibling-count/minorant annex — prerequisite for every option in this section, and for a `snp-*` `--distance` value |
| `--sibling-stats` | Write `<prefix>_siblings.csv`: the family-size distribution, per genome and globally |
| `--sibling-hist` | Print the global family-size histogram (1-4 members) only |
| `--family-overlap` | Write `<prefix>_family_overlap.csv`: for every genome pair, how many variable families both genomes carry a call for |
| `--shannon` | Write `<prefix>_entropy.csv`: per-family Shannon entropy, one row per family, full unsampled scan |
| `--pseudo-alignment` | Write `<prefix>_alignment.fasta`: a SNP-only pseudo-alignment. Requires `--subsample N` |
A family is eligible for a genome pair $(i,j)$ only if both genomes carry exactly one of its observed forms (single-copy, unambiguous).
### `--sibling-stats`
`<prefix>_siblings.csv` — family size = number of distinct central bases observed at a family (1-4).
| Column | Meaning |
|---|---|
| `genome` | genome label, or the literal `global` for the last row |
| `1`, `2`, `3`, `4` | for a genome row: number of families of that size where the genome carries ≥ 1 member. For the `global` row: the actual deduplicated family-size histogram — not the sum of the rows above |
### Family Overlap
`--family-overlap` writes `<prefix>_family_overlap.csv`: header `genome,<label1>,<label2>,...`, one row per genome, cell `[i][j]` = number of variable families (family size ≥ 2) where both genome `i` and genome `j` carry a call. The diagonal is always `0`. Every genome is written, unfiltered by `--exclude-genome`/`--min-shared-family`.
`--min-shared-family N` uses the mean of each genome's own row (excluding the diagonal) against this same matrix as its exclusion statistic. There is no universal value for `N` — inspect `--family-overlap`'s own output to find where the real gap sits in a given genome collection before choosing a threshold.
### `--pseudo-alignment`
`<prefix>_alignment.fasta` — one record per non-excluded genome, one column per variable family (family size ≥ 2). Each site is IUPAC-coded from the genome's presence mask at that family: a single observed form → the plain base; several forms → the matching IUPAC ambiguity code; no form → `-`.
| Option | Default | Description |
|---|---|---|
| `--subsample N` | none (mandatory here) | Target number of families to sample |
| `--free-loss` | off | Treat a genome carrying none of a family's observed members as missing data (`?`) instead of `-` |
| `--no-ambiguity` | off | Treat a genome carrying more than one member of a family as missing data (`?`) instead of an IUPAC ambiguity code |
| `--entropy MU` | off (`1.0` if only `--entropy-sd` is given) | Center of the entropy band to favor when sampling |
| `--entropy-sd SIGMA` | off (`0.5` if only `--entropy` is given) | Width of that band |
## Sampling at scale: `--subsample`, `--shannon`, `--entropy`
`--subsample`, `--free-loss`, `--no-ambiguity`, `--entropy`/`--entropy-sd` are shared by `--pseudo-alignment`, `--sankoff` (and everything it implies: `--tnt`/`--phyg`/`--iqtree`), and a `snp-*` `--distance` value — one draw feeds all of them in a single invocation. `--subsample` is mandatory for `--pseudo-alignment`/`--sankoff`; for a `snp-*` `--distance` value it is optional (omitted means every non-monomorphic family in the index, not an approximation).
Combining `--sankoff` (or `--tnt`/`--phyg`/`--iqtree`) with a `snp-*` `--distance` value in the same command reuses that one draw for both — the distance and the Sankoff calibration/alignment are guaranteed to be computed from the *identical* set of sampled sites, never two independent samples, so the two outputs are directly comparable. This only holds within a single command; running them as two separate `obikmer phylo` invocations draws two independent samples even with the same flags — unless `--session` is used (below).
Every invocation, with or without `--session`, draws its own fresh random sample by default — running the same command twice gives two different (but equally valid) samples, which is useful for measuring sampling variance and is kept that way deliberately. `--session` does not change this: it makes a *specific* sample reusable on request, it does not make sampling itself reproducible from one independent run to the next.
### `--session`: reusing a sample across separate commands
| Option | Default | Description |
|---|---|---|
| `--session DIR` | none | Persist the sample (and, for `--sankoff`, its calibration/alignment) in `DIR` so a later, separate `obikmer phylo` invocation with the exact same selection parameters restores it instead of resampling |
| `--session-force` | off | With `--session DIR`: overwrite its saved parameters and cached sample instead of erroring out when this run's parameters don't match. No effect without `--session` |
`DIR` is created if it doesn't exist. If it already holds a sample built with different `--subsample`/`--free-loss`/`--no-ambiguity`/`--exclude-genome`/`--min-shared-family`/`--entropy`/`--entropy-sd` values than this run, the command exits with an error rather than silently using either the old or the new values — pass `--session-force` to discard the old sample and rebuild under the new parameters, or point `--session` at a different directory to keep both.
A `--sankoff`-family run and a `snp-*` `--distance` run share the same cached sample when pointed at the same `--session DIR` — build it once with either, reuse it from the other, in either order, across separate commands.
If a run using `--session` is interrupted (crash, kill, `Ctrl-C`), the next run against the same `--session DIR` resumes from the last automatic checkpoint (roughly every 8 partitions'-worth of sampling progress) instead of starting over. The resumed run's sample is **not** guaranteed to be identical to what an uninterrupted run would have produced past that checkpoint — each process draws its own independent random sequence, same as any two separate invocations do — but nothing already checkpointed is lost, and no work needs redoing beyond that point.
Without `--subsample`, every variable family (family size ≥ 2) is used. With `--subsample N`, roughly `N` families are kept instead, drawn in proportion to how many candidate families each part of the index actually holds, so the sample stays representative of the whole index. If the index has fewer than `N` candidate families, `--subsample` has no effect.
### `--shannon`: measuring how informative a family is
`<prefix>_entropy.csv` has one row per family visited:
| Column | Meaning |
|---|---|
| `layer` | an internal index-layer identifier — stable within one run, not meaningful across indexes |
| `family_idx` | the family's position within that layer |
| `entropy15` | Shannon entropy (bits) over the 16 possible states (the 15 non-empty subsets of `{A,C,G,T}`), genomes absent from the family excluded from the count |
| `entropy4` | Shannon entropy (bits) reduced to the 4 plain bases, kept alongside `entropy15` for comparison |
| `family_size` | number of distinct central bases observed anywhere in the index for this family (2-4) |
| `n_genomes_present` | how many genomes the entropy was computed over |
Run with `--subsample N --shannon` to get a bounded diagnostic sample instead of a full-index pass — useful for choosing `--entropy`/`--entropy-sd` values before a full run.
### `--entropy MU` / `--entropy-sd SIGMA`: biasing the sample toward informative families
By default, `--subsample` draws families uniformly. With `--entropy`/`--entropy-sd`, each family's chance of being kept is instead weighted by how close its own entropy (`entropy15`) is to `MU`, using a Gaussian curve of width `SIGMA` — no hard cutoff. The filter activates as soon as either flag is given; the other defaults to `1.0`/`0.5`. Combine with `--subsample N` (expect somewhat fewer than `N` families kept in practice) or use alone (a soft filter over the whole index, no size target).
The first `phylo` run on a given index that uses `--entropy`/`--entropy-sd` pays a one-time extra cost (every candidate family's entropy is computed once and saved alongside the index); later runs, even with different `MU`/`SIGMA`, reuse that saved data.
## Sankoff calibration and phylogenetic exports
| Option | Default | Description |
|---|---|---|
| `--sankoff` | off | Calibrate a 16-state parsimony cost matrix and matching pseudo-alignment. Requires `--subsample N` |
| `--sankoff-ratio-ceiling` | `0.5` | Exclude genome pairs whose raw SNP ratio exceeds this value from the base-composition part of the calibration |
| `--free-loss` | off | Recode a family's non-detection as the `?` missing-data symbol instead of an ordinary, costed state, throughout `--sankoff` and every export built from it |
| `--tnt` | off | Also write a TNT script (implies `--sankoff`) |
| `--phyg` | off | Also write PhyG input files (implies `--sankoff`) |
| `--iqtree` | off | Also write an IQ-TREE custom model and alignment (implies `--sankoff`) |
| `--iqtree-min-freq` | `0.001` | With `--iqtree --free-loss`: also treat as missing any state rarer than this in the alignment |
| `--sankoff-cost-scale` | `100` | Integer scaling factor applied to costs before rounding, for TNT/PhyG's integer-only cost commands |
### The 16-state model
Each family is a character with 16 possible states: one per subset of the 4 possible central bases (including the empty subset). Calibration combines a $5 \times 5$ transition matrix over family cardinality (0-4 observed forms) and a $4 \times 4$ base-substitution matrix from unambiguous single-copy loci, both restricted to genome pairs at or below `--sankoff-ratio-ceiling`, into a row-normalized $16 \times 16$ transition probability matrix $P$, converted to a symmetric cost matrix via $\text{cost}(a,b) = -\ln P(a,b)$.
`--sankoff` alone writes the cost matrix, the calibration parameters, and a pseudo-alignment recoded so the empty state uses the symbol `0` (never a gap character). It does not run any external tool.
With `--free-loss`, the empty state is recoded to `?` (TNT/PhyG/IQ-TREE's own missing-data symbol) instead of an ordinary, costed 16th state — `?` rather than `-`, since `-` still carries gap/indel semantics in these tools. `--free-loss` also zeroes the cardinality-transition cost between any two states, not just to/from the empty one: gaining or losing a sibling is priced the same way — for free — as gaining or losing the whole family.
### Exports
All three exports reuse the `--sankoff` calibrated matrix and pseudo-alignment, recoded for the target tool:
- **`--tnt`**: a self-contained TNT script (alignment recoded to TNT's fixed 16-symbol alphabet, integer-scaled cost matrix re-closed to a metric, a default search block).
- **`--phyg`**: a custom cost-matrix file plus a PhyG script reusing the `--sankoff` alignment directly.
- **`--iqtree`**: a custom substitution-model file (exchangeability matrix recovered as $R(a,b) = e^{-\text{cost}(a,b)}$, plus empirical state frequencies) and a matching alignment, for maximum-likelihood inference with real branch lengths. Only states actually occurring in the alignment are kept and compactly renumbered.
TNT and PhyG both write trees with bare numeric leaf labels (`1`, `2`, ..., in the order the genomes appear in `<prefix>_sankoff.fasta`).
## Output files
With `-o/--output PREFIX`, the relevant subset of the files below is written. Without `-o`, only the distance matrix is produced, on stdout. All matrices use genome labels as row/column headers, in index order.
### Distance matrix
| File | Written by | Format | Content |
|---|---|---|---|
| `<prefix>_dist.phy` | always, unless `--csv` | relaxed PHYLIP | the `--distance` matrix |
| `<prefix>_dist.csv` | `--csv` | CSV matrix | the `--distance` matrix, 6 decimals |
| `<prefix>_shared.csv` | `--shared-kmers` | CSV matrix | shared-kmer count per genome pair (integers) |
| `<prefix>_nj.nwk` | `--nj` | Newick | Neighbor-Joining tree |
| `<prefix>_upgma.nwk` | `--upgma` | Newick | UPGMA tree |
CSV matrix layout (`_dist.csv`, `_shared.csv`, `_family_overlap.csv`): header `genome,<label1>,<label2>,...`, one data row per genome, `<label>,<value1>,<value2>,...`.
### Central-position SNP model
| File | Written by | Format | Content |
|---|---|---|---|
| `<prefix>_siblings.csv` | `--sibling-stats` | CSV table | family-size distribution, per genome and global |
| `<prefix>_family_overlap.csv` | `--family-overlap` | CSV matrix | variable families both genomes of a pair carry a call for |
| `<prefix>_entropy.csv` | `--shannon` | CSV table | per-family Shannon entropy, see "Sampling at scale" above |
| `<prefix>_alignment.fasta` | `--pseudo-alignment` | FASTA | SNP-only pseudo-alignment, IUPAC-coded |
### Sankoff calibration and exports
| File | Written by | Format | Content |
|---|---|---|---|
| `<prefix>_sankoff_matrix.csv` | `--sankoff`/`--tnt`/`--phyg`/`--iqtree` | CSV matrix | calibrated 16×16 cost matrix |
| `<prefix>_sankoff_params.yaml` | same flags | YAML | calibration report (raw tallies + derived probabilities) |
| `<prefix>_sankoff.fasta` | same flags | FASTA | Sankoff-recoded pseudo-alignment, header carries an `n_sites` annotation |
| `<prefix>_sankoff.tnt` | `--tnt` | TNT script | ready-to-run parsimony search |
| `<prefix>_sankoff.tcm` | `--phyg` | PhyG TCM | cost matrix in PhyG's own format |
| `<prefix>_sankoff.pg` | `--phyg` | PhyG script | ready-to-run parsimony search |
| `<prefix>_iqtree.model` | `--iqtree` | IQ-TREE model file | custom ML substitution model |
| `<prefix>_iqtree.fasta` | `--iqtree` | FASTA | alignment recoded for that model |
| `<prefix>_iqtree_states.csv` | `--iqtree` | CSV table | maps `_iqtree.model`/`_iqtree.fasta`'s compact state symbols back to `_sankoff_matrix.csv`'s alphabet |
**`_sankoff_matrix.csv`** — header `state,0,A,C,M,G,R,S,V,T,W,Y,H,K,D,B,N`: the 16 symbols are IUPAC codes for the 16 subsets of the 4 possible central bases (bit 0=A, 1=C, 2=G, 3=T), `0` standing for the empty/absent state. One row per source state, one value per destination state, cost $-\ln P(a,b)$, 4 decimals.
**`_sankoff_params.yaml`**:
| Key | Meaning |
|---|---|
| `ratio_ceiling` | the `--sankoff-ratio-ceiling` value used |
| `cardinality_transitions` | 5×5 list of `{from, to, count, probability}`, family cardinality (0-4 observed forms) |
| `composition_transitions` | 4×4 list of `{from, to, count, probability}`, base letters `A/C/G/T`, single-copy substitutions |
**`_sankoff.fasta`** — recoded to match `_sankoff_matrix.csv`'s alphabet: absent state is `0` (or `?` under `--free-loss`). Excluded genomes dropped; columns left monomorphic by that exclusion are re-checked and dropped too.
**`_sankoff.tnt`** (`--tnt`) — `xread` block (alignment recoded to TNT's fixed `0-9A-F` alphabet), an integer-scaled (`--sankoff-cost-scale`) and metric-closed `smatrix`, a default `hold 20; mult; export` search. Run with `printf 'proc <path>;\nquit;\n' | tnt`. Produces `<prefix>_sankoff.tre` (bare numeric leaf labels, order matching `_sankoff.fasta`).
**`_sankoff.tcm`** (`--phyg`) — first line: the 16-symbol alphabet plus a trailing gap symbol (17 total). Each following line: one row of the integer-scaled, metric-closed cost matrix (17 values — the extra gap column/row reuses the cost to/from the empty state `0`).
**`_sankoff.pg`** (`--phyg`) — script: `read(prefasta:..., tcm:...)` against `_sankoff.fasta`/`_sankoff.tcm`, a default 300s/4-instance `search`, `report(...)` writing `<prefix>_sankoff.tre`. Run with `phyg` from the output directory (the script uses relative file names).
**`_iqtree.model`** (`--iqtree`) — lower-triangular exchangeability matrix $R(a,b) = e^{-\text{cost}(a,b)}$ (one row of increasing length per state, whitespace-separated, PAML order), followed by one line of empirical state frequencies. Only states actually occurring in the alignment are kept, compactly renumbered `0..k-1`.
**`_iqtree.fasta`** (`--iqtree`) — alignment recoded to that same compact `0..k-1` alphabet (symbols `0-9A-F`). Under `--free-loss`, non-detection becomes `?` and columns left non-informative once missing calls are ignored are dropped first (required for `+ASC`); with `--iqtree-min-freq` also set (the default), any state rarer than that threshold is folded into the same `?` treatment, and non-informative columns are re-checked and dropped again. Run with:
```
iqtree3 -s <prefix>_iqtree.fasta --seqtype MORPH -m <prefix>_iqtree.model+ASC --prefix <prefix>_iqtree -T AUTO
```
**`_iqtree_states.csv`** (`--iqtree`) — one row per state actually kept in `_iqtree.model`/`_iqtree.fasta` (header `iqtree_symbol,canonical_symbol,frequency`): `iqtree_symbol` is the compact `0-9A-F` symbol as written in those two files, `canonical_symbol` is the matching `_sankoff_matrix.csv` state, `frequency` is that state's empirical frequency at full precision. Under `--free-loss`, absent (`0`/`?`) is never a kept state, so it never appears here — nor does any state `--iqtree-min-freq` folded away for being too rare.
### Rare states and `--iqtree-min-freq`
States that combine 3 or 4 central bases at once (IUPAC `V`/`H`/`K`.../`N`) are inherently rare, and can make `iqtree3` itself numerically unstable ("Numerical underflow for lh-derivative" warnings). With `--free-loss` set, `--iqtree-min-freq` (default `0.001`, one in a thousand) extends the missing-data treatment to any state below this frequency, not just absence. Check `_iqtree_states.csv` to see exactly which states survived and at what frequency; set `--iqtree-min-freq 0` to keep every state that occurs at all. Has no effect without `--free-loss`.
usage-predicates.md
+46
@@ -0,0 +1,46 @@
# Genome predicates and taxonomy paths
Several commands ([[usage-filter|`filter`]], [[usage-select|`select`]], [[usage-dump|`dump`]], [[usage-unitig|`unitig`]]) select or group genomes using the same predicate language over genome metadata (see [[usage-annotate|`annotate`]] for attaching metadata to a genome).
## Predicate syntax
| Form | Meaning |
|---|---|
| `*` or `all` | Matches every genome (case-insensitive) |
| `key=v1\|v2` | Genome's `key` metadata equals one of the listed values |
| `key!=v` | Genome's `key` metadata does not equal `v` |
| `key~path` | Genome's `key` metadata (a taxonomy path) matches `path` (ancestry match) |
| `key!~path` | Genome's `key` metadata does not match `path` |
A genome whose metadata does not contain `key` at all cannot be classified by that predicate and is excluded from the relevant group's quorum count.
Multiple `--ingroup` predicates are combined with AND; multiple `--outgroup` predicates are combined with OR. When both an ingroup and an outgroup predicate would match the same genome, ingroup classification wins.
## Taxonomy paths
A metadata value is treated as a taxonomy path when it starts with the literal prefix `taxonomy:/`; any other value is treated as a plain string and only supports `=`/`!=`.
```
taxonomy:/segment1@rank1/segment2@rank2/...
```
Each segment is a name, optionally annotated with a rank (e.g. `@family`, `@genus`, `@species`); ranks are optional and can be mixed within a path. The `@` character is reserved inside taxonomy paths and cannot appear in segment names or rank labels.
### Path matching (`~` / `!~`)
Matching compares segment names only (ranks are informational, not part of the match), with anchoring controlled by leading/trailing `/`:
| Pattern | Matches |
|---|---|
| `A/B` | anywhere in the path |
| `/A/B` | at the start of the path (prefix) |
| `A/B$` | at the end of the path (suffix) |
| `/A/B$` | the entire path (exact) |
A rank-qualified query, `key@rank=value`, matches only when the path's segment at that specific rank equals `value`.
### Example
```bash
obikmer filter source -o output --ingroup "taxon~/Betulaceae/Betula"
```
usage-query.md
+38
@@ -0,0 +1,38 @@
# query
Query an index with sequences and annotate each query with the kmer matches found.
```bash
obikmer query INDEX INPUTS... [OPTIONS]
```
## Arguments
| Argument | Description |
|---|---|
| `INDEX` | Index directory to query against |
| `INPUTS...` | Input sequence files (FASTA/FASTQ, gzip optional); at least one required |
## Options
| Option | Default | Description |
|---|---|---|
| `--detail` | off | Report per-position, per-genome coverage vectors in the output |
| `--count-missing` | off | Also count query kmers absent from the index |
| `--force-presence` | off | Report presence (0/1) per genome instead of raw counts |
| `--presence-threshold` | `1` | Minimum accumulated count to declare a genome present (implies `--force-presence`) |
| `-z, --findere-z` | derived from the index metadata | Override the Findere z parameter |
| `-T, --threads` | detected core count | Number of worker threads |
| `--chunk-size` | auto-sized (available RAM ÷ threads, clamped to 4–256 MiB) | I/O chunk size, in MiB |
| `--max-open-files` | `threads / 4` (min 1) | Maximum number of input files open simultaneously |
## Output
FASTA on stdout, one record per query, annotated in the OBITools-style header format `>id {"key":value,...}`:
- `kmer_count`: total number of kmers matched
- `kmer_missing`: number of query kmers absent from the index (only with `--count-missing`)
- `kmer_strict_matches`: per-genome match counts
- `coverage`: per-position, per-genome coverage vectors (only with `--detail`)
`--mismatch` is accepted by the CLI but not currently functional; using it produces a warning and is ignored.
usage-select.md
+48
@@ -0,0 +1,48 @@
# select
Project and/or aggregate the genome columns of an index into a new index. Where [[usage-filter|`filter`]] selects rows (kmers), `select` operates on columns (genomes): grouping several genomes into one aggregated column, reordering columns, or dropping some.
```bash
obikmer select SOURCE --output OUTPUT [OPTIONS]
```
## Arguments
| Argument | Description |
|---|---|
| `SOURCE` | Source index directory |
## Options
| Option | Default | Description |
|---|---|---|
| `-o, --output` | — | Output index directory (required) |
| `-f, --force` | off | Overwrite an existing output directory |
| `--group NAME:PRED` | none | Define a named group of genomes by predicate (repeatable; mutually exclusive with `--aggregate-by`) |
| `--group-op NAME:OP` | none | Aggregation operator for a named group |
| `--aggregate-by KEY` | none | Automatically create one group per distinct value of a metadata key (mutually exclusive with `--group`) |
| `--aggregate-op OP` | none | Aggregation operator applied to every auto-generated group |
| `--select COL,...` | all columns | Output columns, in order (group names or genome labels) |
| `--presence-threshold` | `0` | Minimum count for a genome to be considered a carrier (logical operators only) |
| `--dense` | off | Pack the output's presence matrices in the dense format instead of the default sparse one |
| `--force-copy` | off | Copy each layer's unchanged kmer-identity files (mphf/unitigs/evidence/fingerprint) instead of hard-linking them |
## Aggregation operators
`any`, `all`, `none` (logical, evaluated against `--presence-threshold`), `sum`, `min`, `max` (numeric, count index only). If a group's operator is left unspecified, it defaults to `any` when the source is a presence/absence index and `sum` when it stores counts.
A `select` never changes the underlying kmer set — only the per-genome data (counts or presence) is rewritten, so an unaggregated pass-through column (a plain genome label in `--select`) is a cheap copy.
At least one output column must be defined; every name listed in `--select` must resolve to either a defined group or an existing genome label. See [[usage-predicates|Genome predicates and taxonomy paths]] for the predicate syntax used by `--group`.
## Disk usage
`select` always writes to a new output directory — there is no in-place mode. Each layer's kmer-identity files (MPHF, unitigs, evidence, fingerprint) never change under a column projection/aggregation, so they are hard-linked into the output rather than copied: no extra disk is used for them, even on a very large index. Linking falls back to a real copy automatically if it fails (e.g. `SOURCE`/`OUTPUT` on different filesystems). Use `--force-copy` to always copy instead — needed when the output must be able to survive independently of the source on disk (a hard link shares the same underlying data, so overwriting one path outside `select` itself would affect the other).
To replace an index with a selected version of itself, select to a temporary directory and swap it in:
```bash
obikmer select INDEX --output INDEX.tmp --group ... --group-op ... --select ...
rm -rf INDEX
mv INDEX.tmp INDEX
```
usage-superkmer.md
+27
@@ -0,0 +1,27 @@
# superkmer
Extract super-kmers from one or more sequence files and write them to stdout, without building a full index. Useful for inspecting or piping the super-kmer decomposition of a dataset.
```bash
obikmer superkmer [OPTIONS] [INPUTS...]
```
## Arguments
| Argument | Description |
|---|---|
| `INPUTS...` | Input sequence files or directories (FASTA/FASTQ/GenBank, gzip optional). If omitted, reads from stdin. |
## Options
| Option | Default | Description |
|---|---|---|
| `-k, --kmer-size` | `31` | Kmer size (must be odd, in [11, 31]) |
| `-m, --minimizer-size` | `11` | Minimizer size (must be odd, in $[3, k-1]$) |
| `--theta` | `0.7` | Entropy threshold; kmers with a normalized entropy at or below this value are excluded |
| `--level-max` | `6` | Maximum sub-word size used for the entropy score |
| `-p, --partitions` | `256` | Number of partitions (rounded up to the next power of 2) |
| `-T, --threads` | detected core count | Number of worker threads |
| `--max-open-files` | `threads / 4` (min 1) | Maximum number of input files open simultaneously |
Output is written to stdout in the internal scatter format used by `index`; it is primarily intended to be piped into other tools or inspected for debugging.
usage-unitig.md
+19
@@ -0,0 +1,19 @@
# unitig
Dump the unitigs of an index as FASTA. A unitig is a maximal non-branching path through the de Bruijn graph implied by the index's kmers; the concatenation of every unitig reconstructs every stored kmer exactly once.
```bash
obikmer unitig INDEX [OPTIONS]
```
## Arguments
| Argument | Description |
|---|---|
| `INDEX` | Index directory |
## Options
`unitig` accepts the shared [[usage-filter#predicate-options|predicate options]] (`--ingroup`, `--outgroup`, `--min-count`, etc.) to restrict which kmers are included before the unitigs are enumerated.
Output is FASTA on stdout.
usage-utils.md
+26
@@ -0,0 +1,26 @@
# utils
Miscellaneous index maintenance and inspection utilities.
```bash
obikmer utils INDEXES... [OPTIONS]
```
## Arguments
| Argument | Description |
|---|---|
| `INDEXES...` | One or more index directories |
## Options
| Option | Scope | Description |
|---|---|---|
| `--new-label NEW=OLD` | single index only | Rename a genome label |
| `--upgrade-index` | single index only | Add any missing layer metadata files to an older index |
| `--bits-per-kmer` | single index only | Print bits-per-kmer statistics |
| `--stats` | single index only | Print per-genome kmer counts as CSV |
| `--partition-stats` | one or more indexes | Print a partition-size distribution report |
| `--csv FILE` | with `--partition-stats` | Also write raw per-(partition, source) data to FILE as CSV |
At least one operation option must be given. All options except `--partition-stats` require exactly one index directory.