Release 4.0.3: --primary_skani_min_af_both for plasmids (issue #311)

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

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`.

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

drep/VERSION

@@ -1 +1 @@
-4.0.2
+4.0.3

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\

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 "

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]

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