Calibrating a Biosynthetic Gene Cluster Detector: Released Profiles, Added Marker Families and a Reserve-Genome Test Against GECCO

1. Why a Second, Stricter Calibration Was Needed

Our earlier benchmark of a rule-based biosynthetic gene cluster (BGC) detector against GECCO, a conditional-random-field detector (Carroll et al., 2021), reported in Lab Discovery #59, found that the rule-based tool called far fewer clusters but a larger share of them overlapped experimentally curated loci in the MIBiG database (Zdouc et al., 2025). A follow-up on ten new genomes weakened that conclusion: the apparent precision advantage shrank when the rules were applied to genomes they had not been tuned on, and an automated choice between detector settings came down to an arbitrary tie-break. Two lessons followed. Any change to the detector has to be judged on genomes that were not used to make it, and the improvements need to be separable so each can be credited or discarded on its own.

This article reports that calibration. We changed three things in sequence (the source of the protein-family profiles, the set of marker families and core-gene rules, and the cluster-calling parameters), froze the result on a development set, and only then scored it on 9 reserve genomes (Figure 1). Throughout, the quantity we can measure is overlap with curated loci, not true precision: MIBiG is far from exhaustive, so a cluster that overlaps no curated locus may still be a real, uncurated BGC.

========================================================================================================================
  CALIBRATION PROTOCOL (fixed before any reserve genome was scanned)
