Skip to content

# Construction pipeline

All phases after scatter are embarrassingly parallel across partitions.

## Phase 0 — Parameter estimation

The construction parameters p, n, and min_count depend on the kmer frequency spectrum of the dataset. Estimating this spectrum before construction avoids costly re-partitioning if p is badly chosen.

Two approaches are supported:

- External estimation (preferred): run NT-CARD on the input files and pass its histogram output to `obikmer build`. NT-CARD produces a kmer frequency histogram in a single streaming pass using ntHash and a Flajolet-Martin-style estimator; obikmer reads this file and derives p, n, and min_count automatically.
- Internal estimation (future): an `obikmer estimate` subcommand for users who prefer a single-tool workflow. The implementation would combine two components: (1) ntHash, a rolling hash that updates the kmer hash in O(1) per nucleotide by incrementally adding the incoming base and removing the outgoing one — Rust crates exist; (2) a Flajolet-Martin-style streaming estimator that maintains a small table of minimum hash values and infers the frequency histogram from their statistical distribution, as described in the NT-CARD paper (Mohamadi et al. 2017)1.

The histogram gives:

- F0 (number of distinct kmers) → sets p (target ~10M kmers/partition → p = ⌈log₂(F0 / 10M)⌉)
- frequency distribution → sets n (choose n so that fewer than 1% of kmers overflow)
- error valley → suggests min_count (typically the local minimum between the error peak and the coverage peak)

## Phase 1 — Scatter

Single streaming pass over raw input files (FASTA/FASTQ, gzip). FASTQ quality scores are ignored.

Input files are read via `open_nuc_stream`, which opens and decompresses the file, auto-detects the format (FASTA / FASTQ / GenBank), and yields a sequence of `NucPage` buffers. Each `NucPage` is a flat 64 KB buffer of normalised bytes (`ACGT` + `\x00` separators), carrying a k−1 byte overlap from the preceding page so that no k-mer is lost at page boundaries. Per-record identity (sequence id, raw bytes) is not preserved; this is intentional — the scatter phase only needs normalised bases to produce superkmers.

For each read fragment within a page:

1. Ambiguous base filter: cut at any non-ACGT base; discard fragments shorter than k.
2. Entropy filter: scan each fragment with a sliding window of size k. When the kmer \(K_i = S[i \mathinner{..} i+k-1]\) ended by nucleotide \(S[j]\) (with \(j = i+k-1\)) has entropy below threshold \(\theta\), emit the current segment and start a new one (see algorithm below). \(K_i\) belongs to neither segment, and no valid kmer is lost.
3. Length filter: discard any segment shorter than k produced by step 2.
4. Super-kmer extraction: for each clean segment, slide a minimizer window and group consecutive kmers sharing the same canonical minimizer; canonise each super-kmer by lexicographic comparison with its reverse complement (early exit).
5. Partition routing: `hash(canonical_minimizer) → PART` → append super-kmer to `partition/superkmers.bin.gz`.

Segmentation behavior:

When \(K_i\) (ended by \(S[j]\), \(j = i+k-1\)) fails the entropy threshold:

- Current segment \(S[\textit{seg_start} \mathinner{..} j-1]\) is emitted (last valid kmer = \(K_{i-1}\))
- New segment starts at \(S[i+1]\) (first new kmer = \(K_{i+1}\))
- \(K_i\) is excluded: current segment lacks \(S[j]\), new segment lacks \(S[i]\)
- Overlap = \(S[i+1 \mathinner{..} j-1]\) = \(k-2\) nucleotides

Algorithm — Entropy filter: sliding window segmentation

    procedure EntropyFilter(S, N, k, θ):
        seg_start ← 0
        window ← []
        for j ← 0 to N−1:
            window.push(S[j])
            if |window| < k: continue
            i ← j − k + 1
            if entropy(window) ≤ θ:
                emit S[seg_start .. j−1]
                seg_start ← i + 1
                window ← S[i+1 .. j]
            else:
                window.pop_front()
        emit S[seg_start .. N−1]

Writes are sequential and append-only — IO-friendly. Gzip applied at write time. Data volume ≈ raw genome size (2 bits/nt compaction offsets header overhead).

## Phase 2 — Dereplication

Performed independently per partition. Identical super-kmers are consolidated and their COUNT accumulated — analogous to amplicon dereplication in metabarcoding. Uses external bucket sort to stay within RAM bounds:

Pass 1 (streaming): hash the nucleotide payload of each super-kmer, route to one of B bucket files:

    hash(sequence) % B → bucket_i.bin

B ≈ 100 is tunable; RAM needed ≈ partition_size / B.

Pass 2: for each bucket, load into an in-memory `HashMap<sequence, COUNT>`, dereplicate by summing COUNT values, write consolidated super-kmers.

After dereplication: at Nx coverage the partition shrinks by ~x (errors aside). The COUNT field in each super-kmer header = number of times that exact super-kmer sequence was observed across all input reads.

Important: super-kmer COUNT ≠ individual kmer count. A kmer can appear in multiple distinct super-kmers (same partition, different flanking context); its true count = sum of COUNT of all super-kmers containing it. A super-kmer with COUNT=1 may contain only high-abundance kmers, each appearing in many other super-kmers. Abundance filtering therefore cannot be applied at this phase.

## Phase 3 — Per-kmer count aggregation and quorum filtering

For each dereplicated super-kmer, enumerate its kmers and accumulate counts:

    for each super-kmer (sequence, COUNT):
        for each kmer in sequence:
            kmer_counts[canonical(kmer)] += COUNT

Implemented as a three-step pipeline in `count_partition()`:

1. External sort (`kmer_sort::sort_unique_kmers`): read dereplicated superkmers, extract canonical kmer raw `u64` values, sort in RAM-bounded chunks (adaptive: 40% of available RAM ÷ n_threads, min 1 M kmers/chunk), k-way merge with inline dedup → `sorted_unique.bin`. f0 is now known exactly.
2. Provisional MPHF (ptr_hash): built from `sorted_unique.bin` via `new_from_par_iter(f0, ...)`. Stored to `mphf1.bin`; `sorted_unique.bin` deleted immediately.
3. Accumulation pass: re-read dereplicated superkmers; for each kmer, `slot = mphf.index(kmer.raw())`, increment `counts1[slot]` by the superkmer COUNT. Stored in a `PersistentCompactIntVec` (`counts1.bin`).

At the end of this phase, each distinct canonical kmer has its exact total count, and the frequency spectrum (`spectrums/{label}.json`) is written to the index root.

No pre-filter on super-kmer COUNT is possible at phase 2: a super-kmer with COUNT=1 may contain only high-abundance kmers, each present in many other super-kmers across the partition.

## Phase 4 — Super-kmer compaction

The valid kmer set from phase 3 is used as a mask to rewrite the super-kmer files:

    for each dereplicated super-kmer:
        scan kmer by kmer
        kmer not in valid set → break point (terminates current super-kmer)
        kmer in valid set     → extend current super-kmer

Three cases per super-kmer:

- All kmers valid → copied as-is
- No kmer valid → discarded
- Mixed → split into sub-super-kmers at invalid boundaries; each sub-super-kmer inherits the original COUNT

After splitting, re-apply dereplication (bucket sort, phase 2 method) — splitting can produce new identical super-kmers. This re-dereplication is cheap: the volume is already greatly reduced.

Output: a clean super-kmer file where every kmer passes quorum. This file feeds phase 5.

## Phase 5 — Local de Bruijn graph and unitig construction

Within each partition, build a local de Bruijn graph from the valid kmer set and compute its unitigs. All operations are local to the partition — no cross-partition communication.

    valid kmers → HashSet<u64>

    for each kmer K:
        out_degree = |{K[1:]+b | b ∈ {A,C,G,T}} ∩ HashSet|
        in_degree  = |{b+K[:-1] | b ∈ {A,C,G,T}} ∩ HashSet|

    internal node ↔ in_degree=1 AND out_degree=1
    branching / dead-end → unitig start or end

Traverse non-branching paths to assemble unitigs. Kmers whose neighbours fall in other partitions appear as dead ends locally — they terminate the unitig. The result: each kmer appears in exactly one unitig within the partition.

The partition size (controlled by p) must be calibrated so that the HashSet fits in RAM during this phase.

Output: `unitigs.bin` — the permanent evidence structure for the partition. Each kmer in the partition appears at exactly one (unitig_id, offset) location.

Scope of local unitigs: these are unitigs of the partition's local de Bruijn graph, not global unitigs. A kmer whose k-1 successor or predecessor falls in another partition appears as a dead end locally and terminates the unitig. This does not affect correctness of verification but means partition-local unitigs cannot be directly reused for global assembly.

## Phase 6 — MPHF construction and index finalisation

`build_index_layer` is called per partition (in parallel via `build_layers`) with the following parameters sourced from `IndexConfig`:

