bio-genome-annotation-eukaryotic-gene-prediction

What this skill does

## Version Compatibility

Reference examples tested with: BRAKER 3.0+, GALBA 1.0.11+, AUGUSTUS 3.5+, Funannotate 1.8+, HISAT2 2.2.1+, STAR 2.7.11+, BUSCO 5.5+, samtools 1.19+, 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

Reproducibility is engineered, not assumed: run BRAKER/GALBA/Funannotate from their official containers, and pin the **OrthoDB partition version** (v11 vs v12 give different hints -> different models), the repeat library, and the Dfam release. GeneMark (inside BRAKER) historically required an expiring academic `.gm_key`; the requirement has relaxed in recent versions but check the installed version's docs. If code throws an error, introspect the installed tool and adapt rather than retrying.

# Eukaryotic Gene Prediction

**"Predict genes in my eukaryotic genome"** -> Identify protein-coding gene structures from a soft-masked assembly using RNA-seq and/or protein evidence to train and guide a gene finder.
- CLI: `braker.pl --genome=masked.fa --prot_seq=orthodb_clade.fa --rnaseq_sets_ids=SRR... --softmasking` (BRAKER3)

## The Single Most Important Modern Insight -- Training-Set Quality and Evidence Dominate, Not the Algorithm

Gene-prediction accuracy is governed by the quality of the **training set** and the **extrinsic evidence**, not by which gene-finder is chosen. A 2024-era HMM (AUGUSTUS) trained on a clean, evidence-validated gene set beats a fancier model trained on garbage. This reframes every decision:

- BRAKER3's accuracy jump is **not** a better HMM and **not** the TSEBRA combiner (the common misattribution). It is GeneMark-ETP's method of mining a **high-confidence training set** from loci where RNA-seq-assembled transcripts AND protein homology independently agree, then training AUGUSTUS on that (Brůna 2024 *Genome Res* 34:757; Gabriel 2024 *Genome Res* 34:769). BRAKER3 beats BRAKER1/2 because it learns from loci it has two reasons to trust.
- Therefore the **first** question is "what evidence do I have, and is it from the right clade?" - not "which tool?" Soft-masking and OrthoDB-partition choice are upstream of, and more consequential than, the predictor.
- **The #1 silent killer is training on a bad assembly.** A finder trained on a contaminated, fragmented, or repeat-polluted assembly produces confidently wrong models *genome-wide* that are syntactically valid, BUSCO-complete, and undetectable from the GFF3. Decontaminate and check assembly contiguity/BUSCO **before** any self-trainer touches the genome.
- The 2025-2026 deep-learning twist: Tiberius/Helixer can match BRAKER3's evidence-based accuracy with **no evidence at all** on well-represented clades (vertebrates) - but they degrade on under-represented lineages and are isoform/UTR-naive. The field is mid-transition; do not present DL as universally superior.

## Tool Taxonomy

| Pipeline | Citation | Evidence used | When |
|----------|----------|---------------|------|
| BRAKER3 | Gabriel 2024 *Genome Res* | RNA-seq + protein (OrthoDB) | **Default when both exist**; HC-training-set method |
| BRAKER1 | Hoff 2016 *Bioinformatics* | RNA-seq only | spliced-read intron hints train GeneMark-ET |
| BRAKER2 | Brůna 2021 *NAR Genom Bioinform* | protein only (broad OrthoDB) | ProtHint mines hints from a remote protein DB |
| GALBA | Brůna 2023 *BMC Bioinformatics* | protein only (close relatives) | beats BRAKER2 on large vertebrate genomes with good close proteomes |
| GeMoMa | Keilwagen 2018 *BMC Bioinformatics* | reference annotation + intron-position conservation | project an existing close-relative annotation |
| Funannotate | Palmer & Stajich 2020 (Zenodo) | any | **fungal de facto standard**; train->predict(EVM)->update(PASA) |
| MAKER2 / EVM | Holt 2011; Haas 2008 | many tracks | combiner / build-it-yourself route; transparent evidence weighting |
| Helixer / Tiberius | Stiehler 2020; Gabriel 2024 | none (ab initio, DL) | evidence-poor genomes in well-represented clades |
| AUGUSTUS / GeneMark-ES | Stanke 2006; Lomsadze | hints / self-train | components; standalone ab initio is the last resort |

