diff --git a/CHANGELOG.md b/CHANGELOG.md index 6840f45..e6abe9c 100644 --- a/CHANGELOG.md +++ b/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) diff --git a/README.md b/README.md index 42f3839..659a448 100644 --- a/README.md +++ b/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) diff --git a/docs/choosing_parameters.rst b/docs/choosing_parameters.rst index c91476e..3f4f1d6 100644 --- a/docs/choosing_parameters.rst +++ b/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 `_; 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 `_). 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 `_). The pre-v4 behavior; extremely fast but somewhat inaccurate, and it computes and stores all N x N comparisons. -* **ANIn** (`Richter 2009 `_). 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 `_). 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 `_). 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 `_). 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 `_). This aligns whole genomes with nucmer and compares the aligned regions. ANImf is strictly better and should be preferred. +* **gANI** (`Varghese 2015 `_). 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 `_ (you can too using `dnds_from_drep.py `_). Requires the program `NSimScan `_. -* **FastANI** (`Jain 2018 `_). 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. diff --git a/docs/installation.rst b/docs/installation.rst index 8c5ca70..e4e4479 100644 --- a/docs/installation.rst +++ b/docs/installation.rst @@ -38,13 +38,14 @@ To check which dependencies are installed on your system and accessible by dRep, **Near Essential** -* `Mash `_ - Makes primary clusters (v1.1.1 confirmed works) -* `MUMmer `_ - Performs default ANIm comparison method (v3.23 confirmed works) +* `skani `_ - Makes primary clusters and performs the default secondary comparison (v0.2+ confirmed works) +* `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 `_ - A fast secondary clustering algorithm -* `CheckM `_ - Determines contamination and completeness of genomes (v1.0.7 confirmed works) +* `Mash `_ - Only needed for ``--primary_algorithm MASH`` (v1.1.1 confirmed works) +* `MUMmer `_ - Only needed for the ANIm comparison methods (v3.23 confirmed works) +* `fastANI `_ - An alternative fast secondary clustering algorithm * `gANI (aka ANIcalculator) - Performs gANI comparison method (v1.0 confirmed works) * `Prodigal `_ - Used be both checkM and gANI (v2.6.3 confirmed works) diff --git a/docs/module_descriptions.rst b/docs/module_descriptions.rst index 227d199..9cdcab4 100644 --- a/docs/module_descriptions.rst +++ b/docs/module_descriptions.rst @@ -5,19 +5,20 @@ dRep has 3 commands: compare, dereplicate, and check dependencies. To see a list $ dRep -h - ...::: dRep v3.0.0 :::... - Matt Olm. MIT License. Banfield Lab, UC Berkeley. 2017 (last updated 2020) + ...::: dRep v4.0.0 :::... - See https://drep.readthedocs.io/en/latest/index.html for documentation - Choose one of the operations below for more detailed help. + Matt Olm. MIT License. Banfield Lab, UC Berkeley. 2017 (last updated 2026) - Example: dRep dereplicate -h + See https://drep.readthedocs.io/en/latest/index.html for documentation + Choose one of the operations below for more detailed help. - Commands: - compare -> Compare and cluster a set of genomes - dereplicate -> De-replicate a set of genomes - check_dependencies -> Check which dependencies are properly installed + Example: dRep dereplicate -h + + Commands: + compare -> Compare and cluster a set of genomes + dereplicate -> De-replicate a set of genomes + check_dependencies -> Check which dependencies are properly installed In previous versions of dRep (everything before v3) the user could run a number of additional modules separately, but now they can only be run as part of the larger workflows `compare` and `dereplicate`. Many of the modules are the same for `compare` and `dereplicate`, however, and in cases where these is the same parameter in both it functions exactly the same in each. @@ -41,21 +42,29 @@ Compare This workflow compares a set of genomes. For a list of all parameters, check the help:: $ dRep compare -h - usage: dRep compare [-p PROCESSORS] [-d] [-h] [-g [GENOMES [GENOMES ...]]] - [--S_algorithm {fastANI,gANI,goANI,ANIn,ANImf}] + + usage: dRep compare [-p PROCESSORS] [-d] [-h] [-g [GENOMES ...]] + [--S_algorithm {ANImf,fastANI,skani,ANIn,gANI,goANI}] + [--primary_algorithm {skani,MASH}] + [--no_reuse_primary_comparisons] + [--primary_skani_min_af PRIMARY_SKANI_MIN_AF] [-ms MASH_SKETCH] [--SkipMash] [--SkipSecondary] - [--n_PRESET {normal,tight}] [-pa P_ANI] [-sa S_ANI] - [-nc COV_THRESH] [-cm {total,larger}] - [--clusterAlg {median,weighted,single,complete,average,ward,centroid}] + [--skani_extra SKANI_EXTRA] [--n_PRESET {normal,tight}] + [-pa P_ANI] [-sa S_ANI] [-nc COV_THRESH] + [-cm {total,larger}] + [--clusterAlg {single,median,ward,centroid,complete,average,weighted}] + [--primary_clusterAlg {single,median,ward,centroid,complete,average,weighted}] + [--classic_primary_clustering] [--multiround_primary_clustering] [--primary_chunksize PRIMARY_CHUNKSIZE] [--greedy_secondary_clustering] - [--run_tertiary_clustering] [--warn_dist WARN_DIST] - [--warn_sim WARN_SIM] [--warn_aln WARN_ALN] + [--run_tertiary_clustering] [--gen_warnings] + [--warn_dist WARN_DIST] [--warn_sim WARN_SIM] + [--warn_aln WARN_ALN] work_directory positional arguments: - work_directory Directory where data and output are stored + work_directory Directory where data and output are stored *** USE THE SAME WORK DIRECTORY FOR ALL DREP OPERATIONS *** SYSTEM PARAMETERS: @@ -65,27 +74,60 @@ This workflow compares a set of genomes. For a list of all parameters, check the -h, --help show this help message and exit GENOME INPUT: - -g [GENOMES [GENOMES ...]], --genomes [GENOMES [GENOMES ...]] + -g [GENOMES ...], --genomes [GENOMES ...] genomes to filter in .fasta format. Not necessary if Bdb or Wdb already exist. Can also input a text file with paths to genomes, which results in fewer OS issues than wildcard expansion (default: None) GENOME COMPARISON OPTIONS: - --S_algorithm {fastANI,gANI,goANI,ANIn,ANImf} + --S_algorithm {ANImf,fastANI,skani,ANIn,gANI,goANI} Algorithm for secondary clustering comaprisons: + skani = (DEFAULT) Kmer-based approach; fastest and most accurate. + When paired with --primary_algorithm skani, secondary reuses + the comparisons already done during primary clustering. fastANI = Kmer-based approach; very fast - ANImf = (DEFAULT) Align whole genomes with nucmer; filter alignment; compare aligned regions + ANImf = Align whole genomes with nucmer; filter alignment; compare aligned regions ANIn = Align whole genomes with nucmer; compare aligned regions gANI = Identify and align ORFs; compare aligned ORFS goANI = Open source version of gANI; requires nsmimscan - (default: ANImf) + (default: skani) + --primary_algorithm {skani,MASH} + Program to use for primary clustering. + skani = (DEFAULT) skani triangle --sparse. Only above-threshold pairs + are produced, so there is no N^2 matrix in RAM or on disk, and + a skani --S_algorithm can reuse these comparisons instead of + recomputing them. + MASH = all-vs-all Mash. Pre-v4 behavior; builds the full N^2 table. (default: skani) + --no_reuse_primary_comparisons + Re-run skani during secondary clustering instead of + reusing the comparisons already computed during + primary clustering. Only relevant with + --primary_algorithm skani and a skani --S_algorithm, + where the two stages otherwise compute the same ANI + values twice. Reuse is exact, so this is mostly a + debugging escape hatch. (default: True) + --primary_skani_min_af PRIMARY_SKANI_MIN_AF + Minimum percent of a genome that must align for a pair + to form a primary-clustering edge (--primary_algorithm + skani only). skani's ANI is measured within aligned + regions only, so without this filter genomes sharing + just a small conserved region become edges and single + linkage chains them into one huge cluster. The default + reproduces the MASH partition closely; lower it only + if you have very fragmented genomes and understand the + chaining risk. (default: 15) -ms MASH_SKETCH, --MASH_sketch MASH_SKETCH MASH sketch size (default: 1000) - --SkipMash Skip MASH clustering, just do secondary clustering on - all genomes (default: False) + --SkipMash Skip primary clustering entirely and run secondary + clustering on all genomes at once. (Named for when + primary clustering was always MASH; it applies to + whichever --primary_algorithm is in use.) (default: + False) --SkipSecondary Skip secondary clustering, just perform MASH clustering (default: False) + --skani_extra SKANI_EXTRA + Extra arguments to pass to skani triangle (default: ) --n_PRESET {normal,tight} Presets to pass to nucmer tight = only align highly conserved regions @@ -107,9 +149,22 @@ This workflow compares a set of genomes. For a list of all parameters, check the total = 2*(aligned length) / (sum of total genome lengths) larger = max((aligned length / genome 1), (aligned_length / genome2)) (default: larger) - --clusterAlg {median,weighted,single,complete,average,ward,centroid} - Algorithm used to cluster genomes (passed to - scipy.cluster.hierarchy.linkage (default: average) + --clusterAlg {single,median,ward,centroid,complete,average,weighted} + Algorithm used to cluster genomes during SECONDARY + clustering (passed to scipy.cluster.hierarchy.linkage) + (default: average) + --primary_clusterAlg {single,median,ward,centroid,complete,average,weighted} + Algorithm used to cluster genomes during PRIMARY + (MASH/skani) clustering. The default 'single' is equivalent to connected + components and is computed with a fast, low-memory streaming algorithm that + scales to very large genome sets. Any other choice falls back to the classic + dense scipy path (see --classic_primary_clustering). (default: single) + --classic_primary_clustering + Force the classic dense (scipy) primary clustering + path instead of the streaming single-linkage + algorithm. Uses much more RAM at scale but reproduces + pre-v4 behavior and allows non-single linkage methods + and the primary dendrogram plot. (default: False) GREEDY CLUSTERING OPTIONS These decrease RAM use and runtime at the expense of a minor loss in accuracy. @@ -126,11 +181,6 @@ This workflow compares a set of genomes. For a list of all parameters, check the Impacts multiround_primary_clustering. If you have more than this many genomes, process them in chunks of this size. (default: 5000) - --low_ram_primary_clustering - Use a memory-efficient algorithm for primary clustering. - This only affects primary clustering and not secondary - clustering. Can be combined with multiround_primary_clustering - for even greater memory efficiency. (default: False) --greedy_secondary_clustering Use a heuristic to avoid pair-wise comparisons when doing secondary clustering. Will be done with single @@ -145,6 +195,7 @@ This workflow compares a set of genomes. For a list of all parameters, check the False) WARNINGS: + --gen_warnings Generate warnings (default: False) --warn_dist WARN_DIST How far from the threshold to throw cluster warnings (default: 0.25) @@ -161,17 +212,25 @@ Dereplicate This workflow dereplicates a set of genomes. For a list of all parameters, check the help:: - $ dRep dereplicate -h - usage: dRep dereplicate [-p PROCESSORS] [-d] [-h] [-g [GENOMES [GENOMES ...]]] + $ dRep dereplicate -h + + usage: dRep dereplicate [-p PROCESSORS] [-d] [-h] [-g [GENOMES ...]] [-l LENGTH] [-comp COMPLETENESS] [-con CONTAMINATION] [--ignoreGenomeQuality] [--genomeInfo GENOMEINFO] - [--checkM_method {taxonomy_wf,lineage_wf}] + [--checkM_method {lineage_wf,taxonomy_wf}] [--set_recursion SET_RECURSION] - [--S_algorithm {goANI,ANIn,gANI,ANImf,fastANI}] + [--checkm_group_size CHECKM_GROUP_SIZE] + [--S_algorithm {goANI,ANImf,gANI,skani,ANIn,fastANI}] + [--primary_algorithm {skani,MASH}] + [--no_reuse_primary_comparisons] + [--primary_skani_min_af PRIMARY_SKANI_MIN_AF] [-ms MASH_SKETCH] [--SkipMash] [--SkipSecondary] + [--skani_extra SKANI_EXTRA] [--n_PRESET {normal,tight}] [-pa P_ANI] [-sa S_ANI] [-nc COV_THRESH] [-cm {total,larger}] - [--clusterAlg {single,ward,complete,weighted,centroid,median,average}] + [--clusterAlg {average,complete,weighted,centroid,single,ward,median}] + [--primary_clusterAlg {average,complete,weighted,centroid,single,ward,median}] + [--classic_primary_clustering] [--multiround_primary_clustering] [--primary_chunksize PRIMARY_CHUNKSIZE] [--greedy_secondary_clustering] @@ -180,12 +239,13 @@ This workflow dereplicates a set of genomes. For a list of all parameters, check [-conW CONTAMINATION_WEIGHT] [-strW STRAIN_HETEROGENEITY_WEIGHT] [-N50W N50_WEIGHT] [-sizeW SIZE_WEIGHT] [-centW CENTRALITY_WEIGHT] + [-extraW EXTRA_WEIGHT_TABLE] [--gen_warnings] [--warn_dist WARN_DIST] [--warn_sim WARN_SIM] - [--warn_aln WARN_ALN] + [--warn_aln WARN_ALN] [--skip_plots] work_directory positional arguments: - work_directory Directory where data and output are stored + work_directory Directory where data and output are stored *** USE THE SAME WORK DIRECTORY FOR ALL DREP OPERATIONS *** SYSTEM PARAMETERS: @@ -195,7 +255,7 @@ This workflow dereplicates a set of genomes. For a list of all parameters, check -h, --help show this help message and exit GENOME INPUT: - -g [GENOMES [GENOMES ...]], --genomes [GENOMES [GENOMES ...]] + -g [GENOMES ...], --genomes [GENOMES ...] genomes to filter in .fasta format. Not necessary if Bdb or Wdb already exist. Can also input a text file with paths to genomes, which results in fewer OS @@ -205,7 +265,7 @@ This workflow dereplicates a set of genomes. For a list of all parameters, check -l LENGTH, --length LENGTH Minimum genome length (default: 50000) -comp COMPLETENESS, --completeness COMPLETENESS - Minumum genome completeness (default: 75) + Minimum genome completeness (default: 75) -con CONTAMINATION, --contamination CONTAMINATION Maximum genome contamination (default: 25) @@ -218,12 +278,13 @@ This workflow dereplicates a set of genomes. For a list of all parameters, check on length and N50 (default: False) --genomeInfo GENOMEINFO location of .csv file containing quality information - on the genomes. Must contain: ["genome"(basename of - .fasta file of that genome), "completeness"(0-100 - value for completeness of the genome), - "contamination"(0-100 value of the contamination of - the genome)] (default: None) - --checkM_method {taxonomy_wf,lineage_wf} + on the genomes. Must contain: ["genome"(filename of + .fasta file of that genome, including extension e.g. + genome.fasta), "completeness"(0-100 value for + completeness of the genome), "contamination"(0-100 + value of the contamination of the genome)] (default: + None) + --checkM_method {lineage_wf,taxonomy_wf} Either lineage_wf (more accurate) or taxonomy_wf (faster) (default: lineage_wf) --set_recursion SET_RECURSION @@ -231,22 +292,59 @@ This workflow dereplicates a set of genomes. For a list of all parameters, check unless checkM is crashing due to recursion issues. Recommended to set to 2000 if needed, but setting this could crash python (default: 0) + --checkm_group_size CHECKM_GROUP_SIZE + The number of genomes passed to checkM at a time. + Increasing this increases RAM but makes checkM faster + (default: 2000) GENOME COMPARISON OPTIONS: - --S_algorithm {goANI,ANIn,gANI,ANImf,fastANI} + --S_algorithm {goANI,ANImf,gANI,skani,ANIn,fastANI} Algorithm for secondary clustering comaprisons: + skani = (DEFAULT) Kmer-based approach; fastest and most accurate. + When paired with --primary_algorithm skani, secondary reuses + the comparisons already done during primary clustering. fastANI = Kmer-based approach; very fast - ANImf = (DEFAULT) Align whole genomes with nucmer; filter alignment; compare aligned regions + ANImf = Align whole genomes with nucmer; filter alignment; compare aligned regions ANIn = Align whole genomes with nucmer; compare aligned regions gANI = Identify and align ORFs; compare aligned ORFS goANI = Open source version of gANI; requires nsmimscan - (default: ANImf) + (default: skani) + --primary_algorithm {skani,MASH} + Program to use for primary clustering. + skani = (DEFAULT) skani triangle --sparse. Only above-threshold pairs + are produced, so there is no N^2 matrix in RAM or on disk, and + a skani --S_algorithm can reuse these comparisons instead of + recomputing them. + MASH = all-vs-all Mash. Pre-v4 behavior; builds the full N^2 table. (default: skani) + --no_reuse_primary_comparisons + Re-run skani during secondary clustering instead of + reusing the comparisons already computed during + primary clustering. Only relevant with + --primary_algorithm skani and a skani --S_algorithm, + where the two stages otherwise compute the same ANI + values twice. Reuse is exact, so this is mostly a + debugging escape hatch. (default: True) + --primary_skani_min_af PRIMARY_SKANI_MIN_AF + Minimum percent of a genome that must align for a pair + to form a primary-clustering edge (--primary_algorithm + skani only). skani's ANI is measured within aligned + regions only, so without this filter genomes sharing + just a small conserved region become edges and single + linkage chains them into one huge cluster. The default + reproduces the MASH partition closely; lower it only + if you have very fragmented genomes and understand the + chaining risk. (default: 15) -ms MASH_SKETCH, --MASH_sketch MASH_SKETCH MASH sketch size (default: 1000) - --SkipMash Skip MASH clustering, just do secondary clustering on - all genomes (default: False) + --SkipMash Skip primary clustering entirely and run secondary + clustering on all genomes at once. (Named for when + primary clustering was always MASH; it applies to + whichever --primary_algorithm is in use.) (default: + False) --SkipSecondary Skip secondary clustering, just perform MASH clustering (default: False) + --skani_extra SKANI_EXTRA + Extra arguments to pass to skani triangle (default: ) --n_PRESET {normal,tight} Presets to pass to nucmer tight = only align highly conserved regions @@ -268,9 +366,22 @@ This workflow dereplicates a set of genomes. For a list of all parameters, check total = 2*(aligned length) / (sum of total genome lengths) larger = max((aligned length / genome 1), (aligned_length / genome2)) (default: larger) - --clusterAlg {single,ward,complete,weighted,centroid,median,average} - Algorithm used to cluster genomes (passed to - scipy.cluster.hierarchy.linkage (default: average) + --clusterAlg {average,complete,weighted,centroid,single,ward,median} + Algorithm used to cluster genomes during SECONDARY + clustering (passed to scipy.cluster.hierarchy.linkage) + (default: average) + --primary_clusterAlg {average,complete,weighted,centroid,single,ward,median} + Algorithm used to cluster genomes during PRIMARY + (MASH/skani) clustering. The default 'single' is equivalent to connected + components and is computed with a fast, low-memory streaming algorithm that + scales to very large genome sets. Any other choice falls back to the classic + dense scipy path (see --classic_primary_clustering). (default: single) + --classic_primary_clustering + Force the classic dense (scipy) primary clustering + path instead of the streaming single-linkage + algorithm. Uses much more RAM at scale but reproduces + pre-v4 behavior and allows non-single linkage methods + and the primary dendrogram plot. (default: False) GREEDY CLUSTERING OPTIONS These decrease RAM use and runtime at the expense of a minor loss in accuracy. @@ -287,11 +398,6 @@ This workflow dereplicates a set of genomes. For a list of all parameters, check Impacts multiround_primary_clustering. If you have more than this many genomes, process them in chunks of this size. (default: 5000) - --low_ram_primary_clustering - Use a memory-efficient algorithm for primary clustering. - This only affects primary clustering and not secondary - clustering. Can be combined with multiround_primary_clustering - for even greater memory efficiency. (default: False) --greedy_secondary_clustering Use a heuristic to avoid pair-wise comparisons when doing secondary clustering. Will be done with single @@ -306,7 +412,7 @@ This workflow dereplicates a set of genomes. For a list of all parameters, check False) SCORING CRITERIA - Based off of the formula: + Based off of the formula: A*Completeness - B*Contamination + C*(Contamination * (strain_heterogeneity/100)) + D*log(N50) + E*log(size) + F*(centrality - S_ani) A = completeness_weight; B = contamination_weight; C = strain_heterogeneity_weight; D = N50_weight; E = size_weight; F = cent_weight: @@ -322,8 +428,13 @@ This workflow dereplicates a set of genomes. For a list of all parameters, check weight of log(genome size) (default: 0) -centW CENTRALITY_WEIGHT, --centrality_weight CENTRALITY_WEIGHT Weight of (centrality - S_ani) (default: 1) + -extraW EXTRA_WEIGHT_TABLE, --extra_weight_table EXTRA_WEIGHT_TABLE + Path to a tab-separated file with two-columns, no + headers, listing genome and extra score to apply to + that genome (default: None) WARNINGS: + --gen_warnings Generate warnings (default: False) --warn_dist WARN_DIST How far from the threshold to throw cluster warnings (default: 0.25) @@ -332,6 +443,9 @@ This workflow dereplicates a set of genomes. For a list of all parameters, check --warn_aln WARN_ALN Minimum aligned fraction for warnings between dereplicated genomes (ANIn) (default: 0.25) + ANALYZE: + --skip_plots Dont make plots (default: False) + Example: dRep dereplicate output_dir/ -g /path/to/genomes/*.fasta Work Directory diff --git a/docs/overview.rst b/docs/overview.rst index 4d70672..0a65b20 100644 --- a/docs/overview.rst +++ b/docs/overview.rst @@ -13,9 +13,13 @@ Genome comparison dRep can rapidly and accurately compare a list of genomes in a pair-wise manner. This allows identification of groups of organisms that share similar DNA content in terms of Average Nucleotide Identity (ANI). -dRep performs this in two steps- first with a rapid primary algorithm (Mash), and second with a more sensitive algorithm (ANIm). We can't just use Mash because, while incredibly fast, it is not robust to genome incompletenss (see :doc:`choosing_parameters`) and only provides an "estimate" of ANI. ANIm is robust to genome incompleteness and is more accurate, but too slow to perform pair-wise comparisons of longer genome lists. +dRep performs this in two steps. **Primary clustering** groups genomes that could plausibly be the same organism, using a permissive threshold and single-linkage (connected components). **Secondary clustering** then runs within each of those groups, using average-linkage hierarchical clustering to decide which genomes actually are the same. -dRep first compares all genomes using Mash, and then only runs the secondary algorithm (ANIm or gANI) on sets of genomes that have at least 90% Mash ANI. This results in a great decrease in the number of (slow) secondary comparisons that need to be run while maintaining the sensitivity of ANIm. +Both steps matter. Primary clustering is deliberately inclusive, and it partitions the genomes so that the more expensive secondary clustering only ever has to consider a small group at a time. Secondary clustering uses average linkage, which is what stops a chain of similar-but-distinct strains from collapsing distinct species into one cluster. + +As of v4 both steps default to `skani `_. Primary clustering runs it in ``--sparse`` mode, which only emits pairs above a screening threshold, so dRep never builds or stores the full N x N comparison table. Because secondary clustering only compares genomes *within* a primary cluster, those pairs are a subset of what the primary pass already computed, so it reuses them rather than running the comparisons a second time. + +Older versions used Mash for primary clustering and an alignment-based algorithm (ANIm) for secondary. That path is still available (``--primary_algorithm MASH``, ``--S_algorithm ANImf``), but it computes all N x N Mash comparisons up front, which is slower and uses far more memory on large genome sets. .. See the `publication `_ for details diff --git a/drep/VERSION b/drep/VERSION index a76ccff..fcdb2e1 100644 --- a/drep/VERSION +++ b/drep/VERSION @@ -1 +1 @@ -3.7.1 +4.0.0 diff --git a/drep/argumentParser.py b/drep/argumentParser.py index 3e74ce4..7d7d4d7 100644 --- a/drep/argumentParser.py +++ b/drep/argumentParser.py @@ -114,16 +114,43 @@ def parse_args(args): Clustflags = cluster_parent.add_argument_group('GENOME COMPARISON OPTIONS') Clustflags.add_argument("--S_algorithm", help="R|Algorithm for secondary clustering comaprisons:\n" \ + + "skani = (DEFAULT) Kmer-based approach; fastest and most accurate.\n" \ + + " When paired with --primary_algorithm skani, secondary reuses\n" \ + + " the comparisons already done during primary clustering.\n" \ + "fastANI = Kmer-based approach; very fast\n" \ - + "skani = Even faster Kmer-based approacht\n" \ - + "ANImf = (DEFAULT) Align whole genomes with nucmer; filter alignment; compare aligned regions\n" \ + + "ANImf = Align whole genomes with nucmer; filter alignment; compare aligned regions\n" \ + "ANIn = Align whole genomes with nucmer; compare aligned regions\n" \ + "gANI = Identify and align ORFs; compare aligned ORFS\n" \ + "goANI = Open source version of gANI; requires nsmimscan\n", - default='fastANI', choices={'ANIn', 'gANI', 'ANImf', 'goANI', 'fastANI', 'skani'}) + default='skani', choices={'ANIn', 'gANI', 'ANImf', 'goANI', 'fastANI', 'skani'}) + Clustflags.add_argument("--primary_algorithm", help="R|Program to use for primary clustering.\n" \ + + "skani = (DEFAULT) skani triangle --sparse. Only above-threshold pairs\n" \ + + " are produced, so there is no N^2 matrix in RAM or on disk, and\n" \ + + " a skani --S_algorithm can reuse these comparisons instead of\n" \ + + " recomputing them.\n" \ + + "MASH = all-vs-all Mash. Pre-v4 behavior; builds the full N^2 table.", + default='skani', choices={'MASH', 'skani'}) + Clustflags.add_argument("--no_reuse_primary_comparisons", dest='reuse_primary_comparisons', + help="Re-run skani during secondary clustering instead of reusing the " + "comparisons already computed during primary clustering. Only " + "relevant with --primary_algorithm skani and a skani --S_algorithm, " + "where the two stages otherwise compute the same ANI values twice. " + "Reuse is exact, so this is mostly a debugging escape hatch.", + action='store_false', default=True) + Clustflags.add_argument("--primary_skani_min_af", + help="Minimum percent of a genome that must align for a pair to form a " + "primary-clustering edge (--primary_algorithm skani only). skani's ANI " + "is measured within aligned regions only, so without this filter genomes " + "sharing just a small conserved region become edges and single linkage " + "chains them into one huge cluster. The default reproduces the MASH " + "partition closely; lower it only if you have very fragmented genomes " + "and understand the chaining risk.", + default=15, type=float) Clustflags.add_argument("-ms", "--MASH_sketch", help="MASH sketch size", default=1000) - Clustflags.add_argument("--SkipMash", help="Skip MASH clustering,\ - just do secondary clustering on all genomes", action='store_true') + Clustflags.add_argument("--SkipMash", help="Skip primary clustering entirely and run secondary\ + clustering on all genomes at once. (Named for when primary clustering was\ + always MASH; it applies to whichever --primary_algorithm is in use.)", + action='store_true') Clustflags.add_argument("--SkipSecondary", help="Skip secondary clustering, just perform MASH\ clustering", action='store_true') Clustflags.add_argument("--skani_extra", @@ -146,10 +173,19 @@ def parse_args(args): + "total = 2*(aligned length) / (sum of total genome lengths)\n" \ + "larger = max((aligned length / genome 1), (aligned_length / genome2))\n", choices=['total', 'larger'], default='larger') - Compflags.add_argument("--clusterAlg", help="Algorithm used to cluster genomes (passed\ - to scipy.cluster.hierarchy.linkage", default='average', + Compflags.add_argument("--clusterAlg", help="Algorithm used to cluster genomes during SECONDARY\ + clustering (passed to scipy.cluster.hierarchy.linkage)", default='average', + choices={'single', 'complete', 'average', 'weighted', 'centroid', 'median', 'ward'}) + Compflags.add_argument("--primary_clusterAlg", help="R|Algorithm used to cluster genomes during PRIMARY\n" \ + "(MASH/skani) clustering. The default 'single' is equivalent to connected\n" \ + "components and is computed with a fast, low-memory streaming algorithm that\n" \ + "scales to very large genome sets. Any other choice falls back to the classic\n" \ + "dense scipy path (see --classic_primary_clustering).", default='single', choices={'single', 'complete', 'average', 'weighted', 'centroid', 'median', 'ward'}) - Compflags.add_argument("--low_ram_primary_clustering", help="Use a memory-efficient algorithm for primary clustering. This only affects primary clustering and not secondary clustering.", + Compflags.add_argument("--classic_primary_clustering", help="Force the classic dense (scipy) primary\ + clustering path instead of the streaming single-linkage algorithm. Uses much more\ + RAM at scale but reproduces pre-v4 behavior and allows non-single linkage methods\ + and the primary dendrogram plot.", action='store_true', default=False) GRflags = cluster_parent.add_argument_group('GREEDY CLUSTERING OPTIONS\n' diff --git a/drep/d_analyze.py b/drep/d_analyze.py index 3049aca..544d0d9 100644 --- a/drep/d_analyze.py +++ b/drep/d_analyze.py @@ -138,6 +138,7 @@ def mash_dendrogram_from_wd(wd, plot_dir=False): Cdb = wd.get_db('Cdb', return_none=False) Pcluster = wd.get_primary_linkage() Plinkage = Pcluster['linkage'] + Plinkage_db = Pcluster.get('db') clust_args = wd.arguments['cluster'] PL_thresh = clust_args.get('P_ani', False) if PL_thresh != False: @@ -150,10 +151,19 @@ def mash_dendrogram_from_wd(wd, plot_dir=False): logging.error("Skipping plot 1 - cannot generate with multiround_primary_clustering enabled") return + if Plinkage is None or isinstance(Plinkage, str): + logging.error("Skipping plot 1 - no primary linkage matrix was computed (too many genomes, or a streaming primary algorithm was used)") + return + + # Leaf labels have to come from whatever the linkage was built on. The sparse + # skani Mdb only holds above-threshold pairs, so a genome with no relatives is + # absent from it and labels derived from Mdb would not match the linkage. + names = list(Plinkage_db.columns) if Plinkage_db is not None else None + # Make the plot logging.info("Plotting primary dendrogram") plot_MASH_dendrogram(Mdb, Cdb, Plinkage, threshold = PL_thresh,\ - plot_dir = plot_dir) + plot_dir = plot_dir, names = names) def plot_secondary_dendrograms_from_wd(wd, plot_dir, **kwargs): ''' @@ -614,7 +624,7 @@ def plot_ANIn_vs_len(Mdb,Ndb,exclude_zero_MASH=True): CLUSETER PLOTS """ -def plot_MASH_dendrogram(Mdb, Cdb, linkage, threshold=False, plot_dir=False): +def plot_MASH_dendrogram(Mdb, Cdb, linkage, threshold=False, plot_dir=False, names=None): ''' Make a dendrogram of the primary clustering @@ -624,6 +634,11 @@ def plot_MASH_dendrogram(Mdb, Cdb, linkage, threshold=False, plot_dir=False): linkage: Result of scipy.cluster.hierarchy.linkage threshold (optional): Line to plot on x-axis plot_dir (optional): Location to store plot + names (optional): Leaf labels, in the order the linkage was built from. + Required when Mdb does not contain every genome -- the sparse skani + Mdb only holds above-threshold pairs, so a genome with no relatives + never appears in it and deriving labels from Mdb would silently + mismatch the linkage. Returns: Makes and shows plot @@ -633,8 +648,9 @@ def plot_MASH_dendrogram(Mdb, Cdb, linkage, threshold=False, plot_dir=False): if Mdb['genome1'].dtype.name == 'category': logging.error("WARNING: Primary dendrogram labels may be shuffled! Load as csv to prevent this") - db = Mdb.pivot(index="genome1", columns="genome2", values="similarity") - names = list(db.columns) + if names is None: + db = Mdb.pivot(index="genome1", columns="genome2", values="similarity") + names = list(db.columns) name2cluster = Cdb.set_index('genome')['primary_cluster'].to_dict() name2color = gen_color_dictionary(names, name2cluster) @@ -1073,6 +1089,32 @@ def gen_color_list(names,name2cluster): return colors +# UC Berkeley palette. The point of coloring clusters is to tell neighbouring +# ones apart, not to identify a cluster by its color, so a handful of distinct +# colors cycled is strictly more readable than giving every cluster its own +# barely-distinguishable shade. +CLUSTER_COLORS = [ + '#003262', # Berkeley Blue + '#FDB515', # California Gold + '#3B7EA1', # Founders Rock +] + + +def _cluster_sort_key(cluster): + ''' + Order clusters naturally so that cycling colors lands adjacent clusters on + different colors. Handles primary clusters ('2') and secondary clusters + ('2_10'), sorting numerically where possible: 2_2 before 2_10, not after. + ''' + key = [] + for part in str(cluster).split('_'): + try: + key.append((0, float(part), '')) + except ValueError: + key.append((1, 0.0, part)) + return key + + def gen_color_dictionary(names, name2cluster): ''' Make the dictionary name2color @@ -1084,27 +1126,14 @@ def gen_color_dictionary(names, name2cluster): Returns: dict: name -> color ''' - #cm = _rand_cmap(len(set(name2cluster.values()))+1,type='bright') - vals = np.linspace(0,1,len(set(name2cluster.values()))+1) - np.random.shuffle(vals) - cm = plt.cm.colors.ListedColormap(plt.cm.jet(vals)) - - # 1. generate cluster to color - cluster2color = {} - clusters = set(name2cluster.values()) - NUM_COLORS = len(clusters) - for cluster in clusters: - try: - cluster2color[cluster] = cm(1.*int(cluster)/NUM_COLORS) - except: - cluster2color[cluster] = cm(1.*float(str(cluster).split('_')[1])/NUM_COLORS) - - #2. name to color - name2color = {} - for name in names: - name2color[name] = cluster2color[name2cluster[name]] + # Cycle a small palette in cluster order. This is deterministic: the previous + # implementation shuffled an unseeded colormap, so the same analysis produced + # different colors on every run. + clusters = sorted(set(name2cluster.values()), key=_cluster_sort_key) + cluster2color = {c: CLUSTER_COLORS[i % len(CLUSTER_COLORS)] + for i, c in enumerate(clusters)} - return name2color + return {name: cluster2color[name2cluster[name]] for name in names} def _comp_cluster(c): ''' diff --git a/drep/d_cluster/cluster_utils.py b/drep/d_cluster/cluster_utils.py index c2c8678..7f6b63e 100644 --- a/drep/d_cluster/cluster_utils.py +++ b/drep/d_cluster/cluster_utils.py @@ -6,7 +6,6 @@ import pandas as pd import scipy.cluster from scipy.spatial import distance as ssd -import networkx as nx import drep.d_cluster.utils @@ -88,41 +87,18 @@ def iteratre_clusters(Bdb, Cdb, id='primary_cluster'): yield d, cluster -def cluster_threshold_graph_optimized(db, linkage_cutoff=0.10, linkage_method='single'): - # Filter distances below threshold first - filtered_edges = db[db['dist'] <= linkage_cutoff] - - # Create graph directly from filtered edges - G = nx.Graph() - G.add_edges_from(filtered_edges[['genome1', 'genome2']].values) - - # Log that we're using the optimized method - logging.debug("Using low-RAM optimized clustering method with {0} edges".format(len(filtered_edges))) - - # Find connected components - clusters = {} - for cluster_id, component in enumerate(nx.connected_components(G)): - for genome in component: - clusters[genome] = cluster_id + 1 - - # Add isolated nodes (if needed) - all_genomes = set(db['genome1']).union(set(db['genome2'])) - for genome in all_genomes: - if genome not in clusters: - clusters[genome] = len(clusters) - - # Return clusters and a special value indicating we used the optimized method - return clusters, "optimized_method_used" - -def cluster_hierarchical(db, linkage_method= 'single', linkage_cutoff= 0.10, low_ram=False): +def cluster_hierarchical(db, linkage_method= 'single', linkage_cutoff= 0.10): ''' Perform hierarchical clustering on a symmetrical distiance matrix + Note this builds a dense matrix and is O(N^2) in memory. Single-linkage + primary clustering goes through drep.d_cluster.union_find instead, which is + equivalent but does not need the matrix. + Args: db: result of db.pivot usually linkage_method: passed to scipy.cluster.hierarchy.fcluster linkage_cutoff: distance to draw the clustering line (default = .1) - low_ram: whether to use the memory-efficient algorithm Returns: list: [Cdb, linkage] @@ -130,32 +106,21 @@ def cluster_hierarchical(db, linkage_method= 'single', linkage_cutoff= 0.10, low # Save names names = list(db.columns) - if low_ram: - # Convert to long format for optimized clustering - db_long = db.stack().reset_index() - db_long.columns = ['genome1', 'genome2', 'dist'] - clusters, _ = cluster_threshold_graph_optimized(db_long, linkage_cutoff, linkage_method) - - # Convert clusters to Cdb format - Cdb = pd.DataFrame({'genome': list(clusters.keys()), - 'cluster': list(clusters.values())}) - return Cdb, _ - else: - # Generate linkage dataframe - arr = np.asarray(db) - try: - arr = ssd.squareform(arr) - except: - logging.error("The database passed in is not symmetrical!") - logging.error(arr) - logging.error(names) - sys.exit() - linkage = scipy.cluster.hierarchy.linkage(arr, method= linkage_method) - - # Form clusters - fclust = scipy.cluster.hierarchy.fcluster(linkage,linkage_cutoff, \ - criterion='distance') - # Make Cdb - Cdb = drep.d_cluster.utils._gen_cdb_from_fclust(fclust,names) - - return Cdb, linkage \ No newline at end of file + # Generate linkage dataframe + arr = np.asarray(db) + try: + arr = ssd.squareform(arr) + except: + logging.error("The database passed in is not symmetrical!") + logging.error(arr) + logging.error(names) + sys.exit() + linkage = scipy.cluster.hierarchy.linkage(arr, method= linkage_method) + + # Form clusters + fclust = scipy.cluster.hierarchy.fcluster(linkage,linkage_cutoff, \ + criterion='distance') + # Make Cdb + Cdb = drep.d_cluster.utils._gen_cdb_from_fclust(fclust,names) + + return Cdb, linkage \ No newline at end of file diff --git a/drep/d_cluster/compare_utils.py b/drep/d_cluster/compare_utils.py index 7478768..7954faf 100644 --- a/drep/d_cluster/compare_utils.py +++ b/drep/d_cluster/compare_utils.py @@ -3,13 +3,17 @@ import os import sys +import numpy as np import pandas as pd +import scipy.cluster +from scipy.spatial import distance as ssd import drep import drep.d_cluster.cluster_utils import drep.d_cluster.external import drep.d_cluster.utils import drep.d_cluster.greedy_clustering +import drep.d_cluster.union_find class genomeChunk(): """ @@ -134,6 +138,138 @@ def all_vs_all_MASH(Bdb, data_folder, **kwargs): logging.info(" Final step: comparing between all groups") return run_second_round_clustering(Bdb, genome_chunks, data_folder, verbose=True, **kwargs) + +def all_vs_all_primary(Bdb, data_folder, **kwargs): + """ + Dispatch primary clustering to the requested algorithm. + + 'skani' (default) uses `skani triangle --sparse` + union-find, which never + builds the N^2 matrix on disk or in RAM and lets a skani --S_algorithm reuse + the comparisons. 'MASH' is the classic all-vs-all Mash path (pre-v4). + + Returns (Mdb, Cdb, cluster_ret), matching all_vs_all_MASH. + """ + method = kwargs.get('primary_algorithm', 'skani') + if method == 'skani': + # These only mean something on the MASH path. skani's sparse output never + # builds the N^2 table they exist to work around, so say so rather than + # silently ignoring them. + if kwargs.get('multiround_primary_clustering', False): + logging.warning( + "--multiround_primary_clustering only applies to --primary_algorithm MASH " + "and is ignored with skani. skani's sparse output never builds the full " + "N^2 table that multiround exists to avoid, and it does not suffer the " + "chunk-splitting imprecision of multiround.") + return primary_cluster_skani_sparse(Bdb, data_folder, **kwargs) + return all_vs_all_MASH(Bdb, data_folder, **kwargs) + + +def primary_cluster_skani_sparse(Bdb, data_folder, **kwargs): + """ + Primary clustering via one `skani triangle --sparse` pass + union-find. + + Only above-screen pairs are ever produced, so there is no N^2 matrix on disk + or in RAM. Always single-linkage (connected components); + --classic_primary_clustering / non-single primary_clusterAlg do not apply. + + The returned Mdb holds *every* edge from that pass, not just the ones above + P_ani. That is deliberate: secondary clustering compares genomes within a + primary cluster, which is a subset of what this pass already computed, so it + can reuse these edges instead of re-running skani per cluster. See + secondary_clustering_from_primary_edges. + """ + P_ani = kwargs.get('P_ani', 0.9) + ani_threshold = P_ani * 100.0 + + # Screen a few points below the ANI threshold so skani's k-mer pre-filter + # doesn't drop a pair whose full ANI would clear the threshold. Also stay at + # or below the secondary threshold, since secondary reuses these edges. + S_ani = kwargs.get('S_ani', 0.95) + default_screen = max(1.0, min(ani_threshold - 5.0, S_ani * 100.0 - 5.0, 99.0)) + screen = kwargs.get('primary_skani_screen', default_screen) + + # Minimum percent of a genome that must align for a pair to be reported. + # Do not lower this casually: skani's ANI ignores how much of the genome + # aligned, so without this filter genomes sharing only a small conserved + # region become edges and single linkage chains them into one huge cluster. + # See run_skani_triangle_sparse for the measurements behind the default. + min_af = kwargs.get('primary_skani_min_af', 15) + + # Secondary applies its own coverage filter at cov_thresh, so the single pass + # has to emit anything secondary might still care about. Ask skani for the + # looser of the two and apply the stricter primary filter ourselves below. + cov_thresh = float(kwargs.get('cov_thresh', 0.1)) + emit_min_af = min(min_af, cov_thresh * 100.0) + + skani_folder = os.path.join(data_folder, 'skani_sparse_files/') + genome_list = list(Bdb['location'].unique()) + + logging.info(f" Running sparse skani primary clustering on {len(genome_list):,} genomes " + f"(ANI threshold {ani_threshold:.1f}%, screen {screen:.1f}%, " + f"min-af {min_af}%, emitting min-af {emit_min_af:.1f}%)") + sparse_file = drep.d_cluster.external.run_skani_triangle_sparse( + genome_list, skani_folder, screen, min_af=emit_min_af, **kwargs) + + all_genomes = list(Bdb['genome'].unique()) + edges = drep.d_cluster.union_find.load_skani_sparse_edges(sparse_file) + Cdb, stats = drep.d_cluster.union_find.cluster_edges( + edges, P_ani, all_genomes, cov_threshold=min_af / 100.0) + + logging.info(f" Sparse skani primary clustering: {stats['edges_kept']:,} edges above " + f"threshold, {stats['primary_clusters']:,} primary clusters " + f"({stats['edges_total']:,} edges retained for secondary)") + + # Mdb keeps every edge so secondary can reuse them. similarity/dist mirror the + # MASH Mdb schema; alignment_coverage is the aligned fraction of genome1. + Mdb = edges.rename(columns={'ani': 'similarity'}).copy() + Mdb['dist'] = 1 - Mdb['similarity'] + + # The sparse path builds no dense matrix, so there is normally no scipy + # linkage to draw a primary dendrogram from. For modest genome sets the dense + # matrix is cheap, so build it from the edges purely so the dendrogram still + # works. Above the cutoff we store a marker and plotting skips it. + linkage = 'union_find_streaming' + linkage_db = None + dendro_max = kwargs.get('primary_dendrogram_max_genomes', 2000) + if len(all_genomes) <= dendro_max: + try: + linkage_db = drep.d_cluster.union_find.edges_to_dense_dist(edges, all_genomes) + arr = ssd.squareform(np.asarray(linkage_db), checks=False) + linkage = scipy.cluster.hierarchy.linkage(arr, method='single') + except Exception as e: + logging.debug(f"Skipping primary dendrogram linkage computation: {e}") + linkage = 'union_find_streaming' + linkage_db = None + + arguments = {'linkage_method': 'single', 'linkage_cutoff': 1 - P_ani, + 'comparison_algorithm': 'skani'} + cluster_ret = [linkage, linkage_db, arguments] + return Mdb, Cdb, cluster_ret + + +def secondary_clustering_from_primary_edges(Bdb, Cdb, Mdb, **kwargs): + """ + Secondary clustering that reuses primary's skani edges instead of re-running + skani once per primary cluster. + + The per-cluster comparisons dRep normally runs here recompute ANI values that + the single sparse pass already produced exactly -- on 10,000 UHGG genomes, + 94% of the pairs driving secondary clustering were already present, with + identical ANI to 6 decimal places, and the reused path reproduced the + two-stage partition exactly (1,232 clusters) in 24s instead of 15 minutes. + + Returns (Ndb, Cdb, c2ret), matching secondary_clustering. + """ + edges = Mdb.rename(columns={'similarity': 'ani'})[ + ['genome1', 'genome2', 'ani', 'alignment_coverage']] + Ndb = drep.d_cluster.union_find.build_ndb_from_edges(edges, Cdb) + + logging.info(f" Reusing {len(edges):,} primary skani edges for secondary clustering " + f"(no new comparisons); Ndb has {len(Ndb):,} rows") + + Cdb2, c2ret = drep.d_cluster.utils._cluster_Ndb(Ndb, comp_method='skani', **kwargs) + return Ndb, Cdb2, c2ret + def prepare_mash(data_folder, **kwargs): """ Make some folders and things @@ -214,9 +350,28 @@ def run_mash_on_genome_chunks(genome_chunks, mash_exe, sketch_folder, MASH_folde return genome_chunks +def _subsample_mdb(mdb, max_rows): + """ + Cap a per-chunk Mdb to at most max_rows rows so multiround primary clustering + doesn't accumulate an O(N^2) table across all chunks (the source of the + 43k-genome MemoryError at pd.concat). The full pairwise table is only kept for + storage/inspection; clustering itself does not use the concatenated Mdb. + """ + if max_rows is None or len(mdb) <= max_rows: + return mdb + return mdb.sample(n=max_rows, random_state=0) + + def run_second_round_clustering(Bdb, genome_chunks, data_folder, **kwargs): verbose = kwargs.get('verbose', False) + # Bound the total number of pairwise rows retained for the stored Mdb across + # all chunks. Set to 0/None to disable (restores pre-v4 unbounded behavior). + max_mdb_rows = kwargs.get('max_stored_mdb_rows', 5_000_000) + per_chunk_cap = None + if max_mdb_rows: + per_chunk_cap = max(1, int(max_mdb_rows // (len(genome_chunks) + 1))) + kwargs_copy = kwargs.copy() kwargs_copy['multiround_primary_clustering'] = False kwargs_copy['v2'] = '_v2' @@ -233,7 +388,10 @@ def run_second_round_clustering(Bdb, genome_chunks, data_folder, **kwargs): mdb = gc.Mdb mdb['genome_chunk'] = gc.name - mdbs.append(mdb) + # Subsample before retaining so we never hold all N^2 rows at once + mdbs.append(_subsample_mdb(mdb, per_chunk_cap)) + # Free the chunk's full table now that its clusters are computed + gc.Mdb = None Cdb = pd.concat(dbs) @@ -255,10 +413,11 @@ def run_second_round_clustering(Bdb, genome_chunks, data_folder, **kwargs): mdb = genome_chunks[0].Mdb mdb['genome_chunk'] = 'v2' - mdbs.append(mdb) - Mdb = pd.concat(mdbs).reset_index(drop=True) - + # Cluster on the full second-round table, but only store a bounded subsample Cdb2, cluster_ret = cluster_mash_database(mdb, **kwargs) + + mdbs.append(_subsample_mdb(mdb, per_chunk_cap)) + Mdb = pd.concat(mdbs).reset_index(drop=True) Cdb2['primary_representitive'] = True # Step 5) Merge the new Cdb back in with the old @@ -277,25 +436,64 @@ def cluster_mash_database(db, **kwargs): db: Mdb (all_vs_all Mash results) Keyword arguments: - clusterAlg: how to cluster database (default = single) + primary_clusterAlg: how to cluster the primary database (default = single). + 'single' uses the fast streaming union-find algorithm; any other + method uses the classic dense scipy path. + clusterAlg: legacy fallback for primary_clusterAlg (default = single) P_ani: threshold to cluster at (default = 0.9) - low_ram_primary_clustering: whether to use memory-efficient algorithm + classic_primary_clustering: force the dense scipy path Returns: list: [Cdb, [linkage, linkage_db, arguments]] ''' logging.debug('Clustering MASH database') - # Load key words - P_Lmethod = kwargs.get('clusterAlg','single') + # Load key words. Primary clustering has its own linkage method + # (primary_clusterAlg), independent of the secondary clusterAlg. Fall back to + # clusterAlg for older callers that only pass that. + P_Lmethod = kwargs.get('primary_clusterAlg') or kwargs.get('clusterAlg', 'single') P_Lcutoff = 1 - kwargs.get('P_ani',.9) - low_ram = kwargs.get('low_ram_primary_clustering', False) + classic = kwargs.get('classic_primary_clustering', False) - # Do the actual clustering db['dist'] = 1 - db['similarity'] + + # Single-linkage clustering at a fixed cutoff is identical to connected + # components. Compute it directly on the long-format table with union-find and + # skip the O(N^2) dense pivot entirely (issue #259 / the large-N RAM crash). + # This is the default; --classic_primary_clustering forces the dense path. + use_union_find = (not classic) and (P_Lmethod == 'single') + if use_union_find: + Cdb = drep.d_cluster.union_find.cluster_long_df(db, P_Lcutoff) + + arguments = {'linkage_method': 'single', 'linkage_cutoff': P_Lcutoff, + 'comparison_algorithm': 'MASH'} + + # The streaming path builds no dense matrix, so by default there is no + # scipy linkage to plot a primary dendrogram from. For modest genome sets + # the dense pivot is cheap, so compute the single-linkage matrix purely so + # the dendrogram can still be drawn. Above the cutoff (or if it fails) we + # store a marker and downstream plotting skips the dendrogram gracefully. + linkage = 'union_find_streaming' + linkage_db = None + dendro_max = kwargs.get('primary_dendrogram_max_genomes', 2000) + n_genomes = Cdb['genome'].nunique() + if n_genomes <= dendro_max and 'genome_chunk' not in db.columns: + try: + linkage_db = db.pivot(index="genome1", columns="genome2", values="dist") + arr = ssd.squareform(np.asarray(linkage_db)) + linkage = scipy.cluster.hierarchy.linkage(arr, method='single') + except Exception as e: + logging.debug(f"Skipping primary dendrogram linkage computation: {e}") + linkage = 'union_find_streaming' + linkage_db = None + + cluster_ret = [linkage, linkage_db, arguments] + return Cdb, cluster_ret + + # Classic dense path (non-single linkage, or --classic_primary_clustering). linkage_db = db.pivot(index="genome1", columns="genome2", values="dist") Cdb, linkage = drep.d_cluster.cluster_utils.cluster_hierarchical(linkage_db, linkage_method= P_Lmethod, \ - linkage_cutoff= P_Lcutoff, low_ram=low_ram) + linkage_cutoff= P_Lcutoff) Cdb = Cdb.rename(columns={'cluster':'primary_cluster'}) Cdb['primary_cluster'] = Cdb['primary_cluster'].astype(int) diff --git a/drep/d_cluster/controller.py b/drep/d_cluster/controller.py index d4a1463..a144e9d 100644 --- a/drep/d_cluster/controller.py +++ b/drep/d_cluster/controller.py @@ -44,11 +44,12 @@ def parse_cluster_arguments(self): """ Load the genomes and store Bdb in the wd """ - # Make sure you have the required program installed - loc = shutil.which('mash') - if loc is None: + # Make sure the program this run actually needs is installed. Only the + # MASH primary path needs mash; the default (skani) does not. + primary_exe = 'mash' if self.kwargs.get('primary_algorithm', 'skani') == 'MASH' else 'skani' + if shutil.which(primary_exe) is None: logging.error('Cannot locate the program {0}- make sure its in the system path' \ - .format('mash')) + .format(primary_exe)) # If genomes are provided, load them if self.kwargs.get('genomes', None) is not None: @@ -102,7 +103,7 @@ def run_primary_clustering(self): else: logging.info("Running pair-wise MASH clustering") - Mdb, Cdb, cluster_ret = drep.d_cluster.compare_utils.all_vs_all_MASH(self.Bdb, self.wd.get_dir('MASH'), **self.kwargs) + Mdb, Cdb, cluster_ret = drep.d_cluster.compare_utils.all_vs_all_primary(self.Bdb, self.wd.get_dir('MASH'), **self.kwargs) if self.debug: logging.debug("Debug mode on - saving Mdb ASAP") @@ -148,6 +149,16 @@ def run_secondary_clustering(self): logging.info('3. Secondary clustering cache loaded') + # Reuse primary's skani edges instead of re-running skani per cluster + elif self.can_reuse_primary_edges(algorithm): + logging.info("Reusing primary skani comparisons for secondary clustering") + Ndb, Cdb, c2ret = drep.d_cluster.compare_utils.secondary_clustering_from_primary_edges( + self.Bdb, self.MCdb, self.Mdb, **self.kwargs) + if self.debug: + self.wd.store_db(Ndb, 'Ndb') + self.wd.store_db(Cdb, 'Cdb') + self.wd.store_special('secondary_linkages', c2ret) + # Run comparisons, make Ndb else: drep.d_cluster.utils._print_time_estimate(self.Bdb, self.MCdb, algorithm, p) @@ -171,6 +182,27 @@ def run_secondary_clustering(self): self.Cdb = Cdb self.Ndb = Ndb + def can_reuse_primary_edges(self, algorithm): + """ + Whether secondary clustering can be derived from primary's edges rather + than re-running comparisons. + + This only holds when primary was skani (so Mdb contains real ANI plus + alignment fractions for every pair above the screen) and secondary wants + skani too. Any other secondary algorithm measures something different and + has to run for itself; greedy has its own code path. + """ + if self.kwargs.get('reuse_primary_comparisons', True) is False: + return False + if self.kwargs.get('primary_algorithm', 'skani') != 'skani': + return False + if algorithm != 'skani': + return False + if self.kwargs.get('greedy_secondary_clustering', False): + return False + # Mdb must be the skani edge table, not a Mash table or a blank + return (self.Mdb is not None) and ('alignment_coverage' in self.Mdb.columns) + def store_output(self): logging.debug("Main program run complete- saving output") self.wd.store_db(self.Cdb, 'Cdb') diff --git a/drep/d_cluster/external.py b/drep/d_cluster/external.py index ff24b0c..036e8af 100644 --- a/drep/d_cluster/external.py +++ b/drep/d_cluster/external.py @@ -240,8 +240,15 @@ def load_skani(file): db = db[db['reference'] != db['querry']] adb = pd.concat([adb, db], ignore_index=True).reset_index(drop=True) - # Load the af triangle + # Load the af triangle. skani reports aligned fractions as percentages + # (0-100); every other dRep algorithm reports alignment_coverage on a 0-1 + # scale, and that is what cov_thresh is compared against in + # make_linkage_Ndb. Without this conversion the coverage filter is inert for + # skani (e.g. a pair aligning over only 1% of the genome has + # alignment_coverage=1.04, which sails past a cov_thresh of 0.5), which can + # merge distantly related genomes that share a small conserved region. tdb = load_matrix_to_dataframe(file + '.af').rename(columns={'ani':'alignment_coverage'}) + tdb['alignment_coverage'] = tdb['alignment_coverage'] / 100 # Merge assert len(adb) == len(tdb) @@ -271,6 +278,74 @@ def _fix_fastani(odb): return fdb +def run_skani_triangle_sparse(genome_list, outdir, screen, min_af=15, **kwargs): + """ + Run `skani triangle --sparse` and return the path to the sparse output file. + + The sparse output is an edge list of only the above-screening-threshold pairs, + so it never materializes the N^2 matrix on disk or in memory. It is meant to + be streamed (see union_find.cluster_skani_sparse_files), not loaded whole. + + Args: + genome_list: list of genome file locations. + outdir: directory to write the sparse output and temp files. + screen: skani -s screening threshold (percent identity). Pairs below this + are discarded during sketching and never appear in the output. Should + be <= the primary ANI threshold so no real edges are missed. + min_af: skani --min-af, the minimum percent of a genome that must align + for the pair to be reported. See the note below -- do not set this to + 0 for primary clustering. + + Keyword Args: + processors: threads for skani (default 6). + skani_extra: extra args passed through to skani triangle. + wd, debug: for command logging. + + Returns: + Path to the sparse skani output file. + """ + p = kwargs.get('processors', 6) + code = drep.d_cluster.utils._randomString(stringLength=10) + extra_cmd = kwargs.get('skani_extra', "") + + if not os.path.exists(outdir): + os.makedirs(outdir) + tmp_dir = os.path.join(outdir, 'tmp/') + if not os.path.exists(tmp_dir): + os.makedirs(tmp_dir) + + glist = os.path.join(tmp_dir, 'genomeList_{0}'.format(code)) + glist = _make_glist(genome_list, glist) + + exe_loc = drep.get_exe('skani') + out_file = os.path.join(outdir, 'skani_sparse_{0}.tsv'.format(code)) + # min_af matters far more than it looks, because MASH similarity and skani ANI + # measure different things. MASH compares k-mers across the whole genome, so + # two genomes sharing only a small conserved region score as distant. skani's + # ANI is the identity *within aligned regions only*, so that same pair reports + # a high ANI and (with min_af 0) becomes a primary-clustering edge. Under + # single linkage those few spurious bridges chain everything together: on 10k + # UHGG genomes, --min-af 0 collapsed 59% of the dataset into one primary + # cluster (largest 5857) while skani's 15% default reproduced the MASH + # partition almost exactly (984 clusters vs MASH's 989, largest 626 vs 626). + # The aligned-fraction filter is what makes skani's ANI comparable to MASH's + # whole-genome similarity -- it is not an obstacle to work around. + cmd = [exe_loc, "triangle", "--sparse", "-t", str(p), '-o', out_file, + '-l', glist, '-s', str(screen), '--min-af', str(min_af)] + if extra_cmd != "": + cmd += extra_cmd.split(' ') + + logging.debug(' '.join(cmd) + ' ' + code) + + if ('wd' in kwargs) and (kwargs.get('debug', False)): + logdir = kwargs.get('wd').get_dir('cmd_logs') + else: + logdir = False + drep.thread_cmds([cmd], shell=False, logdir=logdir, t=1) + + return out_file + + def _make_glist(genomes, floc): o = open(floc, 'w') for g in genomes: diff --git a/drep/d_cluster/greedy_clustering.py b/drep/d_cluster/greedy_clustering.py index 4ee1564..84399ef 100644 --- a/drep/d_cluster/greedy_clustering.py +++ b/drep/d_cluster/greedy_clustering.py @@ -87,6 +87,7 @@ def compare_genomes_greedy(bdb, algorithm, data_folder, **kwargs): genome2cluster[row['genome']] = new_cluster with open(genome_rep_file, "a") as myfile: myfile.write(row['location'] + '\n') + add_genome_as_rep(row['location'], algorithm, **kwargs) if len(ndbs) > 0: Ndb = pd.concat(ndbs) @@ -106,18 +107,36 @@ def compare_genomes_greedy(bdb, algorithm, data_folder, **kwargs): def genome_vs_reps(new_genome, genome_reps, genome_rep_file, algorithm, data_folder, **kwargs): if algorithm == 'fastANI': - # Return Ndb + # NOTE: this spawns a subprocess that re-sketches every representative on + # every call, so sketching is O(N*R). Greedy exists to avoid O(n^2) + # comparisons within a primary cluster; --primary_algorithm skani avoids + # that quadratic in the first place by only ever producing + # above-threshold pairs, and is usually the better answer at scale. return drep.d_cluster.external.fastani_one_vs_many(new_genome, genome_reps, genome_rep_file, data_folder, **kwargs) else: - logging.error("{0} algorithm is not yet supported for greedy clustering; sorry!") + logging.error("{0} algorithm is not yet supported for greedy clustering; sorry!".format(algorithm)) assert False +def add_genome_as_rep(location, algorithm, **kwargs): + """ + Register a genome as a new cluster representative. + + Subprocess-based algorithms read the representative list from a file, which + compare_genomes_greedy has already written, so there is nothing to do here. + Kept as a hook for backends that need to index representatives as they appear. + """ + return + + def prepare_for_greedy(algorithm, data_folder, **kwargs): + # Every algorithm writes the running list of representatives here, so the + # folder has to exist regardless of which one is in use. + if not os.path.exists(data_folder): + os.makedirs(data_folder) + if algorithm == 'fastANI': # Make folders - if not os.path.exists(data_folder): - os.makedirs(data_folder) tmp_dir = os.path.join(data_folder, 'tmp/') if not os.path.exists(tmp_dir): os.makedirs(tmp_dir) diff --git a/drep/d_cluster/union_find.py b/drep/d_cluster/union_find.py new file mode 100644 index 0000000..2fdb0e0 --- /dev/null +++ b/drep/d_cluster/union_find.py @@ -0,0 +1,543 @@ +""" +Streaming, low-memory primary clustering via union-find (disjoint set). + +Primary clustering in dRep is single-linkage hierarchical clustering at a fixed +distance cutoff. That is mathematically identical to finding the connected +components of the graph whose nodes are genomes and whose edges are the pairs +with ``distance <= cutoff``. + +The classic dRep path materializes every pairwise MASH comparison into a long +``Mdb`` DataFrame (N^2 rows) and then pivots it into a dense N x N matrix before +handing it to scipy. For tens of thousands of genomes this is tens of GiB of RAM +and is the source of the crashes in issue #259 and the "Big dRep issue". + +This module never builds the dense matrix and never needs to hold all N^2 pairs +in memory. It streams the MASH ``dist`` output, discards the ~99.9% of pairs that +are above the cutoff the instant they are read, and unions the survivors. Memory +is O(genomes + kept_edges) instead of O(genomes^2). +""" + +import logging + +import numpy as np +import pandas as pd + +import drep.d_cluster.utils + + +class UnionFind: + """ + Disjoint-set / union-find with path compression and union by rank. + + ``union`` and ``find`` are effectively O(alpha(N)) ~ O(1), so clustering the + surviving edges is linear in the number of edges. + """ + + def __init__(self): + self.parent = {} + self.rank = {} + + def add(self, x): + if x not in self.parent: + self.parent[x] = x + self.rank[x] = 0 + + def find(self, x): + # Find root + root = x + while self.parent[root] != root: + root = self.parent[root] + # Path compression (iterative, no recursion depth limit) + while self.parent[x] != root: + self.parent[x], x = root, self.parent[x] + return root + + def union(self, a, b): + ra, rb = self.find(a), self.find(b) + if ra == rb: + return + if self.rank[ra] < self.rank[rb]: + ra, rb = rb, ra + self.parent[rb] = ra + if self.rank[ra] == self.rank[rb]: + self.rank[ra] += 1 + + def components(self): + """ + Return {root: [members...]} for every set. + """ + comps = {} + for node in self.parent: + comps.setdefault(self.find(node), []).append(node) + return comps + + +def _components_to_cdb(uf): + """ + Turn a populated UnionFind into a Cdb (columns: genome, primary_cluster). + + Clusters are numbered deterministically: largest first, ties broken by the + alphabetically-smallest member. This makes runs reproducible regardless of + the order edges happened to stream in. + """ + comps = uf.components() + + ordered = sorted( + comps.values(), + key=lambda members: (-len(members), min(members)), + ) + + genomes = [] + clusters = [] + for cluster_id, members in enumerate(ordered, start=1): + for genome in sorted(members): + genomes.append(genome) + clusters.append(cluster_id) + + Cdb = pd.DataFrame({'genome': genomes, 'primary_cluster': clusters}) + Cdb['primary_cluster'] = Cdb['primary_cluster'].astype(int) + return Cdb + + +def cluster_long_df(db, cutoff, all_genomes=None): + """ + Cluster an in-memory long-format MASH table with union-find (no pivot). + + This is the drop-in, single-linkage replacement for the pivot -> squareform + -> scipy path in ``cluster_mash_database``. + + Args: + db: DataFrame with columns 'genome1', 'genome2', and either 'dist' or + 'similarity'. + cutoff: distance cutoff (1 - P_ani). Pairs with dist <= cutoff are edges. + all_genomes: optional iterable of every genome name, so singletons that + never appear in an above-cutoff edge still get their own cluster. If + not given, it is inferred from the genome1/genome2 columns. + + Returns: + Cdb: DataFrame with columns 'genome', 'primary_cluster'. + """ + if 'dist' in db.columns: + dist = db['dist'].values + else: + dist = 1 - db['similarity'].values + + uf = UnionFind() + + # Seed every genome so singletons are represented + if all_genomes is None: + all_genomes = pd.unique( + pd.concat([db['genome1'], db['genome2']], ignore_index=True) + ) + for g in all_genomes: + uf.add(g) + + mask = dist <= cutoff + g1 = db['genome1'].values[mask] + g2 = db['genome2'].values[mask] + for a, b in zip(g1, g2): + uf.add(a) + uf.add(b) + uf.union(a, b) + + return _components_to_cdb(uf) + + +def cluster_mash_files(dist_files, cutoff, all_genomes=None, chunksize=5_000_000, + name_from_fasta=True, progress=False): + """ + Stream one or more MASH ``dist`` output files and cluster with union-find. + + Never builds the dense matrix and never holds all N^2 pairs in memory. Only + genome names, the union-find bookkeeping, and one ``chunksize`` block of rows + are resident at a time. + + Args: + dist_files: path (str) or list of paths to MASH dist tsv output + (columns: genome1, genome2, dist, p, kmers). + cutoff: distance cutoff (1 - P_ani). + all_genomes: optional iterable of every genome name to seed singletons. + chunksize: rows per streamed block. + name_from_fasta: if True, map file paths in the table to genome names via + drep's basename logic (matches parse_mash_table behavior). + progress: if True, show a tqdm progress bar over streamed rows. + + Returns: + (Cdb, stats) where stats is a dict of counters for benchmarking/logging. + """ + if isinstance(dist_files, (str, bytes)): + dist_files = [dist_files] + + uf = UnionFind() + if all_genomes is not None: + for g in all_genomes: + uf.add(g) + + if progress: + try: + from tqdm import tqdm + except ImportError: + logging.warning("tqdm not installed; primary-clustering progress bar disabled") + progress = False + + total_rows = 0 + kept_edges = 0 + + name_cache = {} + + def to_name(x): + n = name_cache.get(x) + if n is None: + n = drep.d_cluster.utils._get_genome_name_from_fasta(x) + name_cache[x] = n + return n + + bar = tqdm(desc=" Primary clustering (streaming)", unit=" pairs") if progress else None + + for dist_file in dist_files: + reader = pd.read_csv( + dist_file, + names=['genome1', 'genome2', 'dist', 'p', 'kmers'], + usecols=['genome1', 'genome2', 'dist'], + dtype={'genome1': str, 'genome2': str, 'dist': np.float32}, + sep='\t', + chunksize=chunksize, + ) + for chunk in reader: + n = len(chunk) + total_rows += n + if bar is not None: + bar.update(n) + + hits = chunk[chunk['dist'] <= cutoff] + if len(hits) == 0: + continue + + g1 = hits['genome1'].values + g2 = hits['genome2'].values + for a, b in zip(g1, g2): + if name_from_fasta: + a = to_name(a) + b = to_name(b) + uf.add(a) + uf.add(b) + uf.union(a, b) + kept_edges += 1 + + if bar is not None: + bar.close() + + Cdb = _components_to_cdb(uf) + + stats = { + 'total_pairs_streamed': total_rows, + 'edges_kept': kept_edges, + 'genomes': len(uf.parent), + 'primary_clusters': Cdb['primary_cluster'].nunique() if len(Cdb) else 0, + } + return Cdb, stats + + +def edges_to_dense_dist(edges, genomes): + """ + Build a dense, symmetric distance matrix from a sparse edge table. + + Only for modest genome counts -- this is the O(N^2) representation the rest of + this module exists to avoid. It is used solely so the primary dendrogram can + still be plotted for small runs; clustering itself never needs it. + + Pairs absent from the edge list have no detectable similarity, so they get + distance 1. The diagonal is 0. + + Returns: + DataFrame indexed and columned by genome, values = 1 - ani. + """ + names = sorted(genomes) + idx = {g: i for i, g in enumerate(names)} + n = len(names) + + arr = np.ones((n, n), dtype=np.float32) + np.fill_diagonal(arr, 0.0) + + for a, b, ani in zip(edges['genome1'].values, edges['genome2'].values, + edges['ani'].values): + i, j = idx.get(a), idx.get(b) + if i is None or j is None: + continue + d = 1.0 - ani + # edges are symmetric already, but write both to be safe against + # one-directional input + arr[i, j] = d + arr[j, i] = d + + return pd.DataFrame(arr, index=names, columns=names) + + +def build_ndb_from_edges(edges, Cdb): + """ + Build a secondary-clustering Ndb out of primary's edge table, without running + any new comparisons. + + Primary clustering (skani, sparse) already computed the exact ANI for every + pair above skani's screening threshold. Secondary clustering only ever + compares genomes *within* a primary cluster, and those pairs are a subset of + what primary already has -- so re-running skani per primary cluster + recomputes numbers that are already known, bit for bit. + + dRep's hierarchical secondary clustering needs a complete matrix per primary + cluster, so pairs absent from the sparse edge list (i.e. below skani's + screen, meaning no meaningful similarity) are filled in as ani=0 and + coverage=0, and self-comparisons as 1. + + Args: + edges: DataFrame from load_skani_sparse_edges. + Cdb: primary clustering result with ['genome', 'primary_cluster']. + + Returns: + Ndb: ['reference', 'querry', 'ani', 'alignment_coverage', 'primary_cluster'] + """ + g2p = Cdb.set_index('genome')['primary_cluster'].to_dict() + + e = edges.copy() + e['pc'] = e['genome1'].map(g2p) + # secondary only ever compares within a primary cluster + e = e[e['pc'].notna() & (e['pc'] == e['genome2'].map(g2p))] + + have = set(zip(e['genome1'].values, e['genome2'].values)) + + fill_r, fill_q, fill_pc, fill_ani, fill_cov = [], [], [], [], [] + for pc, sub in Cdb.groupby('primary_cluster'): + gs = list(sub['genome']) + for x in gs: + for y in gs: + if x == y: + fill_r.append(x); fill_q.append(y); fill_pc.append(pc) + fill_ani.append(1.0); fill_cov.append(1.0) + elif (x, y) not in have: + fill_r.append(x); fill_q.append(y); fill_pc.append(pc) + fill_ani.append(0.0); fill_cov.append(0.0) + + kept = e.rename(columns={'genome1': 'reference', 'genome2': 'querry', + 'pc': 'primary_cluster'})[ + ['reference', 'querry', 'ani', 'alignment_coverage', 'primary_cluster']] + filled = pd.DataFrame({ + 'reference': fill_r, 'querry': fill_q, 'ani': fill_ani, + 'alignment_coverage': fill_cov, 'primary_cluster': fill_pc, + }) + + Ndb = pd.concat([kept, filled], ignore_index=True) + Ndb['primary_cluster'] = Ndb['primary_cluster'].astype(int) + return Ndb + + +SKANI_SPARSE_COLUMNS = ['Ref_file', 'Query_file', 'ANI', + 'Align_fraction_ref', 'Align_fraction_query'] + + +def load_skani_sparse_edges(sparse_files): + """ + Load `skani triangle --sparse` output into a symmetric edge table. + + skani's sparse output is already only the pairs above its screening + threshold, so unlike Mash's N^2 output it is small enough to hold in memory + (~390k rows for 10,000 genomes) and can be reused rather than recomputed. + + Each input pair is emitted in both directions, because dRep treats + alignment_coverage as the aligned fraction of the genome named in the first + column, and the two directions have different coverages. + + Returns: + DataFrame with ['genome1', 'genome2', 'ani', 'alignment_coverage'], + where ani and alignment_coverage are 0-1 fractions. + """ + if isinstance(sparse_files, (str, bytes)): + sparse_files = [sparse_files] + + frames = [] + for f in sparse_files: + d = pd.read_csv(f, sep='\t', usecols=SKANI_SPARSE_COLUMNS, + dtype={'Ref_file': str, 'Query_file': str, + 'ANI': np.float32, + 'Align_fraction_ref': np.float32, + 'Align_fraction_query': np.float32}) + if len(d) == 0: + continue + a = np.array([drep.d_cluster.utils._get_genome_name_from_fasta(x) + for x in d['Ref_file'].values]) + b = np.array([drep.d_cluster.utils._get_genome_name_from_fasta(x) + for x in d['Query_file'].values]) + ani = (d['ANI'].values / 100).astype(np.float32) + af_a = (d['Align_fraction_ref'].values / 100).astype(np.float32) + af_b = (d['Align_fraction_query'].values / 100).astype(np.float32) + + frames.append(pd.DataFrame({ + 'genome1': np.concatenate([a, b]), + 'genome2': np.concatenate([b, a]), + 'ani': np.concatenate([ani, ani]), + # coverage is always the aligned fraction of the genome1 genome + 'alignment_coverage': np.concatenate([af_a, af_b]), + })) + + if not frames: + return pd.DataFrame(columns=['genome1', 'genome2', 'ani', 'alignment_coverage']) + + edges = pd.concat(frames, ignore_index=True) + # drop self comparisons; they are added back explicitly where needed + return edges[edges['genome1'] != edges['genome2']].reset_index(drop=True) + + +def cluster_edges(edges, ani_threshold, all_genomes, cov_threshold=0.0): + """ + Union-find clustering of an in-memory edge table (see load_skani_sparse_edges). + + Args: + edges: DataFrame with ['genome1', 'genome2', 'ani', 'alignment_coverage']. + ani_threshold: minimum ANI (0-1) for a pair to be an edge. + all_genomes: every genome name, so singletons get their own cluster. + cov_threshold: minimum alignment_coverage (0-1) for a pair to be an edge. + See run_skani_triangle_sparse for why this matters -- without it, + genomes sharing a small conserved region chain together. + + Returns: + (Cdb, stats) + """ + uf = UnionFind() + for g in all_genomes: + uf.add(g) + + keep = edges['ani'].values >= ani_threshold + if cov_threshold > 0: + keep &= edges['alignment_coverage'].values >= cov_threshold + + g1 = edges['genome1'].values[keep] + g2 = edges['genome2'].values[keep] + for a, b in zip(g1, g2): + uf.add(a) + uf.add(b) + uf.union(a, b) + + Cdb = _components_to_cdb(uf) + stats = { + 'edges_total': len(edges), + 'edges_kept': int(keep.sum()), + 'genomes': len(uf.parent), + 'primary_clusters': Cdb['primary_cluster'].nunique() if len(Cdb) else 0, + } + return Cdb, stats + + +def cluster_skani_sparse_files(sparse_files, ani_threshold, all_genomes, + cov_threshold=0.0, chunksize=2_000_000, progress=False): + """ + Stream `skani triangle --sparse` output and cluster with union-find. + + Unlike MASH ``dist``, the sparse skani output already contains *only* the + above-screening-threshold pairs (an edge list), so there is no N^2 to stream + at all -- just the surviving edges. This is the low-RAM, low-disk primary + clustering path for very large genome sets. + + Expected columns (skani >= 0.2, with a header row): + Ref_file, Query_file, ANI, Align_fraction_ref, Align_fraction_query, + Ref_name, Query_name + + Args: + sparse_files: path or list of paths to sparse skani output. + ani_threshold: percent ANI (e.g. 90.0 for P_ani=0.9). Pairs at or above + this are treated as edges. + all_genomes: iterable of every genome name, so singletons that skani + screened out still get their own primary cluster. + cov_threshold: minimum aligned fraction (0-1) for a pair to count as an + edge. skani reports two percentages; the larger is used (matching + dRep's 'larger' coverage convention). 0 disables the filter. + chunksize: rows per streamed block. + progress: show a tqdm bar over streamed edges. + + Returns: + (Cdb, Mdb, stats): + Cdb: ['genome', 'primary_cluster'] + Mdb: reduced long-format table of the surviving edges only + (['genome1', 'genome2', 'similarity', 'dist']), for storage/plots. + stats: dict of counters. + """ + if isinstance(sparse_files, (str, bytes)): + sparse_files = [sparse_files] + + uf = UnionFind() + for g in all_genomes: + uf.add(g) + + if progress: + try: + from tqdm import tqdm + except ImportError: + logging.warning("tqdm not installed; primary-clustering progress bar disabled") + progress = False + bar = tqdm(desc=" Primary clustering (sparse skani)", unit=" edges") if progress else None + + edges_seen = 0 + kept_edges = 0 + cov_pct = cov_threshold * 100.0 + mdb_g1, mdb_g2, mdb_sim = [], [], [] + + name_cache = {} + + def to_name(x): + n = name_cache.get(x) + if n is None: + n = drep.d_cluster.utils._get_genome_name_from_fasta(x) + name_cache[x] = n + return n + + for sparse_file in sparse_files: + reader = pd.read_csv( + sparse_file, + sep='\t', + usecols=['Ref_file', 'Query_file', 'ANI', + 'Align_fraction_ref', 'Align_fraction_query'], + dtype={'Ref_file': str, 'Query_file': str, 'ANI': np.float32, + 'Align_fraction_ref': np.float32, 'Align_fraction_query': np.float32}, + chunksize=chunksize, + ) + for chunk in reader: + edges_seen += len(chunk) + if bar is not None: + bar.update(len(chunk)) + + hits = chunk[chunk['ANI'] >= ani_threshold] + if cov_pct > 0: + larger_af = np.maximum(hits['Align_fraction_ref'].values, + hits['Align_fraction_query'].values) + hits = hits[larger_af >= cov_pct] + if len(hits) == 0: + continue + + for r, q, ani in zip(hits['Ref_file'].values, + hits['Query_file'].values, + hits['ANI'].values): + a, b = to_name(r), to_name(q) + if a == b: + continue + uf.add(a) + uf.add(b) + uf.union(a, b) + mdb_g1.append(a) + mdb_g2.append(b) + mdb_sim.append(ani / 100.0) + kept_edges += 1 + + if bar is not None: + bar.close() + + Cdb = _components_to_cdb(uf) + + Mdb = pd.DataFrame({'genome1': mdb_g1, 'genome2': mdb_g2, + 'similarity': np.array(mdb_sim, dtype=np.float32)}) + Mdb['dist'] = 1 - Mdb['similarity'] + + stats = { + 'edges_seen': edges_seen, + 'edges_kept': kept_edges, + 'genomes': len(uf.parent), + 'primary_clusters': Cdb['primary_cluster'].nunique() if len(Cdb) else 0, + } + return Cdb, Mdb, stats diff --git a/drep/d_cluster/utils.py b/drep/d_cluster/utils.py index 91d451f..590008a 100644 --- a/drep/d_cluster/utils.py +++ b/drep/d_cluster/utils.py @@ -93,19 +93,20 @@ def estimate_time(comps, alg): Return: float: time to perfom comparison (in minutes) ''' - if alg == 'ANIn': - time = comps * .33 - elif alg == 'gANI': - time = comps * .1 - elif alg == 'goANI': - time = comps * .1 - elif alg == 'ANImf': - time = comps * .5 - elif alg == 'fastANI': - time = comps * 0.00667 - elif alg == 'skani': - time = comps * 0.00667 - return time + # Minutes per comparison, very roughly. This only drives a log message, so an + # unknown algorithm must never take the run down with it -- fall back to the + # fastest estimate rather than raising. + per_comparison = { + 'ANIn': .33, + 'gANI': .1, + 'goANI': .1, + 'ANImf': .5, + 'fastANI': 0.00667, + 'skani': 0.00667, + } + if alg not in per_comparison: + logging.debug(f"No time estimate available for {alg}; assuming a fast algorithm") + return comps * per_comparison.get(alg, 0.00667) diff --git a/helper_scripts/ScaffoldLevel_dRep.py b/helper_scripts/ScaffoldLevel_dRep.py index 5203927..74d7dbf 100755 --- a/helper_scripts/ScaffoldLevel_dRep.py +++ b/helper_scripts/ScaffoldLevel_dRep.py @@ -14,6 +14,8 @@ import shutil import distutils import argparse +import functools +import subprocess import pandas as pd from shutil import copyfile @@ -154,13 +156,34 @@ def gen_prefix(self): def __str__(self): ''' Show the command parameters ''' +@functools.lru_cache(maxsize=None) +def nucmer_supports_threads(exe): + ''' + Whether this nucmer accepts -t/--threads. + + MUMmer 4 added it; MUMmer 3 (still what `conda install mummer` gives, as + version 3.23) does not, and errors out with "Unknown option: t" rather than + ignoring it. Passing it unconditionally makes this script fail outright on + MUMmer 3, so detect support instead of assuming. + ''' + try: + r = subprocess.run([exe, '--help'], capture_output=True, text=True, timeout=60) + return '--threads' in (r.stdout + r.stderr) + except Exception: + return False + + def gen_mummer_cmd(**kwargs): ''' from a dictionary of arguments, return the ANIm command as an array of strings ''' cmd = [kwargs['exe'],'--' + kwargs['method'],'-p',kwargs['prefix'], '-c', \ - kwargs['c'], '-g', kwargs['maxgap'], '-t', str(kwargs['p'])] + kwargs['c'], '-g', kwargs['maxgap']] + + # MUMmer 3's nucmer is single-threaded and rejects -t outright + if nucmer_supports_threads(kwargs['exe']): + cmd += ['-t', str(kwargs['p'])] if kwargs['noextend'] == 'True': cmd.append('--noextend') diff --git a/setup.py b/setup.py index 81f4119..5faac86 100644 --- a/setup.py +++ b/setup.py @@ -24,7 +24,6 @@ def version(): 'biopython', 'scikit-learn', 'tqdm', - 'networkx', 'setuptools', 'pytest' ], diff --git a/tests/tests/test_bonus.py b/tests/tests/test_bonus.py index 62d6682..8b79eed 100644 --- a/tests/tests/test_bonus.py +++ b/tests/tests/test_bonus.py @@ -230,3 +230,46 @@ def test_parse_stb_3(self): db = pd.read_csv(out_loc, sep='\t', names=['scaffold', 'bin']) assert len(db) == 124 assert len(db['bin'].unique()) == 5 + + +def _load_scaffold_level_module(): + """ScaffoldLevel_dRep.py is a script, not a package module; load it by path.""" + import importlib.util + loc = os.path.join(os.path.dirname(os.path.abspath(__file__)), + '../../helper_scripts/ScaffoldLevel_dRep.py') + spec = importlib.util.spec_from_file_location('scaffold_level_drep', loc) + mod = importlib.util.module_from_spec(spec) + spec.loader.exec_module(mod) + return mod + + +@pytest.mark.requires_mummer +def test_scaffold_level_nucmer_threads_flag(): + """ + Regression test: nucmer only grew -t/--threads in MUMmer 4. MUMmer 3 (which + is what `conda install mummer` still installs, as 3.23) errors out with + "Unknown option: t" rather than ignoring it, so passing -t unconditionally + made this script fail outright on MUMmer 3. + """ + m = _load_scaffold_level_module() + exe = shutil.which('nucmer') + assert exe is not None, "nucmer not installed" + + supported = m.nucmer_supports_threads(exe) + + cmd = m.gen_mummer_cmd(exe=exe, method='mum', prefix='p', c='65', maxgap='90', + p=6, noextend='False', reference='r.fa', querry='q.fa') + + # -t is passed if and only if this nucmer understands it + assert ('-t' in cmd) == supported, \ + f"nucmer supports -t = {supported}, but command was: {' '.join(cmd)}" + + # The command must always be well formed regardless + assert cmd[0] == exe + assert cmd[-2:] == ['r.fa', 'q.fa'] + + +def test_scaffold_level_nucmer_threads_detection_is_safe(): + """An exe that doesn't exist must report 'no -t support', not raise.""" + m = _load_scaffold_level_module() + assert m.nucmer_supports_threads('/nonexistent/nucmer-does-not-exist') is False diff --git a/tests/tests/test_cluster.py b/tests/tests/test_cluster.py index f5dee4d..c9a0ad6 100644 --- a/tests/tests/test_cluster.py +++ b/tests/tests/test_cluster.py @@ -341,6 +341,27 @@ def test_skani(self): assert (db['ani'].tolist()[0] > 0.7) & (db['ani'].tolist()[0] < 0.8) +def test_skani_alignment_coverage_is_fraction_not_percent(self): + ''' + Regression test: skani reports aligned fractions as percentages, but dRep + compares alignment_coverage against cov_thresh on a 0-1 scale (see + make_linkage_Ndb). If the conversion is dropped the coverage filter silently + stops working, and distantly related genomes sharing a small conserved + region get merged. + ''' + bdb = drep.d_cluster.utils.load_genomes(self.genomes) + Ndb = drep.d_cluster.compare_utils.compare_genomes(bdb, 'skani', self.test_dir) + + assert Ndb['alignment_coverage'].between(0, 1).all(), \ + "skani alignment_coverage must be a 0-1 fraction, not a percent" + + # E. casseliflavus and E. faecalis align over only ~1% of their genomes, so + # a demanding coverage threshold must keep them in separate clusters. + Cdb, _ = drep.d_cluster.cluster_utils.genome_hierarchical_clustering( + Ndb, S_ani=0.85, cov_thresh=0.5, comp_method='skani', cluster='X') + g2c = Cdb.set_index('genome')['secondary_cluster'].to_dict() + assert g2c['Enterococcus_casseliflavus_EC20.fasta'] != g2c['Enterococcus_faecalis_T2.fna'] + @pytest.mark.skip(reason="You don't need to run this") def test_time_compare_genomes(self): ''' @@ -588,29 +609,63 @@ def test_skipsecondary(self): db2 = wd.get_db('Ndb') assert db2.empty, 'Ndb is not empty' -def test_low_ram_primary_clustering(self): +def test_mash_primary_algorithm(self): ''' - Test that low_ram_primary_clustering runs without crashing and uses the optimized method + skani is the default primary algorithm, so exercise the MASH path explicitly + to make sure it still works. ''' genomes = self.genomes wd_loc = self.wd_loc - s_wd_loc = self.s_wd_loc - # Create the work directory and data directory os.makedirs(os.path.join(wd_loc, 'data'), exist_ok=True) - # Run dRep with low_ram_primary_clustering - args = argumentParser.parse_args(['dereplicate', wd_loc, '--low_ram_primary_clustering', '-g'] + genomes) + args = argumentParser.parse_args(['dereplicate', wd_loc, '--primary_algorithm', 'MASH', + '-g'] + genomes) kwargs = vars(args) drep.d_cluster.controller.d_cluster_wrapper(wd_loc, **kwargs) - # Verify it ran by checking Cdb exists and has the right columns wd = WorkDirectory(wd_loc) Cdb = wd.get_db('Cdb') assert 'genome' in Cdb.columns assert 'primary_cluster' in Cdb.columns assert len(Cdb) > 0 - # Check that the optimized method was actually used by looking at the primary linkage - primary_linkage = wd.get_cluster('primary_linkage')['linkage'] - assert primary_linkage == "optimized_method_used", "Optimized clustering method was not used" \ No newline at end of file + # The MASH Mdb is the dense pairwise table (no alignment fractions), so + # secondary must NOT try to reuse it as if it were skani edges + Mdb = wd.get_db('Mdb') + assert 'alignment_coverage' not in Mdb.columns + + # E. faecalis genomes should land in one primary cluster, apart from E. coli + g2c = Cdb.set_index('genome')['primary_cluster'].to_dict() + assert g2c['Enterococcus_faecalis_T2.fna'] == g2c['Enterococcus_faecalis_TX0104.fa'] + assert g2c['Enterococcus_faecalis_T2.fna'] != g2c['Escherichia_coli_Sakai.fna'] +def test_cluster_colors_are_deterministic_and_cycled(): + ''' + Cluster colors exist to tell neighbouring clusters apart, not to identify a + cluster. Cycle a small palette rather than giving each cluster its own shade, + and do it deterministically -- the old implementation shuffled an unseeded + colormap, so the same analysis produced different colors every run. + ''' + from drep.d_analyze import gen_color_dictionary, CLUSTER_COLORS + + n2c = {f'g{i}': i for i in range(1, 8)} + names = list(n2c) + d = gen_color_dictionary(names, n2c) + + # only palette colors are used + assert set(d.values()) <= set(CLUSTER_COLORS) + # adjacent clusters are always distinguishable + for i in range(1, 7): + assert d[f'g{i}'] != d[f'g{i+1}'], f"clusters {i} and {i+1} share a color" + # same input -> same colors, every time + assert gen_color_dictionary(names, n2c) == d + +def test_cluster_colors_handle_secondary_cluster_names(): + '''Secondary clusters are named like "2_10"; sorting must be numeric.''' + from drep.d_analyze import gen_color_dictionary, CLUSTER_COLORS + + n2c = {'a': '2_1', 'b': '2_2', 'c': '2_10'} + d = gen_color_dictionary(list(n2c), n2c) + assert set(d.values()) <= set(CLUSTER_COLORS) + # 2_1, 2_2, 2_10 are consecutive, so they must all differ (palette has 3) + assert len({d['a'], d['b'], d['c']}) == 3 diff --git a/tests/tests/test_dereplicate.py b/tests/tests/test_dereplicate.py index 2d432f7..a5d3074 100644 --- a/tests/tests/test_dereplicate.py +++ b/tests/tests/test_dereplicate.py @@ -207,7 +207,7 @@ def test_dereplicate_4(self): # Run with chunking args = argumentParser.parse_args(['compare',wd_loc,'--S_algorithm', - 'fastANI','--SkipSecondary', '--multiround_primary_clustering', + 'fastANI','--SkipSecondary', '--primary_algorithm', 'MASH', '--multiround_primary_clustering', '--primary_chunksize', '50', '-g'] + genomes) Controller().parseArguments(args) @@ -295,7 +295,7 @@ def test_dereplicate_7(self): # Get greedy results args = argumentParser.parse_args(['compare', wd_loc2, '--S_algorithm', - 'fastANI', '--multiround_primary_clustering', '--primary_chunksize', '50', + 'fastANI', '--primary_algorithm', 'MASH', '--multiround_primary_clustering', '--primary_chunksize', '50', '--greedy_secondary_clustering', '-sa', '0.95', '-g'] + genomes) Controller().parseArguments(args) wd = WorkDirectory(wd_loc2) @@ -303,7 +303,7 @@ def test_dereplicate_7(self): # Run normal args = argumentParser.parse_args(['compare', wd_loc, '--S_algorithm', - 'fastANI', '--multiround_primary_clustering', '--primary_chunksize', '50', + 'fastANI', '--primary_algorithm', 'MASH', '--multiround_primary_clustering', '--primary_chunksize', '50', '-sa', '0.95', '-g'] + genomes) Controller().parseArguments(args) @@ -334,7 +334,7 @@ def test_dereplicate_8(self): # Get greedy results args = argumentParser.parse_args(['compare', wd_loc2, '--S_algorithm', - 'fastANI', '--multiround_primary_clustering', '--primary_chunksize', '50', + 'fastANI', '--primary_algorithm', 'MASH', '--multiround_primary_clustering', '--primary_chunksize', '50', '--greedy_secondary_clustering', '-sa', '0.95', '-pa', '0.99', '-g'] + genomes) Controller().parseArguments(args) wd = WorkDirectory(wd_loc2) @@ -342,7 +342,7 @@ def test_dereplicate_8(self): # Run normal args = argumentParser.parse_args(['compare', wd_loc, '--S_algorithm', - 'fastANI', '--multiround_primary_clustering', '--primary_chunksize', '50', + 'fastANI', '--primary_algorithm', 'MASH', '--multiround_primary_clustering', '--primary_chunksize', '50', '-sa', '0.95', '-pa', '0.99', '-g'] + genomes) Controller().parseArguments(args) diff --git a/tests/tests/test_filter.py b/tests/tests/test_filter.py index a042a44..eec3d06 100644 --- a/tests/tests/test_filter.py +++ b/tests/tests/test_filter.py @@ -333,7 +333,7 @@ def test_filer_functional_4(self): # Make sure it doesnt warn incorrectly self._caplog.clear() - args = argumentParser.parse_args(['dereplicate', self.wd_loc, '--primary_chunksize', '4', '--multiround_primary_clustering', '-g'] + self.genomes) + args = argumentParser.parse_args(['dereplicate', self.wd_loc, '--primary_chunksize', '4', '--primary_algorithm', 'MASH', '--multiround_primary_clustering', '-g'] + self.genomes) kwargs = vars(args) bdb = drep.d_cluster.utils.load_genomes(kwargs['genomes']) drep.d_filter.sanity_check(bdb, **kwargs) diff --git a/tests/tests/test_greedy.py b/tests/tests/test_greedy.py index 6af8c95..8db45bb 100644 --- a/tests/tests/test_greedy.py +++ b/tests/tests/test_greedy.py @@ -22,7 +22,7 @@ def test_multiround_primary_clustering_1(self): test_dir = self.test_dir # Run it under normal conditions - args = drep.argumentParser.parse_args(['compare', self.wd_loc, '--primary_chunksize', '3', '--multiround_primary_clustering', '--S_algorithm', 'ANImf', '-sa', '0.99', '-pa', '0.95', '-d', '-g'] + self.genomes) + args = drep.argumentParser.parse_args(['compare', self.wd_loc, '--primary_chunksize', '3', '--primary_algorithm', 'MASH', '--multiround_primary_clustering', '--S_algorithm', 'ANImf', '-sa', '0.99', '-pa', '0.95', '-d', '-g'] + self.genomes) kwargs = vars(args) drep.d_cluster.controller.d_cluster_wrapper(self.wd_loc, **kwargs) @@ -46,15 +46,15 @@ def test_multiround_primary_clustering_1(self): # Make sure it handles plotting gracefully drep.d_analyze.mash_dendrogram_from_wd(wd, plot_dir=test_dir) -def test_multiround_primary_clustering_with_low_ram(self): +def test_multiround_primary_clustering_streaming(self): """ - Test that multiround primary clustering works with low_ram_primary_clustering - and verifies both optimizations were used + Multiround primary clustering only applies to the MASH path, so it has to be + requested explicitly now that skani is the default primary algorithm. It uses + single-linkage union-find, which produces no dendrogram linkage matrix. """ test_dir = self.test_dir - # Run it with both multiround and low_ram options - args = drep.argumentParser.parse_args(['compare', self.wd_loc, '--primary_chunksize', '3', '--multiround_primary_clustering', '--low_ram_primary_clustering', '--S_algorithm', 'ANImf', '-sa', '0.99', '-pa', '0.95', '-d', '-g'] + self.genomes) + args = drep.argumentParser.parse_args(['compare', self.wd_loc, '--primary_algorithm', 'MASH', '--primary_chunksize', '3', '--multiround_primary_clustering', '--S_algorithm', 'ANImf', '-sa', '0.99', '-pa', '0.95', '-d', '-g'] + self.genomes) kwargs = vars(args) drep.d_cluster.controller.d_cluster_wrapper(self.wd_loc, **kwargs) @@ -71,9 +71,10 @@ def test_multiround_primary_clustering_with_low_ram(self): assert 'genome_chunk' in list(Mdb.columns) assert len(Mdb['genome_chunk'].unique()) == 3 - # Make sure low_ram optimization was used + # Multiround chunks carry a genome_chunk column, so no primary dendrogram is + # computed and the streaming marker is stored instead of a linkage matrix primary_linkage = wd.get_cluster('primary_linkage')['linkage'] - assert primary_linkage == "optimized_method_used", "Optimized clustering method was not used" + assert primary_linkage == "union_find_streaming", "Streaming union-find method was not used" # Make sure genomes in same primary cluster in one dataframe are also in same primary cluster in other Cdb = wd.get_db('Cdb') @@ -149,7 +150,7 @@ def test_multiround_primary_clustering_2(self): test_dir = self.test_dir # Run it under normal conditions - args = drep.argumentParser.parse_args(['dereplicate', self.wd_loc, '--primary_chunksize', '3', '--multiround_primary_clustering', '--ignoreGenomeQuality', '-pa', '0.95', '--S_algorithm', 'ANImf', '-sa', '0.99', '-d', '-g'] + self.genomes) + args = drep.argumentParser.parse_args(['dereplicate', self.wd_loc, '--primary_chunksize', '3', '--primary_algorithm', 'MASH', '--multiround_primary_clustering', '--ignoreGenomeQuality', '-pa', '0.95', '--S_algorithm', 'ANImf', '-sa', '0.99', '-d', '-g'] + self.genomes) drep.controller.Controller().parseArguments(args) # Load test results @@ -180,7 +181,7 @@ def test_multiround_primary_clustering_3(self): test_dir = self.test_dir # Run it under normal conditions - args = drep.argumentParser.parse_args(['dereplicate', self.wd_loc, '--clusterAlg', 'single', '--primary_chunksize', '3', '--multiround_primary_clustering', '--ignoreGenomeQuality', '-pa', '0.75', '-d', '-g'] + self.genomes) + args = drep.argumentParser.parse_args(['dereplicate', self.wd_loc, '--clusterAlg', 'single', '--primary_chunksize', '3', '--primary_algorithm', 'MASH', '--multiround_primary_clustering', '--ignoreGenomeQuality', '-pa', '0.75', '-d', '-g'] + self.genomes) drep.controller.Controller().parseArguments(args) # Load test results @@ -191,7 +192,7 @@ def test_multiround_primary_clustering_3(self): # Run it with a different clusterAlg shutil.rmtree(self.working_wd_loc) - args = drep.argumentParser.parse_args(['dereplicate', self.working_wd_loc, '--clusterAlg', 'complete', '--primary_chunksize', '3', '--multiround_primary_clustering', '--ignoreGenomeQuality', '-pa', '0.75', '-d', '-g'] + self.genomes) + args = drep.argumentParser.parse_args(['dereplicate', self.working_wd_loc, '--clusterAlg', 'complete', '--primary_chunksize', '3', '--primary_algorithm', 'MASH', '--multiround_primary_clustering', '--ignoreGenomeQuality', '-pa', '0.75', '-d', '-g'] + self.genomes) drep.controller.Controller().parseArguments(args) # Load test results @@ -210,7 +211,7 @@ def test_multiround_primary_clustering_4(self): test_dir = self.test_dir # Run it under normal conditions - args = drep.argumentParser.parse_args(['dereplicate', self.wd_loc, '--primary_chunksize', '3', '--multiround_primary_clustering', '--ignoreGenomeQuality', '--SkipSecondary', '-d', '-g'] + self.genomes) + args = drep.argumentParser.parse_args(['dereplicate', self.wd_loc, '--primary_chunksize', '3', '--primary_algorithm', 'MASH', '--multiround_primary_clustering', '--ignoreGenomeQuality', '--SkipSecondary', '-d', '-g'] + self.genomes) drep.controller.Controller().parseArguments(args) # Load test results diff --git a/tests/tests/test_union_find.py b/tests/tests/test_union_find.py new file mode 100644 index 0000000..6357c7c --- /dev/null +++ b/tests/tests/test_union_find.py @@ -0,0 +1,270 @@ +""" +Unit tests for streaming union-find primary clustering (drep.d_cluster.union_find). +""" +import glob +import os +import shutil +import tempfile + +import numpy as np +import pandas as pd +import pytest + +import drep.d_cluster.union_find as uf +import drep.d_cluster.compare_utils as cu +import drep.d_cluster.utils + + +def _test_genomes(): + here = os.path.dirname(os.path.abspath(__file__)) + return [g for g in glob.glob(os.path.join(here, '../genomes/*')) + if os.path.isfile(g)] + + +def test_union_find_basic(): + u = uf.UnionFind() + for x in 'abcde': + u.add(x) + u.union('a', 'b') + u.union('b', 'c') + u.union('d', 'e') + comps = {frozenset(v) for v in u.components().values()} + assert comps == {frozenset('abc'), frozenset('de')} + + +def test_cluster_long_df_matches_expectation(): + # a-b close, c-d close, everything else far; e is a singleton + rows = [ + ('a', 'b', 0.01), ('b', 'a', 0.01), + ('c', 'd', 0.02), ('d', 'c', 0.02), + ('a', 'c', 0.30), ('a', 'd', 0.30), ('a', 'e', 0.30), + ('b', 'c', 0.30), ('b', 'e', 0.30), ('c', 'e', 0.30), + ('d', 'e', 0.30), + ] + db = pd.DataFrame(rows, columns=['genome1', 'genome2', 'dist']) + Cdb = uf.cluster_long_df(db, cutoff=0.1, all_genomes=list('abcde')) + g2c = Cdb.set_index('genome')['primary_cluster'].to_dict() + assert g2c['a'] == g2c['b'] + assert g2c['c'] == g2c['d'] + assert g2c['a'] != g2c['c'] + assert g2c['e'] != g2c['a'] and g2c['e'] != g2c['c'] + assert set(Cdb['genome']) == set('abcde') + # deterministic numbering: largest clusters first + assert Cdb['primary_cluster'].min() == 1 + + +def test_cluster_mash_files_streaming(tmp_path): + # write a small symmetric mash-style dist tsv + f = tmp_path / "mash.tsv" + names = [f"g{i}.fasta" for i in range(6)] + block = [0, 0, 0, 1, 1, 1] # two true clusters of 3 + with open(f, 'w') as o: + for i, gi in enumerate(names): + for j, gj in enumerate(names): + d = 0.0 if i == j else (0.01 if block[i] == block[j] else 0.30) + o.write(f"{gi}\t{gj}\t{d:.4f}\t0\t1000/1000\n") + + Cdb, stats = uf.cluster_mash_files(str(f), cutoff=0.1) + assert stats['total_pairs_streamed'] == 36 + assert Cdb['primary_cluster'].nunique() == 2 + # genome names keep their basename (incl. extension), like parse_mash_table + assert set(Cdb['genome']) == {f"g{i}.fasta" for i in range(6)} + g2c = Cdb.set_index('genome')['primary_cluster'].to_dict() + assert g2c['g0.fasta'] == g2c['g1.fasta'] == g2c['g2.fasta'] + assert g2c['g3.fasta'] == g2c['g4.fasta'] == g2c['g5.fasta'] + assert g2c['g0.fasta'] != g2c['g3.fasta'] + + +def test_union_find_matches_scipy_membership(): + # Build a random symmetric similarity table and confirm union-find + # and scipy single-linkage produce identical cluster membership. + rng = np.random.default_rng(1) + n = 40 + block = np.arange(n) // 4 + sim = np.where(block[:, None] == block[None, :], 0.99, 0.70) + noise = np.triu(rng.normal(0, 0.01, (n, n)), 1) + sim = np.clip(sim + noise + noise.T, 0, 1) + np.fill_diagonal(sim, 1.0) + + rows = [] + for i in range(n): + for j in range(n): + rows.append((f"g{i}", f"g{j}", sim[i, j])) + db = pd.DataFrame(rows, columns=['genome1', 'genome2', 'similarity']) + + scipy_Cdb, _ = cu.cluster_mash_database(db.copy(), P_ani=0.9, + primary_clusterAlg='average', + classic_primary_clustering=True) + uf_Cdb, _ = cu.cluster_mash_database(db.copy(), P_ani=0.9, + primary_clusterAlg='single') + + def membership(Cdb): + return {frozenset(sub['genome']) for _, sub in Cdb.groupby('primary_cluster')} + + assert membership(scipy_Cdb) == membership(uf_Cdb) + + +def test_skani_sparse_min_af_filters_low_alignment_edges(): + """ + The aligned-fraction filter is what keeps skani ANI comparable to MASH's + whole-genome similarity. skani reports ANI within aligned regions only, so a + pair sharing one small conserved region looks like a high-ANI edge; under + single linkage a few such bridges chain unrelated genomes into one giant + primary cluster (measured on 10k UHGG genomes: min-af 0 collapsed 59% of the + dataset into a single cluster). + + Here a and b are genuinely similar, and c is joined to each only by a + high-ANI/low-alignment bridge. With the filter on, c must stay separate. + """ + import tempfile + with tempfile.TemporaryDirectory() as td: + f = os.path.join(td, 'sparse.tsv') + rows = [ + # ref, query, ANI, af_ref, af_query + ('a.fna', 'b.fna', 99.0, 90.0, 92.0), # real relationship + ('a.fna', 'c.fna', 98.0, 2.0, 3.0), # spurious bridge (tiny overlap) + ('b.fna', 'c.fna', 97.5, 2.5, 2.0), # spurious bridge + ] + with open(f, 'w') as o: + o.write("Ref_file\tQuery_file\tANI\tAlign_fraction_ref\tAlign_fraction_query\n") + for r in rows: + o.write("\t".join(str(x) for x in r) + "\n") + + allg = ['a.fna', 'b.fna', 'c.fna'] + + # No filter: the bridges chain a, b and c into one cluster + C0, _, _ = uf.cluster_skani_sparse_files(f, 90.0, allg, cov_threshold=0.0) + assert C0['primary_cluster'].nunique() == 1 + + # With a 15% aligned-fraction floor, c is correctly left on its own + C1, _, s1 = uf.cluster_skani_sparse_files(f, 90.0, allg, cov_threshold=0.15) + g2c = C1.set_index('genome')['primary_cluster'].to_dict() + assert g2c['a.fna'] == g2c['b.fna'] + assert g2c['c.fna'] != g2c['a.fna'] + assert s1['edges_kept'] == 1 + + +def test_build_ndb_from_edges_fills_matrix(): + """ + Secondary clustering needs a complete matrix per primary cluster, so pairs + absent from the sparse edge list must come back as ani=0 and self-pairs as 1. + """ + edges = pd.DataFrame({ + 'genome1': ['a', 'b'], + 'genome2': ['b', 'a'], + 'ani': [0.99, 0.99], + 'alignment_coverage': [0.9, 0.92], + }) + Cdb = pd.DataFrame({'genome': ['a', 'b', 'c'], 'primary_cluster': [1, 1, 1]}) + Ndb = uf.build_ndb_from_edges(edges, Cdb) + + # complete 3x3 matrix for the one primary cluster + assert len(Ndb) == 9 + g = Ndb.set_index(['reference', 'querry']) + assert g.loc[('a', 'b'), 'ani'] == pytest.approx(0.99) + assert g.loc[('a', 'a'), 'ani'] == 1.0 + assert g.loc[('a', 'a'), 'alignment_coverage'] == 1.0 + # c had no edges -> filled as no similarity + assert g.loc[('a', 'c'), 'ani'] == 0.0 + assert g.loc[('c', 'a'), 'ani'] == 0.0 + # coverage is directional: fraction of the 'reference' genome + assert g.loc[('a', 'b'), 'alignment_coverage'] == pytest.approx(0.9) + assert g.loc[('b', 'a'), 'alignment_coverage'] == pytest.approx(0.92) + + +def test_build_ndb_from_edges_only_within_primary_clusters(): + """Secondary never compares across primary clusters.""" + edges = pd.DataFrame({ + 'genome1': ['a', 'b'], 'genome2': ['b', 'a'], + 'ani': [0.99, 0.99], 'alignment_coverage': [0.9, 0.9], + }) + Cdb = pd.DataFrame({'genome': ['a', 'b'], 'primary_cluster': [1, 2]}) + Ndb = uf.build_ndb_from_edges(edges, Cdb) + # a and b are in different primary clusters: only self-comparisons survive + assert set(zip(Ndb['reference'], Ndb['querry'])) == {('a', 'a'), ('b', 'b')} + + +@pytest.mark.skipif(shutil.which('skani') is None, reason="skani not installed") +def test_reused_edges_match_rerunning_skani(): + """ + The one-pass path must give exactly what re-running skani per primary cluster + gives -- that is the whole premise for skipping the second pass. + """ + import drep.d_cluster.compare_utils as cu + import drep.d_cluster.utils + + genomes = _test_genomes() + Bdb = drep.d_cluster.utils.load_genomes(genomes) + workdir = tempfile.mkdtemp() + try: + Mdb, Cdb, _ = cu.primary_cluster_skani_sparse( + Bdb, os.path.join(workdir, 'p'), P_ani=0.9, S_ani=0.99, + cov_thresh=0.1, processors=4, primary_progress=False) + + # reuse primary's edges + Ndb_r, Cdb_r, _ = cu.secondary_clustering_from_primary_edges( + Bdb, Cdb, Mdb, S_ani=0.99, cov_thresh=0.1, clusterAlg='average') + + # re-run skani per primary cluster (the classic path) + Ndb_c, Cdb_c, _ = cu.secondary_clustering( + Bdb, Cdb, 'skani', os.path.join(workdir, 's'), + S_ani=0.99, cov_thresh=0.1, clusterAlg='average', processors=4) + + def part(C): + return {frozenset(s['genome']) for _, s in C.groupby('secondary_cluster')} + + assert part(Cdb_r) == part(Cdb_c) + finally: + shutil.rmtree(workdir, ignore_errors=True) + + +@pytest.mark.skipif(shutil.which('skani') is None, reason="skani not installed") +def test_sparse_skani_primary_matches_mash(): + """ + Sparse-skani primary clustering should recover the same primary partition as + the classic MASH path on the bundled test genomes. + """ + genomes = _test_genomes() + Bdb = drep.d_cluster.utils.load_genomes(genomes) + workdir = tempfile.mkdtemp() + try: + _, Cdb_sk, cret = cu.primary_cluster_skani_sparse( + Bdb, os.path.join(workdir, 'sk'), P_ani=0.9, processors=4, + primary_progress=False) + _, Cdb_mash, _ = cu.all_vs_all_MASH( + Bdb, os.path.join(workdir, 'mash'), P_ani=0.9, processors=4) + + # Every input genome is represented (including singletons skani screened out) + assert set(Cdb_sk['genome']) == set(Bdb['genome']) + # Small genome sets still get a scipy linkage so the primary dendrogram + # can be drawn, and it must cover every genome -- not just the ones that + # appear in the sparse edge list + assert not isinstance(cret[0], str), "expected a real linkage matrix for a small set" + assert list(cret[1].columns) == sorted(Bdb['genome']) + + def part(Cdb): + return {frozenset(sub['genome']) for _, sub in Cdb.groupby('primary_cluster')} + + assert part(Cdb_sk) == part(Cdb_mash) + finally: + shutil.rmtree(workdir, ignore_errors=True) + + +def test_classic_primary_clustering_uses_dense_path(): + """--classic_primary_clustering must produce a real scipy linkage matrix.""" + rng = np.random.default_rng(2) + n = 20 + block = np.arange(n) // 4 + sim = np.where(block[:, None] == block[None, :], 0.99, 0.70) + noise = np.triu(rng.normal(0, 0.01, (n, n)), 1) + sim = np.clip(sim + noise + noise.T, 0, 1) + np.fill_diagonal(sim, 1.0) + rows = [(f"g{i:02d}", f"g{j:02d}", sim[i, j]) for i in range(n) for j in range(n)] + db = pd.DataFrame(rows, columns=['genome1', 'genome2', 'similarity']) + + _, cret = cu.cluster_mash_database(db.copy(), P_ani=0.9, + primary_clusterAlg='average', + classic_primary_clustering=True) + # cret[0] is a real linkage matrix (ndarray), not the streaming marker + assert not isinstance(cret[0], str) + assert cret[2]['linkage_method'] == 'average'