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