OrthoDB rule: BRAKER2/3 want a **clade partition** (pick the smallest partition that still contains the target clade - Vertebrata for a fish, not all Metazoa). GALBA wants a few **close-relative proteomes**, not a broad clade.

## Decision Tree by Scenario

| Evidence / scenario | Recommended | Why |
|---------------------|-------------|-----|
| RNA-seq + protein DB | BRAKER3 | state-of-the-art at all genome sizes |
| RNA-seq only | BRAKER1 | intron hints from spliced reads |
| Protein only, close relatives | GALBA | miniprot-aligned close proteomes train AUGUSTUS |
| Protein only, broad/distant | BRAKER2 | ProtHint mines a remote OrthoDB clade |
| No evidence, represented clade | Tiberius or Helixer | DL ab initio now matches evidence-based on vertebrates |
| Reference annotation of a close relative | GeMoMa | homology + intron-position projection |
| Fungus (any evidence) | Funannotate | tiny-intron-aware; bundled EVM + PASA update |
| Need isoforms + UTRs | add PASA / Iso-Seq update step | predictors emit one CDS-only model per locus |
| Genome not yet masked | -> repeat-annotation (soft-mask first) | mandatory prerequisite |
| Want to assess the result | -> annotation-qc | BUSCO genome-vs-proteome, OMArk, sanity metrics |

## Soft-Masking Is Mandatory (Prerequisite)

Run repeat-annotation to **soft-mask** (repeats -> lowercase) before prediction. Unmasked TEs contain ORFs and pseudo-splice-sites; the predictor calls thousands of spurious genes inside repeats (catastrophic in plants where >80% of the genome can be TE) and TE domains pollute training. Pass `--softmasking` so BRAKER honors lowercase as a soft penalty (a real gene can still span a repeat). **Hard-masking (repeats -> N) destroys sequence and truncates real repeat-overlapping genes - avoid it for prediction.** But over-aggressive masking with an uncurated library deletes real multi-copy families (NLR/R-genes, zinc-fingers): filter the repeat library against a protein DB and confirm conserved families survive.

## BRAKER3 (RNA-seq + Protein)

```bash
# Align RNA-seq with a splice-aware aligner; output sorted BAM
hisat2-build masked.fasta idx
hisat2 -x idx -1 R1.fq.gz -2 R2.fq.gz --dta -p 16 | samtools sort -@4 -o rnaseq.bam
samtools index rnaseq.bam

# BRAKER3: protein = an OrthoDB clade partition (smallest that contains the clade)
braker.pl --genome=masked.fasta --prot_seq=Vertebrata.fa \
    --bam=rnaseq.bam --softmasking --threads=16 --species=my_species \
    --gff3 --workingdir=braker3_out
```

`--rnaseq_sets_ids=SRR...,SRR... --rnaseq_sets_dirs=/fastq/` auto-downloads/aligns reads in place of `--bam`. Outputs: `braker.gtf`/`braker.gff3` (TSEBRA-combined), `braker.codingseq`, `braker.aa`.

## GALBA (Protein-Only, Close Relatives)

```bash
galba.pl --genome=masked.fasta --prot_seq=close_relatives.faa \
    --species=my_species --threads=16
```

Prefer BRAKER3 whenever RNA-seq exists - intron evidence substantially improves splice-site accuracy.

## Funannotate (Fungi)

```bash
funannotate mask -i assembly.fa -o masked.fa --cpus 16
funannotate train -i masked.fa -o out -l R1.fq -r R2.fq --species "Genus species"
funannotate predict -i masked.fa -o out -s "Genus species" \
    --transcript_evidence transcripts.fa --protein_evidence proteins.fa
funannotate update -i out --cpus 16     # PASA adds UTRs and isoforms
```

## Gene-Model Sanity Statistics with Python

**Goal:** Compute the triage panel that reveals annotation health where gene count and BUSCO cannot - the isoform ratio, mono-exonic fraction, and protein-length distribution.

**Approach:** Load the GFF3 into gffutils; compute the mRNA:gene ratio (1.00 = isoform-naive), the single-exon fraction, and CDS-length stats; flag clade-anomalous values.

```