V4 streaming primary

CHANGELOG.md

@@ -4,6 +4,139 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](http://keepachangelog.com/)
 and this project (attempts to) adhere to [Semantic Versioning](http://semver.org/).
 
+## [4.0.0] - 2026-07-15
+
+dRep v4 makes genome clustering scale. The headline change is that primary and
+secondary clustering now both default to skani, and run as a **single pass** over
+the data instead of comparing every genome twice.
+
+**This release changes results.** Read "Breaking changes" below before upgrading
+an existing analysis. If you need the old behavior, `--primary_algorithm MASH
+--S_algorithm fastANI --primary_clusterAlg average` gets close, but the skani
+alignment-coverage fix (see Fixed) cannot be turned off, and it was a genuine
+bug.
+
+### Breaking changes
+
+- **skani is now the default for both clustering steps** (`--primary_algorithm`
+  defaults to `skani`, was MASH; `--S_algorithm` defaults to `skani`, was
+  fastANI). skani must be installed. Mash is now only needed for
+  `--primary_algorithm MASH`, and the startup dependency check no longer demands
+  it otherwise.
+- **Primary clustering now defaults to single linkage.** Previously `--clusterAlg`
+  (default `average`) drove *both* clustering steps. Primary now has its own
+  `--primary_clusterAlg`, defaulting to `single`. Single linkage is the right
+  choice for a deliberately inclusive pre-filter, and it is what makes the
+  low-memory algorithm possible. `--clusterAlg` still controls secondary
+  clustering and still defaults to `average`.
+- **`Mdb.csv` means something different under `--primary_algorithm skani`.** It is
+  now a *sparse* table of real skani ANI values plus alignment coverage, holding
+  only pairs above skani's screening threshold — roughly 779k rows for 10,000
+  genomes, versus 100M rows of dense Mash distances. Genomes with no
+  above-threshold pairs do not appear in it at all. Anything parsing `Mdb.csv`
+  needs to account for this.
+- `--S_algorithm skani` results change; see the coverage fix under Fixed.
+- `--low_ram_primary_clustering` was removed (see Removed).
+
+### Added
+
+- `--primary_algorithm {skani,MASH}` — choose the primary clustering program.
+- `--primary_clusterAlg` — linkage method for primary clustering, independent of
+  the secondary `--clusterAlg`.
+- `--classic_primary_clustering` — force the pre-v4 dense scipy primary path.
+- `--primary_skani_min_af` — minimum aligned fraction for a pair to form a
+  primary-clustering edge (skani only, default 15).
+- `--no_reuse_primary_comparisons` — re-run skani during secondary clustering
+  rather than reusing primary's comparisons. A debugging escape hatch; reuse is
+  exact.
+
+### Changed
+
+- **Primary clustering no longer builds the N x N distance matrix.** Single-linkage
+  clustering at a fixed threshold is identical to finding connected components,
+  so it is now computed directly with union-find, removing the pivot-then-unpivot
+  RAM spike (issue #259). Memory is O(genomes + edges) instead of O(genomes^2).
+  On a synthetic 8,000-genome set, peak RAM for this step dropped from 7.42 GB to
+  0.55 GB; the old path grew quadratically while the new one stays flat. Together
+  with `--primary_algorithm skani`, this addresses the out-of-memory crashes
+  reported when clustering tens of thousands of genomes.
+- **Secondary clustering reuses primary's comparisons.** With skani for both
+  steps, dRep previously sketched every genome twice and computed the same ANI
+  values twice: once across all genomes, then again within each primary cluster.
+  Since secondary only compares genomes *within* a primary cluster, those pairs
+  are a subset of what primary already computed. On 10,000 UHGG genomes, 94% of
+  the pairs driving secondary clustering were already present with identical ANI
+  to 6 decimal places. dRep now runs skani once and derives both steps from it.
+  Secondary clustering went from 15 minutes to 24 seconds, producing an identical
+  partition; whole-pipeline `dereplicate` went from 22.4 to 13.8 minutes.
+- Primary clustering with skani writes far less to disk: 82 MB vs 15 GB of Mash
+  output for 10,000 genomes.
+- `--multiround_primary_clustering` and `--primary_chunksize` now warn that they
+  only apply to `--primary_algorithm MASH`. skani's sparse output never builds
+  the N x N table that multiround exists to avoid, and has none of multiround's
+  chunk-splitting imprecision.
+- `--SkipMash` help text clarified: it skips primary clustering whatever the
+  primary algorithm is. The name is historical.
+
+### Removed
+
+- `--low_ram_primary_clustering`. Union-find is now the default for single-linkage
+  primary clustering, so the flag had become a no-op.
+- **networkx is no longer a dependency.** It was only used by the connected-components
+  path behind `--low_ram_primary_clustering`.
+
+### Fixed
+
+- **skani alignment coverage was a percent, not a fraction (results-affecting).**
+  skani reports aligned fractions as 0-100, but `load_skani` only divided ANI by
+  100 and passed the aligned fraction through untouched. Every other algorithm
+  reports `alignment_coverage` on a 0-1 scale, which is the scale `cov_thresh` is
+  compared against, so the coverage filter was effectively inert for
+  `--S_algorithm skani`: a pair aligning over 1% of the genome had
+  `alignment_coverage=1.04` and sailed past a `cov_thresh` of 0.5. On the bundled
+  test genomes this merged *E. casseliflavus* with *E. faecalis* — two different
+  species — into one secondary cluster. `--S_algorithm skani` users should expect
+  different (more conservative) clusters as a result. Note this fixes the
+  *units*, not how low-coverage pairs are handled once measured; see Known
+  limitations.
+- `ScaffoldLevel_dRep.py` crashed on MUMmer 3 with "nucmer failed with exit code
+  1". It passed `-t` (threads) unconditionally, but that option only exists in
+  MUMmer 4, and MUMmer 3 rejects it rather than ignoring it. `conda install
+  mummer` still installs 3.23. The script now detects whether nucmer supports
+  `-t` and only passes it if so; MUMmer 4 keeps its threading.
+- The primary dendrogram is still produced for modest genome sets under the
+  streaming/sparse paths, which build no linkage matrix of their own. Above
+  `--primary_dendrogram_max_genomes` (2000) it is skipped, as multiround already
+  did.
+
+### Known limitations
+
+- **Alignment coverage still only filters pairs, not clusters.** `cov_thresh` is
+  applied by setting a low-coverage pair's ANI to 0 before hierarchical
+  clustering. That stops the pair itself from pulling two genomes together, but
+  average linkage can still route around it: a genome joins a cluster on the
+  strength of its *other* relationships and ends up grouped with a member it
+  never had adequate coverage with. Concretely, with `cov_thresh=0.5`, a genome
+  aligning over only 1% of another still lands in the same secondary cluster as
+  it — via their mutual neighbours — and so one of the two is discarded as
+  redundant. Three genomes is enough to trigger this.
+
+  This is a real and long-standing behavior, not a regression, and it is separate
+  from the skani units bug fixed above. Fixing it properly means validating
+  cluster membership after clustering (e.g. rejecting a genome that lacks
+  sufficient coverage with the cluster) rather than adjusting distances, and it
+  has to hold for hierarchical, greedy, and multiround paths alike. Deferred.
+
+### Validation
+
+v4 clustering was validated against 10,000 real genomes from the UHGG catalogue
+(24 GB, 1,248 species, 9,599 MAGs + 401 isolates). dRep independently recovered
+1,232 secondary clusters at `-sa 0.95`, against UHGG's own 1,248 species
+assignments: 2.4% of clusters spanned more than one UHGG species, and 2.9% of
+UHGG species were split across clusters. The one-pass path reproduced the
+two-pass result exactly — same primary clusters, same secondary clusters, same
+representative genomes.
+
 ## [3.7.1] - 2026-06-30
 - Fix crash when fewer than 2 genomes remain after filtering (issue #300)
 - Fix argument list bug (issue #288)

README.md

@@ -14,6 +14,19 @@ Publication is available at
 Open source pre-print publication is available at
 [bioRxiv](https://doi.org/10.1101/108142)
 
+## ⚡ New in v4
+
+dRep v4 uses [skani](https://github.com/bluenote-1577/skani) for **both** primary and secondary genome comparisons by default, replacing v3's default of MASH (primary) + fastANI (secondary). skani is much faster than that pair, and it *streams* its comparisons instead of building an all-vs-all matrix in memory — so `dereplicate` runs far quicker and its memory footprint grows roughly linearly with genome count rather than quadratically.
+
+**Whole-pipeline `dRep dereplicate` on 10,000 genomes** — identical inputs and settings:
+
+| | v3 defaults (MASH → fastANI) | v4 defaults (skani) |
+|---|---|---|
+| Wall-clock time | 6 h 15 min | **14.5 min** (~26× faster) |
+| Peak memory (RSS) | ~13 GB | ~7 GB |
+
+*Benchmarked on an Apple M1 Pro with `-p 10`; the memory advantage widens further at larger genome counts.*
+
 ## Installation with pip
 ```
 $ pip install drep
@@ -38,13 +51,14 @@ $ dRep check_dependencies
 
 ## Dependencies
 ### Near Essential
-* [Mash](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0997-x>) - Makes primary clusters (v1.1.1 confirmed works)
-* [MUMmer](http://mummer.sourceforge.net/) - Performs default ANIm comparison method (v3.23 confirmed works)
+* [skani](https://github.com/bluenote-1577/skani) - Makes primary clusters and performs the default secondary comparison (v0.2+ confirmed works)
+* [CheckM](http://ecogenomics.github.io/CheckM/) - Determines contamination and completeness of genomes (v1.0.7 confirmed works). Only needed for `dereplicate`; skip it with `--genomeInfo` or `--ignoreGenomeQuality`
 
 ### Optional
 
-* [fastANI](https://github.com/ParBLiSS/FastANI) - A fast secondary clustering algorithm
-* [CheckM](http://ecogenomics.github.io/CheckM/)_ - Determines contamination and completeness of genomes (v1.0.7 confirmed works)
+* [Mash](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0997-x>) - Only needed for `--primary_algorithm MASH` (v1.1.1 confirmed works)
+* [MUMmer](http://mummer.sourceforge.net/) - Only needed for the ANIm comparison methods (v3.23 confirmed works)
+* [fastANI](https://github.com/ParBLiSS/FastANI) - An alternative fast secondary clustering algorithm
 * [gANI (aka ANIcalculator)](https://ani.jgi-psf.org/html/download.php?) - Performs gANI comparison method (v1.0 confirmed works)
 * [Prodigal](http://prodigal.ornl.gov/) - Used be both checkM and gANI (v2.6.3 confirmed works)
 * [NSimScan](https://pubmed.ncbi.nlm.nih.gov/27153714/) - Only needed for goANI algorithm (open source version of gANI)

docs/choosing_parameters.rst

@@ -154,15 +154,22 @@ dRep can use any method of linkage listed at the following webpage by using the
 7. Overview of genome comparison algorithms
 ----------------------------------------------
 
-**Primary clustering** is always performed with `Mash <https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0997-x>`_; an extremely fast but somewhat inaccurate algorithm.
+**Primary clustering** groups genomes that could plausibly be "the same", so that the more accurate secondary algorithm only has to run within those groups. Two programs are supported:
 
-There are several supported **secondary clustering algorithms**. These calculate the accurate Average Nucleotide Identity (ANI) between genomes that is used to cluster genomes into secondary clusters. The following algorithms are currently supported as of version 3:
+* **skani** (DEFAULT as of v4) (`Shaw 2023 <https://doi.org/10.1038/s41592-023-02018-3>`_). Run as ``skani triangle --sparse``, which only ever emits pairs above its screening threshold. It therefore never builds the full N x N table, which is what made older versions of dRep run out of memory on large genome sets. It is also more accurate than Mash near the clustering threshold.
+* **MASH** (`Ondov 2016 <https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0997-x>`_). The pre-v4 behavior; extremely fast but somewhat inaccurate, and it computes and stores all N x N comparisons.
 
-* **ANIn** (`Richter 2009 <https://www.ncbi.nlm.nih.gov/pubmed/19855009>`_). This aligns whole genomes with nucmer and compares the aligned regions.
-* **ANImf** (DEFAULT). This is the same as ANIn, but filters the alignments such that each region of genome 1 and only align to a single region of genome 2. This takes slightly more time, but is much more accurate on genomes with repeat regions
-* **gANI** (`Varghese 2015 <https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4538840/>`_). This aligns genes (ORFs) called by Prodigal instead of aligning whole genomes. This algorithm is a bit faster than ANIm-based algorithms, but only aligns coding regions.
+There are several supported **secondary clustering algorithms**. These calculate the accurate Average Nucleotide Identity (ANI) between genomes that is used to cluster genomes into secondary clusters. The following algorithms are currently supported:
+
+* **skani** (DEFAULT as of v4) (`Shaw 2023 <https://doi.org/10.1038/s41592-023-02018-3>`_). Fast and accurate, including on incomplete genomes. When paired with ``--primary_algorithm skani`` (the default), secondary clustering reuses the comparisons already computed during primary clustering rather than recomputing them, which makes the secondary stage nearly free.
+* **FastANI** (`Jain 2018 <https://doi.org/10.1038/s41467-018-07641-9>`_). A really fast Mash-based algorithm that can also handle incomplete genomes. Seems to be just as accurate as alignment-based algorithms. Was the default in v3.
+* **ANImf**. This is the same as ANIn, but filters the alignments such that each region of genome 1 can only align to a single region of genome 2. This takes slightly more time, but is much more accurate on genomes with repeat regions. Was the default in earlier versions.
+* **ANIn** (`Richter 2009 <https://www.ncbi.nlm.nih.gov/pubmed/19855009>`_). This aligns whole genomes with nucmer and compares the aligned regions. ANImf is strictly better and should be preferred.
+* **gANI** (`Varghese 2015 <https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4538840/>`_). This aligns genes (ORFs) called by Prodigal instead of aligning whole genomes. This algorithm is a bit faster than ANIm-based algorithms, but only aligns coding regions. Requires the ANIcalculator program.
 * **goANI**. This is my own open-source implementation of gANI, which is not open source (and for which the authors would not share the source code when asked). I wrote this algorithm so that I could calculate dN/dS between aligned genes for `this study <https://msystems.asm.org/content/5/1/e00731-19>`_ (you can too using `dnds_from_drep.py <https://github.com/MrOlm/bacterialEvolutionMetrics>`_). Requires the program `NSimScan <https://pubmed.ncbi.nlm.nih.gov/27153714/>`_.
-* **FastANI** (`Jain 2018 <https://doi.org/10.1038/s41467-018-07641-9>`_). A really fast Mash-based algorithm that can also handle incomplete genomes. Seems to be just as accurate as alignment-based algorithms. **Should probably be the default algorithm when you care about runtime.***
+
+.. note::
+  **A note on skani and alignment coverage.** Mash compares k-mers across the whole genome, so two genomes that share only a small conserved region look distant. skani instead reports the identity *within the aligned regions only*, so that same pair can report a high ANI. This is why the skani primary path requires a minimum aligned fraction (``--primary_skani_min_af``, default 15%) before a pair counts as a primary-clustering edge. Without it, a handful of genomes sharing small conserved regions chain unrelated organisms together under single linkage. Lower it only if you have very fragmented genomes and understand that risk.
 
 .. note::
   None of these algorithms are perfect, especially in repeat-prone genomes. Regions of the genome which are not homologous can align to each other and artificially decrease ANI. In fact, when a genome is compared to itself, the algorithms often reports values <100% for this reason.

docs/installation.rst

@@ -38,13 +38,14 @@ To check which dependencies are installed on your system and accessible by dRep,
 
 **Near Essential**
 
-* `Mash <https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0997-x>`_ - Makes primary clusters (v1.1.1 confirmed works)
-* `MUMmer <http://mummer.sourceforge.net/>`_ - Performs default ANIm comparison method (v3.23 confirmed works)
+* `skani <https://github.com/bluenote-1577/skani>`_ - Makes primary clusters and performs the default secondary comparison (v0.2+ confirmed works)
+* `CheckM <http://ecogenomics.github.io/CheckM/>`_ - Determines contamination and completeness of genomes (v1.0.7 confirmed works). Only needed for ``dereplicate``; you can skip it with ``--genomeInfo`` or ``--ignoreGenomeQuality``
 
 **Recommended**
 
-* `fastANI <https://github.com/ParBLiSS/FastANI>`_ - A fast secondary clustering algorithm
-* `CheckM <http://ecogenomics.github.io/CheckM/>`_ - Determines contamination and completeness of genomes (v1.0.7 confirmed works)
+* `Mash <https://genomebiology.biomedcentral.com/articles/10.1186/s13059-016-0997-x>`_ - Only needed for ``--primary_algorithm MASH`` (v1.1.1 confirmed works)
+* `MUMmer <http://mummer.sourceforge.net/>`_ - Only needed for the ANIm comparison methods (v3.23 confirmed works)
+* `fastANI <https://github.com/ParBLiSS/FastANI>`_ - An alternative fast secondary clustering algorithm
 * `gANI (aka ANIcalculator) - Performs gANI comparison method (v1.0 confirmed works)
 * `Prodigal <https://github.com/hyattpd/Prodigal>`_ - Used be both checkM and gANI (v2.6.3 confirmed works)