bio-genome-annotation-annotation-qc

What this skill does

## Version Compatibility

Reference examples tested with: BUSCO 5.5+, OMArk 0.3+, CheckM2 1.0+, compleasm 0.2.6+, gffutils 0.12+, matplotlib 3.8+.

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

BUSCO results depend on the **lineage dataset** (e.g. `vertebrata_odb10` vs the shallow `eukaryota_odb10`) and the OrthoDB version - record them; a 99% on a shallow ~255-gene set and a 99% on a deep ~5,500-gene clade set are very different claims. If code throws an error, introspect the installed tool and adapt rather than retrying.

# Annotation QC

**"Is my genome annotation any good?"** -> Measure conserved-gene completeness, proteome consistency/contamination, and gene-set sanity, and decide whether the limiting factor is the assembly or the predictor.
- CLI: `busco -i proteins.faa -m proteins -l <lineage>_odb10` (proteome) and `busco -i genome.fa -m genome -l <lineage>_odb10` (assembly), `omark`, `checkm2 predict` (prokaryotes)

## The Single Most Important Modern Insight -- BUSCO Measures the Easy 10%; the Diagnostic Is Genome-vs-Proteome

BUSCO completeness measures recovery of the ~1,000-5,500 most conserved single-copy orthologs - the most conserved, highly expressed, intron-stable genes any half-broken pipeline will find. A 98%-complete BUSCO confirms the housekeeping core is intact; it says **nothing** about whether the other ~20,000 models are chimeric, fragmented, frame-shifted, fused, or hallucinated from TEs. BUSCO is a smoke detector in one room of a burning house. Three load-bearing moves:

1. **Run assembly-BUSCO vs proteome-BUSCO on the same assembly - the diagnostic almost nobody runs.** `-m genome` does its own gene-finding; `-m proteins` scores the delivered proteome. If assembly is 98% and proteome is 85%, the genes are physically present and the *predictor* missed them - a training/evidence/masking problem, fixable without touching the assembly. If both are low, the genes aren't in the assembly - stop annotating and fix the assembly. This single fork redirects more wasted effort than any other check. (Always run BUSCO on the *delivered proteome*, `-m proteins`, not genome mode reported as if it described the annotation.)
2. **Gene count is a vanity metric.** The *same* assembly annotated by two pipelines routinely differs 20-40% in gene count, and a higher count is as likely to mean spurious TE ORFs and split models as more real genes. The diagnostic signal is in the **mono-exonic fraction, protein-length distribution, mean exons/gene, and mRNA:gene ratio** - not the headline count.
3. **Complement BUSCO with OMArk.** BUSCO's single-copy-ortholog lens is blind to over-prediction, chimeras, and contamination; OMArk assesses the *whole* proteome for completeness AND consistency (are genes consistent with their homologs?) and detects contaminant species. Use both.

## Tool Taxonomy

| Tool | Citation | Measures | Domain |
|------|----------|----------|--------|
| BUSCO | Manni 2021 *Mol Biol Evol* | conserved single-copy ortholog recovery (C/D/F/M) | euk + prok + viral; proteome/genome/transcriptome modes |
| OMArk | Nevers 2025 *Nat Biotechnol* | proteome completeness + consistency + contamination | eukaryotic proteomes; catches over-prediction BUSCO misses |
| compleasm | Huang 2023 | faster miniprot-based BUSCO-compatible completeness | euk; genome mode |
| CheckM2 | Chklovski 2023 *Nat Methods* | completeness + contamination (ML, lineage-agnostic) | bacteria/archaea, isolates + MAGs |
| Gene-set sanity panel | (gffutils) | gene count, mono-exonic %, protein length, mRNA:gene, coding density | any; the panel BUSCO can't provide |
| LAI | Ou 2018 *NAR* | LTR-RT assembly resolution (repeat-space contiguity) | LTR-rich genomes; assembly-side QC |

## Decision Tree by Scenario

| Scenario | QC to run | Why |
|----------|-----------|-----|
| Prokaryotic genome/MAG | CheckM2 (completeness/contamination) + coding density + tRNA/rRNA counts | CheckM2 is the field-standard completeness call; BUSCO is conservative on prokaryotes |
| Eukaryotic annotation | BUSCO `-m proteins` AND `-m genome` + OMArk + gene-set sanity panel | the genome-vs-proteome fork + over-prediction/contamination |
| Suspect high gene count | mono-exonic fraction + protein-length + BUSCO-Duplicated | distinguish haplotigs / unmasked TEs / split models |
| High BUSCO-Duplicated | check synteny + clade ploidy | uncollapsed haplotigs (purge) vs real WGD (keep) |
| Transferred annotation | BUSCO on the lifted set vs reference + ORF integrity | quantify silently lost conserved genes |
| Low proteome-BUSCO, high genome-BUSCO | fix evidence/masking/training | the genes are present; the predictor missed them |
| Both BUSCO low | -> genome-assembly/assembly-qc (purge_dups, more data) | the limit is the assembly |

## BUSCO and How to Read It

```bash
busco -i proteins.faa -m proteins -l vertebrata_odb10 -o busco_prot -c 16   # the delivered proteome
busco -i genome.fa    -m genome   -l vertebrata_odb10 -o busco_genome -c 16 # the assembly
```

Reported as `C:[S,D],F,M`: **Complete** (full-length match), **Duplicated** (subset of Complete, found ≥2x), **Fragmented** (partial), **Missing**. Use the **deepest applicable clade dataset**, not the shallow `eukaryota_odb10` (a 99% on a 255-gene set is trivially easy and not comparable to a 99% on a 5,500-gene clade set). **High Duplicated is the most-misread signal** - three causes, opposite responses: uncollapsed haplotigs (genome-wide, heterozygosity-scaled, non-syntenic -> purge_dups *before* annotating), real recent WGD (syntenic, ploidy-consistent, Ks duplication peak -> keep), or split models from a fragmented assembly. A clean haploid annotation runs ~1-3% Duplicated; >5-8% with no known WGD screams purge_dups first. compleasm is a faster miniprot-based alternative that often reports higher completeness than BUSCO's metaeuk path.

## OMArk (Proteome Consistency and Contamination)

```bash
omamer search --db LUCA.h5 --query proteins.faa --out proteins.omamer
omark -f proteins.omamer -d LUCA.h5 -o omark_out
```

OMArk catches what BUSCO's single-copy lens misses: it classifies the *whole* proteome as consistent / inconsistent (a gene whose structure conflicts with its homologs - a chimera or fragment) / unknown, and detects **contaminant species** mixed into the proteome. A high "inconsistent" fraction flags over-prediction or fused/split models even when BUSCO is green.

## CheckM2 (Prokaryotic Completeness/Contamination)

```bash
checkm2 predict --input genome.fna --output-directory checkm2_out --threads 16
```

For bacteria/archaea, CheckM2's lineage-agnostic ML model is the standard completeness/contamination call (handles reduced/novel lineages where marker sets fail). Gate annotation on it: **contamination >5%** mixes two organisms' genes into a chimeric set; **completeness <90% with contamination >5-10%** makes gene count, coding density, and hypothetical fraction uninterpretable - fix the assembly/binning first.

## Gene-Set Sanity Panel with Python

**Goal:** Compute the triage panel that reveals annotation health where gene count and BUSCO cannot.

**Approach:** Load the GFF3 into gffutils; compute the mRNA:gene ratio (1.00 = isoform/UTR-naive), the mono-exonic fraction, and the protein-length distribution; flag clade-anomalous values.

```python
import gffutils

MONOEXONIC_FLAG = 0.30   # >30% single-exon in a vertebrate suggests unmasked TEs/pseudogenes/fragments (calibrate per clade; fungi run higher)

def sanity_panel(gff_file):
    db = gffutils.create_db(gff_file, ':memory:', merge_strategy='merge')
    genes = list(db.features_of_type('gene'))
    mrnas = list(db.features_of_type(['mRNA', 'transcript']))
    exon_counts = [len(list(db.children(tx, featuretype='exon'))) for tx in mrnas]
    mono_frac = sum(1 for e in exon_counts if e == 1) /