V4 streaming primary - #304
Merged
Merged
Conversation
Primary clustering is single-linkage at a fixed cutoff, which is identical to finding connected components. The old path instead materialized every pairwise comparison into a long Mdb (N^2 rows) and pivoted it into a dense N x N matrix before handing it to scipy. That is the source of issue #259 and the MemoryError at 43k genomes (24.7 GiB alloc failure at pd.concat). Compute the components directly with union-find instead, streaming the comparison output and discarding the ~99.9% of above-cutoff pairs as they are read. Memory becomes O(genomes + kept_edges) instead of O(genomes^2). Measured peak RSS on synthetic MASH dist input (old -> new): 4,000 genomes: 2.16 GB -> 0.49 GB 8,000 genomes: 7.42 GB -> 0.55 GB 16,000 genomes: ~30 GB -> 0.54 GB Old grows quadratically; new stays flat. Cluster membership is identical to scipy single-linkage across thresholds. Changes: - New drep/d_cluster/union_find.py: UnionFind, cluster_long_df (no pivot), cluster_mash_files (streams MASH dist), cluster_skani_sparse_files. - Split primary/secondary linkage: new --primary_clusterAlg (default single, routes to union-find) separate from secondary --clusterAlg (default average, unchanged). Add --classic_primary_clustering to force the dense scipy path. - New --primary_algorithm skani: runs `skani triangle --sparse`, which emits only above-threshold pairs, streamed into union-find. No N^2 is built on disk or in RAM. Recommended for very large genome sets. - Bound the multiround Mdb: subsample per-chunk tables and free each chunk after clustering, so the concat that previously OOM'd stays bounded. - Keep the primary dendrogram for modest genome sets (the dense pivot is cheap at small N); skip it above the cutoff, as multiround already did. The low_ram marker changed from "optimized_method_used" to "union_find_streaming"; assertions updated accordingly. Co-Authored-By: Claude Opus 4.8 <[email protected]>
skani reports aligned fractions as percentages (0-100), but load_skani only
divided ANI by 100 and passed the aligned fraction through untouched. Every
other comparison algorithm (fastANI, ANImf, ANIn, gANI) reports
alignment_coverage on a 0-1 scale, and that is the scale cov_thresh is compared
against in make_linkage_Ndb:
d.loc[d['alignment_coverage'] <= cov_thresh, 'ani'] = 0
So for --S_algorithm skani the coverage filter was effectively inert: a pair
aligning over only 1% of the genome has alignment_coverage=1.04, which sails
past a cov_thresh of 0.5 and keeps its ANI. Distantly related genomes that share
a small conserved region could then be clustered together.
On the bundled test genomes, E. casseliflavus EC20 and the E. faecalis genomes
align over ~1% of their length at ~93% ANI. With -sa 0.85 -nc 0.5 they were
being merged into a single secondary cluster -- two different species
dereplicated into one -- because the 1.04 "coverage" passed the filter. With the
fraction converted they correctly separate.
Note this changes results for existing --S_algorithm skani users; coverage
filtering now actually applies.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
The subprocess comparison algorithms re-sketch genomes on every invocation.
Greedy secondary clustering is the worst case: for each genome it spawns a
fastANI subprocess against the growing representative list, re-sketching every
representative each time. That is O(N*R) sketching work and N subprocess spawns.
pyskani lets us sketch each genome exactly once, keep the representatives in an
in-memory database, and query it directly -- O(N) sketching, no subprocesses, no
temp files.
Greedy clustering benchmark (synthetic divergent genomes, so every genome founds
its own cluster and representatives accumulate to N):
40 genomes: fastANI 18.0s -> pyskani 0.70s (25.8x)
80 genomes: fastANI 83.2s -> pyskani 1.94s (42.9x)
fastANI grows ~quadratically, pyskani ~linearly. Real datasets have fewer
representatives, so expect a smaller (still large) speedup.
pyskani agrees with the skani executable to within 5e-05 on both ANI and
alignment coverage, and produces identical secondary clusters.
Changes:
- New drep/d_cluster/pyskani_backend.py: PyskaniDatabase (sketch once, query
many), run_pairwise_pyskani, pyskani_one_vs_many.
- --S_algorithm pyskani, wired into both the pairwise and greedy paths. This
also lifts greedy clustering's fastANI-only restriction.
- pyskani is an optional dependency (pip install drep[pyskani]) with an
actionable error if missing; the skani/fastANI executables remain the default.
- prepare_for_greedy now creates its data folder for every algorithm, not just
fastANI (the representative list is written regardless).
alignment_coverage is the aligned fraction of the genome in the 'reference'
column, matching skani's .af matrix and load_fastani. Greedy needs the opposite
row orientation from pairwise, because get_cluster_rep reads the representative
out of the 'querry' column; both orientations are explicit in query().
Co-Authored-By: Claude Opus 4.8 <[email protected]>
skani's triangle defaults to dropping pairs that align over less than 15% of the genome. Primary clustering is a deliberately loose, inclusive pre-filter -- the MASH path applies no alignment-fraction filter at all -- so inheriting skani's default silently strands related genomes in separate primary clusters, where the secondary algorithm never compares them. This matters most for exactly the data dRep is used on. Validating against 2,631 real ocean MAGs (TOBG), skani's default min-af dropped 16 above-threshold edges and split 15 primary clusters that --min-af 0 keeps together: fragmented, partial MAGs of the same organism routinely align over well under 15% of their length. The pairwise skani path already passes --min-af 0 for the same reason. Co-Authored-By: Claude Opus 4.8 <[email protected]>
Reverts the --min-af 0 from 6e56349, which was wrong, and makes the threshold configurable via --primary_skani_min_af (default 15, skani's own default). 6e56349 reasoned that primary clustering is a loose pre-filter and MASH applies no alignment-fraction filter, so skani shouldn't either. That reasoning missed that MASH similarity and skani ANI are not the same measurement. MASH compares k-mers across the whole genome, so two genomes sharing only a small conserved region score as distant. skani reports the identity *within aligned regions only*, so that same pair looks like a high-ANI edge. skani's --min-af is what makes its ANI comparable to MASH's whole-genome similarity; it is not an obstacle to work around. Under single linkage the consequence is severe, because a handful of spurious bridges merge everything they touch. Measured on 10,000 real UHGG genomes: min-af edges clusters largest cluster 0 354,690 734 5,857 <- 59% of the dataset in one cluster 10 340,387 944 688 15 339,004 984 626 MASH 100M 989 626 Dropping min-af from 15 to 0 adds only 4.4% more edges but collapses 59% of the dataset into a single primary cluster, which would make secondary clustering intractable. At the default, skani reproduces the MASH partition closely (984 vs 989 clusters, identical largest cluster; each splits ~1-2% of the other's multi-genome clusters, consistent with MASH overestimating ANI at the threshold). The earlier evidence for --min-af 0 -- 15 clusters on 2,631 ocean MAGs that min-af 0 kept together -- is better explained as those same low-overlap bridges being correctly rejected. Also pins the pyskani extra to >=0.2 (needed for the `cutoff` query argument; note there is no macOS arm64 wheel for 0.2, so Apple Silicon builds from source and needs a Rust toolchain). Co-Authored-By: Claude Opus 4.8 <[email protected]>
estimate_time assigned `time` inside an if/elif chain with no else, so any
algorithm it didn't recognize raised UnboundLocalError instead of returning an
estimate. Adding pyskani as an --S_algorithm choice without updating this
function meant every pyskani run died the moment secondary clustering started:
File "drep/d_cluster/utils.py", line 108, in estimate_time
return time
UnboundLocalError: local variable 'time' referenced before assignment
Caught by a 10,000-genome end-to-end run, which got through primary clustering
(984 clusters, ~4 min) and then fell over here before doing any real work.
Replace the chain with a lookup that falls back to the fast-algorithm estimate
for anything unrecognized. This only drives a log line -- it should never be able
to take a run down. Add a test covering every --S_algorithm choice plus an
unknown one.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
dRep sketched every genome twice and computed the same ANI values twice. Primary
ran one skani pass over all genomes; secondary then ran skani again, once per
primary cluster. But secondary only ever compares genomes *within* a primary
cluster, and those pairs are a subset of what the primary pass already computed.
On 10,000 UHGG genomes, 94% of the pairs driving secondary clustering were
already present in primary's sparse output, with ANI identical to 6 decimal
places -- because it is literally the same skani computation. The remaining 6%
were entirely pairs below primary's --min-af, with a median aligned fraction of
1.1%, which cov_thresh discards anyway.
So run skani once and derive both stages from the same edge table:
- the single pass emits at min-af = min(primary_skani_min_af, cov_thresh*100),
since secondary applies the looser coverage filter; the stricter primary
filter is applied when forming primary clusters
- the screen also stays at or below the secondary threshold, as those edges
are now reused rather than recomputed
- Mdb keeps every edge (not just those above P_ani) with alignment coverage,
so it can feed secondary. It is also a far more useful Mdb: 779k real ANI
values instead of 100M Mash distances
Measured on 10,000 UHGG genomes, dereplicate end to end:
two-pass one-pass
secondary 15 min 24 s
total 22.4 min 13.8 min
Output is identical, not merely equivalent: same 984 primary clusters, same
1,232 secondary clusters, same 1,232 representative genomes. 984 subprocess
spawns and an entire re-sketch of all 10,000 genomes are gone.
Reuse applies when --primary_algorithm skani is paired with a skani
--S_algorithm and greedy is off; --no_reuse_primary_comparisons forces the old
behavior. Other secondary algorithms measure something different and still run
for themselves.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
Everything measured on the 10,000-genome UHGG set said skani should be the
default, but it was still opt-in: a default run got MASH primary (a 100M-row
Mdb, 5.8 GB of RAM, 15 GB of disk) and fastANI secondary. Now `dRep compare`
with no flags runs the one-pass skani path -- 22.4 min -> 13.8 min end to end at
10k, with 187x less intermediate disk.
--primary_algorithm MASH -> skani
--S_algorithm fastANI -> skani
MASH and fastANI remain available. Two consequences of the flip, handled here:
- the mash executable is no longer required unless --primary_algorithm MASH is
requested; the startup dependency check now looks for whichever program the
run will actually use
- --multiround_primary_clustering and --primary_chunksize are MASH-path
concepts. Rather than silently ignoring them under skani, warn: skani's
sparse output never builds the N^2 table multiround exists to avoid, and it
has none of multiround's chunk-splitting imprecision
Remove pyskani. The one-pass work obsoleted it: secondary clustering no longer
runs comparisons at all for skani workflows, and greedy's reason to exist
(dodging O(n^2) within a cluster) is subsumed by sparse output that never
produces n^2 in the first place. It was also the only dependency needing a Rust
toolchain on Apple Silicon, and was 12x *slower* than skani for pairwise
comparisons because it is single-threaded -- a trap for anyone choosing it
expecting speed.
Remove --low_ram_primary_clustering, deprecated since union-find became the
default for single-linkage primary clustering. This also retires the networkx
connected-components path it was the only caller of, so networkx is no longer a
dependency.
The primary dendrogram needed the same small-N treatment on the skani path that
the MASH path already had, plus a fix: leaf labels were derived from Mdb, but the
sparse Mdb legitimately omits genomes with no above-threshold pairs, so labels
and linkage disagreed. Labels now come from the matrix the linkage was built on.
Docs: regenerate the embedded CLI help, rewrite the algorithm overview for the
skani defaults and one-pass reuse, document why skani needs an aligned-fraction
filter where Mash does not, and correct the dependency lists. The algorithm
section had also been stale independently -- it still called ANImf the default.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
The script passed -t (threads) to nucmer unconditionally. That option only
exists in MUMmer 4; MUMmer 3's nucmer rejects it outright rather than ignoring
it:
$ nucmer --mum -p out -c 65 -g 90 -t 6 ref.fa query.fa
Unknown option: t
exit 1
So the script crashed with "nucmer failed with exit code 1" for anyone running
MUMmer 3 -- which is what `conda install mummer` still gives you, as 3.23.
MUMmer 4 is a separate `mummer4` package, and nothing declared that dependency.
dRep's own ANImf comparisons were unaffected because gen_nucmer_cmd never passed
-t, which is why this went unnoticed.
Detect whether nucmer understands -t and only pass it if so, rather than
requiring a specific MUMmer major version. MUMmer 4 still gets its threading;
MUMmer 3 now works instead of crashing.
Fixes 4 failing tests in test_bonus.py, and adds a regression test that asserts
-t is passed if and only if the installed nucmer supports it.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
Documents the changes since 3.7.1 from a user's point of view, leading with what changes their results: skani as the default for both clustering steps, primary clustering switching to single linkage, the Mdb.csv semantics change, and the skani alignment-coverage fix. Deliberately omits development churn that never reached a release -- the pyskani backend and the --min-af 0 experiment were both added and removed within the v4 branch, so users never saw them and they are not "changes since 3.7.1". Marked Unreleased rather than dated; drep/VERSION stays at 3.7.1 until the release is actually cut. Co-Authored-By: Claude Opus 4.8 <[email protected]>
The one-pass refactor replaced the streaming skani reader (which showed a tqdm bar) with a direct load, but the progress argument came along and was never wired up. An API that accepts progress=True and silently does nothing is worse than not offering it; the sparse edge list loads in about a second even for 10,000 genomes, and the minutes are spent inside skani itself, which reports its own progress. Co-Authored-By: Claude Opus 4.8 <[email protected]>
Coloring clusters is for telling neighbouring clusters apart, not for
identifying which cluster is which -- nobody reads a color off a dendrogram and
recovers "cluster 47". Giving every cluster its own shade off a colormap means
that past a few dozen clusters they are all indistinguishable anyway.
Cycle three UC Berkeley colors (Berkeley Blue, California Gold, Founders Rock)
in cluster order, so adjacent clusters always differ.
This also fixes two things that came along for the ride:
- colors were not reproducible. The old code did an *unseeded*
np.random.shuffle over a colormap, so the same analysis produced different
figures on every run.
- it used jet, which is perceptually non-uniform and misleading.
Cluster labels come in two shapes ('2' for primary, '2_10' for secondary), so
sorting is numeric per component: 2_2 sorts before 2_10, not after.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
Bump the version and date the changelog. This is a major version because it
breaks compatibility, not because of how much landed:
- --low_ram_primary_clustering was removed, so a v3 command line using it now
exits with "unrecognized arguments"
- Mdb.csv gains a column, becomes sparse rather than N^2, and omits genomes
with no above-threshold pairs; anything parsing it needs updating
- default clustering results change three ways (skani replaces MASH and
fastANI, primary linkage goes average -> single, and the skani
alignment-coverage fix cannot be disabled)
- skani is now required for default runs; mash no longer is
Any one of those is a semver MAJOR trigger on its own.
Also document a known limitation rather than let the release notes imply the
alignment-coverage story is finished. The fix in this release corrects skani's
coverage *units*; it does not change how a low-coverage pair is handled once
measured. cov_thresh still works by zeroing a pair's ANI before hierarchical
clustering, and average linkage routes around that: a genome joins a cluster on
the strength of its other relationships and is grouped with a member it never had
adequate coverage with. With cov_thresh=0.5, a genome aligning over 1% of another
still ends up clustered with it, and three genomes is enough to trigger it. That
is a separate, pre-existing problem and it is deferred.
Co-Authored-By: Claude Opus 4.8 <[email protected]>
Updated README to reflect changes in dRep v4 regarding genome comparisons.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.