I’ve been building an automated pipeline over GTDB r232 to find placeholder clades with no existing nomenclatural claim and characterise them to SeqCode type-material standard. I’m a software engineer with very little specific knowledge in this field.
Two things in that process each cost me a full re-run before I worked out what was going on. I haven’t seen either written down anywhere obvious, so I’m posting them here in case they save somebody else the time.
The rp16 marker set fails as a unit on MAGs
I started my phylogeny with the sixteen syntenic ribosomal proteins from Hug et al. 2016, which seemed like the safe choice. Across 251 quality-passing genomes, eight of them carried none of the sixteen at all.
Those same eight genomes carried between 26 and 36 other ribosomal proteins, and CheckM2 scored every one of them at 90% complete or better. The explanation turned out to be simple once I looked: all sixteen sit in the S10-spc-alpha operon, so a single contig break takes out the entire marker set in one go. CheckM2 is model-based rather than marker-based, so nothing flags it.
That’s worth knowing because a marker set whose members all fail together isn’t really giving you sixteen independent observations. I ended up switching to every ribosomal protein KEGG defines, filtered to those present in at least 90% of genomes and single-copy in at least 95% of the ones carrying them. That left 41 markers, spread across the genome rather than concentrated in one locus.
KOfam HMMs can be much longer than the protein you’re aligning
The KOfam model for ribosomal protein L2 (K02886) has 515 match states. The bacterial protein is around 275 residues. KOfam models span all domains of life, so a bacterium ends up roughly half gaps against them purely by construction.
I didn’t anticipate this, and on a 41-marker concatenation it put median genome occupancy at 0.44. Every occupancy floor I’d consider reasonable then dropped all 251 genomes, which is a fairly clear signal that something is wrong with your filter rather than your data.
The fix is to trim the model positions that the whole taxon set leaves empty before you judge any individual genome, and only then do the real trimming on the genomes that survive. That took my alignment from 11,865 columns to 5,844 and moved median occupancy from 0.44 to 0.89.
Related
A KOfam per-KO threshold can land inside a continuous score distribution rather than between two separate ones. When that happens you get a clean-looking presence/absence split between two clades whose proteins actually differ by a few percent of a bit score.
Two of my four headline pathway findings died this way. In the worst case, twelve genomes cleared a 334-bit cutoff at scores between 334.7 and 346.0, while the sister clade missed it at 298.9 to 313.2. That reads as a presence and an absence and it very clearly isn’t one. I now report the best below-threshold hit for whichever side is recorded as lacking the function, calculated per lineage rather than pooled, because pooling near and distant relatives hid this case the first time.
Questions
Rule 23d gives a Candidatus name priority on validation, and withdraws that protection if the name isn’t validated by 1 January 2027.
Pallen et al. 2022 in IJSEM, which coined roughly 65,000 names, was peer reviewed. The r214 follow-up that added another 23,000 or so exists only as a Research Square preprint (10.21203/rs.3.rs-4235597/v1), and Rule 24a excludes preprints from effective publication.
So do the r214 names hold Rule 23d priority at all? I’ve been treating them as blocking, which at least fails closed, but it puts something like 19,000 clades in a bucket they may not belong in. If anyone here knows the answer, even a one-line reply would be genuinely useful to me.
What came out of the screen
For what it’s worth, the pipeline produced five unnamed phylum-level candidates with no claim against them. I characterised them over 296 genomes: Bakta annotation, KEGG module reconstruction, AAI by reciprocal best hits, branch-length placement against the reference tree, and a 41-marker concatenated tree of 240 taxa and 5,277 columns under LG+F+G4 with 1000 ultrafast bootstrap and 1000 SH-aLRT replicates.
All five come out monophyletic and well supported, and the protein tree independently reproduces the sister relationships I’d read off the GTDB reference tree, which was reassuring.
I’d rather be straightforward about my motivation here. I would like to describe one of these under the SeqCode, and part of the reason is getting a name of my own choosing into the record.
I’m also aware of what I’m not in a position to judge, mainly which rank any of these actually supports and whether the one metabolic result that survived my own checks means what it looks like it means. If a critical read would interest anyone, I’d value that a great deal, and I’m happy to put the numbers and the holes in them somewhere public.