feat(sort): per-dataset minimizer cutoff exponent

[ivan-aksamentov]

PR: per-dataset minimizer cutoff exponent

feat/minimizer-cutoff -> feat/multiref

Adds an optional per-dataset minimizerIndex.cutoff exponent to pathogen.json, allowing datasets to retain a larger fraction of k-mers in the minimizer index. Builds on feat/multiref (multi-reference scoring). Test with CVA16 dataset using branch cva16-testing-sample-with-cutoff in nextclade_data.

• Depends on: feat/multiref in both nextclade and nextclade_data

Related

• Meeting notes: multi-reference minimizer -- original design discussion
• nextclade_data PR #404 -- multi-ref PR (base for this work)
• nextclade_data PR #412 -- initial CVA16 dataset upload
• nextclade PR #1649 -- sort performance optimization (relevant to the runtime cost of wider cutoffs)

Motivation

The default cutoff 1<<28 retains ~6.25% of k-mers. For some diverse pathogens, a wider cutoff (e.g. 1<<31, ~50%) might improve discrimination by retaining more k-mers per reference. This PR adds the mechanism; whether any dataset benefits from it is a separate evaluation.

What changes

• nextclade: add cutoff: Option<i64> to MinimizerIndexConfig in pathogen.json schema, update generated schema files
• nextclade_data: add cutoff_from_exponent() to resolve the exponent to a raw threshold, plumb per-dataset cutoffs through make_ref_search_index() and rebuild, set global params.cutoff to max(all per-dataset cutoffs)

Design: global vs per-dataset cutoff

The minimizer index has a single params.cutoff that the client uses to hash query k-mers. This must be the max of all per-dataset cutoffs because the client hashes queries once before knowing which dataset matches. Per-dataset cutoff is enforced on the index side: each dataset stores only hashes below its own cutoff. Extra query k-mers generated by the wider global cutoff don't match lower-cutoff datasets (verified: zero hit-count differences).

Deviations from earlier findings

• Results not comparable to earlier cutoff experiments: the earlier finding "cutoff 28 -> 31 reduced unclassified from 540 -> 79" was measured against the old single-ref baseline using unclassified count as the metric. This PR measures cutoff 31 against the new multi-ref baseline on 36 classified genomes using mean score (-4%). Different baselines, metrics, and sample sets. A rerun with the original data and metric would be needed for a direct comparison.
• Per-dataset cutoff has global cost: any single dataset raising its cutoff forces all queries to hash more k-mers for all datasets (+67% sort runtime across all 105 datasets for cutoff 31), plus one borderline SARS-CoV-2 reassignment. The cost is not scoped to the dataset that requested the wider cutoff.

Verification

Tests: 11/11 Python tests pass (7 multi-ref + 4 cutoff-specific). Rust unit tests pass (cutoff field parses and round-trips correctly).

CLI: cutoff isolation experiment

Three-way comparison to isolate the effect of per-dataset cutoff from multi-ref. All three variants use the feat/multiref nextclade binary (which reads expectedMinimizerHits); only the minimizer index differs.

Variant	nextclade binary	nextclade_data branch	CVA16 pathogen.json	Global params.cutoff
A: single-ref	feat/multiref	master	single reference, no minimizerIndex	1<<28
B: multi-ref	feat/multiref	cva16-testing-sample	7 additional refs, no cutoff	1<<28
C: multi-ref + cutoff 31	feat/multiref	cva16-testing-sample-with-cutoff	7 additional refs, "cutoff": 31	1<<31

Setup: build binary and three indexes

# from a parent directory containing both repos
cd nextclade
OUT="$PWD/tmp/cutoff-test"
mkdir -p "$OUT"

# -- build the feat/multiref binary (used for all three variants) --
git checkout feat/multiref
./docker/dev b nextclade
cp .build/docker/debug/nextclade "$OUT/nextclade-multiref"

# -- variant A: single-ref 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 "$OUT/index_single_ref.json"

# -- variant B: multi-ref index from cva16-testing-sample (default cutoff) --
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 "$OUT/index_multiref.json"

# -- variant C: multi-ref + cutoff 31 from cva16-testing-sample-with-cutoff --
git checkout cva16-testing-sample-with-cutoff
# pathogen.json already has minimizerIndex.cutoff = 31 on this branch
./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 "$OUT/index_multiref_cutoff31.json"

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

cd ../nextclade

Test 1: CVA16 score improvement (36 full-length genomes, ~7400 bp)

Commands

cd nextclade
for variant in single_ref multiref multiref_cutoff31; do
  ./docker/dev "$OUT/nextclade-multiref" sort \
    -m "$OUT/index_${variant}.json" \
    -r "$OUT/sort_cva16_${variant}.tsv" \
    "$OUT/cva16_examples.fasta"
done

All 36 sequences classified as enpen/enterovirus/cva16 in all three variants.

Metric	A: single-ref	B: multi-ref	C: multi-ref + cutoff 31	B vs A	C vs B
Mean score	0.352	0.625	0.599	+0.273 (+77%)	-0.026 (-4%)
Min score	0.281	0.501	0.464	+0.220 (+78%)	-0.037 (-7%)
Max score	0.420	0.987	0.928	+0.567 (+135%)	-0.059 (-6%)

