From 95c4520f0af4c5f9beb871d561279b35713b6f2e Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 8 Oct 2026 22:37:18 +0000 Subject: [PATCH] Add --primary_skani_min_af_both for plasmids; release 4.0.3 (issue #311) A primary-clustering edge only needed --primary_skani_min_af in one direction. A small plasmid that is mostly a shared transposon or resistance cassette aligns ~100% to every larger plasmid carrying it, so a few of them chain thousands of unrelated plasmids into one primary cluster, and secondary clustering then builds every pairwise comparison in it (741M Ndb rows and 350 GB for 27k Acinetobacter plasmids). The new flag requires the aligned fraction in both directions. Secondary clustering only groups pairs that align in both directions anyway, so with --primary_skani_min_af <= 100 x --cov_thresh and -pa <= -sa this splits primary clusters without changing the final clusters. On simulated plasmids sharing a transposon it turned 1 primary cluster into 9, with identical final clusters and winners and 9x fewer Ndb rows. It is opt-in, since the either-direction 15% default was tuned to reproduce the MASH partition on UHGG. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_013UpcPFwkh6QMusnD6ApThL --- CHANGELOG.md | 23 +++++++++++++++++++++ docs/choosing_parameters.rst | 2 ++ docs/module_descriptions.rst | 30 +++++++++++++++++++++++++++ drep/VERSION | 2 +- drep/argumentParser.py | 12 +++++++++++ drep/d_cluster/compare_utils.py | 7 +++++-- drep/d_cluster/union_find.py | 33 ++++++++++++++++++++++++++++-- tests/tests/test_union_find.py | 36 +++++++++++++++++++++++++++++++++ 8 files changed, 140 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 912b9b6..bad1b8e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,29 @@ 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.3] - 2026-10-08 + +**Plasmid datasets no longer collapse into one giant primary cluster.** Add +`--primary_skani_min_af_both` when dereplicating plasmids. + +### Added +- `--primary_skani_min_af_both` requires `--primary_skani_min_af` in both + directions of a pair, rather than either, for it to form a primary-clustering + edge (issue #311). A small plasmid that is mostly a shared transposon or + resistance cassette aligns ~100% to every larger plasmid carrying it, so with + the default either-direction check a few of them chain thousands of unrelated + plasmids into one primary cluster, and secondary clustering then has to build + every pairwise comparison in that cluster. Secondary clustering only groups + pairs that align in both directions anyway, so this splits primary clusters + without changing which genomes can end up in the same secondary cluster. On + simulated plasmids sharing a transposon it turned 1 primary cluster into 9, + with identical final clusters and winners and 9x fewer Ndb rows + +### Changed +- Releases are now published to PyPI (with a matching GitHub release) + automatically when a version bump is merged; 4.0.2 was the first. 4.0.1 was + never uploaded to PyPI, but its changes are all in 4.0.2 + ## [4.0.2] - 2026-09-30 ### Fixed diff --git a/docs/choosing_parameters.rst b/docs/choosing_parameters.rst index 3f4f1d6..72d30f8 100644 --- a/docs/choosing_parameters.rst +++ b/docs/choosing_parameters.rst @@ -181,6 +181,8 @@ dRep was developed for the use-case of bacterial dereplication, and there are so A major thing to be aware of is primary clustering. As described in :ref:`5. Importance of genome completeness`, genomes need to be >50% complete for primary clustering to work. If you're comparing entities in which you cannot assess completeness or in which you want to compare genomes that share only a limited number of genes (e.g. phage or plasmids), this a problem. The easiest way to handle it is to avoid primary cluster altogether with the parameter `--SkipMash`, or lower the primary clustering threshold with `-pa`. +For plasmids, also use ``--primary_skani_min_af_both``. Plasmids share mobile elements (transposons, resistance cassettes, replication modules), and a small plasmid that is mostly one of those elements aligns ~100% to every larger plasmid carrying it. By default a pair only needs ``--primary_skani_min_af`` in one direction to become a primary-clustering edge, so a few small plasmids like that can chain thousands of unrelated plasmids into one giant primary cluster, which makes secondary clustering extremely slow and memory-hungry. ``--primary_skani_min_af_both`` requires the aligned fraction in both directions instead. Secondary clustering only groups pairs that align in both directions anyway, so as long as ``--primary_skani_min_af`` is at most 100 x ``--cov_thresh`` this splits primary clusters without changing which plasmids can end up in the same secondary cluster. ``--skani_extra=--small-genomes`` is also worth adding for plasmids. + Also consider the effect of alignment coverage (:ref:`2. Minimum alignment coverage`) on hierarchical clustering (:ref:`6. Oddities of hierarchical clustering`). If your working with entities that are especially mosaic, like phage, this can be a bigger problem than with bacteria. Genome filtering and scoring is also a major factor. If your genomes can't be assessed by checkM, you can turn off quality filtering and the use of completeness and contamination when picking genomes with the flag `--ignoreGenomeQuality`. diff --git a/docs/module_descriptions.rst b/docs/module_descriptions.rst index dbdb6fa..2afbf0d 100644 --- a/docs/module_descriptions.rst +++ b/docs/module_descriptions.rst @@ -48,6 +48,7 @@ This workflow compares a set of genomes. For a list of all parameters, check the [--primary_algorithm {skani,MASH}] [--no_reuse_primary_comparisons] [--primary_skani_min_af PRIMARY_SKANI_MIN_AF] + [--primary_skani_min_af_both] [-ms MASH_SKETCH] [--SkipMash] [--SkipSecondary] [--skani_extra SKANI_EXTRA] [--n_PRESET {normal,tight}] [-pa P_ANI] [-sa S_ANI] [-nc COV_THRESH] @@ -117,6 +118,20 @@ This workflow compares a set of genomes. For a list of all parameters, check the reproduces the MASH partition closely; lower it only if you have very fragmented genomes and understand the chaining risk. (default: 15) + --primary_skani_min_af_both + Require --primary_skani_min_af in both directions of a + pair, not just one, for it to form a primary- + clustering edge (--primary_algorithm skani only). + Recommended for plasmids and other small genomes: a + small genome that is mostly a shared element (e.g. a + transposon) aligns ~100% to every larger genome + carrying it, and single linkage then chains unrelated + genomes into one huge primary cluster. Secondary + clustering only groups pairs that align in both + directions, so setting --primary_skani_min_af to at + most 100 x --cov_thresh does not change which genomes + can end up in the same secondary cluster. (default: + False) -ms MASH_SKETCH, --MASH_sketch MASH_SKETCH MASH sketch size (default: 1000) --SkipMash Skip primary clustering entirely and run secondary @@ -224,6 +239,7 @@ This workflow dereplicates a set of genomes. For a list of all parameters, check [--primary_algorithm {skani,MASH}] [--no_reuse_primary_comparisons] [--primary_skani_min_af PRIMARY_SKANI_MIN_AF] + [--primary_skani_min_af_both] [-ms MASH_SKETCH] [--SkipMash] [--SkipSecondary] [--skani_extra SKANI_EXTRA] [--n_PRESET {normal,tight}] [-pa P_ANI] [-sa S_ANI] @@ -336,6 +352,20 @@ This workflow dereplicates a set of genomes. For a list of all parameters, check reproduces the MASH partition closely; lower it only if you have very fragmented genomes and understand the chaining risk. (default: 15) + --primary_skani_min_af_both + Require --primary_skani_min_af in both directions of a + pair, not just one, for it to form a primary- + clustering edge (--primary_algorithm skani only). + Recommended for plasmids and other small genomes: a + small genome that is mostly a shared element (e.g. a + transposon) aligns ~100% to every larger genome + carrying it, and single linkage then chains unrelated + genomes into one huge primary cluster. Secondary + clustering only groups pairs that align in both + directions, so setting --primary_skani_min_af to at + most 100 x --cov_thresh does not change which genomes + can end up in the same secondary cluster. (default: + False) -ms MASH_SKETCH, --MASH_sketch MASH_SKETCH MASH sketch size (default: 1000) --SkipMash Skip primary clustering entirely and run secondary diff --git a/drep/VERSION b/drep/VERSION index 4d54dad..c4e41f9 100644 --- a/drep/VERSION +++ b/drep/VERSION @@ -1 +1 @@ -4.0.2 +4.0.3 diff --git a/drep/argumentParser.py b/drep/argumentParser.py index 80a178c..09235c8 100644 --- a/drep/argumentParser.py +++ b/drep/argumentParser.py @@ -147,6 +147,18 @@ def parse_args(args): "partition closely; lower it only if you have very fragmented genomes " "and understand the chaining risk.", default=15, type=float) + Clustflags.add_argument("--primary_skani_min_af_both", + help="Require --primary_skani_min_af in both directions of a pair, " + "not just one, for it to form a primary-clustering edge " + "(--primary_algorithm skani only). Recommended for plasmids and " + "other small genomes: a small genome that is mostly a shared " + "element (e.g. a transposon) aligns ~100%% to every larger genome " + "carrying it, and single linkage then chains unrelated genomes " + "into one huge primary cluster. Secondary clustering only groups " + "pairs that align in both directions, so setting " + "--primary_skani_min_af to at most 100 x --cov_thresh does not " + "change which genomes can end up in the same secondary cluster.", + action='store_true', default=False) Clustflags.add_argument("-ms", "--MASH_sketch", help="MASH sketch size", default=1000) Clustflags.add_argument("--SkipMash", help="Skip primary clustering entirely and run secondary\ clustering on all genomes at once. (Named for when primary clustering was\ diff --git a/drep/d_cluster/compare_utils.py b/drep/d_cluster/compare_utils.py index 7954faf..086a65d 100644 --- a/drep/d_cluster/compare_utils.py +++ b/drep/d_cluster/compare_utils.py @@ -194,6 +194,8 @@ def primary_cluster_skani_sparse(Bdb, data_folder, **kwargs): # 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) + # Require min_af in both directions rather than either (issue #311) + min_af_both = kwargs.get('primary_skani_min_af_both', False) # 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 @@ -206,14 +208,15 @@ def primary_cluster_skani_sparse(Bdb, data_folder, **kwargs): logging.info(f" Running sparse skani primary clustering on {len(genome_list):,} genomes " f"(ANI threshold {ani_threshold:.1f}%, screen {screen:.1f}%, " - f"min-af {min_af}%, emitting min-af {emit_min_af:.1f}%)") + f"min-af {min_af}%{' in both directions' if min_af_both else ''}, " + f"emitting min-af {emit_min_af:.1f}%)") sparse_file = drep.d_cluster.external.run_skani_triangle_sparse( genome_list, skani_folder, screen, min_af=emit_min_af, **kwargs) all_genomes = list(Bdb['genome'].unique()) edges = drep.d_cluster.union_find.load_skani_sparse_edges(sparse_file) Cdb, stats = drep.d_cluster.union_find.cluster_edges( - edges, P_ani, all_genomes, cov_threshold=min_af / 100.0) + edges, P_ani, all_genomes, cov_threshold=min_af / 100.0, cov_both=min_af_both) logging.info(f" Sparse skani primary clustering: {stats['edges_kept']:,} edges above " f"threshold, {stats['primary_clusters']:,} primary clusters " diff --git a/drep/d_cluster/union_find.py b/drep/d_cluster/union_find.py index 2fdb0e0..6f43201 100644 --- a/drep/d_cluster/union_find.py +++ b/drep/d_cluster/union_find.py @@ -386,7 +386,25 @@ def load_skani_sparse_edges(sparse_files): return edges[edges['genome1'] != edges['genome2']].reset_index(drop=True) -def cluster_edges(edges, ani_threshold, all_genomes, cov_threshold=0.0): +def _reverse_coverage(edges): + """ + For each row of an edge table, the alignment_coverage of the same pair in + the other direction (genome2's aligned fraction), or NaN if that row is + missing. + """ + g1 = edges['genome1'].values + g2 = edges['genome2'].values + names = pd.Index(pd.unique(np.concatenate([g1, g2]))) + n = len(names) + a = names.get_indexer(g1).astype(np.int64) + b = names.get_indexer(g2).astype(np.int64) + + cov = pd.Series(edges['alignment_coverage'].values, index=a * n + b) + cov = cov[~cov.index.duplicated(keep='last')] + return cov.reindex(b * n + a).values + + +def cluster_edges(edges, ani_threshold, all_genomes, cov_threshold=0.0, cov_both=False): """ Union-find clustering of an in-memory edge table (see load_skani_sparse_edges). @@ -397,6 +415,13 @@ def cluster_edges(edges, ani_threshold, all_genomes, cov_threshold=0.0): 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. + cov_both: require cov_threshold in both directions rather than either. + With either, a small genome that is mostly one shared element (e.g. + a plasmid that is mostly a transposon) aligns ~100% to every larger + genome carrying that element and chains them all together (issue + #311). Secondary clustering only joins pairs that align in both + directions anyway, so this splits primary clusters without changing + which genomes secondary clustering can put together. Returns: (Cdb, stats) @@ -407,7 +432,11 @@ def cluster_edges(edges, ani_threshold, all_genomes, cov_threshold=0.0): keep = edges['ani'].values >= ani_threshold if cov_threshold > 0: - keep &= edges['alignment_coverage'].values >= cov_threshold + cov = edges['alignment_coverage'].values + if cov_both: + # A missing other direction gives NaN, which fails the comparison + cov = np.minimum(cov, _reverse_coverage(edges)) + keep &= cov >= cov_threshold g1 = edges['genome1'].values[keep] g2 = edges['genome2'].values[keep] diff --git a/tests/tests/test_union_find.py b/tests/tests/test_union_find.py index 6357c7c..3628891 100644 --- a/tests/tests/test_union_find.py +++ b/tests/tests/test_union_find.py @@ -144,6 +144,42 @@ def test_skani_sparse_min_af_filters_low_alignment_edges(): assert s1['edges_kept'] == 1 +def test_cov_both_stops_small_genomes_chaining_clusters(): + """ + A small plasmid that is mostly a shared transposon aligns ~100% (in its own + direction) to every large plasmid carrying it, so with the aligned fraction + required in either direction it chains unrelated plasmids together (issue + #311). Requiring it in both directions keeps them apart. + """ + def pair(a, b, ani, af_a, af_b): + return [(a, b, ani, af_a), (b, a, ani, af_b)] + + rows = (pair('A1', 'A2', 0.995, 0.95, 0.96) # same plasmid family + + pair('B1', 'B2', 0.995, 0.94, 0.95) # another family + + pair('s', 'A1', 0.99, 0.95, 0.08) # small plasmid = transposon + + pair('s', 'B1', 0.99, 0.95, 0.10)) + edges = pd.DataFrame(rows, columns=['genome1', 'genome2', 'ani', 'alignment_coverage']) + allg = ['A1', 'A2', 'B1', 'B2', 's'] + + # Either direction (the default): s bridges the two families + C0, _ = uf.cluster_edges(edges, 0.90, allg, cov_threshold=0.5) + assert C0['primary_cluster'].nunique() == 1 + + # Both directions: two families plus s on its own + C1, s1 = uf.cluster_edges(edges, 0.90, allg, cov_threshold=0.5, cov_both=True) + g2c = C1.set_index('genome')['primary_cluster'].to_dict() + assert g2c['A1'] == g2c['A2'] + assert g2c['B1'] == g2c['B2'] + assert len({g2c['A1'], g2c['B1'], g2c['s']}) == 3 + assert s1['edges_kept'] == 4 + + # A pair with only one direction present can't pass a both-directions check + one_way = edges[~((edges['genome1'] == 'A2') & (edges['genome2'] == 'A1'))] + C2, _ = uf.cluster_edges(one_way, 0.90, allg, cov_threshold=0.5, cov_both=True) + g2c = C2.set_index('genome')['primary_cluster'].to_dict() + assert g2c['A1'] != g2c['A2'] + + def test_build_ndb_from_edges_fills_matrix(): """ Secondary clustering needs a complete matrix per primary cluster, so pairs