feat(sort): multi-reference minimizer scoring for dataset suggestion

[ivan-aksamentov]

PR: multi-reference minimizer scoring for dataset suggestion

feat/multiref -> master

Paired with nextstrain/nextclade_data branch feat/multiref (PR #404), which produces the index this change consumes. Test with CVA16 dataset using branch cva16-testing-sample in nextclade_data.

Related

• Meeting notes: multi-reference minimizer -- original design discussion
• nextclade_data PR #404 -- data-side multi-ref PR
• nextclade_data PR #412 -- initial CVA16 dataset upload
• nextclade_data PR #442 -- prior minimizerReferences refactoring (superseded)
• nextclade PR #1639 -- prior fix to dataset suggestion scoring for short references
• nextclade PR #1409 -- minimizer threshold tuning (0.3 -> 0.1)
• nextclade issue #1554 -- sort without pre-calculated index

What this enables

A single dataset can now be built from several reference sequences and still be suggested correctly by nextclade sort. This matters for genetically diverse pathogens (for example enteroviruses), where one reference is not close enough to every circulating sequence for reliable auto-detection, and where splitting the pathogen into several datasets only creates ambiguity about which one to analyze a sequence with.

The problem

Dataset suggestion sketches a query into a set of minimizers and, for each candidate dataset, counts how many of the query's minimizers appear in that dataset's stored set. The count is turned into a score by dividing by the dataset's stored minimizer count, nMinimizers.

When a dataset is built from several references, the index stores the union of the references' minimizers. A query only ever resembles one of those references, so its hit count stays roughly constant while nMinimizers grows with every reference added. Dividing by the union count therefore drives the score down as a dataset is given better coverage, which is exactly backwards: the datasets best able to recognize a sequence are penalized for it.

The fix

Divide by the expected number of hits from a single reference instead of by the union size. For a reference of length $L$, k-mer length $k$, and acceptance cutoff $c$ over a $b$-bit hash space, the expected count is

$$E = (L - k),\frac{c}{2^{b}}$$

which is the number of scanned k-mers times the fraction retained by the cutoff. $E$ is invariant to how many references the dataset merges, so a query that matches one reference scores the same whether the dataset carries one reference or ten. Dividing by the union size would instead penalize a dataset purely for describing more of its own diversity.

The index now stores $E$ as a new optional expectedMinimizerHits field, computed per reference at index-build time. A separate field is preferred over repurposing nMinimizers because:

• nMinimizers is an integer -- rounding the fractional $E$ loses precision that can affect score ordering between close datasets
• Older clients read nMinimizers as the denominator; a new field lets them keep working unchanged

For single-reference datasets $E \approx$ nMinimizers, so the fallback to the integer field is safe and keeps older indexes scoring correctly. nMinimizers is deprecated but kept required for older clients.

The index is also validated when loaded: a non-positive k, cutoff, or reference length, and a non-finite or non-positive denominator, are rejected with an error rather than silently producing infinite or negative scores.

What changes in each repository

• nextclade (this PR): consume expectedMinimizerHits as the score denominator with a nMinimizers fallback, deprecate nMinimizers, validate the index on load, regenerate the index JSON schema, and recognize the minimizerIndex config in the pathogen.json schema so dataset builds no longer flag it as an unknown property.
• nextclade_data (feat/multiref): author minimizerIndex.references (union of several FASTA files per dataset), compute and store expectedMinimizerHits as the denominator, and guard degenerate inputs during the build.

Deviations from agreed plan

• Code change, not data-only: the agreed plan was a data-side-only fix. The implementation adds a new expectedMinimizerHits field (f64) to the Nextclade index schema because nMinimizers is i64 -- rounding the fractional analytic denominator loses precision that affects score ordering between close datasets. Older clients fall back to nMinimizers (same value, rounded).
• +77% CVA16 improvement is confounded: the before/after comparison bundles two effects: (1) the denominator fix (expectedMinimizerHits replacing union-count nMinimizers) and (2) seven additional CVA16 references in the test index (denominator drops ~2020 -> ~462 due to shorter per-reference expectation). The test does not isolate the multi-reference contribution alone.

Backward compatibility

The two changes are backward compatible and can be reviewed independently, but they must reach production together: until the data side emits expectedMinimizerHits, this consumer falls back to nMinimizers, which on a multi-reference index is the union count and therefore the deflated denominator. No runtime warning distinguishes a healthy single-reference index from a deflated multi-reference one, so the pairing is enforced by coordinating the merge, not by the code.

Verification

The scoring and validation paths are covered by serde round-trip tests over both the old and new index formats and by rejection tests for each validation rule. Full unit suite and clippy pass.

CLI: end-to-end nextclade sort comparison

Before/after comparison using the CVA16 dataset (Coxsackievirus A16, lineages B1a/B1b/B1c/C/D/E/F) -- the motivating use case for multi-reference scoring.

Setup: build before/after binaries and indexes

Four artifacts needed. Assumes nextclade and nextclade_data repos are checked out side by side.

Artifact	Repo	Branch	Role
nextclade-master	nextclade	master	Before binary (scores by nMinimizers)
nextclade-multiref	nextclade	feat/multiref	After binary (scores by expectedMinimizerHits)
minimizer_index_before.json	nextclade_data	master	Before index (single-ref)
minimizer_index_after.json	nextclade_data	cva16-testing-sample	After index (multi-ref, 7 additional CVA16 refs)

# -- nextclade: build "before" binary from master --
cd nextclade
OUT="$PWD/tmp/multiref-test"
mkdir -p "$OUT"
git checkout master
./docker/dev b nextclade
cp .build/docker/debug/nextclade "$OUT/nextclade-master"

# -- nextclade: build "after" binary from feat/multiref --
git checkout feat/multiref
./docker/dev b nextclade
cp .build/docker/debug/nextclade "$OUT/nextclade-multiref"

# -- nextclade_data: rebuild "before" index from master --
cd ../nextclade_data
git checkout master
./scripts/rebuild --input-dir data/ --output-dir data_output/ \
  --no-pull --allow-dirty \
  --nextclade-schemas-dir=../nextclade/packages/nextclade-schemas
cp data_output/minimizer_index.json "../nextclade/$OUT/minimizer_index_before.json"

# -- nextclade_data: rebuild "after" index from cva16-testing-sample --
# (test branch with CVA16 multi-ref dataset configured on top of feat/multiref)
git checkout cva16-testing-sample
./scripts/rebuild --input-dir data/ --output-dir data_output/ \
  --no-pull --allow-dirty \
  --nextclade-schemas-dir=../nextclade/packages/nextclade-schemas
cp data_output/minimizer_index.json "../nextclade/$OUT/minimizer_index_after.json"

# -- copy query files --
cp data/enpen/enterovirus/cva16/sequences.fasta "../nextclade/$OUT/cva16_examples.fasta"
cp tests/minimizers/data/queries.fasta "../nextclade/$OUT/queries.fasta"

cd ../nextclade

Test 1: CVA16 example sequences (36 full-length genomes, ~7300-7400 bp)

Commands

# from nextclade repo root
./docker/dev $OUT/nextclade-master sort \
  -m $OUT/minimizer_index_before.json \
  -r $OUT/sort_cva16_before.tsv \
  $OUT/cva16_examples.fasta

./docker/dev $OUT/nextclade-multiref sort \
  -m $OUT/minimizer_index_after.json \
  -r $OUT/sort_cva16_after.tsv \
  $OUT/cva16_examples.fasta

All 36 sequences correctly classified as enpen/enterovirus/cva16 in both runs. Scores improved across the board:

Metric	Before (single-ref)	After (multi-ref)	Delta
Mean score	0.352	0.625	+0.273 (+77%)
Min score	0.281 (OQ091684)	0.501 (PX448982)	+0.220 (+78%)
Max score	0.420 (OQ091677)	0.987 (OP562190)	+0.567 (+135%)

Per-sequence scores (36 sequences, sorted by delta)

seqName           before   after    delta
OQ091677            0.420   0.595   +0.175
PX448982            0.316   0.501   +0.185
KY792582            0.370   0.571   +0.201
PX449037            0.313   0.520   +0.207
PP585416            0.344   0.550   +0.207
PX448882            0.323   0.536   +0.214
PX448985            0.350   0.565   +0.215
PX448865            0.337   0.554   +0.216
PX449023            0.316   0.534   +0.218
MN541006            0.332   0.553   +0.222
PX449025            0.314   0.538   +0.224
PX448962            0.351   0.582   +0.231
PX448850            0.324   0.556   +0.232
KY792581            0.354   0.589   +0.235
MT641429            0.341   0.588   +0.248
PX448978            0.367   0.614   +0.248
OP562200            0.394   0.648   +0.254
OQ091684            0.281   0.536   +0.255
MG571837            0.333   0.594   +0.261
PQ182742            0.308   0.569   +0.261
PQ182744            0.366   0.639   +0.273
MW036460            0.356   0.633   +0.277
MW999270            0.362   0.676   +0.281
MF189180            0.352   0.635   +0.283
MH361019            0.350   0.636   +0.286
MT212012            0.379   0.671   +0.292
MW999256            0.360   0.650   +0.291
OQ791560            0.356   0.658   +0.302
MW999270            0.362   0.676   +0.313
MW030435            0.352   0.676   +0.323
MT212004            0.364   0.688   +0.324
MW030431            0.366   0.699   +0.333
MT212000            0.370   0.712   +0.342
KY425537            0.398   0.766   +0.368
KP751558            0.418   0.846   +0.428
OP562190            0.375   0.987   +0.612

OP562190 scores 0.987 -- it is one of the lineage references in the multi-ref set.

Test 2: Regression check -- mixed-pathogen queries (177 sequences)

9 CVA16 fragments (120 bp), 28 HBV, 12 HIV, 128 SARS-CoV-2.

Commands

./docker/dev $OUT/nextclade-master sort \
  -m $OUT/minimizer_index_before.json \
  -r $OUT/sort_before.tsv \
  $OUT/queries.fasta

./docker/dev $OUT/nextclade-multiref sort \
  -m $OUT/minimizer_index_after.json \
  -r $OUT/sort_after.tsv \
  $OUT/queries.fasta

Dataset	Before	After
nextstrain/sars-cov-2/wuhan-hu-1/orfs	95	95
nextstrain/sars-cov-2/BA.2.86	32	32
community/neherlab/hiv-1/hxb2	12	12
(unclassified)	38	38

0 assignment changes across all 177 sequences. 0 hit-count differences for the 138 classified sequences. The 9 CVA16-named queries are 120 bp fragments -- too short for minimizer matching (k=17). The 28 HBV sequences have no dataset in the index.

Conclusions

1. Multi-ref: +77% mean score improvement for CVA16 (0.352 -> 0.625) with zero regressions on other datasets
2. expectedMinimizerHits correctly compensates for union size inflation -- verified by zero hit-count differences for non-CVA16 datasets
3. Short query fragments (<500 bp) are below the practical threshold for minimizer-based sorting regardless of multi-ref support

Work items

• Divide the suggestion score by expectedMinimizerHits, falling back to nMinimizers when absent
• Add optional expectedMinimizerHits to the minimizer index schema (regenerated JSON and YAML)
• Deprecate nMinimizers while keeping it required for older clients
• Validate the minimizer index on load (reject non-positive k, cutoff, length, and non-finite or non-positive denominator)
• Add serde round-trip and validation-rejection tests; add the approx dev-dependency for float assertions
• Add a typed minimizerIndex field (references) to the pathogen.json config so the generated schema documents it instead of leaving it in the untyped catch-all

[github-actions]

Preview: https://nextstrain--nextclade--pr-1771.previews.neherlab.click

(ci)