File structure checks

[epbrenner]

The problem:

It has been a known issue that sometimes, pangenomes will fail to complete through data_processing.R. This has prevented Klebsiella from running in prior attempts, and Abhirupa isolated the issue to a subset when the pangenome generation was run in batches and smaller pangenomes merged. Some genomes in one particular subset appeared to cause a failure, but we had no further diagnostic information. Because we have also seen inflation in feature matrix sizes as more genomes have been added, I originally sought to implement QC filtering to exclude contaminated or incomplete genomes. However, in my current testing using these QC filters, I have found another, smaller subset of genomes that do not complete in Panaroo.

Details:

prepareGenomes(user_bacs = "Staphylococcus epidermidis", max_checkm_contam = 0) returns 39 genomes with AST data, a very small set. Panaroo fails to generate a pangenome after 1 hour of running. Logs stop printing when invalid genes are being removed from GFFs, a necessary step for non-Prokka GFFs. Despite this freeze, CPUs are locked at 400% use, indicating 4 cores are fully occupied during this entire time. Disk I/O is static, with no read and a constant small write, but no output is observed. RAM use is static and low. No hardware approaches its limits during the freeze. This occurs with our older Panaroo 1.5.3 container, and a much more recent 1.7.0 container, and after purging all Docker images and caches and rebuilding everything freshly, and when purging the package and installing everything freshly. Pangenomes run and complete in 5-10min for other bugs with 400 genomes or fewer.

Approach:

In this PR, I am planning on subsetting this batch of 39 genomes down to identify an even smaller set that does not complete. Once I identify a few candidate GFFs that may be responsible, I will try to find atypical patterns in the way their gene headers are formatted that might lead to a CPU deadlock.

[epbrenner]

An addendum to the problem.

Positron:

prepareGenomes(user_bacs = "Staphylococcus epidermidis")
Querying cache for 36 taxon IDs.
Collected 872 distinct genome IDs (cache).
Initial summary: targets=872 | AMR genomes=38 | genome_data genomes=847
QC summary: kept=847 | dropped=26
Final summary:
targets : 872
bac_data : 872
genome_data : 847
amr_phenotype: 38
Post-filter distinct genomes: 3

.updateBVBRCdata(max_age_days = 0)

prepareGenomes(user_bacs = "Staphylococcus epidermidis")
Querying cache for 93 taxon IDs.
Collected 2211 distinct genome IDs (cache).
Initial summary: targets=2211 | AMR genomes=150 | genome_data genomes=2044
QC summary: kept=2044 | dropped=181
Final summary:
targets : 2211
bac_data : 2211
genome_data : 2044
amr_phenotype: 150
Post-filter distinct genomes: 79

RStudio:

prepareGenomes(user_bacs = "Staphylococcus epidermidis")
Querying cache for 90 taxon IDs.
Collected 926 distinct genome IDs (cache).
Initial summary: targets=926 | AMR genomes=62 | genome_data genomes=901
QC summary: kept=901 | dropped=26
Final summary:
targets : 926
bac_data : 926
genome_data : 901
amr_phenotype: 62
Post-filter distinct genomes: 27

.updateBVBRCdata(max_age_days = 0)

prepareGenomes(user_bacs = "Staphylococcus epidermidis")
Querying cache for 94 taxon IDs.
Collected 2212 distinct genome IDs (cache).
Initial summary: targets=2212 | AMR genomes=150 | genome_data genomes=2045
QC summary: kept=2045 | dropped=181
Final summary:
targets : 2212
bac_data : 2212
genome_data : 2045
amr_phenotype: 150
Post-filter distinct genomes: 79

[jananiravi]

My initial testing with Claude looking for the source of the inconsistency issue

Testing prepareGenomes(user_bacs = "Staphylococcus epidermidis") in two different IDEs (Positron and RStudio), before/after forcing a cache refresh:

The point here is that the taxon/genome counts are not reproducible run-to-run or machine-to-machine unless you force a full cache refresh first. Two fresh sessions, same species query, no code changes — and one returns 872 targets/3 post-filter genomes, the other 926/27, with no warning that the local BV-BRC cache is partially stale. Only after .updateBVBRCdata(max_age_days = 0) do both environments converge to the same numbers (2211/2212 targets, 79 post-filter genomes both times).

This is a distinct, arguably more fundamental bug than the Panaroo hang itself: whatever cache-staleness logic decides how many taxon IDs/genomes to pull is silently non-deterministic depending on what's already sitting in the local cache, and there's no signal to the user that their result set is incomplete. It also matters for the Panaroo investigation directly — the "39 genomes" from the original PR description were almost certainly generated against a stale/partial cache, so the actual offending-genome subset may not be stable across your own re-runs either.

[jananiravi]

@epbrenner, I only retrieved 17 (of your 39) Sep genomes.

[epbrenner]

Providing a status update while I'm working on this: I'm now testing if inconsistency in the ordering of entries in the metadata download plus our initial parsing into tabular format might be leading to variation in the counts. Just auditing the freshly downloaded cache file yields different numbers of rows each time the download is repeated (e.g., 830,312 vs 836,173).

[jananiravi]

Quick update on the GFFs of the 17 Sepi I downloaded:

On the GFFs themselves — one real structural anomaly found among the 17:

1282.3129 is the clear structural outlier of the 17. Quantified across all genomes:

Genome	CDS count	CDS with locus_tag=
1282.3129	2563	0
1282.2719	2490	2402
904315.3	2483	2452
904316.3	2385	2361
(all others)	~2200–2500	~97–99%

Every other genome has locus_tag on 96-99% of its CDS features (the small gaps are normal — pseudogenes/truncated genes). 1282.3129 has zero — not one of its 2563 CDS records carries a locus_tag attribute at all. Combined with the two other things I found on it earlier (no ##sequence-region headers, and self-referential placeholder contig IDs accn|1282.3129.con.0001 instead of a real NCBI accession like every other genome uses), this genome looks like it came through an older/different BV-BRC annotation pathway that never populated the standard GFF3 metadata fields. 1282.2719 is otherwise normal (2402/2490 locus_tags) — its only quirk is the missing sequence-region headers, which is probably cosmetic.

locus_tag is exactly the kind of field a Prokka-oriented tool like Panaroo would lean on as a stable per-feature key when doing non-Prokka GFF cleanup (you had mentioned this: "a necessary step for non-Prokka GFFs"). A genome where that key is entirely absent is a strong candidate for hitting a degenerate fallback code path inside Panaroo's gene-validity filter.

@epbrenner, maybe we can compare with 'good' Campy or S. argenteus? Should we look for this error in bad Klebs, too?

[jananiravi]

The S. argenteus run finished completely — full pangenome output written (gene_presence_absence_roary.csv, final_graph.gml, etc.), total wall time ~2 minutes for 73 genomes. That's right in your normal "5-10 min for ≤400 genomes" range.

From C:
New lead, and it's a good one: assembly fragmentation. Argenteus contig counts range 10-30 per genome. Staph epidermidis genomes range 4-144, with most sitting in the 56-144 range (only 1282.2719 is low, at 4). "Refinding genes" runs immediately after "trimming contig ends" — its entire job is reconciling genes split across contig boundaries, so a batch of mostly-fragmented draft assemblies is exactly the kind of input that step is built to struggle with.

Notably, this isn't "more total searches" — Staph epi's ~10,460/9,641 searches over 17/16 genomes is actually a similar per-genome rate to argenteus's 37,278 over 73 genomes (~600 vs ~510 per genome). The difference is entirely in per-item rate: argenteus held steady at ~3 items/sec throughout; Staph epi collapsed toward zero after item ~8-12. So whatever's wrong makes individual searches pathologically expensive, not more numerous — consistent with fragmentation forcing expensive boundary-reconciliation work per gene, or (alternate theory still on the table) IS-element/transposase-driven near-identical-sequence ambiguity making individual alignment searches degenerate.

Good next test: take the least fragmented Staph epi genomes (starting with 1282.2719 at 4 contigs) as a new mini-batch, and separately try adding one of the highly-fragmented ones (say 904315.3 at 144 contigs) to a slice of the clean argenteus set, to see which side of that experiment reproduces the stall.

Update:
Big update — Test 2 broke through. Its log now shows genuinely new content (not the same frozen line) — it's past "refinding genes" entirely and into a later validity-check pass, flagging 1282.2719's genes for internal stop codons ("invalid gene! ... Has stop: True"). CPU also dropped from ~608% to ~300%, consistent with moving out of the heavy parallel search phase into more serial post-processing.

So: the batch without 904315.3 sailed past the exact step where every batch containing 904315.3 has permanently frozen (Test 1 still stuck at item 11, now ~15 minutes in). This is about as clean a confirmation as we're going to get — 904315.3 (144 contigs, the most fragmented Staph epidermidis genome in the set) looks like the real trigger, not 1282.3129's missing locus_tags, not a generic environment slowdown.

[jananiravi]

Intermediate update from our conversation

Neither metric (#contigs or #pseudogenes) singles out 904315.3 cleanly by simple magnitude:

Pseudogene count doesn't correlate with contig count, and 904315.3 isn't the max:
36 contigs → 22 pseudogenes
56 contigs → 47 pseudogenes
...
87 contigs → 59 pseudogenes (904320.3 — the MOST pseudogenes, and it's in the successful Test 2 batch!)
90 contigs → 26 pseudogenes
144 contigs → 18 pseudogenes (904315.3 — the suspect, but mid-range on pseudogenes)
904315.3 is the clear outlier only on raw contig count (144, next-highest is 1282.3129 at 103, then a steep drop to ~90 and below). It is not the pseudogene outlier — 904320.3 has more than 3x its pseudogene count and sailed through fine in Test 2.

Boundary-truncation check (genes within 50bp of a contig edge, i.e. plausibly split/truncated by fragmentation): only 53 of 2501 features (2.1%) in 904315.3 are edge-adjacent. That's a small fraction, not an overwhelming majority — so it's not that most of its genes are boundary-split either.

So your pseudogene theory isn't confirmed by raw counts, and my fragmentation theory is weaker than I made it sound too — 904315.3 is really only unusual in one respect (144 contigs, a genuine outlier vs. everyone else's ≤103), not in pseudogene load or boundary-truncation rate. This matches his broader intuition that "more contigs → faster" generally, since 904320.3 (87 contigs, most pseudogenes) is fine — 904315.3 may be an outlier for a reason specific to that extreme contig count rather than fragmentation-in-general or pseudogene load-in-general.

⚠️ Correction — Test 1 actually finished! It broke through "refinding genes" (the "11it" display was just stale tqdm), continued through gene-family collapsing, and reached "writing output..." — the container already exited. So it wasn't a permanent hang after all, just a severe multi-minute-to-tens-of-minutes slowdown before resolving. That matches Evan's own experience exactly (his job "eventually finished... but took >5 hours rather than 5 minutes") — this is degraded performance, not an infinite loop. Good catch to keep watching rather than accept my earlier "permanently stuck" framing.

[jananiravi]

@epbrenner @AbhirupaGhosh | Claude and Evan agree with each other!

Root cause found: it's Panaroo's refind step, not a specific "bad" genome.

Confirmed with two independent tests on different genome sets:

• Our 17-genome S. epidermidis batch: --refind-mode off → completes in ~1-2 min (vs. permanently degraded / never finishing with refinding on).
• Evan's 39-genome batch: --refind-mode off → 70 seconds (vs. >5 hours with refinding on).

What we ruled out:

• Missing locus_tag (1282.3129, 0/2563 CDS annotated) — removing it alone didn't fix anything; a clean 73-genome S. argenteus control batch had 5 genomes with the identical zero-locus_tag pattern and ran fine (~2 min, 37k refind searches, no stall).
• Contig-count fragmentation (904315.3, 144 contigs, the outlier of the batch) — looked promising initially (adding it alone to a fast argenteus batch reproduced a slow window), but 904320.3 (87 contigs, more pseudogenes than 904315.3) ran through the same step fine elsewhere, so raw contig count isn't the determinant either.

What we found: every batch we tested that degraded still contained 1282.2719. With --refind-mode off, Panaroo's --remove-invalid-genes step flagged 25 genes in the entire 17-genome run — all 25 from 1282.2719 alone — each with an internal premature stop codon (pseudogene-like truncations). No other genome in the batch contributed any. That matches your original instinct about pseudogenes: refinding tries to relocate/reconcile these truncated, ambiguous ORFs against the rest of the pangenome, and that reconciliation is where the cost explodes, not a generic "many contigs" or "missing metadata" property.

Suggested next steps:

1. Short-term workaround: run with --refind-mode off (or strict) for affected batches — near-zero cost, per both our test and Evan's.
2. Longer-term: add a QC filter to .apply_metadata_qc() that caps the number of Panaroo-flagged/BV-BRC-annotated pseudogene features per genome, alongside the existing CheckM completeness/contamination gate — a genome like 1282.2719 with dozens of internal-stop genes should probably be excluded before it ever reaches Panaroo, independent of contig count.
3. Worth Evan double-checking whether his 39-genome bad set has a similar single-genome concentration of invalid/pseudogene calls, to confirm this generalizes beyond our 17-genome case.