========================================================================================================================
  DEVELOPMENT SET  : 15 genomes (76.3 Mb) examined in the first two benchmark rounds
                     9 bacterial genomes with curated loci | 2 fungal genomes (outside both tools' regime) | 4 BGC-free controls
        |
        |  step 1  swap seed-built profiles for released Pfam profiles        (measured on development set)
        |  step 2  add 14 marker families + new core-gene rules (revision 3)    (families shortlisted on development loci)
        |  step 3  grid search: 288 configurations (rule revision x window x seed score x gap x padding)
        |          objective: <= 1 cluster on development negative controls, then balanced score of
        |          locus recall and MIBiG-overlap share, then boundary overlap (IoU), then fewer clusters
        v
  FREEZE           : winning configuration written to disk  (V3 w=6 seed>=1.2 gap=5000 pad=1000)
        |
        v
  RESERVE SET      : 9 genomes (42.9 Mb) never examined before: 5 bacterial genomes with curated loci
                     + 4 BGC-free controls (3 bacteria, 1 fission yeast); GECCO run at its default cut-off
        |
        v
  SCORING          : loci overlapped, clusters called, boundary IoU, paired comparison with GECCO, negative-control calls
========================================================================================================================

2. Step One: Released Profiles Close Most of the Domain-Level Gap

The detector recognises biosynthetic enzymes by scanning every predicted protein against profile hidden Markov models (HMMs) of protein families from Pfam (Mistry et al., 2021), searched with the forward algorithm described by Eddy (2011). In the first version those profiles were rebuilt from the small Pfam seed alignments, some containing only 5–8 sequences, using a simplified prior. We exported the released Pfam profiles for the same families and loaded them into the identical search engine. On the development genomes the fraction of GECCO’s domain calls that the detector also made rose from 0.661 to 0.891, with precision against GECCO essentially unchanged (0.971 versus 0.962); on the second-round bacterial genomes it rose from 0.728 to 0.895 (Table 1). Locus recall changed far less (13 to 14 of 28 on the development genomes), so most of the domain-level gap did not translate into missed clusters, which are called from several genes at once.

The profile source was the only variable in this step, but it also changed the Pfam release (the released profiles shipped with GECCO come from an earlier Pfam version than the seed alignments we downloaded), so some of the effect could reflect that difference. One family, the SPASM domain, is absent from the released set and keeps its seed-built profile. A related variable, the domain cut-off rule (Pfam curated gathering threshold, a fraction of it, or GECCO’s p-value filter), changed domain agreement slightly but locus recall by at most one locus in either direction.

3. Step Two: Fourteen Added Marker Families and New Core-Gene Rules

Several curated loci were missed because they are built around families the detector did not model. To choose additions without guessing, we asked which Pfam domains that GECCO detects are over-represented in genes lying inside the 73 curated loci of the development genomes, relative to all other genes in the same genomes, and kept families that were both strongly enriched and biologically sensible for specialised metabolism (Table 2). The shortlist includes the terminal reductase domain of nonribosomal peptide synthetases, the MbtH-like partner protein of those enzymes, a bacterial terpene synthase family, the starter-unit transacylase of fungal polyketide synthases, pathway-regulator and siderophore-enzyme domains, and the melanin pathway cofactor. Because the enrichment was computed on the same genomes used for development, it is in-sample evidence; the reserve set tests whether it generalises.

Three new core-gene definitions use these families: an adenylation domain fused to a terminal reductase counts as an NRPS gene, a stand-alone adenylation gene beside an MbtH-like gene counts as an NRPS-like cluster seed, and the bacterial terpene synthase and the fungal starter-unit transacylase fused to a ketosynthase each define a core gene. The safeguards from the first version are unchanged: a lone ketosynthase or a lone adenylation gene still cannot seed a cluster. The MbtH family illustrates what these small marker proteins look like in the seed alignment (Figure 2).

========================================================================================================================
  FOCAL ALIGNMENT: MbtH-LIKE PROTEIN (PF03621) PFAM SEED, 60 MOST CONSERVED COLUMNS (112 sequences; columns 1-60 of 81)
  MbtH is a small partner protein of nonribosomal peptide synthetases and is used in revision 3 to recognise stand-alone adenylation genes
========================================================================================================================
  Column index    : .........1.........2.........3.........4.........5.........6   (tick = every 10 columns)
  D2B244_STRRD/1  : MSNPFEN-ADGTYLALINDEGQYSLWPSWA--EIPAGWTVAF---GE---DTRQACLDHI
  D1ACX0_THECD/1  : MSNPFED-PDGTFLVLVNDEGQYSLWPSFA--EVPAGWTVAK---SA---DTRQACLEFV
  Q9Z388_STRCO/1  : MTNPFED-ADGRYLVLVNDEGQHSLWPSFV--DVPAGWTVAL---GE---SDREACLEYV
  W5WFX2_9PSEU/1  : MTNPFEN-NDGRYLVLVNEEGQHSLWPSFV--DVPAGWTQVH---GE---ESRQSCLDYI
  Q2JA74_FRACC/3  : MTNPFED-ESASYFVLVNDEGQYSLWPVFA--EVPGGWSIVH---GE---DSRKACLDFI
  K0JVK5_SACES/1  : MTNPFDD-EAGTFLVLVNDEGQHSLWPQFA--EVPPGWTAAH---GP---DSRQSCVEYV
  Seed consensus  : MTNPFDD-EDGTFLVLVNDEGQYSLWPAFA--EVPAGWTVVF---GP---GSRQACLDYI
  Conservation    : :.:::..  .. . *::*.:.*.*:** ..   :*.**  .    :      :  *: :.   (* >=90% , : >=60% , . >=40% identity)
  Entropy (bits)  : ▄▅▃▄▃▃▄▁▆▆▅▆▃▆▁▂▄▁▅▂▄▂▃▁▃▁▁▆▅▆▂▁▅▄▁▆▁▁▅▆▄▆▃▂▁▃▆▅▁▁▆▆▃▆▅▁▃▅▃▃   (taller = more variable)
========================================================================================================================

4. Step Three: Parameter Search on the Development Set Only

The cluster caller has five tunable settings: how many genes on each side contribute support to a gene (window), the support score a core gene needs to seed a cluster, the largest gap between genes during extension, and the padding added to each cluster end. We scored 288 combinations of rule revision and settings on the development set. The selection rule was fixed beforehand: allow at most 1 cluster in the BGC-free development genomes, then maximise a balanced score (the harmonic mean of locus recall and the fraction of predicted clusters overlapping a curated locus), then the mean overlap of cluster boundaries with curated boundaries, then prefer fewer clusters. The winner was revision 3 with a window of 6 genes, a seed score of 1.2, a gap of 5,000 bp and 1,000 bp of padding. The configuration was written to disk before any reserve genome was scanned.

Two caveats apply to this choice. The winning window and seed score sit at the edge of the grid (the widest window and the lowest seed score tried), so a better setting may lie outside it and the optimum is not established. Many neighbouring settings scored identically (Table 3), which suggests the objective is flat in this region and the choice among them carries little weight.

5. Reserve-Set Results

The reserve set contains 5 bacterial genomes with 27 curated loci (Xenorhabdus nematophila, Streptomyces scabiei, Streptomyces avermitilis, Saccharopolyspora erythraea and Salinispora arenicola) and four BGC-free controls (Buchnera aphidicola, Chlamydia trachomatis, Treponema pallidum and the fission yeast Schizosaccharomyces pombe). On these genomes the published detector overlapped 9 of 27 curated loci (19–52% 95% interval). Released profiles raised that to 12, the revision-3 rules to 15, and the frozen parameters to 19 of 27 (70.4%, interval 52–84%), against 24 of 27 for GECCO (72–96%) (Table 4, Figure 3).

The gain came with more clusters: the frozen configuration called 105 clusters against 66 for the published detector, and GECCO called 178. The fraction of predicted clusters that overlap a curated locus was 18.1% for the frozen configuration and 12.9% for GECCO, a smaller difference than the one seen in the first benchmark, and one that cannot be interpreted as a precision difference because most clusters in every list are uncurated. Mean overlap between predicted and curated boundaries (computed over the loci each tool found) was 0.592 for the frozen configuration, 0.746 for the published detector and 0.550 for GECCO; the published detector’s higher value reflects that it only found the easiest loci, and the frozen configuration’s boundaries are looser than before.

========================================================================================================================
  RESERVE-SET LADDER: 5 BACTERIAL GENOMES, 27 CURATED LOCI (scale: loci bar = fraction of 27; cluster bar = fraction of 178)
========================================================================================================================
  Published detector (seed profiles) loci 9/27   #############                            |
                                     clusters 66    ===============                          | share 13.6%
  + released Pfam profiles           loci 12/27  ##################                       |
                                     clusters 73    ================                         | share 16.4%
  + revision-3 rules, default params loci 15/27  ######################                   |
                                     clusters 84    ===================                      | share 17.9%
  + frozen parameters                loci 19/27  ############################             |
                                     clusters 105   ========================                 | share 18.1%
  GECCO v0.11.0 (default 0.8)        loci 24/27  ####################################     |
                                     clusters 178   ======================================== | share 12.9%
------------------------------------------------------------------------------------------------------------------------
  share = fraction of predicted clusters overlapping any curated locus (a lower bound on precision; MIBiG is incomplete)
  Negative controls (4 BGC-free genomes): published 1 | frozen 1 | GECCO 18 clusters
========================================================================================================================

For comparison, the same configurations on the development genomes, where the settings were chosen and the numbers are therefore optimistic, are in Table 5; the improvement there was of similar size (30 to 46 of 64 loci), and GECCO overlapped 52.

========================================================================================================================
  RESERVE GENOMES: LOCI OVERLAPPED / CLUSTERS CALLED (published | released profiles | revision 3 | frozen | GECCO)
========================================================================================================================
  FN667742.1   Xenorhabdus nematophila ATCC 1  4.43 Mb  loci  6 |   3/10   3/11   3/13   4/13  |   5/13
  FN554889.1   Streptomyces scabiei 87.22     10.15 Mb  loci  6 |   3/15   4/15   5/18   5/22  |   5/47
  BA000030.4   Streptomyces avermitilis MA-46  9.03 Mb  loci  5 |   1/18   1/18   1/20   2/26  |   4/49
  AM420293.1   Saccharopolyspora erythraea NR  8.21 Mb  loci  5 |   1/10   1/13   2/15   4/22  |   5/42
  CP000850.1   Salinispora arenicola CNS-205   5.79 Mb  loci  5 |   1/13   3/16   4/18   4/22  |   5/27
  NC_002528.1  Buchnera aphidicola APS         0.64 Mb  loci  0 |    0/0    0/0    0/0    0/0  |    0/1
  NC_000117.1  Chlamydia trachomatis D/UW-3    1.04 Mb  loci  0 |    0/0    0/0    0/0    0/0  |    0/1
  NC_000919.1  Treponema pallidum Nichols      1.14 Mb  loci  0 |    0/1    0/1    0/1    0/1  |    0/1
  NC_003421.2  Schizosaccharomyces pombe chr   2.45 Mb  loci  0 |    0/0    0/0    0/0    0/0  |   0/15
------------------------------------------------------------------------------------------------------------------------
  TOTAL (bacterial, with loci)                              loci 27 |   9/66  12/73  15/84 19/105  | 24/178
========================================================================================================================

6. Paired Comparison With GECCO and Which Loci Remain Missed

Because both tools were scored on the same 27 loci, they can be compared pairwise. 19 loci were overlapped by both, 0 by the frozen configuration alone, 5 by GECCO alone, and 3 by neither. The exact McNemar test on the discordant loci gives p = 0.063, so with this few loci a real difference in sensitivity can be neither established nor excluded at the conventional 0.05 level; the interval for GECCO (72–96%) overlaps that of the frozen configuration. The 5 loci found only by GECCO are of the classes 2 other, 1 terpene, 2 ribosomal, which are the classes the detector models least well: ribosomal peptides, terpenes and loci annotated ‘other’ (Table 7). Within nonribosomal peptide and polyketide classes, which the detector was designed around, the two tools agreed on every reserve locus (Table 8).

7. Negative Controls and Known Failure Modes

In the 4 reserve controls the frozen configuration called 1 cluster (in T. pallidum) against 18 for GECCO, which included 15 in the fission yeast chromosome. In the 4 development controls the counts were 1 and 14. The detector’s one reserve call is a 334-residue adenylation gene with a 78-residue carrier-protein gene about 1.9 kb away, the pattern the NRPS-like rule is designed to catch but which can also arise in primary metabolism. No curated cluster is known in these genomes, but absence of a curated cluster is not proof that none exists.

The lower seed score also admits single-gene calls that are probably housekeeping. Applied to E. coli K-12 with the frozen configuration, the detector reports the enterobactin locus plus two RiPP-class clusters: one built on a conserved YcaO-family gene beside a radical-SAM gene, the other on a fused radical-SAM/SPASM gene (the detector does not search for precursor peptides, so neither call has that supporting evidence). We treat both as likely false positives. They were found after the configuration was frozen, so we have not changed the rules in response; doing so would require a new reserve set. The same loosening explains the larger cluster counts and the looser boundaries noted above.

Discussion & Key Takeaways