perf(mash): select tophits with one groupby, de-duplicate PLSDB handling (-316 lines)

[sanjaynagi-eit]

Problem

process_mash_tsv exists in two near-identical ~90 line copies (Plass and Assembly). Each one, for every contig, filters and sorts the entire mash dataframe twice inside a python loop:

for contig in contigs:
    hit_df = mash_df.loc[mash_df["contig"] == contig].sort_values("mash_distance").reset_index(drop=True)
    hits = len(hit_df["mash_distance"])
    if hits > 0:
        tmp_df = mash_df.loc[mash_df["contig"] == contig].sort_values("mash_distance").reset_index(drop=True).loc[0]

Cost grows as contigs x hits. The 68-column PLSDB column list is also duplicated verbatim in both copies.

Change

One stable sort plus groupby("contig").first(), in a shared mash_tophits_df(); the PLSDB read moves to a shared load_plsdb_metadata(). Both classes now call the same two helpers. Net -316 lines.

contigs	hits	before	after
5	15	0.006 s	0.004 s
25	75	0.019 s	0.004 s
100	300	0.072 s	0.004 s
400	1200	0.283 s	0.006 s

Equivalence

Identical output whenever a contig's best hit is unique (verified across all sizes above, plus the empty-mash-file path).

Ties: mash distance is a function of matching hash count, so exact ties do occur. The old code took whichever row pandas' unstable quicksort happened to place first — arbitrary. The stable sort here always takes the first such row in mash's output order. At the handful of hits per contig that mash dist -d 0.1 -v 0.1 actually returns, numpy sorts stably anyway, so this reproduces the old result and merely makes it deterministic instead of implementation-defined. There is a test pinning it.

One optimisation deliberately not taken

The PLSDB table is read whole (~240 MB of dataframe) to serve a merge that keeps at most one row per contig, so chunked reading looks like an easy memory win — 237 MB -> 174 MB at 20k-row chunks, -> 102 MB at 2k.

It changes user output. pandas infers dtypes per chunk, and several PLSDB columns are only mixed-type when the whole file is seen. ASSEMBLY_coverage is object over the full table and writes 120; within a small chunk it is float64 and writes 120.0. That frame goes verbatim into <prefix>_summary.tsv.

Not worth pinning 68 dtypes to the current database release for 0.03 s and ~100 MB of transient RSS, so the whole-file read stays and the reason is recorded as a comment next to the column list.

Tests

New tests/test_run_mash.py: lowest-distance selection, a row per contig including misses, the empty-mash path, deterministic tie-breaking, and that extra contigs in the mash output do not add rows. All existing tests pass (75 non-slow).

[sanjaynagi-eit]

Updated — I found and fixed a real output regression in my own patch while diffing end-to-end runs against upstream.

Building the tophits table with a left merge against 1..contig_count introduces NaN for contigs with no hit, and pandas silently promotes the integer columns to float. mash reports mash_pval = 0 for a strong hit, and that came back out of <prefix>_summary.tsv as 0.0 instead of 0.

My micro-benchmark had not caught it because its synthetic p-values were all floats. The rows are now assembled from the selected hits instead of merged, which preserves the original dtypes exactly — and is faster than the merge form (400 contigs x 1200 hits: 0.286 s -> 0.021 s, previously 0.006 s measured on the merge version but with the wrong dtypes).

Two regression tests added: test_integer_mash_pval_is_not_promoted_to_float and test_contig_column_stays_integral.

Same class of bug as the chunked-PLSDB-read issue described in the PR body — worth flagging that a left merge and a chunked parse both change dtypes in ways that reach user output.