Cutoff 31 (C) scores lower than default cutoff (B): mean -4%, min -7%. The wider cutoff increases both numerator (more k-mer hits) and denominator (expectedMinimizerHits), but the denominator grows faster because it is computed analytically as (L - k) * cutoff / 2^32 while actual hits are reduced by non-ACGT k-mers and hash collisions.

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

Commands

cd nextclade
for variant in single_ref multiref multiref_cutoff31; do
  ./docker/dev "$OUT/nextclade-multiref" sort \
    -m "$OUT/index_${variant}.json" \
    -r "$OUT/sort_mixed_${variant}.tsv" \
    "$OUT/queries.fasta"
done

Dataset assignment counts (9 CVA16 fragments, 28 HBV, 12 HIV, 128 SARS-CoV-2):

Dataset	A: single-ref	B: multi-ref	C: multi-ref + cutoff 31
nextstrain/sars-cov-2/wuhan-hu-1/orfs	95	95	96
nextstrain/sars-cov-2/BA.2.86	32	32	31
community/neherlab/hiv-1/hxb2	12	12	12
(unclassified)	38	38	38

A vs B: zero assignment changes. Multi-ref alone causes no regressions.

B vs C: one SARS-CoV-2 sequence (IMS-10004-CVDP-E64B5426) moved from BA.2.86 (score 0.978) to wuhan-hu-1/orfs (score 0.967). This is a borderline case between two valid SARS-CoV-2 datasets -- not a false positive, but a score-ordering change caused by the denominator shift. Both datasets use the default cutoff (1<<28), so their stored minimizers are identical in B and C; only the expectedMinimizerHits denominator differs because nMinimizers (the rounded expectation) changed from round((L-k)*cutoff_28/2^32) to round((L-k)*cutoff_28/2^32) -- the same formula, but the index was rebuilt and the literal stored nMinimizers value shifted by 1-2 due to floating-point rounding. Hit counts are identical (verified: 0 hit-count differences for same-assignment sequences in B vs C).

Index-level comparison

Field (CVA16 entry)	A: single-ref	B: multi-ref	C: multi-ref + cutoff 31
nMinimizers	478	462	3698
expectedMinimizerHits	(absent)	462.19	3697.50
length	7413	7412	7412
Global params.cutoff	1<<28	1<<28	1<<31

Field (SARS-CoV-2 BA.2.86)	A	B	C
nMinimizers	1839	1868	1868
expectedMinimizerHits	(absent)	1867.88	1867.88

The SARS-CoV-2 expectedMinimizerHits is identical in B and C (both use default cutoff for SC2). The score difference comes from the nMinimizers change between A and B (478 -> 462 for CVA16 is irrelevant; 1839 -> 1868 for SC2 matters because A uses nMinimizers as denominator while B and C use expectedMinimizerHits).

Runtime

Variant	Sort time (177 queries)
A: single-ref	0.61s
B: multi-ref	0.61s
C: multi-ref + cutoff 31	1.02s (+67%)

The runtime increase is caused by the global params.cutoff rising from 1<<28 to 1<<31, forcing the client to hash 8x more k-mers per query for all 105 datasets, not just CVA16.

Conclusions

1. Cutoff 31 hurts CVA16 scores vs default cutoff (mean -0.026, min -0.037)
2. Cutoff 31 causes one borderline SARS-CoV-2 reassignment (between two valid SC2 datasets)
3. Cutoff 31 increases sort runtime +67% for all queries due to global cutoff spillover
4. Multi-ref alone (variant B) provides all the CVA16 score improvement with zero regressions and no runtime cost
5. The per-dataset cutoff mechanism works correctly (set membership enforced, zero hit-count changes for non-CVA16 datasets), but has no beneficial application for CVA16

Work items

• Add cutoff: Option<i64> to MinimizerIndexConfig in pathogen.json schema
• Add cutoff_from_exponent() with range validation (exponent in [1, 32])
• Plumb per-dataset cutoffs through make_ref_search_index() and rebuild
• Set global params.cutoff to max of all per-dataset cutoffs
• Add tests for cutoff exponent resolution, superset invariant, and rejection of invalid exponents

[github-actions]

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

(ci)

[rneher]

Does the sort code actually discard hashes that are above the cut-off, or does it look up everything? But even if it looks up everything, it could be that the larger span of values in the map makes the look-up less efficient.

[ivan-aksamentov]

@rneher Query hashes above the cutoff are discarded before lookup [src].

I believe the slowdown comes from a different mechanism: the global params.cutoff rises from 1<<28 to 1<<31, so ~8x more query k-mers pass the filter, producing ~8x more BTreeMap::get lookups [src] per query -- across all 105 datasets, not just CVA16.

Needs to be profiled to be sure what causes it.

I could not figure out how to do the per-dataset cutoff, so this sets the global to max(28, all_dataset_cutoffs). I don't think with current algo per-dataset cutoff is possible (prove me wrong)

Overall it's still very fast, so runtime perf might not have much practical consideration. The decision point is whether it improves anything or not scientifically. This requires some more relevant data and some testing.

[ivan-aksamentov]

@rneher It was .unique() and some other silly things. Report + proposed improvements: #1773