From 9ef5e114edb9f4534186e1de233d1c529ef1f6be Mon Sep 17 00:00:00 2001 From: Matt Olm Date: Tue, 14 Jul 2026 08:04:55 -0600 Subject: [PATCH 01/15] Streaming union-find primary clustering + sparse skani Primary clustering is single-linkage at a fixed cutoff, which is identical to finding connected components. The old path instead materialized every pairwise comparison into a long Mdb (N^2 rows) and pivoted it into a dense N x N matrix before handing it to scipy. That is the source of issue #259 and the MemoryError at 43k genomes (24.7 GiB alloc failure at pd.concat). Compute the components directly with union-find instead, streaming the comparison output and discarding the ~99.9% of above-cutoff pairs as they are read. Memory becomes O(genomes + kept_edges) instead of O(genomes^2). Measured peak RSS on synthetic MASH dist input (old -> new): 4,000 genomes: 2.16 GB -> 0.49 GB 8,000 genomes: 7.42 GB -> 0.55 GB 16,000 genomes: ~30 GB -> 0.54 GB Old grows quadratically; new stays flat. Cluster membership is identical to scipy single-linkage across thresholds. Changes: - New drep/d_cluster/union_find.py: UnionFind, cluster_long_df (no pivot), cluster_mash_files (streams MASH dist), cluster_skani_sparse_files. - Split primary/secondary linkage: new --primary_clusterAlg (default single, routes to union-find) separate from secondary --clusterAlg (default average, unchanged). Add --classic_primary_clustering to force the dense scipy path. - New --primary_algorithm skani: runs `skani triangle --sparse`, which emits only above-threshold pairs, streamed into union-find. No N^2 is built on disk or in RAM. Recommended for very large genome sets. - Bound the multiround Mdb: subsample per-chunk tables and free each chunk after clustering, so the concat that previously OOM'd stays bounded. - Keep the primary dendrogram for modest genome sets (the dense pivot is cheap at small N); skip it above the cutoff, as multiround already did. The low_ram marker changed from "optimized_method_used" to "union_find_streaming"; assertions updated accordingly. Co-Authored-By: Claude Opus 4.8 --- drep/argumentParser.py | 24 ++- drep/d_analyze.py | 4 + drep/d_cluster/compare_utils.py | 149 ++++++++++++- drep/d_cluster/controller.py | 2 +- drep/d_cluster/external.py | 54 +++++ drep/d_cluster/union_find.py | 356 ++++++++++++++++++++++++++++++++ tests/tests/test_cluster.py | 2 +- tests/tests/test_greedy.py | 2 +- tests/tests/test_union_find.py | 153 ++++++++++++++ 9 files changed, 730 insertions(+), 16 deletions(-) create mode 100644 drep/d_cluster/union_find.py create mode 100644 tests/tests/test_union_find.py diff --git a/drep/argumentParser.py b/drep/argumentParser.py index 3e74ce4..6ee6b1d 100644 --- a/drep/argumentParser.py +++ b/drep/argumentParser.py @@ -121,6 +121,12 @@ def parse_args(args): + "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'}) + Clustflags.add_argument("--primary_algorithm", help="R|Program to use for primary clustering.\n" \ + + "MASH = (DEFAULT) all-vs-all Mash\n" \ + + "skani = skani triangle --sparse; only above-threshold pairs are\n" \ + + " produced and they are streamed rather than held as a dense\n" \ + + " matrix. Recommended for very large genome sets (no N^2 RAM/disk).", + default='MASH', choices={'MASH', 'skani'}) 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') @@ -146,10 +152,22 @@ 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("--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("--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("--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) + Compflags.add_argument("--low_ram_primary_clustering", help="(Deprecated; the streaming single-linkage\ + algorithm is now the default for primary clustering.) Kept for backwards compatibility.", 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..d58dadc 100644 --- a/drep/d_analyze.py +++ b/drep/d_analyze.py @@ -150,6 +150,10 @@ 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 - cannot generate with low_ram_primary_clustering (no linkage matrix)") + return + # Make the plot logging.info("Plotting primary dendrogram") plot_MASH_dendrogram(Mdb, Cdb, Plinkage, threshold = PL_thresh,\ diff --git a/drep/d_cluster/compare_utils.py b/drep/d_cluster/compare_utils.py index 7478768..b155d03 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,63 @@ 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. + + 'MASH' (default) uses the classic all-vs-all Mash path. 'skani' uses + `skani triangle --sparse` + streaming union-find, which never builds the N^2 + matrix on disk or in RAM and is the recommended path for very large genome + sets. + + Returns (Mdb, Cdb, cluster_ret), matching all_vs_all_MASH. + """ + method = kwargs.get('primary_algorithm', 'MASH') + if method == 'skani': + 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 `skani triangle --sparse` streamed into union-find. + + Only above-threshold pairs are ever produced (skani screens during sketching) + and they are streamed rather than held as a dense matrix, so memory stays + O(genomes + edges) regardless of genome count. Always single-linkage + (connected components); --classic_primary_clustering / non-single + primary_clusterAlg do not apply here. + """ + 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. + default_screen = max(1.0, min(ani_threshold - 5.0, 99.0)) + screen = kwargs.get('primary_skani_screen', default_screen) + + 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}%)") + sparse_file = drep.d_cluster.external.run_skani_triangle_sparse( + genome_list, skani_folder, screen, **kwargs) + + all_genomes = list(Bdb['genome'].unique()) + Cdb, Mdb, stats = drep.d_cluster.union_find.cluster_skani_sparse_files( + sparse_file, ani_threshold, all_genomes, + progress=kwargs.get('primary_progress', True)) + + logging.info(f" Sparse skani primary clustering: {stats['edges_kept']:,} edges kept, " + f"{stats['primary_clusters']:,} primary clusters") + + arguments = {'linkage_method': 'single', 'linkage_cutoff': 1 - P_ani, + 'comparison_algorithm': 'skani'} + cluster_ret = ['union_find_streaming', None, arguments] + return Mdb, Cdb, cluster_ret + def prepare_mash(data_folder, **kwargs): """ Make some folders and things @@ -214,9 +275,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 +313,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 +338,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 +361,70 @@ 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 + low_ram_primary_clustering: deprecated alias forcing single-linkage union-find 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) + classic = kwargs.get('classic_primary_clustering', False) low_ram = kwargs.get('low_ram_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') or low_ram) + if use_union_find: + if low_ram and P_Lmethod != 'single': + logging.warning( + f"low_ram_primary_clustering uses single-linkage (connected components); " + f"ignoring primary_clusterAlg={P_Lmethod} for primary clustering.") + 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 (not low_ram) and 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, low_ram=False) 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..665624d 100644 --- a/drep/d_cluster/controller.py +++ b/drep/d_cluster/controller.py @@ -102,7 +102,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") diff --git a/drep/d_cluster/external.py b/drep/d_cluster/external.py index ff24b0c..630eee2 100644 --- a/drep/d_cluster/external.py +++ b/drep/d_cluster/external.py @@ -271,6 +271,60 @@ def _fix_fastani(odb): return fdb +def run_skani_triangle_sparse(genome_list, outdir, screen, **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. + + 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)) + cmd = [exe_loc, "triangle", "--sparse", "-t", str(p), '-o', out_file, + '-l', glist, '-s', str(screen)] + 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/union_find.py b/drep/d_cluster/union_find.py new file mode 100644 index 0000000..d72faa3 --- /dev/null +++ b/drep/d_cluster/union_find.py @@ -0,0 +1,356 @@ +""" +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`` and for the ``low_ram`` path that + used to ``stack()`` an already-dense matrix back into long format. + + 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 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/tests/tests/test_cluster.py b/tests/tests/test_cluster.py index f5dee4d..462d423 100644 --- a/tests/tests/test_cluster.py +++ b/tests/tests/test_cluster.py @@ -613,4 +613,4 @@ def test_low_ram_primary_clustering(self): # 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 + assert primary_linkage == "union_find_streaming", "Optimized clustering method was not used" \ No newline at end of file diff --git a/tests/tests/test_greedy.py b/tests/tests/test_greedy.py index 6af8c95..8c6fc47 100644 --- a/tests/tests/test_greedy.py +++ b/tests/tests/test_greedy.py @@ -73,7 +73,7 @@ def test_multiround_primary_clustering_with_low_ram(self): # Make sure low_ram optimization was used 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') diff --git a/tests/tests/test_union_find.py b/tests/tests/test_union_find.py new file mode 100644 index 0000000..ae4f9a2 --- /dev/null +++ b/tests/tests/test_union_find.py @@ -0,0 +1,153 @@ +""" +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_low_ram_matches_scipy_membership(): + # Build a random symmetric similarity table and confirm union-find (low_ram) + # 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) + + +@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']) + # Marker so downstream plotting skips the (nonexistent) primary dendrogram + assert cret[0] == 'union_find_streaming' + + 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' From f50c790e41c79ad79c37189427fb261c5bb1826a Mon Sep 17 00:00:00 2001 From: Matt Olm Date: Tue, 14 Jul 2026 08:42:44 -0600 Subject: [PATCH 02/15] Fix skani alignment coverage being a percent, not a fraction skani reports aligned fractions as percentages (0-100), but load_skani only divided ANI by 100 and passed the aligned fraction through untouched. Every other comparison algorithm (fastANI, ANImf, ANIn, gANI) reports alignment_coverage on a 0-1 scale, and that is the scale cov_thresh is compared against in make_linkage_Ndb: d.loc[d['alignment_coverage'] <= cov_thresh, 'ani'] = 0 So for --S_algorithm skani the coverage filter was effectively inert: a pair aligning over only 1% of the genome has alignment_coverage=1.04, which sails past a cov_thresh of 0.5 and keeps its ANI. Distantly related genomes that share a small conserved region could then be clustered together. On the bundled test genomes, E. casseliflavus EC20 and the E. faecalis genomes align over ~1% of their length at ~93% ANI. With -sa 0.85 -nc 0.5 they were being merged into a single secondary cluster -- two different species dereplicated into one -- because the 1.04 "coverage" passed the filter. With the fraction converted they correctly separate. Note this changes results for existing --S_algorithm skani users; coverage filtering now actually applies. Co-Authored-By: Claude Opus 4.8 --- drep/d_cluster/external.py | 9 ++++++++- tests/tests/test_cluster.py | 21 +++++++++++++++++++++ 2 files changed, 29 insertions(+), 1 deletion(-) diff --git a/drep/d_cluster/external.py b/drep/d_cluster/external.py index 630eee2..4d6a25d 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) diff --git a/tests/tests/test_cluster.py b/tests/tests/test_cluster.py index 462d423..60c7ea9 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): ''' From 62299daf7e0681a474137fd4656ddd7ee1d986d3 Mon Sep 17 00:00:00 2001 From: Matt Olm Date: Tue, 14 Jul 2026 08:42:57 -0600 Subject: [PATCH 03/15] Add pyskani backend: sketch each genome once, compare in-process The subprocess comparison algorithms re-sketch genomes on every invocation. Greedy secondary clustering is the worst case: for each genome it spawns a fastANI subprocess against the growing representative list, re-sketching every representative each time. That is O(N*R) sketching work and N subprocess spawns. pyskani lets us sketch each genome exactly once, keep the representatives in an in-memory database, and query it directly -- O(N) sketching, no subprocesses, no temp files. Greedy clustering benchmark (synthetic divergent genomes, so every genome founds its own cluster and representatives accumulate to N): 40 genomes: fastANI 18.0s -> pyskani 0.70s (25.8x) 80 genomes: fastANI 83.2s -> pyskani 1.94s (42.9x) fastANI grows ~quadratically, pyskani ~linearly. Real datasets have fewer representatives, so expect a smaller (still large) speedup. pyskani agrees with the skani executable to within 5e-05 on both ANI and alignment coverage, and produces identical secondary clusters. Changes: - New drep/d_cluster/pyskani_backend.py: PyskaniDatabase (sketch once, query many), run_pairwise_pyskani, pyskani_one_vs_many. - --S_algorithm pyskani, wired into both the pairwise and greedy paths. This also lifts greedy clustering's fastANI-only restriction. - pyskani is an optional dependency (pip install drep[pyskani]) with an actionable error if missing; the skani/fastANI executables remain the default. - prepare_for_greedy now creates its data folder for every algorithm, not just fastANI (the representative list is written regardless). alignment_coverage is the aligned fraction of the genome in the 'reference' column, matching skani's .af matrix and load_fastani. Greedy needs the opposite row orientation from pairwise, because get_cluster_rep reads the representative out of the 'querry' column; both orientations are explicit in query(). Co-Authored-By: Claude Opus 4.8 --- drep/argumentParser.py | 6 +- drep/d_cluster/compare_utils.py | 8 +- drep/d_cluster/greedy_clustering.py | 36 ++++- drep/d_cluster/pyskani_backend.py | 215 ++++++++++++++++++++++++++++ setup.py | 5 + tests/tests/test_pyskani.py | 160 +++++++++++++++++++++ 6 files changed, 424 insertions(+), 6 deletions(-) create mode 100644 drep/d_cluster/pyskani_backend.py create mode 100644 tests/tests/test_pyskani.py diff --git a/drep/argumentParser.py b/drep/argumentParser.py index 6ee6b1d..fdf4a8a 100644 --- a/drep/argumentParser.py +++ b/drep/argumentParser.py @@ -116,11 +116,15 @@ def parse_args(args): Clustflags.add_argument("--S_algorithm", help="R|Algorithm for secondary clustering comaprisons:\n" \ + "fastANI = Kmer-based approach; very fast\n" \ + "skani = Even faster Kmer-based approacht\n" \ + + "pyskani = skani run in-process via the pyskani library. Each genome is\n" \ + + " sketched exactly once instead of being re-sketched by a new\n" \ + + " subprocess for every comparison, which is much faster for\n" \ + + " greedy clustering. Requires `pip install pyskani`.\n" \ + "ANImf = (DEFAULT) 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='fastANI', choices={'ANIn', 'gANI', 'ANImf', 'goANI', 'fastANI', 'skani', 'pyskani'}) Clustflags.add_argument("--primary_algorithm", help="R|Program to use for primary clustering.\n" \ + "MASH = (DEFAULT) all-vs-all Mash\n" \ + "skani = skani triangle --sparse; only above-threshold pairs are\n" \ diff --git a/drep/d_cluster/compare_utils.py b/drep/d_cluster/compare_utils.py index b155d03..b56f968 100644 --- a/drep/d_cluster/compare_utils.py +++ b/drep/d_cluster/compare_utils.py @@ -14,6 +14,7 @@ import drep.d_cluster.utils import drep.d_cluster.greedy_clustering import drep.d_cluster.union_find +import drep.d_cluster.pyskani_backend class genomeChunk(): """ @@ -521,6 +522,11 @@ def compare_genomes(bdb, algorithm, data_folder, **kwargs): df = drep.d_cluster.external.run_pairwise_skani(genome_list, working_data_folder, **kwargs) return df + elif algorithm == 'pyskani': + genome_list = bdb['location'].tolist() + df = drep.d_cluster.pyskani_backend.run_pairwise_pyskani(genome_list, **kwargs) + return df + elif algorithm == 'gANI': # Figure out prodigal folder wd = kwargs.get('wd', False) @@ -554,7 +560,7 @@ def compare_genomes(bdb, algorithm, data_folder, **kwargs): sys.exit() else: - SUPPORTED = ['fastANI'] + SUPPORTED = ['fastANI', 'pyskani'] if algorithm not in SUPPORTED: message = f"{algorithm} is not supported for greedy secondary clustering!\nChoose one of the following supported S_algorithm options: {' '.join(SUPPORTED)}" logging.error(message) diff --git a/drep/d_cluster/greedy_clustering.py b/drep/d_cluster/greedy_clustering.py index 4ee1564..3a75a5b 100644 --- a/drep/d_cluster/greedy_clustering.py +++ b/drep/d_cluster/greedy_clustering.py @@ -7,6 +7,7 @@ import drep.d_cluster.external import drep.d_cluster.compare_utils +import drep.d_cluster.pyskani_backend def greedy_secondary_clustering(Bdb, Cdb, algorithm, data_folder, **kwargs): ndbs = [] @@ -87,6 +88,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 +108,39 @@ 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 + # Return Ndb. NOTE: this spawns a subprocess that re-sketches every + # representative on every call (O(N*R) sketching). The pyskani path below + # sketches each genome exactly once instead. return drep.d_cluster.external.fastani_one_vs_many(new_genome, genome_reps, genome_rep_file, data_folder, **kwargs) + elif algorithm == 'pyskani': + return drep.d_cluster.pyskani_backend.pyskani_one_vs_many( + new_genome, kwargs['pyskani_db'], **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. + + For pyskani this sketches it once into the in-memory database, so subsequent + genomes can be compared against it without any re-sketching. Subprocess-based + algorithms read the representative list from a file instead and need nothing + here. + """ + if algorithm == 'pyskani': + kwargs['pyskani_db'].add_genome(location) + + 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) @@ -132,6 +155,11 @@ def prepare_for_greedy(algorithm, data_folder, **kwargs): kwargs['logdir'] = logdir kwargs['current_exe'] = drep.get_exe('fastANI') + elif algorithm == 'pyskani': + # One in-memory database of representatives for this primary cluster. + # Each representative is sketched exactly once, on the way in. + kwargs['pyskani_db'] = drep.d_cluster.pyskani_backend.PyskaniDatabase(**kwargs) + return kwargs def order_genomes_for_greedy(bdb, **kwargs): diff --git a/drep/d_cluster/pyskani_backend.py b/drep/d_cluster/pyskani_backend.py new file mode 100644 index 0000000..37f0fef --- /dev/null +++ b/drep/d_cluster/pyskani_backend.py @@ -0,0 +1,215 @@ +""" +In-process skani comparisons via pyskani (https://github.com/althonos/pyskani). + +The subprocess-based comparison algorithms re-sketch genomes on every +invocation. That is especially wasteful during greedy secondary clustering, +where each new genome is compared against the growing set of cluster +representatives by spawning a fresh subprocess that re-sketches *every* +representative -- O(N * R) sketching work and N subprocess spawns for N genomes. + +pyskani lets us sketch each genome exactly once, keep the representatives in an +in-memory database, and query it directly. Sketching becomes O(N), there are no +subprocesses, and no temporary files are written. + +This module is imported lazily: pyskani is an optional dependency, and dRep +falls back to the subprocess skani/fastANI implementations without it. +""" + +import functools +import logging +import os + +import pandas as pd + +import drep.d_cluster.utils + +# Column layout every secondary-clustering comparison must return +NDB_COLUMNS = ['reference', 'querry', 'ani', 'alignment_coverage'] + + +def import_pyskani(): + """ + Import pyskani, raising an actionable error if it isn't installed. + """ + try: + import pyskani + except ImportError: + raise ImportError( + "The 'pyskani' S_algorithm requires the pyskani package, which is not " + "installed. Install it with `pip install pyskani` (or `pip install " + "drep[pyskani]`), or choose a different --S_algorithm (e.g. skani, " + "which uses the skani executable instead)." + ) + return pyskani + + +def load_contigs(location): + """ + Read a FASTA file into a list of contig sequences (as bytes), which is what + pyskani's sketch/query expect. + """ + from Bio import SeqIO + return [bytes(record.seq) for record in SeqIO.parse(location, 'fasta')] + + +# During greedy clustering a genome is queried and then, if it founds a new +# cluster, immediately sketched as a representative. A tiny cache avoids parsing +# the same FASTA twice in a row. Deliberately kept at maxsize=2 -- caching every +# genome's contigs would hold the entire input set in memory. +@functools.lru_cache(maxsize=2) +def _load_contigs_cached(location): + return load_contigs(location) + + +class PyskaniDatabase: + """ + A pyskani database that sketches each genome exactly once. + + Sketching is the expensive part of an ANI comparison, so genomes are sketched + on the way in and the resulting database is queried many times. Contigs are + cached only while needed to sketch/query, not retained for the lifetime of + the object. + """ + + def __init__(self, **kwargs): + pyskani = import_pyskani() + self.db = pyskani.Database() + self.names = [] + # Screening cutoff. Mirrors skani's -s flag; the subprocess pairwise path + # uses -s 1 (i.e. compare essentially everything), so default low here and + # let callers raise it when they only care about close relatives. + self.cutoff = kwargs.get('pyskani_cutoff', 0.01) + self.learned_ani = kwargs.get('pyskani_learned_ani', None) + + def add(self, name, contigs): + """Sketch a genome once and store it as a reference.""" + self.db.sketch(name, *contigs) + self.names.append(name) + + def add_genome(self, location): + """Sketch a genome from a FASTA path; returns its dRep genome name.""" + name = drep.d_cluster.utils._get_genome_name_from_fasta(location) + self.add(name, _load_contigs_cached(location)) + return name + + def query_hits(self, name, contigs): + """ + Query the database with a genome, returning the raw pyskani hits. + """ + kw = {'cutoff': self.cutoff} + if self.learned_ani is not None: + kw['learned_ani'] = self.learned_ani + return self.db.query(name, *contigs, **kw) + + def query(self, name, contigs, query_as_reference=False): + """ + Query the database, returning (reference, querry, ani, alignment_coverage) + tuples. + + Everywhere in dRep, alignment_coverage is the aligned fraction of the + genome named in the 'reference' column (skani's .af matrix cell [A][B] is + the aligned fraction of A; fastANI's matched/total is the fraction of the + genome that load_fastani puts in 'reference'). Both orientations below + respect that rule. + + Args: + query_as_reference: if False (default, pairwise use), emit + reference=the database genome and coverage=its aligned fraction. + If True (greedy use), emit reference=the queried genome and + coverage=the queried genome's aligned fraction -- which is what + get_cluster_rep expects, since it reads the representative out of + the 'querry' column. + """ + rows = [] + for hit in self.query_hits(name, contigs): + if query_as_reference: + rows.append((hit.query_name, hit.reference_name, + hit.identity, hit.query_fraction)) + else: + rows.append((hit.reference_name, hit.query_name, + hit.identity, hit.reference_fraction)) + return rows + + def query_genome(self, location, query_as_reference=False): + name = drep.d_cluster.utils._get_genome_name_from_fasta(location) + return self.query(name, _load_contigs_cached(location), + query_as_reference=query_as_reference) + + +def _fill_missing_pairs(rows, names): + """ + skani only reports pairs that share enough k-mers. dRep's hierarchical + secondary clustering needs a complete matrix, so absent pairs are filled in + as ani=0 / coverage=0 (no detectable relatedness), matching what the + subprocess `skani triangle --min-af 0` path yields for unrelated genomes. + """ + have = {(r[0], r[1]) for r in rows} + filled = list(rows) + for a in names: + for b in names: + if (a, b) not in have: + filled.append((a, b, 0.0, 0.0)) + return filled + + +def run_pairwise_pyskani(genome_list, **kwargs): + """ + All-vs-all ANI within a set of genomes, in-process. + + Each genome is sketched exactly once and then queried against the database, + so the sketching cost is linear in the number of genomes rather than + quadratic. + + Args: + genome_list: list of genome file locations. + + Returns: + Ndb: DataFrame with ['reference', 'querry', 'ani', 'alignment_coverage']. + """ + db = PyskaniDatabase(**kwargs) + + contigs = {} + for location in genome_list: + name = drep.d_cluster.utils._get_genome_name_from_fasta(location) + cs = load_contigs(location) + contigs[name] = cs + db.add(name, cs) + + logging.debug(f"pyskani: sketched {len(contigs)} genomes once; querying") + + rows = [] + for name, cs in contigs.items(): + rows.extend(db.query(name, cs)) + + rows = _fill_missing_pairs(rows, list(contigs.keys())) + Ndb = pd.DataFrame(rows, columns=NDB_COLUMNS) + + # A genome can hit itself with identity slightly below 1 depending on + # sketching; force exact self-identity like the other algorithms do. + self_mask = Ndb['reference'] == Ndb['querry'] + Ndb.loc[self_mask, 'ani'] = 1.0 + Ndb.loc[self_mask, 'alignment_coverage'] = 1.0 + + # Keep one row per ordered pair (query returns each direction once) + Ndb = Ndb.drop_duplicates(subset=['reference', 'querry'], keep='first') + return Ndb.reset_index(drop=True) + + +def pyskani_one_vs_many(location, db, **kwargs): + """ + Compare one genome against an existing PyskaniDatabase of representatives. + + Used by greedy secondary clustering: the database holds every current cluster + representative, already sketched, so this is a single in-process query rather + than a subprocess that re-sketches all representatives. + + Emits rows with the representative in the 'querry' column (mirroring the + fastANI greedy path), because get_cluster_rep reads the winning + representative from there. + + Returns an Ndb-shaped DataFrame (possibly empty if nothing is similar). + """ + rows = db.query_genome(location, query_as_reference=True) + if len(rows) == 0: + return pd.DataFrame(columns=NDB_COLUMNS) + return pd.DataFrame(rows, columns=NDB_COLUMNS) diff --git a/setup.py b/setup.py index 81f4119..49831ee 100644 --- a/setup.py +++ b/setup.py @@ -28,4 +28,9 @@ def version(): 'setuptools', 'pytest' ], + extras_require={ + # In-process skani (--S_algorithm pyskani). Optional: dRep falls back + # to the skani/fastANI executables when it isn't installed. + 'pyskani': ['pyskani'], + }, zip_safe=False) diff --git a/tests/tests/test_pyskani.py b/tests/tests/test_pyskani.py new file mode 100644 index 0000000..f8702e4 --- /dev/null +++ b/tests/tests/test_pyskani.py @@ -0,0 +1,160 @@ +""" +Tests for the in-process pyskani backend (drep.d_cluster.pyskani_backend). + +pyskani is an optional dependency, so every test here skips cleanly when it +isn't installed. +""" +import glob +import os +import shutil +import tempfile + +import pandas as pd +import pytest + +import drep.d_cluster.compare_utils as cu +import drep.d_cluster.cluster_utils as clu +import drep.d_cluster.external as ext +import drep.d_cluster.greedy_clustering as gc +import drep.d_cluster.utils +import drep.d_filter + +def _has_pyskani(): + try: + import pyskani # noqa: F401 + return True + except ImportError: + return False + + +requires_pyskani = pytest.mark.skipif(not _has_pyskani(), reason="pyskani not installed") +requires_skani = pytest.mark.skipif(shutil.which('skani') is None, reason="skani not installed") +requires_fastani = pytest.mark.skipif(shutil.which('fastANI') is None, reason="fastANI not installed") + + +def _test_genomes(): + here = os.path.dirname(os.path.abspath(__file__)) + return sorted([g for g in glob.glob(os.path.join(here, '../genomes/*')) + if os.path.isfile(g)]) + + +def _partition(Cdb, col='secondary_cluster'): + return {frozenset(sub['genome']) for _, sub in Cdb.groupby(col)} + + +@requires_pyskani +def test_pyskani_ndb_shape(): + """run_pairwise_pyskani returns a complete, well-formed Ndb.""" + import drep.d_cluster.pyskani_backend as pb + genomes = _test_genomes() + Ndb = pb.run_pairwise_pyskani(genomes) + + assert list(Ndb.columns) == ['reference', 'querry', 'ani', 'alignment_coverage'] + # Every ordered pair present (dRep's hierarchical clustering needs a full matrix) + assert len(Ndb) == len(genomes) ** 2 + # ANI and coverage are on a 0-1 scale + assert Ndb['ani'].between(0, 1).all() + assert Ndb['alignment_coverage'].between(0, 1).all() + # Self comparisons are exactly 1 + selfs = Ndb[Ndb['reference'] == Ndb['querry']] + assert (selfs['ani'] == 1).all() + assert (selfs['alignment_coverage'] == 1).all() + + +@requires_pyskani +@requires_skani +def test_pyskani_agrees_with_subprocess_skani(): + """ + pyskani and the skani executable should report the same ANI/coverage, and + produce the same secondary clusters. + """ + import drep.d_cluster.pyskani_backend as pb + genomes = _test_genomes() + workdir = tempfile.mkdtemp() + try: + sub = ext.run_pairwise_skani(genomes, os.path.join(workdir, 'skani/'), processors=4) + py = pb.run_pairwise_pyskani(genomes) + + m = pd.merge(sub, py, on=['reference', 'querry'], suffixes=('_sub', '_py')) + # Only close pairs matter for clustering; skani's min-af filter drops + # distant pairs from pyskani's output (they're filled in as ani=0). + close = m[m['ani_sub'] >= 0.95] + assert len(close) > 0 + assert (close['ani_sub'] - close['ani_py']).abs().max() < 0.001 + assert (close['alignment_coverage_sub'] - close['alignment_coverage_py']).abs().max() < 0.001 + + for sa, nc in [(0.99, 0.1), (0.95, 0.1)]: + c1, _ = clu.genome_hierarchical_clustering(sub, S_ani=sa, cov_thresh=nc, + comp_method='skani', cluster='X') + c2, _ = clu.genome_hierarchical_clustering(py, S_ani=sa, cov_thresh=nc, + comp_method='pyskani', cluster='X') + assert _partition(c1) == _partition(c2), f"clusters differ at S_ani={sa}" + finally: + shutil.rmtree(workdir, ignore_errors=True) + + +@requires_pyskani +@requires_fastani +def test_greedy_pyskani_matches_greedy_fastani(): + """Greedy clustering should give the same answer via pyskani as via fastANI.""" + genomes = _test_genomes() + bdb = drep.d_cluster.utils.load_genomes(genomes) + bdb = drep.d_filter._add_lengthN50(bdb, bdb) + + workdir = tempfile.mkdtemp() + try: + parts = {} + for alg in ['fastANI', 'pyskani']: + d = os.path.join(workdir, alg + '/') + os.makedirs(d, exist_ok=True) + Ndb, Cdb, _ = gc.compare_genomes_greedy( + bdb, alg, d, S_ani=0.95, cov_thresh=0.1, cluster='P1', processors=4) + parts[alg] = _partition(Cdb) + assert parts['fastANI'] == parts['pyskani'] + finally: + shutil.rmtree(workdir, ignore_errors=True) + + +@requires_pyskani +def test_greedy_pyskani_sketches_each_genome_once(): + """ + The whole point of the pyskani greedy path: every genome is sketched exactly + once, no matter how many representatives accumulate. + """ + import drep.d_cluster.pyskani_backend as pb + genomes = _test_genomes() + bdb = drep.d_cluster.utils.load_genomes(genomes) + bdb = drep.d_filter._add_lengthN50(bdb, bdb) + + sketch_calls = [] + orig = pb.PyskaniDatabase.add + + def counting_add(self, name, contigs): + sketch_calls.append(name) + return orig(self, name, contigs) + + pb.PyskaniDatabase.add = counting_add + workdir = tempfile.mkdtemp() + try: + gc.compare_genomes_greedy(bdb, 'pyskani', os.path.join(workdir, 'g/'), + S_ani=0.95, cov_thresh=0.1, cluster='P1') + # One sketch per representative, and never the same genome twice + assert len(sketch_calls) == len(set(sketch_calls)), "a genome was sketched more than once" + assert len(sketch_calls) <= len(genomes) + finally: + pb.PyskaniDatabase.add = orig + shutil.rmtree(workdir, ignore_errors=True) + + +@requires_pyskani +def test_compare_genomes_dispatches_pyskani(): + """--S_algorithm pyskani is reachable through the normal dispatch path.""" + genomes = _test_genomes() + bdb = drep.d_cluster.utils.load_genomes(genomes) + workdir = tempfile.mkdtemp() + try: + Ndb = cu.compare_genomes(bdb, 'pyskani', workdir) + assert len(Ndb) == len(genomes) ** 2 + assert 'ani' in Ndb.columns + finally: + shutil.rmtree(workdir, ignore_errors=True) From 6e563494d8c9fb1f5a0a2afd02156429aec1903d Mon Sep 17 00:00:00 2001 From: Matt Olm Date: Tue, 14 Jul 2026 09:29:14 -0600 Subject: [PATCH 04/15] Pass --min-af 0 to sparse skani primary clustering skani's triangle defaults to dropping pairs that align over less than 15% of the genome. Primary clustering is a deliberately loose, inclusive pre-filter -- the MASH path applies no alignment-fraction filter at all -- so inheriting skani's default silently strands related genomes in separate primary clusters, where the secondary algorithm never compares them. This matters most for exactly the data dRep is used on. Validating against 2,631 real ocean MAGs (TOBG), skani's default min-af dropped 16 above-threshold edges and split 15 primary clusters that --min-af 0 keeps together: fragmented, partial MAGs of the same organism routinely align over well under 15% of their length. The pairwise skani path already passes --min-af 0 for the same reason. Co-Authored-By: Claude Opus 4.8 --- drep/d_cluster/external.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/drep/d_cluster/external.py b/drep/d_cluster/external.py index 4d6a25d..1786e12 100644 --- a/drep/d_cluster/external.py +++ b/drep/d_cluster/external.py @@ -316,8 +316,15 @@ def run_skani_triangle_sparse(genome_list, outdir, screen, **kwargs): exe_loc = drep.get_exe('skani') out_file = os.path.join(outdir, 'skani_sparse_{0}.tsv'.format(code)) + # --min-af 0 is essential here. skani defaults to dropping pairs that align + # over <15% of the genome, but primary clustering is a deliberately loose, + # inclusive pre-filter -- the MASH path applies no alignment-fraction filter + # at all. Fragmented/partial MAGs of the same organism routinely align over + # less than 15%, and dropping those pairs would strand related genomes in + # separate primary clusters, where they are never compared by the secondary + # algorithm. (The pairwise skani path passes --min-af 0 for the same reason.) cmd = [exe_loc, "triangle", "--sparse", "-t", str(p), '-o', out_file, - '-l', glist, '-s', str(screen)] + '-l', glist, '-s', str(screen), '--min-af', '0'] if extra_cmd != "": cmd += extra_cmd.split(' ') From 30ef829726f4eb133f11dbb07c77fe0d103c5c2b Mon Sep 17 00:00:00 2001 From: Matt Olm Date: Tue, 14 Jul 2026 11:07:29 -0600 Subject: [PATCH 05/15] Keep skani's aligned-fraction filter for primary clustering Reverts the --min-af 0 from 6e56349, which was wrong, and makes the threshold configurable via --primary_skani_min_af (default 15, skani's own default). 6e56349 reasoned that primary clustering is a loose pre-filter and MASH applies no alignment-fraction filter, so skani shouldn't either. That reasoning missed that MASH similarity and skani ANI are not the same measurement. MASH compares k-mers across the whole genome, so two genomes sharing only a small conserved region score as distant. skani reports the identity *within aligned regions only*, so that same pair looks like a high-ANI edge. skani's --min-af is what makes its ANI comparable to MASH's whole-genome similarity; it is not an obstacle to work around. Under single linkage the consequence is severe, because a handful of spurious bridges merge everything they touch. Measured on 10,000 real UHGG genomes: min-af edges clusters largest cluster 0 354,690 734 5,857 <- 59% of the dataset in one cluster 10 340,387 944 688 15 339,004 984 626 MASH 100M 989 626 Dropping min-af from 15 to 0 adds only 4.4% more edges but collapses 59% of the dataset into a single primary cluster, which would make secondary clustering intractable. At the default, skani reproduces the MASH partition closely (984 vs 989 clusters, identical largest cluster; each splits ~1-2% of the other's multi-genome clusters, consistent with MASH overestimating ANI at the threshold). The earlier evidence for --min-af 0 -- 15 clusters on 2,631 ocean MAGs that min-af 0 kept together -- is better explained as those same low-overlap bridges being correctly rejected. Also pins the pyskani extra to >=0.2 (needed for the `cutoff` query argument; note there is no macOS arm64 wheel for 0.2, so Apple Silicon builds from source and needs a Rust toolchain). Co-Authored-By: Claude Opus 4.8 --- drep/argumentParser.py | 9 ++++++++ drep/d_cluster/compare_utils.py | 11 +++++++-- drep/d_cluster/external.py | 25 +++++++++++++-------- setup.py | 5 ++++- tests/tests/test_union_find.py | 40 +++++++++++++++++++++++++++++++++ 5 files changed, 78 insertions(+), 12 deletions(-) diff --git a/drep/argumentParser.py b/drep/argumentParser.py index fdf4a8a..acebd5d 100644 --- a/drep/argumentParser.py +++ b/drep/argumentParser.py @@ -131,6 +131,15 @@ def parse_args(args): + " produced and they are streamed rather than held as a dense\n" \ + " matrix. Recommended for very large genome sets (no N^2 RAM/disk).", default='MASH', choices={'MASH', 'skani'}) + 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') diff --git a/drep/d_cluster/compare_utils.py b/drep/d_cluster/compare_utils.py index b56f968..f34cefa 100644 --- a/drep/d_cluster/compare_utils.py +++ b/drep/d_cluster/compare_utils.py @@ -175,13 +175,20 @@ def primary_cluster_skani_sparse(Bdb, data_folder, **kwargs): default_screen = max(1.0, min(ani_threshold - 5.0, 99.0)) screen = kwargs.get('primary_skani_screen', default_screen) + # Minimum percent of a genome that must align for a pair to count as an edge. + # 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) + 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"(ANI threshold {ani_threshold:.1f}%, screen {screen:.1f}%, min-af {min_af}%)") sparse_file = drep.d_cluster.external.run_skani_triangle_sparse( - genome_list, skani_folder, screen, **kwargs) + genome_list, skani_folder, screen, min_af=min_af, **kwargs) all_genomes = list(Bdb['genome'].unique()) Cdb, Mdb, stats = drep.d_cluster.union_find.cluster_skani_sparse_files( diff --git a/drep/d_cluster/external.py b/drep/d_cluster/external.py index 1786e12..036e8af 100644 --- a/drep/d_cluster/external.py +++ b/drep/d_cluster/external.py @@ -278,7 +278,7 @@ def _fix_fastani(odb): return fdb -def run_skani_triangle_sparse(genome_list, outdir, screen, **kwargs): +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. @@ -292,6 +292,9 @@ def run_skani_triangle_sparse(genome_list, outdir, screen, **kwargs): 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). @@ -316,15 +319,19 @@ def run_skani_triangle_sparse(genome_list, outdir, screen, **kwargs): exe_loc = drep.get_exe('skani') out_file = os.path.join(outdir, 'skani_sparse_{0}.tsv'.format(code)) - # --min-af 0 is essential here. skani defaults to dropping pairs that align - # over <15% of the genome, but primary clustering is a deliberately loose, - # inclusive pre-filter -- the MASH path applies no alignment-fraction filter - # at all. Fragmented/partial MAGs of the same organism routinely align over - # less than 15%, and dropping those pairs would strand related genomes in - # separate primary clusters, where they are never compared by the secondary - # algorithm. (The pairwise skani path passes --min-af 0 for the same reason.) + # 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', '0'] + '-l', glist, '-s', str(screen), '--min-af', str(min_af)] if extra_cmd != "": cmd += extra_cmd.split(' ') diff --git a/setup.py b/setup.py index 49831ee..dd25acd 100644 --- a/setup.py +++ b/setup.py @@ -31,6 +31,9 @@ def version(): extras_require={ # In-process skani (--S_algorithm pyskani). Optional: dRep falls back # to the skani/fastANI executables when it isn't installed. - 'pyskani': ['pyskani'], + # >=0.2 is required for the `cutoff` query argument. Note there is no + # macOS arm64 wheel for 0.2 yet, so Apple Silicon builds it from source + # and needs a Rust toolchain (conda install -c conda-forge rust). + 'pyskani': ['pyskani>=0.2'], }, zip_safe=False) diff --git a/tests/tests/test_union_find.py b/tests/tests/test_union_find.py index ae4f9a2..4605b24 100644 --- a/tests/tests/test_union_find.py +++ b/tests/tests/test_union_find.py @@ -104,6 +104,46 @@ def membership(Cdb): 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 + + @pytest.mark.skipif(shutil.which('skani') is None, reason="skani not installed") def test_sparse_skani_primary_matches_mash(): """ From 5ed4d6a1580c09c79e0cc3ead2fe53e6197826eb Mon Sep 17 00:00:00 2001 From: Matt Olm Date: Tue, 14 Jul 2026 13:48:22 -0600 Subject: [PATCH 06/15] Don't let estimate_time kill a run on an unknown S_algorithm estimate_time assigned `time` inside an if/elif chain with no else, so any algorithm it didn't recognize raised UnboundLocalError instead of returning an estimate. Adding pyskani as an --S_algorithm choice without updating this function meant every pyskani run died the moment secondary clustering started: File "drep/d_cluster/utils.py", line 108, in estimate_time return time UnboundLocalError: local variable 'time' referenced before assignment Caught by a 10,000-genome end-to-end run, which got through primary clustering (984 clusters, ~4 min) and then fell over here before doing any real work. Replace the chain with a lookup that falls back to the fast-algorithm estimate for anything unrecognized. This only drives a log line -- it should never be able to take a run down. Add a test covering every --S_algorithm choice plus an unknown one. Co-Authored-By: Claude Opus 4.8 --- drep/d_cluster/utils.py | 29 ++++++++++++++++------------- tests/tests/test_pyskani.py | 17 +++++++++++++++++ 2 files changed, 33 insertions(+), 13 deletions(-) diff --git a/drep/d_cluster/utils.py b/drep/d_cluster/utils.py index 91d451f..771182b 100644 --- a/drep/d_cluster/utils.py +++ b/drep/d_cluster/utils.py @@ -93,19 +93,22 @@ 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, + # in-process; sketches each genome once, so at least as fast as skani + 'pyskani': 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/tests/tests/test_pyskani.py b/tests/tests/test_pyskani.py index f8702e4..1d2d311 100644 --- a/tests/tests/test_pyskani.py +++ b/tests/tests/test_pyskani.py @@ -146,6 +146,23 @@ def counting_add(self, name, contigs): shutil.rmtree(workdir, ignore_errors=True) +def test_estimate_time_handles_every_S_algorithm(): + """ + estimate_time only drives a log line, but it used to raise UnboundLocalError + for any algorithm it didn't know about -- which killed a 10k-genome run at + the start of secondary clustering. Every --S_algorithm choice must work, and + unknown ones must not raise. + """ + from drep.d_cluster.utils import estimate_time + + for alg in ['ANIn', 'gANI', 'goANI', 'ANImf', 'fastANI', 'skani', 'pyskani']: + t = estimate_time(100, alg) + assert t > 0, f"{alg} gave {t!r}" + + # An unrecognized algorithm must degrade gracefully, not raise + assert estimate_time(100, 'some_future_algorithm') > 0 + + @requires_pyskani def test_compare_genomes_dispatches_pyskani(): """--S_algorithm pyskani is reachable through the normal dispatch path.""" From 71e2328ef6b5f513beb02621313e615329766bdc Mon Sep 17 00:00:00 2001 From: Matt Olm Date: Tue, 14 Jul 2026 14:46:20 -0600 Subject: [PATCH 07/15] One-pass skani: reuse primary's comparisons for secondary clustering dRep sketched every genome twice and computed the same ANI values twice. Primary ran one skani pass over all genomes; secondary then ran skani again, once per primary cluster. But secondary only ever compares genomes *within* a primary cluster, and those pairs are a subset of what the primary pass already computed. On 10,000 UHGG genomes, 94% of the pairs driving secondary clustering were already present in primary's sparse output, with ANI identical to 6 decimal places -- because it is literally the same skani computation. The remaining 6% were entirely pairs below primary's --min-af, with a median aligned fraction of 1.1%, which cov_thresh discards anyway. So run skani once and derive both stages from the same edge table: - the single pass emits at min-af = min(primary_skani_min_af, cov_thresh*100), since secondary applies the looser coverage filter; the stricter primary filter is applied when forming primary clusters - the screen also stays at or below the secondary threshold, as those edges are now reused rather than recomputed - Mdb keeps every edge (not just those above P_ani) with alignment coverage, so it can feed secondary. It is also a far more useful Mdb: 779k real ANI values instead of 100M Mash distances Measured on 10,000 UHGG genomes, dereplicate end to end: two-pass one-pass secondary 15 min 24 s total 22.4 min 13.8 min Output is identical, not merely equivalent: same 984 primary clusters, same 1,232 secondary clusters, same 1,232 representative genomes. 984 subprocess spawns and an entire re-sketch of all 10,000 genomes are gone. Reuse applies when --primary_algorithm skani is paired with a skani --S_algorithm and greedy is off; --no_reuse_primary_comparisons forces the old behavior. Other secondary algorithms measure something different and still run for themselves. Co-Authored-By: Claude Opus 4.8 --- drep/argumentParser.py | 7 ++ drep/d_cluster/compare_utils.py | 75 ++++++++++++---- drep/d_cluster/controller.py | 31 +++++++ drep/d_cluster/union_find.py | 153 ++++++++++++++++++++++++++++++++ tests/tests/test_union_find.py | 74 +++++++++++++++ 5 files changed, 324 insertions(+), 16 deletions(-) diff --git a/drep/argumentParser.py b/drep/argumentParser.py index acebd5d..704056e 100644 --- a/drep/argumentParser.py +++ b/drep/argumentParser.py @@ -131,6 +131,13 @@ def parse_args(args): + " produced and they are streamed rather than held as a dense\n" \ + " matrix. Recommended for very large genome sets (no N^2 RAM/disk).", default='MASH', 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 " diff --git a/drep/d_cluster/compare_utils.py b/drep/d_cluster/compare_utils.py index f34cefa..216bbe9 100644 --- a/drep/d_cluster/compare_utils.py +++ b/drep/d_cluster/compare_utils.py @@ -159,50 +159,93 @@ def all_vs_all_primary(Bdb, data_folder, **kwargs): def primary_cluster_skani_sparse(Bdb, data_folder, **kwargs): """ - Primary clustering via `skani triangle --sparse` streamed into union-find. + Primary clustering via one `skani triangle --sparse` pass + union-find. - Only above-threshold pairs are ever produced (skani screens during sketching) - and they are streamed rather than held as a dense matrix, so memory stays - O(genomes + edges) regardless of genome count. Always single-linkage - (connected components); --classic_primary_clustering / non-single - primary_clusterAlg do not apply here. + 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. - default_screen = max(1.0, min(ani_threshold - 5.0, 99.0)) + # 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 count as an edge. + # 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}%, min-af {min_af}%)") + 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=min_af, **kwargs) + genome_list, skani_folder, screen, min_af=emit_min_af, **kwargs) all_genomes = list(Bdb['genome'].unique()) - Cdb, Mdb, stats = drep.d_cluster.union_find.cluster_skani_sparse_files( - sparse_file, ani_threshold, all_genomes, - progress=kwargs.get('primary_progress', True)) + 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)") - logging.info(f" Sparse skani primary clustering: {stats['edges_kept']:,} edges kept, " - f"{stats['primary_clusters']:,} primary clusters") + # 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'] arguments = {'linkage_method': 'single', 'linkage_cutoff': 1 - P_ani, 'comparison_algorithm': 'skani'} cluster_ret = ['union_find_streaming', None, 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 diff --git a/drep/d_cluster/controller.py b/drep/d_cluster/controller.py index 665624d..6b291f6 100644 --- a/drep/d_cluster/controller.py +++ b/drep/d_cluster/controller.py @@ -148,6 +148,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 +181,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', 'MASH') != 'skani': + return False + if algorithm not in ('skani', 'pyskani'): + 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/union_find.py b/drep/d_cluster/union_find.py index d72faa3..98564f0 100644 --- a/drep/d_cluster/union_find.py +++ b/drep/d_cluster/union_find.py @@ -239,6 +239,159 @@ def to_name(x): return Cdb, stats +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, progress=False): + """ + 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): """ diff --git a/tests/tests/test_union_find.py b/tests/tests/test_union_find.py index 4605b24..3f61a55 100644 --- a/tests/tests/test_union_find.py +++ b/tests/tests/test_union_find.py @@ -144,6 +144,80 @@ def test_skani_sparse_min_af_filters_low_alignment_edges(): 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(): """ From 783821da291f4b568930111c328fc13ac55767c3 Mon Sep 17 00:00:00 2001 From: Matt Olm Date: Tue, 14 Jul 2026 16:51:21 -0600 Subject: [PATCH 08/15] Make skani the default; remove pyskani and --low_ram_primary_clustering Everything measured on the 10,000-genome UHGG set said skani should be the default, but it was still opt-in: a default run got MASH primary (a 100M-row Mdb, 5.8 GB of RAM, 15 GB of disk) and fastANI secondary. Now `dRep compare` with no flags runs the one-pass skani path -- 22.4 min -> 13.8 min end to end at 10k, with 187x less intermediate disk. --primary_algorithm MASH -> skani --S_algorithm fastANI -> skani MASH and fastANI remain available. Two consequences of the flip, handled here: - the mash executable is no longer required unless --primary_algorithm MASH is requested; the startup dependency check now looks for whichever program the run will actually use - --multiround_primary_clustering and --primary_chunksize are MASH-path concepts. Rather than silently ignoring them under skani, warn: skani's sparse output never builds the N^2 table multiround exists to avoid, and it has none of multiround's chunk-splitting imprecision Remove pyskani. The one-pass work obsoleted it: secondary clustering no longer runs comparisons at all for skani workflows, and greedy's reason to exist (dodging O(n^2) within a cluster) is subsumed by sparse output that never produces n^2 in the first place. It was also the only dependency needing a Rust toolchain on Apple Silicon, and was 12x *slower* than skani for pairwise comparisons because it is single-threaded -- a trap for anyone choosing it expecting speed. Remove --low_ram_primary_clustering, deprecated since union-find became the default for single-linkage primary clustering. This also retires the networkx connected-components path it was the only caller of, so networkx is no longer a dependency. The primary dendrogram needed the same small-N treatment on the skani path that the MASH path already had, plus a fix: leaf labels were derived from Mdb, but the sparse Mdb legitimately omits genomes with no above-threshold pairs, so labels and linkage disagreed. Labels now come from the matrix the linkage was built on. Docs: regenerate the embedded CLI help, rewrite the algorithm overview for the skani defaults and one-pass reuse, document why skani needs an aligned-fraction filter where Mash does not, and correct the dependency lists. The algorithm section had also been stale independently -- it still called ANImf the default. Co-Authored-By: Claude Opus 4.8 --- README.md | 9 +- docs/choosing_parameters.rst | 19 ++- docs/installation.rst | 9 +- docs/module_descriptions.rst | 230 +++++++++++++++++++++------- docs/overview.rst | 8 +- drep/argumentParser.py | 32 ++-- drep/d_analyze.py | 22 ++- drep/d_cluster/cluster_utils.py | 81 +++------- drep/d_cluster/compare_utils.py | 57 ++++--- drep/d_cluster/controller.py | 13 +- drep/d_cluster/greedy_clustering.py | 27 ++-- drep/d_cluster/pyskani_backend.py | 215 -------------------------- drep/d_cluster/union_find.py | 38 ++++- drep/d_cluster/utils.py | 2 - setup.py | 9 -- tests/tests/test_cluster.py | 24 +-- tests/tests/test_dereplicate.py | 10 +- tests/tests/test_filter.py | 2 +- tests/tests/test_greedy.py | 23 +-- tests/tests/test_pyskani.py | 177 --------------------- tests/tests/test_union_find.py | 11 +- 21 files changed, 382 insertions(+), 636 deletions(-) delete mode 100644 drep/d_cluster/pyskani_backend.py delete mode 100644 tests/tests/test_pyskani.py diff --git a/README.md b/README.md index 42f3839..18fbda0 100644 --- a/README.md +++ b/README.md @@ -38,13 +38,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..26a29bc 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 v3.7.1 :::... - 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,goANI,fastANI,gANI,skani,ANIn}] + [--primary_algorithm {MASH,skani}] + [--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 {average,weighted,single,median,centroid,complete,ward}] + [--primary_clusterAlg {average,weighted,single,median,centroid,complete,ward}] + [--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,goANI,fastANI,gANI,skani,ANIn} 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 {MASH,skani} + 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 {average,weighted,single,median,centroid,complete,ward} + Algorithm used to cluster genomes during SECONDARY + clustering (passed to scipy.cluster.hierarchy.linkage) + (default: average) + --primary_clusterAlg {average,weighted,single,median,centroid,complete,ward} + 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}] [--set_recursion SET_RECURSION] - [--S_algorithm {goANI,ANIn,gANI,ANImf,fastANI}] + [--checkm_group_size CHECKM_GROUP_SIZE] + [--S_algorithm {ANIn,skani,gANI,fastANI,ANImf,goANI}] + [--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 {weighted,ward,complete,single,average,median,centroid}] + [--primary_clusterAlg {weighted,ward,complete,single,average,median,centroid}] + [--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,11 +278,12 @@ 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) + 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 {taxonomy_wf,lineage_wf} Either lineage_wf (more accurate) or taxonomy_wf (faster) (default: lineage_wf) @@ -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 {ANIn,skani,gANI,fastANI,ANImf,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 @@ -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 {weighted,ward,complete,single,average,median,centroid} + Algorithm used to cluster genomes during SECONDARY + clustering (passed to scipy.cluster.hierarchy.linkage) + (default: average) + --primary_clusterAlg {weighted,ward,complete,single,average,median,centroid} + 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/argumentParser.py b/drep/argumentParser.py index 704056e..7d7d4d7 100644 --- a/drep/argumentParser.py +++ b/drep/argumentParser.py @@ -114,23 +114,22 @@ 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" \ - + "pyskani = skani run in-process via the pyskani library. Each genome is\n" \ - + " sketched exactly once instead of being re-sketched by a new\n" \ - + " subprocess for every comparison, which is much faster for\n" \ - + " greedy clustering. Requires `pip install pyskani`.\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', 'pyskani'}) + default='skani', choices={'ANIn', 'gANI', 'ANImf', 'goANI', 'fastANI', 'skani'}) Clustflags.add_argument("--primary_algorithm", help="R|Program to use for primary clustering.\n" \ - + "MASH = (DEFAULT) all-vs-all Mash\n" \ - + "skani = skani triangle --sparse; only above-threshold pairs are\n" \ - + " produced and they are streamed rather than held as a dense\n" \ - + " matrix. Recommended for very large genome sets (no N^2 RAM/disk).", - default='MASH', choices={'MASH', 'skani'}) + + "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 " @@ -148,8 +147,10 @@ def parse_args(args): "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", @@ -186,9 +187,6 @@ def parse_args(args): RAM at scale but reproduces pre-v4 behavior and allows non-single linkage methods\ and the primary dendrogram plot.", action='store_true', default=False) - Compflags.add_argument("--low_ram_primary_clustering", help="(Deprecated; the streaming single-linkage\ - algorithm is now the default for primary clustering.) Kept for backwards compatibility.", - action='store_true', default=False) GRflags = cluster_parent.add_argument_group('GREEDY CLUSTERING OPTIONS\n' 'These decrease RAM use and runtime at the expense of a minor loss in ' diff --git a/drep/d_analyze.py b/drep/d_analyze.py index d58dadc..454b539 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: @@ -151,13 +152,18 @@ def mash_dendrogram_from_wd(wd, plot_dir=False): return if Plinkage is None or isinstance(Plinkage, str): - logging.error("Skipping plot 1 - cannot generate with low_ram_primary_clustering (no linkage matrix)") + 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): ''' @@ -618,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 @@ -628,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 @@ -637,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) 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 216bbe9..7954faf 100644 --- a/drep/d_cluster/compare_utils.py +++ b/drep/d_cluster/compare_utils.py @@ -14,7 +14,6 @@ import drep.d_cluster.utils import drep.d_cluster.greedy_clustering import drep.d_cluster.union_find -import drep.d_cluster.pyskani_backend class genomeChunk(): """ @@ -144,15 +143,23 @@ def all_vs_all_primary(Bdb, data_folder, **kwargs): """ Dispatch primary clustering to the requested algorithm. - 'MASH' (default) uses the classic all-vs-all Mash path. 'skani' uses - `skani triangle --sparse` + streaming union-find, which never builds the N^2 - matrix on disk or in RAM and is the recommended path for very large genome - sets. + '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', '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) @@ -217,9 +224,26 @@ def primary_cluster_skani_sparse(Bdb, data_folder, **kwargs): 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 = ['union_find_streaming', None, arguments] + cluster_ret = [linkage, linkage_db, arguments] return Mdb, Cdb, cluster_ret @@ -418,7 +442,6 @@ def cluster_mash_database(db, **kwargs): clusterAlg: legacy fallback for primary_clusterAlg (default = single) P_ani: threshold to cluster at (default = 0.9) classic_primary_clustering: force the dense scipy path - low_ram_primary_clustering: deprecated alias forcing single-linkage union-find Returns: list: [Cdb, [linkage, linkage_db, arguments]] @@ -431,7 +454,6 @@ def cluster_mash_database(db, **kwargs): P_Lmethod = kwargs.get('primary_clusterAlg') or kwargs.get('clusterAlg', 'single') P_Lcutoff = 1 - kwargs.get('P_ani',.9) classic = kwargs.get('classic_primary_clustering', False) - low_ram = kwargs.get('low_ram_primary_clustering', False) db['dist'] = 1 - db['similarity'] @@ -439,12 +461,8 @@ def cluster_mash_database(db, **kwargs): # 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') or low_ram) + use_union_find = (not classic) and (P_Lmethod == 'single') if use_union_find: - if low_ram and P_Lmethod != 'single': - logging.warning( - f"low_ram_primary_clustering uses single-linkage (connected components); " - f"ignoring primary_clusterAlg={P_Lmethod} for primary clustering.") Cdb = drep.d_cluster.union_find.cluster_long_df(db, P_Lcutoff) arguments = {'linkage_method': 'single', 'linkage_cutoff': P_Lcutoff, @@ -459,7 +477,7 @@ def cluster_mash_database(db, **kwargs): linkage_db = None dendro_max = kwargs.get('primary_dendrogram_max_genomes', 2000) n_genomes = Cdb['genome'].nunique() - if (not low_ram) and n_genomes <= dendro_max and 'genome_chunk' not in db.columns: + 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)) @@ -475,7 +493,7 @@ def cluster_mash_database(db, **kwargs): # 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=False) + linkage_cutoff= P_Lcutoff) Cdb = Cdb.rename(columns={'cluster':'primary_cluster'}) Cdb['primary_cluster'] = Cdb['primary_cluster'].astype(int) @@ -572,11 +590,6 @@ def compare_genomes(bdb, algorithm, data_folder, **kwargs): df = drep.d_cluster.external.run_pairwise_skani(genome_list, working_data_folder, **kwargs) return df - elif algorithm == 'pyskani': - genome_list = bdb['location'].tolist() - df = drep.d_cluster.pyskani_backend.run_pairwise_pyskani(genome_list, **kwargs) - return df - elif algorithm == 'gANI': # Figure out prodigal folder wd = kwargs.get('wd', False) @@ -610,7 +623,7 @@ def compare_genomes(bdb, algorithm, data_folder, **kwargs): sys.exit() else: - SUPPORTED = ['fastANI', 'pyskani'] + SUPPORTED = ['fastANI'] if algorithm not in SUPPORTED: message = f"{algorithm} is not supported for greedy secondary clustering!\nChoose one of the following supported S_algorithm options: {' '.join(SUPPORTED)}" logging.error(message) diff --git a/drep/d_cluster/controller.py b/drep/d_cluster/controller.py index 6b291f6..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: @@ -193,9 +194,9 @@ def can_reuse_primary_edges(self, algorithm): """ if self.kwargs.get('reuse_primary_comparisons', True) is False: return False - if self.kwargs.get('primary_algorithm', 'MASH') != 'skani': + if self.kwargs.get('primary_algorithm', 'skani') != 'skani': return False - if algorithm not in ('skani', 'pyskani'): + if algorithm != 'skani': return False if self.kwargs.get('greedy_secondary_clustering', False): return False diff --git a/drep/d_cluster/greedy_clustering.py b/drep/d_cluster/greedy_clustering.py index 3a75a5b..84399ef 100644 --- a/drep/d_cluster/greedy_clustering.py +++ b/drep/d_cluster/greedy_clustering.py @@ -7,7 +7,6 @@ import drep.d_cluster.external import drep.d_cluster.compare_utils -import drep.d_cluster.pyskani_backend def greedy_secondary_clustering(Bdb, Cdb, algorithm, data_folder, **kwargs): ndbs = [] @@ -108,13 +107,12 @@ 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 (O(N*R) sketching). The pyskani path below - # sketches each genome exactly once instead. + # 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) - elif algorithm == 'pyskani': - return drep.d_cluster.pyskani_backend.pyskani_one_vs_many( - new_genome, kwargs['pyskani_db'], **kwargs) else: logging.error("{0} algorithm is not yet supported for greedy clustering; sorry!".format(algorithm)) assert False @@ -124,13 +122,11 @@ def add_genome_as_rep(location, algorithm, **kwargs): """ Register a genome as a new cluster representative. - For pyskani this sketches it once into the in-memory database, so subsequent - genomes can be compared against it without any re-sketching. Subprocess-based - algorithms read the representative list from a file instead and need nothing - here. + 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. """ - if algorithm == 'pyskani': - kwargs['pyskani_db'].add_genome(location) + return def prepare_for_greedy(algorithm, data_folder, **kwargs): @@ -155,11 +151,6 @@ def prepare_for_greedy(algorithm, data_folder, **kwargs): kwargs['logdir'] = logdir kwargs['current_exe'] = drep.get_exe('fastANI') - elif algorithm == 'pyskani': - # One in-memory database of representatives for this primary cluster. - # Each representative is sketched exactly once, on the way in. - kwargs['pyskani_db'] = drep.d_cluster.pyskani_backend.PyskaniDatabase(**kwargs) - return kwargs def order_genomes_for_greedy(bdb, **kwargs): diff --git a/drep/d_cluster/pyskani_backend.py b/drep/d_cluster/pyskani_backend.py deleted file mode 100644 index 37f0fef..0000000 --- a/drep/d_cluster/pyskani_backend.py +++ /dev/null @@ -1,215 +0,0 @@ -""" -In-process skani comparisons via pyskani (https://github.com/althonos/pyskani). - -The subprocess-based comparison algorithms re-sketch genomes on every -invocation. That is especially wasteful during greedy secondary clustering, -where each new genome is compared against the growing set of cluster -representatives by spawning a fresh subprocess that re-sketches *every* -representative -- O(N * R) sketching work and N subprocess spawns for N genomes. - -pyskani lets us sketch each genome exactly once, keep the representatives in an -in-memory database, and query it directly. Sketching becomes O(N), there are no -subprocesses, and no temporary files are written. - -This module is imported lazily: pyskani is an optional dependency, and dRep -falls back to the subprocess skani/fastANI implementations without it. -""" - -import functools -import logging -import os - -import pandas as pd - -import drep.d_cluster.utils - -# Column layout every secondary-clustering comparison must return -NDB_COLUMNS = ['reference', 'querry', 'ani', 'alignment_coverage'] - - -def import_pyskani(): - """ - Import pyskani, raising an actionable error if it isn't installed. - """ - try: - import pyskani - except ImportError: - raise ImportError( - "The 'pyskani' S_algorithm requires the pyskani package, which is not " - "installed. Install it with `pip install pyskani` (or `pip install " - "drep[pyskani]`), or choose a different --S_algorithm (e.g. skani, " - "which uses the skani executable instead)." - ) - return pyskani - - -def load_contigs(location): - """ - Read a FASTA file into a list of contig sequences (as bytes), which is what - pyskani's sketch/query expect. - """ - from Bio import SeqIO - return [bytes(record.seq) for record in SeqIO.parse(location, 'fasta')] - - -# During greedy clustering a genome is queried and then, if it founds a new -# cluster, immediately sketched as a representative. A tiny cache avoids parsing -# the same FASTA twice in a row. Deliberately kept at maxsize=2 -- caching every -# genome's contigs would hold the entire input set in memory. -@functools.lru_cache(maxsize=2) -def _load_contigs_cached(location): - return load_contigs(location) - - -class PyskaniDatabase: - """ - A pyskani database that sketches each genome exactly once. - - Sketching is the expensive part of an ANI comparison, so genomes are sketched - on the way in and the resulting database is queried many times. Contigs are - cached only while needed to sketch/query, not retained for the lifetime of - the object. - """ - - def __init__(self, **kwargs): - pyskani = import_pyskani() - self.db = pyskani.Database() - self.names = [] - # Screening cutoff. Mirrors skani's -s flag; the subprocess pairwise path - # uses -s 1 (i.e. compare essentially everything), so default low here and - # let callers raise it when they only care about close relatives. - self.cutoff = kwargs.get('pyskani_cutoff', 0.01) - self.learned_ani = kwargs.get('pyskani_learned_ani', None) - - def add(self, name, contigs): - """Sketch a genome once and store it as a reference.""" - self.db.sketch(name, *contigs) - self.names.append(name) - - def add_genome(self, location): - """Sketch a genome from a FASTA path; returns its dRep genome name.""" - name = drep.d_cluster.utils._get_genome_name_from_fasta(location) - self.add(name, _load_contigs_cached(location)) - return name - - def query_hits(self, name, contigs): - """ - Query the database with a genome, returning the raw pyskani hits. - """ - kw = {'cutoff': self.cutoff} - if self.learned_ani is not None: - kw['learned_ani'] = self.learned_ani - return self.db.query(name, *contigs, **kw) - - def query(self, name, contigs, query_as_reference=False): - """ - Query the database, returning (reference, querry, ani, alignment_coverage) - tuples. - - Everywhere in dRep, alignment_coverage is the aligned fraction of the - genome named in the 'reference' column (skani's .af matrix cell [A][B] is - the aligned fraction of A; fastANI's matched/total is the fraction of the - genome that load_fastani puts in 'reference'). Both orientations below - respect that rule. - - Args: - query_as_reference: if False (default, pairwise use), emit - reference=the database genome and coverage=its aligned fraction. - If True (greedy use), emit reference=the queried genome and - coverage=the queried genome's aligned fraction -- which is what - get_cluster_rep expects, since it reads the representative out of - the 'querry' column. - """ - rows = [] - for hit in self.query_hits(name, contigs): - if query_as_reference: - rows.append((hit.query_name, hit.reference_name, - hit.identity, hit.query_fraction)) - else: - rows.append((hit.reference_name, hit.query_name, - hit.identity, hit.reference_fraction)) - return rows - - def query_genome(self, location, query_as_reference=False): - name = drep.d_cluster.utils._get_genome_name_from_fasta(location) - return self.query(name, _load_contigs_cached(location), - query_as_reference=query_as_reference) - - -def _fill_missing_pairs(rows, names): - """ - skani only reports pairs that share enough k-mers. dRep's hierarchical - secondary clustering needs a complete matrix, so absent pairs are filled in - as ani=0 / coverage=0 (no detectable relatedness), matching what the - subprocess `skani triangle --min-af 0` path yields for unrelated genomes. - """ - have = {(r[0], r[1]) for r in rows} - filled = list(rows) - for a in names: - for b in names: - if (a, b) not in have: - filled.append((a, b, 0.0, 0.0)) - return filled - - -def run_pairwise_pyskani(genome_list, **kwargs): - """ - All-vs-all ANI within a set of genomes, in-process. - - Each genome is sketched exactly once and then queried against the database, - so the sketching cost is linear in the number of genomes rather than - quadratic. - - Args: - genome_list: list of genome file locations. - - Returns: - Ndb: DataFrame with ['reference', 'querry', 'ani', 'alignment_coverage']. - """ - db = PyskaniDatabase(**kwargs) - - contigs = {} - for location in genome_list: - name = drep.d_cluster.utils._get_genome_name_from_fasta(location) - cs = load_contigs(location) - contigs[name] = cs - db.add(name, cs) - - logging.debug(f"pyskani: sketched {len(contigs)} genomes once; querying") - - rows = [] - for name, cs in contigs.items(): - rows.extend(db.query(name, cs)) - - rows = _fill_missing_pairs(rows, list(contigs.keys())) - Ndb = pd.DataFrame(rows, columns=NDB_COLUMNS) - - # A genome can hit itself with identity slightly below 1 depending on - # sketching; force exact self-identity like the other algorithms do. - self_mask = Ndb['reference'] == Ndb['querry'] - Ndb.loc[self_mask, 'ani'] = 1.0 - Ndb.loc[self_mask, 'alignment_coverage'] = 1.0 - - # Keep one row per ordered pair (query returns each direction once) - Ndb = Ndb.drop_duplicates(subset=['reference', 'querry'], keep='first') - return Ndb.reset_index(drop=True) - - -def pyskani_one_vs_many(location, db, **kwargs): - """ - Compare one genome against an existing PyskaniDatabase of representatives. - - Used by greedy secondary clustering: the database holds every current cluster - representative, already sketched, so this is a single in-process query rather - than a subprocess that re-sketches all representatives. - - Emits rows with the representative in the 'querry' column (mirroring the - fastANI greedy path), because get_cluster_rep reads the winning - representative from there. - - Returns an Ndb-shaped DataFrame (possibly empty if nothing is similar). - """ - rows = db.query_genome(location, query_as_reference=True) - if len(rows) == 0: - return pd.DataFrame(columns=NDB_COLUMNS) - return pd.DataFrame(rows, columns=NDB_COLUMNS) diff --git a/drep/d_cluster/union_find.py b/drep/d_cluster/union_find.py index 98564f0..c2baa8b 100644 --- a/drep/d_cluster/union_find.py +++ b/drep/d_cluster/union_find.py @@ -104,8 +104,7 @@ 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`` and for the ``low_ram`` path that - used to ``stack()`` an already-dense matrix back into long format. + -> scipy path in ``cluster_mash_database``. Args: db: DataFrame with columns 'genome1', 'genome2', and either 'dist' or @@ -239,6 +238,41 @@ def to_name(x): 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 diff --git a/drep/d_cluster/utils.py b/drep/d_cluster/utils.py index 771182b..590008a 100644 --- a/drep/d_cluster/utils.py +++ b/drep/d_cluster/utils.py @@ -103,8 +103,6 @@ def estimate_time(comps, alg): 'ANImf': .5, 'fastANI': 0.00667, 'skani': 0.00667, - # in-process; sketches each genome once, so at least as fast as skani - 'pyskani': 0.00667, } if alg not in per_comparison: logging.debug(f"No time estimate available for {alg}; assuming a fast algorithm") diff --git a/setup.py b/setup.py index dd25acd..5faac86 100644 --- a/setup.py +++ b/setup.py @@ -24,16 +24,7 @@ def version(): 'biopython', 'scikit-learn', 'tqdm', - 'networkx', 'setuptools', 'pytest' ], - extras_require={ - # In-process skani (--S_algorithm pyskani). Optional: dRep falls back - # to the skani/fastANI executables when it isn't installed. - # >=0.2 is required for the `cutoff` query argument. Note there is no - # macOS arm64 wheel for 0.2 yet, so Apple Silicon builds it from source - # and needs a Rust toolchain (conda install -c conda-forge rust). - 'pyskani': ['pyskani>=0.2'], - }, zip_safe=False) diff --git a/tests/tests/test_cluster.py b/tests/tests/test_cluster.py index 60c7ea9..de9c15e 100644 --- a/tests/tests/test_cluster.py +++ b/tests/tests/test_cluster.py @@ -609,29 +609,33 @@ 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 == "union_find_streaming", "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'] \ No newline at end of file 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 8c6fc47..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,7 +71,8 @@ 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 == "union_find_streaming", "Streaming union-find method was not used" @@ -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_pyskani.py b/tests/tests/test_pyskani.py deleted file mode 100644 index 1d2d311..0000000 --- a/tests/tests/test_pyskani.py +++ /dev/null @@ -1,177 +0,0 @@ -""" -Tests for the in-process pyskani backend (drep.d_cluster.pyskani_backend). - -pyskani is an optional dependency, so every test here skips cleanly when it -isn't installed. -""" -import glob -import os -import shutil -import tempfile - -import pandas as pd -import pytest - -import drep.d_cluster.compare_utils as cu -import drep.d_cluster.cluster_utils as clu -import drep.d_cluster.external as ext -import drep.d_cluster.greedy_clustering as gc -import drep.d_cluster.utils -import drep.d_filter - -def _has_pyskani(): - try: - import pyskani # noqa: F401 - return True - except ImportError: - return False - - -requires_pyskani = pytest.mark.skipif(not _has_pyskani(), reason="pyskani not installed") -requires_skani = pytest.mark.skipif(shutil.which('skani') is None, reason="skani not installed") -requires_fastani = pytest.mark.skipif(shutil.which('fastANI') is None, reason="fastANI not installed") - - -def _test_genomes(): - here = os.path.dirname(os.path.abspath(__file__)) - return sorted([g for g in glob.glob(os.path.join(here, '../genomes/*')) - if os.path.isfile(g)]) - - -def _partition(Cdb, col='secondary_cluster'): - return {frozenset(sub['genome']) for _, sub in Cdb.groupby(col)} - - -@requires_pyskani -def test_pyskani_ndb_shape(): - """run_pairwise_pyskani returns a complete, well-formed Ndb.""" - import drep.d_cluster.pyskani_backend as pb - genomes = _test_genomes() - Ndb = pb.run_pairwise_pyskani(genomes) - - assert list(Ndb.columns) == ['reference', 'querry', 'ani', 'alignment_coverage'] - # Every ordered pair present (dRep's hierarchical clustering needs a full matrix) - assert len(Ndb) == len(genomes) ** 2 - # ANI and coverage are on a 0-1 scale - assert Ndb['ani'].between(0, 1).all() - assert Ndb['alignment_coverage'].between(0, 1).all() - # Self comparisons are exactly 1 - selfs = Ndb[Ndb['reference'] == Ndb['querry']] - assert (selfs['ani'] == 1).all() - assert (selfs['alignment_coverage'] == 1).all() - - -@requires_pyskani -@requires_skani -def test_pyskani_agrees_with_subprocess_skani(): - """ - pyskani and the skani executable should report the same ANI/coverage, and - produce the same secondary clusters. - """ - import drep.d_cluster.pyskani_backend as pb - genomes = _test_genomes() - workdir = tempfile.mkdtemp() - try: - sub = ext.run_pairwise_skani(genomes, os.path.join(workdir, 'skani/'), processors=4) - py = pb.run_pairwise_pyskani(genomes) - - m = pd.merge(sub, py, on=['reference', 'querry'], suffixes=('_sub', '_py')) - # Only close pairs matter for clustering; skani's min-af filter drops - # distant pairs from pyskani's output (they're filled in as ani=0). - close = m[m['ani_sub'] >= 0.95] - assert len(close) > 0 - assert (close['ani_sub'] - close['ani_py']).abs().max() < 0.001 - assert (close['alignment_coverage_sub'] - close['alignment_coverage_py']).abs().max() < 0.001 - - for sa, nc in [(0.99, 0.1), (0.95, 0.1)]: - c1, _ = clu.genome_hierarchical_clustering(sub, S_ani=sa, cov_thresh=nc, - comp_method='skani', cluster='X') - c2, _ = clu.genome_hierarchical_clustering(py, S_ani=sa, cov_thresh=nc, - comp_method='pyskani', cluster='X') - assert _partition(c1) == _partition(c2), f"clusters differ at S_ani={sa}" - finally: - shutil.rmtree(workdir, ignore_errors=True) - - -@requires_pyskani -@requires_fastani -def test_greedy_pyskani_matches_greedy_fastani(): - """Greedy clustering should give the same answer via pyskani as via fastANI.""" - genomes = _test_genomes() - bdb = drep.d_cluster.utils.load_genomes(genomes) - bdb = drep.d_filter._add_lengthN50(bdb, bdb) - - workdir = tempfile.mkdtemp() - try: - parts = {} - for alg in ['fastANI', 'pyskani']: - d = os.path.join(workdir, alg + '/') - os.makedirs(d, exist_ok=True) - Ndb, Cdb, _ = gc.compare_genomes_greedy( - bdb, alg, d, S_ani=0.95, cov_thresh=0.1, cluster='P1', processors=4) - parts[alg] = _partition(Cdb) - assert parts['fastANI'] == parts['pyskani'] - finally: - shutil.rmtree(workdir, ignore_errors=True) - - -@requires_pyskani -def test_greedy_pyskani_sketches_each_genome_once(): - """ - The whole point of the pyskani greedy path: every genome is sketched exactly - once, no matter how many representatives accumulate. - """ - import drep.d_cluster.pyskani_backend as pb - genomes = _test_genomes() - bdb = drep.d_cluster.utils.load_genomes(genomes) - bdb = drep.d_filter._add_lengthN50(bdb, bdb) - - sketch_calls = [] - orig = pb.PyskaniDatabase.add - - def counting_add(self, name, contigs): - sketch_calls.append(name) - return orig(self, name, contigs) - - pb.PyskaniDatabase.add = counting_add - workdir = tempfile.mkdtemp() - try: - gc.compare_genomes_greedy(bdb, 'pyskani', os.path.join(workdir, 'g/'), - S_ani=0.95, cov_thresh=0.1, cluster='P1') - # One sketch per representative, and never the same genome twice - assert len(sketch_calls) == len(set(sketch_calls)), "a genome was sketched more than once" - assert len(sketch_calls) <= len(genomes) - finally: - pb.PyskaniDatabase.add = orig - shutil.rmtree(workdir, ignore_errors=True) - - -def test_estimate_time_handles_every_S_algorithm(): - """ - estimate_time only drives a log line, but it used to raise UnboundLocalError - for any algorithm it didn't know about -- which killed a 10k-genome run at - the start of secondary clustering. Every --S_algorithm choice must work, and - unknown ones must not raise. - """ - from drep.d_cluster.utils import estimate_time - - for alg in ['ANIn', 'gANI', 'goANI', 'ANImf', 'fastANI', 'skani', 'pyskani']: - t = estimate_time(100, alg) - assert t > 0, f"{alg} gave {t!r}" - - # An unrecognized algorithm must degrade gracefully, not raise - assert estimate_time(100, 'some_future_algorithm') > 0 - - -@requires_pyskani -def test_compare_genomes_dispatches_pyskani(): - """--S_algorithm pyskani is reachable through the normal dispatch path.""" - genomes = _test_genomes() - bdb = drep.d_cluster.utils.load_genomes(genomes) - workdir = tempfile.mkdtemp() - try: - Ndb = cu.compare_genomes(bdb, 'pyskani', workdir) - assert len(Ndb) == len(genomes) ** 2 - assert 'ani' in Ndb.columns - finally: - shutil.rmtree(workdir, ignore_errors=True) diff --git a/tests/tests/test_union_find.py b/tests/tests/test_union_find.py index 3f61a55..6357c7c 100644 --- a/tests/tests/test_union_find.py +++ b/tests/tests/test_union_find.py @@ -75,8 +75,8 @@ def test_cluster_mash_files_streaming(tmp_path): assert g2c['g0.fasta'] != g2c['g3.fasta'] -def test_low_ram_matches_scipy_membership(): - # Build a random symmetric similarity table and confirm union-find (low_ram) +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 @@ -236,8 +236,11 @@ def test_sparse_skani_primary_matches_mash(): # Every input genome is represented (including singletons skani screened out) assert set(Cdb_sk['genome']) == set(Bdb['genome']) - # Marker so downstream plotting skips the (nonexistent) primary dendrogram - assert cret[0] == 'union_find_streaming' + # 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')} From 626db0c447d90dee197390a1b509fc3cdbb22c33 Mon Sep 17 00:00:00 2001 From: Matt Olm Date: Wed, 15 Jul 2026 12:00:49 -0600 Subject: [PATCH 09/15] Fix ScaffoldLevel_dRep.py crashing on MUMmer 3 The script passed -t (threads) to nucmer unconditionally. That option only exists in MUMmer 4; MUMmer 3's nucmer rejects it outright rather than ignoring it: $ nucmer --mum -p out -c 65 -g 90 -t 6 ref.fa query.fa Unknown option: t exit 1 So the script crashed with "nucmer failed with exit code 1" for anyone running MUMmer 3 -- which is what `conda install mummer` still gives you, as 3.23. MUMmer 4 is a separate `mummer4` package, and nothing declared that dependency. dRep's own ANImf comparisons were unaffected because gen_nucmer_cmd never passed -t, which is why this went unnoticed. Detect whether nucmer understands -t and only pass it if so, rather than requiring a specific MUMmer major version. MUMmer 4 still gets its threading; MUMmer 3 now works instead of crashing. Fixes 4 failing tests in test_bonus.py, and adds a regression test that asserts -t is passed if and only if the installed nucmer supports it. Co-Authored-By: Claude Opus 4.8 --- helper_scripts/ScaffoldLevel_dRep.py | 25 +++++++++++++++- tests/tests/test_bonus.py | 43 ++++++++++++++++++++++++++++ 2 files changed, 67 insertions(+), 1 deletion(-) 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/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 From 1650d073b6f8b1d688000e1c186370c59dd7d504 Mon Sep 17 00:00:00 2001 From: Matt Olm Date: Wed, 15 Jul 2026 12:03:45 -0600 Subject: [PATCH 10/15] Add v4.0.0 release notes Documents the changes since 3.7.1 from a user's point of view, leading with what changes their results: skani as the default for both clustering steps, primary clustering switching to single linkage, the Mdb.csv semantics change, and the skani alignment-coverage fix. Deliberately omits development churn that never reached a release -- the pyskani backend and the --min-af 0 experiment were both added and removed within the v4 branch, so users never saw them and they are not "changes since 3.7.1". Marked Unreleased rather than dated; drep/VERSION stays at 3.7.1 until the release is actually cut. Co-Authored-By: Claude Opus 4.8 --- CHANGELOG.md | 111 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 111 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6840f45..4b9ca0a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,117 @@ 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] - Unreleased + +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 +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. Coverage filtering now actually applies. +- `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. + +### 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) From 034fb1872a09411a80b4572d706a25fe34fc579e Mon Sep 17 00:00:00 2001 From: Matt Olm Date: Wed, 15 Jul 2026 14:58:05 -0600 Subject: [PATCH 11/15] Drop unused progress arg from load_skani_sparse_edges The one-pass refactor replaced the streaming skani reader (which showed a tqdm bar) with a direct load, but the progress argument came along and was never wired up. An API that accepts progress=True and silently does nothing is worse than not offering it; the sparse edge list loads in about a second even for 10,000 genomes, and the minutes are spent inside skani itself, which reports its own progress. Co-Authored-By: Claude Opus 4.8 --- drep/d_cluster/union_find.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/drep/d_cluster/union_find.py b/drep/d_cluster/union_find.py index c2baa8b..2fdb0e0 100644 --- a/drep/d_cluster/union_find.py +++ b/drep/d_cluster/union_find.py @@ -334,7 +334,7 @@ def build_ndb_from_edges(edges, Cdb): 'Align_fraction_ref', 'Align_fraction_query'] -def load_skani_sparse_edges(sparse_files, progress=False): +def load_skani_sparse_edges(sparse_files): """ Load `skani triangle --sparse` output into a symmetric edge table. From bb51af7f9b464d6e87ca44e7d5f55ce5129da54d Mon Sep 17 00:00:00 2001 From: Matt Olm Date: Wed, 15 Jul 2026 16:25:42 -0600 Subject: [PATCH 12/15] Cycle a small palette for cluster colors instead of one shade each Coloring clusters is for telling neighbouring clusters apart, not for identifying which cluster is which -- nobody reads a color off a dendrogram and recovers "cluster 47". Giving every cluster its own shade off a colormap means that past a few dozen clusters they are all indistinguishable anyway. Cycle three UC Berkeley colors (Berkeley Blue, California Gold, Founders Rock) in cluster order, so adjacent clusters always differ. This also fixes two things that came along for the ride: - colors were not reproducible. The old code did an *unseeded* np.random.shuffle over a colormap, so the same analysis produced different figures on every run. - it used jet, which is perceptually non-uniform and misleading. Cluster labels come in two shapes ('2' for primary, '2_10' for secondary), so sorting is numeric per component: 2_2 sorts before 2_10, not after. Co-Authored-By: Claude Opus 4.8 --- drep/d_analyze.py | 53 +++++++++++++++++++++++-------------- tests/tests/test_cluster.py | 32 +++++++++++++++++++++- 2 files changed, 64 insertions(+), 21 deletions(-) diff --git a/drep/d_analyze.py b/drep/d_analyze.py index 454b539..544d0d9 100644 --- a/drep/d_analyze.py +++ b/drep/d_analyze.py @@ -1089,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 @@ -1100,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/tests/tests/test_cluster.py b/tests/tests/test_cluster.py index de9c15e..c9a0ad6 100644 --- a/tests/tests/test_cluster.py +++ b/tests/tests/test_cluster.py @@ -638,4 +638,34 @@ def test_mash_primary_algorithm(self): # 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'] \ No newline at end of file + 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 From 569c311feea592c6a4a56c1d7058e8c304ade8db Mon Sep 17 00:00:00 2001 From: Matt Olm Date: Wed, 15 Jul 2026 16:45:46 -0600 Subject: [PATCH 13/15] Release 4.0.0 Bump the version and date the changelog. This is a major version because it breaks compatibility, not because of how much landed: - --low_ram_primary_clustering was removed, so a v3 command line using it now exits with "unrecognized arguments" - Mdb.csv gains a column, becomes sparse rather than N^2, and omits genomes with no above-threshold pairs; anything parsing it needs updating - default clustering results change three ways (skani replaces MASH and fastANI, primary linkage goes average -> single, and the skani alignment-coverage fix cannot be disabled) - skani is now required for default runs; mash no longer is Any one of those is a semver MAJOR trigger on its own. Also document a known limitation rather than let the release notes imply the alignment-coverage story is finished. The fix in this release corrects skani's coverage *units*; it does not change how a low-coverage pair is handled once measured. cov_thresh still works by zeroing a pair's ANI before hierarchical clustering, and average linkage routes around that: a genome joins a cluster on the strength of its other relationships and is grouped with a member it never had adequate coverage with. With cov_thresh=0.5, a genome aligning over 1% of another still ends up clustered with it, and three genomes is enough to trigger it. That is a separate, pre-existing problem and it is deferred. Co-Authored-By: Claude Opus 4.8 --- CHANGELOG.md | 28 +++++++++++++++++++++++++--- docs/module_descriptions.rst | 34 +++++++++++++++++----------------- drep/VERSION | 2 +- 3 files changed, 43 insertions(+), 21 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4b9ca0a..e6abe9c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,7 +4,7 @@ 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] - Unreleased +## [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 @@ -13,7 +13,8 @@ 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 -coverage fix (see Fixed) cannot be turned off, and it was a genuine bug. +alignment-coverage fix (see Fixed) cannot be turned off, and it was a genuine +bug. ### Breaking changes @@ -94,7 +95,10 @@ coverage fix (see Fixed) cannot be turned off, and it was a genuine bug. `--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. Coverage filtering now actually applies. + 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 @@ -105,6 +109,24 @@ coverage fix (see Fixed) cannot be turned off, and it was a genuine bug. `--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 diff --git a/docs/module_descriptions.rst b/docs/module_descriptions.rst index 26a29bc..9cdcab4 100644 --- a/docs/module_descriptions.rst +++ b/docs/module_descriptions.rst @@ -6,7 +6,7 @@ dRep has 3 commands: compare, dereplicate, and check dependencies. To see a list $ dRep -h - ...::: dRep v3.7.1 :::... + ...::: dRep v4.0.0 :::... Matt Olm. MIT License. Banfield Lab, UC Berkeley. 2017 (last updated 2026) @@ -44,16 +44,16 @@ This workflow compares a set of genomes. For a list of all parameters, check the $ dRep compare -h usage: dRep compare [-p PROCESSORS] [-d] [-h] [-g [GENOMES ...]] - [--S_algorithm {ANImf,goANI,fastANI,gANI,skani,ANIn}] - [--primary_algorithm {MASH,skani}] + [--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] [--skani_extra SKANI_EXTRA] [--n_PRESET {normal,tight}] [-pa P_ANI] [-sa S_ANI] [-nc COV_THRESH] [-cm {total,larger}] - [--clusterAlg {average,weighted,single,median,centroid,complete,ward}] - [--primary_clusterAlg {average,weighted,single,median,centroid,complete,ward}] + [--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] @@ -81,7 +81,7 @@ This workflow compares a set of genomes. For a list of all parameters, check the issues than wildcard expansion (default: None) GENOME COMPARISON OPTIONS: - --S_algorithm {ANImf,goANI,fastANI,gANI,skani,ANIn} + --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 @@ -92,7 +92,7 @@ This workflow compares a set of genomes. For a list of all parameters, check the gANI = Identify and align ORFs; compare aligned ORFS goANI = Open source version of gANI; requires nsmimscan (default: skani) - --primary_algorithm {MASH,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 @@ -149,11 +149,11 @@ 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 {average,weighted,single,median,centroid,complete,ward} + --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 {average,weighted,single,median,centroid,complete,ward} + --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 @@ -217,10 +217,10 @@ This workflow dereplicates a set of genomes. For a list of all parameters, check 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] [--checkm_group_size CHECKM_GROUP_SIZE] - [--S_algorithm {ANIn,skani,gANI,fastANI,ANImf,goANI}] + [--S_algorithm {goANI,ANImf,gANI,skani,ANIn,fastANI}] [--primary_algorithm {skani,MASH}] [--no_reuse_primary_comparisons] [--primary_skani_min_af PRIMARY_SKANI_MIN_AF] @@ -228,8 +228,8 @@ This workflow dereplicates a set of genomes. For a list of all parameters, check [--skani_extra SKANI_EXTRA] [--n_PRESET {normal,tight}] [-pa P_ANI] [-sa S_ANI] [-nc COV_THRESH] [-cm {total,larger}] - [--clusterAlg {weighted,ward,complete,single,average,median,centroid}] - [--primary_clusterAlg {weighted,ward,complete,single,average,median,centroid}] + [--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] @@ -284,7 +284,7 @@ This workflow dereplicates a set of genomes. For a list of all parameters, check completeness of the genome), "contamination"(0-100 value of the contamination of the genome)] (default: None) - --checkM_method {taxonomy_wf,lineage_wf} + --checkM_method {lineage_wf,taxonomy_wf} Either lineage_wf (more accurate) or taxonomy_wf (faster) (default: lineage_wf) --set_recursion SET_RECURSION @@ -298,7 +298,7 @@ This workflow dereplicates a set of genomes. For a list of all parameters, check (default: 2000) GENOME COMPARISON OPTIONS: - --S_algorithm {ANIn,skani,gANI,fastANI,ANImf,goANI} + --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 @@ -366,11 +366,11 @@ 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 {weighted,ward,complete,single,average,median,centroid} + --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 {weighted,ward,complete,single,average,median,centroid} + --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 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 From b888eb4c57542baa63cdd7b354ac207d8afb6839 Mon Sep 17 00:00:00 2001 From: Matt Olm Date: Thu, 16 Jul 2026 16:39:02 -0600 Subject: [PATCH 14/15] update readme --- README.md | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/README.md b/README.md index 18fbda0..9d56200 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 (the O(N²) scaling that previously limited large runs — see [#259](https://github.com/MrOlm/drep/issues/259)). + +**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 From 9eebddd4d696b09a343de554592b83d20f0bdd53 Mon Sep 17 00:00:00 2001 From: Matt Olm Date: Thu, 16 Jul 2026 16:39:54 -0600 Subject: [PATCH 15/15] Revise dRep v4 section in README Updated README to reflect changes in dRep v4 regarding genome comparisons. --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 9d56200..659a448 100644 --- a/README.md +++ b/README.md @@ -16,7 +16,7 @@ Open source pre-print publication is available at ## ⚡ 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 (the O(N²) scaling that previously limited large runs — see [#259](https://github.com/MrOlm/drep/issues/259)). +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: