bio-genome-annotation-prokaryotic-annotation

What this skill does

## Version Compatibility

Reference examples tested with: Bakta 1.9+, Prokka 1.14.6 (legacy), tRNAscan-SE 2.0+, CheckM2 1.0+, gffutils 0.12+.

Before using code patterns, verify installed versions match. If versions differ:
- CLI: `<tool> --version` then `<tool> --help` to confirm flags
- Python: `pip show <package>` then `help(module.function)` to check signatures

Annotation content tracks the **database version**, not just the binary: a Bakta full DB (schema v5+, ~30 GB zipped) and a Bakta light DB give different functional calls, and two runs months apart can differ purely from DB updates. Record the Bakta DB version in methods. Prokka's bundled databases are frozen (~2019-2021). If code throws an error, introspect the installed tool and adapt rather than retrying.

# Prokaryotic Genome Annotation

**"Annotate my bacterial genome"** -> Call protein-coding genes, tRNAs/rRNAs/ncRNAs, and other features, then assign function by database identity, and emit submission-ready files.
- CLI: `bakta --db db/ assembly.fa` (default), `prokka --outdir out assembly.fa` (legacy), NCBI PGAP (submission/RefSeq-grade)

## The Single Most Important Modern Insight -- Gene Calling Is Near-Solved; Function Is Not, and Wrong Labels Propagate

For a typical isolate, Prodigal/Pyrodigal recovers >95-99% of true coding genes. The unsolved problems are the **start codon** (translation initiation site), the **small/overlapping/recoded ORFs**, and above all the **function**. Three consequences a postdoc must internalize:

1. **A confidently wrong product name is worse than "hypothetical protein."** Names assigned by loose homology transfer across distant lineages and *self-amplify* through databases - there is no mechanism to retract a correction once it spreads, so annotation accuracy has gone *down*, not up, as sequencing scaled (Salzberg 2019 *Genome Biol* 20:92). Bakta's design counters this: gene-level identity only via exact (MD5/UniRef100) match, falling back to cluster level, then "hypothetical" rather than guessing. 25-50% hypothetical is healthy for a non-model organism; near 0% means over-confident transfer.

2. **Comparing gene counts across tools/versions/DBs is invalid.** Tool agreement is "highly dependent on the organism of study" and biased toward model organisms (Dimonaco 2022 *Bioinformatics* 38:1198) - there is no universal best tool. Differences between two annotations are mostly tool artifacts, not biology. For any collection (pangenomics), **re-annotate every assembly from FASTA with one pipeline + one pinned DB**; merging published annotations inflates the accessory genome ~10× (Tonkin-Hill 2020 *Genome Biol* 21:180).

3. **Annotation completeness ≠ assembly completeness.** A fragmented/contaminated assembly produces truncated partial CDS at every contig break, inflated counts, and missing rRNA operons. Run CheckM2 *before* trusting any annotation QC number.

## Tool Taxonomy

| Tool | Maintained | Database approach | Output | When |
|------|-----------|-------------------|--------|------|
| Bakta | Yes (active) | Curated, **versioned**, alignment-free (UniRef + AMRFinderPlus + expert systems) | GFF3/GBFF/EMBL/FASTA/TSV/JSON + plot | **Default for new work**; reproducible, MAG-aware |
| Prokka | Frozen ~2021 | BLAST hierarchy vs frozen UniProt/RefSeq + HMM | GFF3/GBK/FAA + tbl | Legacy only; pangenome pipelines (Roary) expect Prokka GFF |
| NCBI PGAP | Yes (NCBI) | RefSeq protein-family models + ProSplign | ASN.1/GenBank, submission-ready | **GenBank/RefSeq submission**; best pseudogene/frameshift/selenoprotein handling |
| DFAST | Yes (DDBJ) | DFAST reference DBs + swappable refs | INSDC files | DDBJ submission; fast, flexible reference swap |
| RAST / BV-BRC | Yes | SEED subsystems | GenBank/GFF + subsystems | Subsystem/metabolic framing; web; not for direct INSDC submission |

Prokka uses Aragorn (tRNA) + barrnap (rRNA); Bakta uses tRNAscan-SE 2.0 (tRNA, domain-specific) + Infernal/Rfam (rRNA/ncRNA) + PILER-CR (CRISPR arrays) + DeepSig (signal peptides). Bakta is deliberately conservative on naming, so it reports "hypothetical" where Prokka confidently (sometimes wrongly) names a gene - Bakta looking "less annotated" is it being honest.

## Decision Tree by Scenario

| Scenario | Recommended | Why |
|----------|-------------|-----|
| Routine isolate, reproducible | Bakta, full DB | standardized, versioned, current functional calls |
| Submitting to GenBank/RefSeq | PGAP | RefSeq re-annotates with PGAP regardless; best pseudogene handling |
| Submitting to DDBJ | DFAST | INSDC-ready, DDBJ-aligned |
| MAG / metagenome bin | Bakta `--meta` + CheckM2 gate | anonymous-mode calling; QC before trust |
| Mycoplasma / Mollicutes | Bakta `--translation-table 4` | UGA=Trp; table 11 splits genes at every UGA |
| Archaeon | Bakta + verify tRNAscan-SE archaeal model; PGAP if N-termini matter | leaderless mRNAs degrade Prodigal TIS |
| Inherited Prokka pangenome pipeline | pin Prokka version, or re-call all in Bakta | tool consistency vs current biology |
| Subsystem/metabolic view | RAST / BV-BRC | SEED subsystem categories |
| Genome not yet assembled / poor QC | -> genome-assembly/assembly-qc | fix assembly before annotating |
| AMR for clinical report | -> run AMRFinderPlus/CARD-RGI directly | `--organism` context, point mutations, partials matter |

## Bakta (Default)

```bash
# Database (record the version): full ~30 GB for publishable annotation; light for triage/CI
bakta_db download --output db/ --type full

bakta --db db/ --output bakta_out --prefix ecoli_k12 \
    --genus Escherichia --species coli --strain K-12 \
    --locus-tag ECK12 --gram - --complete --threads 16 \
    assembly.fasta
```

Key flags: `--translation-table {11,4,25}` (default 11), `--gram {+,-,?}` (gates DeepSig signal-peptide calls; default `?`), `--complete` (all sequences are finished replicons; enables oriC detection - do NOT use on draft contigs), `--meta` (metagenome/MAG mode), `--compliant` (enforce INSDC structure), `--keep-contig-headers`, `--proteins <faa>` (trusted-protein transfer). Set `--genus`/`--species` from a GTDB-Tk classification, not a guess. Primary outputs: `.gff3`, `.gbff`, `.faa`, `.ffn`, `.fna`, `.tsv`, plus `.hypotheticals.tsv` and `.inference.tsv` (open the inference column to ask *why* a product was assigned).

## Prokka (Legacy)

```bash
prokka --outdir prokka_out --prefix my_genome --locustag MYORG \
    --genus Escherichia --species coli --cpus 8 --rfam assembly.fasta
```

Use only for tool-chain consistency with an existing Prokka-based pangenome workflow, and pin the version. Its bundled databases are frozen: a gene family characterized after ~2019 is "hypothetical" in Prokka but named by current Bakta/PGAP, so the *same gene* flips core/accessory purely on DB vintage.

## Coding Density and CDS Extraction with Python

**Goal:** Load Bakta/Prokka GFF3 into a queryable database and compute coding density, the first sanity number.

**Approach:** Build a gffutils database, sum CDS lengths, divide by genome length; flag values outside the expected band.

```python
import gffutils

CODING_DENSITY_LOW = 0.85   # <0.85 in a free-living bacterium suggests wrong table, fragmentation, or heavy pseudogenization
CODING_DENSITY_HIGH = 0.93  # >0.93 suggests ORF over-calling (spurious short hypotheticals)

def coding_density(gff_file, genome_length):
    db = gffutils.create_db(gff_file, ':memory:', merge_strategy='merge')
    coding_bp = sum(c.end - c.start + 1 for c in db.features_of_type('CDS'))
    density = coding_bp / genome_length
    if density < CODING_DENSITY_LOW or density > CODING_DENSITY_HIGH:
        print(f'WARNING: coding density {density:.1%} outside expected 85-93%')
    return density
```

## Hard Cases the Caller Gets Wrong

- **Wrong genetic code is silent and looks like fragmentation.** A Mycoplasma run under table 11 yields anomalously low coding density + many short "hypothetical" fragments because every internal UGA split a gene. Confirm