- `block_bits` — from `IndexConfig::block_bits`; controls the `.idx` block size (2^block_bits unitig chunks per block) for exact evidence
- `evidence` — `EvidenceKind::Exact` or `EvidenceKind::Approx { b, z }`; propagated unchanged from `IndexConfig::evidence`
- `min_ab` / `max_ab` — abundance bounds applied before graph construction
- `with_counts` — whether to store kmer counts alongside set membership

Abundance filtering: when `min_ab > 1` or `max_ab.is_some()`, the provisional `mphf1.bin` and `counts1.bin` produced in phase 3 are memory-mapped. Each canonical kmer is accepted only if its count in `counts1` satisfies the bounds. If either file is absent, filtering is skipped (all kmers accepted).

    for each kmer in dereplicated super-kmer:
        ab = counts1[mphf1.index(kmer.raw())]
        if ab < min_ab || ab > max_ab: skip
        graph.push(kmer)

Graph build and unitig write:

The surviving kmers are fed into `GraphDeBruijn`, which computes degrees and yields unitigs. Unitigs are written to `layer_0/unitigs.bin` via a `UnitigFileWriter`.

MPHF and evidence build:

`Layer::build` (membership-only) or `Layer::<PersistentCompactIntMatrix>::build` (with counts) is called next. Internally, `MphfLayer::build` performs two passes:

1. Pass 1 (parallel): build `unitigs.bin.idx` (block size = 2^`block_bits`) then construct the MPHF from all canonical kmers in `unitigs.bin`; store to `mphf.bin`.
2. Pass 2 (sequential): for each kmer in `unitigs.bin`, compute its slot and write `evidence.bin` (`chunk_id: 25 bits | rank: 7 bits` packed into a `u32`); also invoke the payload callback (`fill_slot`) to populate `counts/` if `with_counts`.

After `Layer::build` completes, `layer_meta.json` records `EvidenceKind::Exact`.

Approximate evidence override:

If `evidence` is `EvidenceKind::Approx { b, z }`, `build_approx_evidence` is called immediately after `Layer::build`. It overwrites the exact evidence bundle with `fingerprint.bin` (b-bit hash per slot) and rewrites `layer_meta.json` with `EvidenceKind::Approx { b, z }`. No `.idx` file is needed at query time in this mode.

    // Exact path  →  evidence.bin + unitigs.bin.idx + layer_meta.json(Exact)
    // Approx path →  fingerprint.bin + layer_meta.json(Approx{b,z})
    //               (evidence.bin left on disk but not used)

Partition metadata:

After all layer files are written, `PartitionMeta { n_layers: 1 }` is serialised to `index/meta.json` inside the partition directory. This file is required by `LayeredMap::open` for subsequent merge operations.

File layout per partition after phase 6:

    part_XXXXX/
      index/
        meta.json                ← PartitionMeta { n_layers: 1 }
        layer_0/
          unitigs.bin            ← permanent evidence (all modes)
          unitigs.bin.idx        ← block index (exact mode only)
          mphf.bin               ← MPHF
          evidence.bin           ← exact evidence (exact mode)
          fingerprint.bin        ← b-bit fingerprints (approx mode)
          layer_meta.json        ← EvidenceKind tag
          counts/                ← PersistentCompactIntMatrix (with_counts only)

Cleanup: unless `--keep-intermediate` is set, `remove_build_artifacts` deletes `dereplicated.skmer.zst`, `mphf1.bin`, and `counts1.bin` after all partitions are indexed.

See obilayeredmap and MPHF selection for data structure details.

Query path (exact evidence):

    query kmer q
      → canonical_minimizer(q) → hash → PART → part_XXXXX/
      → MPHF(q) → slot s
      → evidence[s] = (chunk_id, rank)
      → read k nucleotides at rank in unitigs[chunk_id] → compare with q
      → match   : return payload[s]   ← exact hit
      → no match: kmer absent         ← MPHF collision on absent kmer

Query path (approximate evidence):

    query kmer q
      → MPHF(q) → slot s
      → fingerprint[s] matches seq_hash(q)?
      → yes : probable hit (FP rate = 1/2^b per kmer, 1/2^(b·z) per z-window)
      → no  : kmer absent

1. Mohamadi, H., Khan, H. & Birol, I. (2017). ntCard: A streaming algorithm for cardinality estimation in genomics data. Bioinformatics (Oxford, England), 33, 1324--1330. ↩