Version Compatibility
Reference examples tested with: BRAKER3 3.0+, BUSCO 5.5+, Bakta 1.9+, Infernal 1.1+, InterProScan 5.66+, Prokka 1.14+, RepeatMasker 4.1+, RepeatModeler 2.0.4+ (-threads replaced -pa in 2.0.4), eggNOG-mapper 2.1+, pandas 2.2+, tRNAscan-SE 2.0+
Before using code patterns, verify installed versions match. If versions differ:
- Python:
pip show <package> then help(module.function) to check signatures
- CLI:
<tool> --version then <tool> --help to confirm flags
If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.
Genome Annotation Pipeline
"Annotate my genome assembly" -> Orchestrate prokaryotic (Bakta) or eukaryotic (BRAKER3) gene prediction, repeat masking (RepeatMasker), functional annotation (eggNOG-mapper, InterProScan), and ncRNA annotation (Infernal).
This is a workflow skill: it owns the chaining decisions and hand-offs, not the internals of any one step.
The governing principle
A gene set is ~95% right and 100% confident; its trustworthiness is decided at four seams, not inside the gene-finder.
- Pro- vs eukaryotic is THE fork, committed from taxonomy up front, and it fixes the genetic-code table. Prokaryote -> Bakta one-step (verify the genetic-code TABLE from GTDB-Tk classification, never guess — a Mycoplasma under table 11 splits every gene at internal UGA). Eukaryote -> multi-step RepeatMask -> BRAKER3 -> functional -> ncRNA. There is no general-purpose eukaryote annotator: alternative genetic codes, trans-splicing, and polycistronic transcription break standard pipelines.
- Annotate ONLY a decontaminated, QC-passed assembly. CheckM2 before prokaryotic annotation is non-negotiable: contamination >5% mixes two organisms' genes into one chimeric set; a gene-finder trained on a contaminated/fragmented assembly produces confidently-wrong models genome-wide that are invisible in the GFF3. Annotation quality is bounded above by assembly quality.
- The evidence set is a committed input, not an afterthought. Eukaryotic: RNA-seq BAM + protein (OrthoDB) evidence drives BRAKER3's high-confidence training-set mining (the real advance — learning from loci where transcripts AND homology agree). Committing RNA-seq (ideally Iso-Seq for isoforms+UTRs) is decided at project design; without it the annotation is one-isoform, CDS-only, UTR-less and silently poisons AS/3'-tag/APA analyses.
- Tool + DB version + date is a reproducibility commitment. Bakta's DB is versioned (record it); Prokka's is frozen ~2019 (a post-2019 gene is "hypothetical" in Prokka, "named" in Bakta — accessory-vs-core flips on tool vintage alone). For any comparison, re-annotate everyone with ONE pipeline + ONE DB version from FASTA.
Made-once commitments
| Commitment | Consequence inherited downstream |
|---|
| Pro- vs eukaryotic path + genetic-code table (from taxonomy) | The whole tool chain; a wrong code table splits genes at recoded stops |
| Decontaminated, QC-passed assembly (CheckM2 gate) | Chimeric gene set / genome-wide corrupt training if skipped; annotation quality is bounded by assembly quality |
| Evidence set (RNA-seq + protein) | BRAKER3 training quality; without RNA-seq the annotation is isoform-naive, UTR-less |
| Tool + DB version | Named-vs-hypothetical and accessory-vs-core flip on tool vintage; re-annotate all with one version for comparison |
Pipeline Overview
Assembled contigs
|
v
[0. Assembly QC] ----------> QUAST, BUSCO (confirm assembly quality)
|
+----- Prokaryotic? -----> Path A: Bakta (one-step annotation)
| |
| v
| Annotated genome (GFF3, GenBank, FASTA)
|
+----- Eukaryotic? ------> Path B: Multi-step pipeline
|
v
[1. Repeat Masking] ----> RepeatModeler + RepeatMasker
|
v
[2. Gene Prediction] ---> BRAKER3 (RNA-seq + protein evidence)
|
v
[3. Functional Annotation] -> eggNOG-mapper + InterProScan
|
v
[4. ncRNA Annotation] ---> Infernal + tRNAscan-SE
|
v
Annotated genome (GFF3, proteins, functional tables)
Path A: Prokaryotic Annotation (Bakta)
Bakta provides comprehensive one-step annotation for bacteria and archaea. Preferred over Prokka for new projects.
Database Setup
bakta_db download --output /path/to/bakta_db --type full
Run Bakta
bakta \
--db /path/to/bakta_db \
--output bakta_out \
--prefix my_genome \
--locus-tag MYORG \
--genus Escherichia --species "coli" \
--strain K12 \
--gram - \
--translation-table 11 \
--threads 8 \
assembly.fasta
Prokaryotic QC Checkpoint
import subprocess
import json
def validate_prokaryotic_annotation(bakta_dir, prefix, expected_cds_range=(500, 8000)):
'''
QC gates for prokaryotic annotation.
- CDS count in expected range for genome size
- tRNA count >= 20 (typical minimum for free-living bacteria)
- rRNA operons detected
'''
gff_file = f'{bakta_dir}/{prefix}.gff3'
feature_counts = {'CDS': 0, 'tRNA': 0, 'rRNA': 0, 'tmRNA': 0, 'ncRNA': 0}
with open(gff_file) as f:
for line in f:
if line.startswith('#'):
continue
fields = line.strip().split('\t')
if len(fields) >= 3 and fields[2] in feature_counts:
feature_counts[fields[2]] += 1
qc_pass = True
if not (expected_cds_range[0] <= feature_counts['CDS'] <= expected_cds_range[1]):
print(f'WARNING: CDS count {feature_counts["CDS"]} outside expected range ')
qc_pass =
feature_counts[] < :
()
qc_pass =
()
qc_pass, feature_counts
Path B: Eukaryotic Annotation
Step 1: Repeat Masking
BuildDatabase -name mygenome assembly.fasta
RepeatModeler -database mygenome -threads 8 -LTRStruct
RepeatMasker \
-lib mygenome-families.fa \
-pa 8 \
-xsmall \
-gff \
-dir repeat_out \
assembly.fasta
Repeat Masking QC Checkpoint
def check_repeat_content(repeatmasker_tbl, taxon='vertebrate'):
'''
Verify repeat content is within expected range for taxon.
Typical ranges:
- Vertebrate: 30-60%
- Insect: 15-45%
- Plant: 20-85%
- Fungus: 3-20%
'''
expected_ranges = {
'vertebrate': (30, 60), 'insect': (15, 45),
'plant': (20, 85), 'fungus': (3, 20)
}
low, high = expected_ranges.get(taxon, (5, 80))
with open(repeatmasker_tbl) as f:
for line in f:
if 'total interspersed' in line.lower():
pct = float(line.strip().split()[-1].replace('%', ''))
break
qc_pass = low <= pct <= high
if not qc_pass:
print(f'WARNING: Repeat content {pct:.1f}% outside expected range ({low}-{high}%) for {taxon}')
return qc_pass, pct
Step 2: Gene Prediction with BRAKER3
braker.pl \
--genome=repeat_out/assembly.fasta.masked \
--bam=rnaseq_sorted.bam \
--prot_seq=proteins.fa \
--softmasking \
--threads 8 \
--species=my_species \
--gff3 \
--workingdir=braker_out
braker.pl \
--genome=repeat_out/assembly.fasta.masked \
--bam=rnaseq_sorted.bam \
--softmasking \
--threads 8 \
--species=my_species \
--gff3
braker.pl \
--genome=repeat_out/assembly.fasta.masked \
--prot_seq=orthodb_proteins.fa \
--softmasking \
--threads 8 \
--species=my_species \
--gff3
Gene Prediction QC Checkpoint
busco \
-i braker_out/braker.aa \
-l <clade>_odb10 \
-o busco_annotation \
-m proteins \
--cpu 8
def check_gene_prediction(braker_gff, busco_summary, expected_genes_range=(15000, 35000)):
'''
QC gates after gene prediction.
- Gene count within expected range for genome
- BUSCO completeness > 90%
- Mean exons per gene > 1 (spliced genes expected in eukaryotes)
'''
gene_count = 0
exon_count = 0
with open(braker_gff) as f:
for line in f:
if line.startswith('#'):
continue
feature = line.strip().split('\t')[2] if len(line.strip().split('\t')) >= 3 else ''
if feature == 'gene':
gene_count += 1
elif feature == 'exon':
exon_count += 1
mean_exons = exon_count / gene_count if gene_count > 0 else 0
with open(busco_summary) as f:
for line in f:
if line.strip().startswith('C:'):
completeness = float(line.strip().split('C:')[1].split('%')[0])
issues = []
(expected_genes_range[] <= gene_count <= expected_genes_range[]):
issues.append()
completeness < :
issues.append()
mean_exons < :
issues.append()
()
(issues) == , issues
Step 3: Functional Annotation
emapper.py \
-i braker_out/braker.aa \
--output eggnog_results \
--cpu 8 \
-m diamond \
--tax_scope auto \
--go_evidence non-electronic \
--target_orthologs all \
--seed_ortholog_evalue 1e-5 \
--override
interproscan.sh \
-i braker_out/braker.aa \
-b interpro_results \
-f tsv,gff3 \
-goterms \
-pa \
-cpu 8
Functional Annotation QC Checkpoint
import pandas as pd
def check_functional_annotation(eggnog_annotations, total_genes):
'''
QC gate: > 60% of genes should have functional assignment.
Below 50% suggests database issues or highly divergent organism.
'''
cols = ['query', 'seed_ortholog', 'evalue', 'score', 'eggNOG_OGs', 'max_annot_lvl',
'COG_category', 'Description', 'Preferred_name', 'GOs', 'EC', 'KEGG_ko']
df = pd.read_csv(eggnog_annotations, sep='\t', comment='#', header=None)
df.columns = (cols + [f'c{i}' for i in range(len(df.columns) - len(cols))])[:len(df.columns)]
annotated = len(df[df['Description'] != '-'])
pct_annotated = annotated / total_genes * 100
has_go = len(df[df['GOs'] != '-'])
has_kegg = len(df[df['KEGG_ko'] != '-'])
print(f'Annotated: {annotated}/{total_genes} ({pct_annotated:.1f}%)')
print(f'With GO terms: {has_go}, With KEGG: ')
pct_annotated < :
()
pct_annotated >=
Step 4: ncRNA Annotation
tRNAscan-SE \
-E \
--thread 8 \
-o trna_results.txt \
--gff trna.gff \
assembly.fasta
cmpress Rfam.cm
cmscan \
--cpu 8 \
--cut_ga --rfam --nohmmonly \
--tblout rfam_results.tbl \
--fmt 2 \
--clanin Rfam.clanin \
Rfam.cm \
assembly.fasta
grep -v ' = ' rfam_results.tbl > rfam_results.deoverlapped.tbl
Merging Annotations
def merge_annotations(braker_gff, trna_gff, rfam_tbl, eggnog_tsv, output_gff):
'''Merge gene predictions, ncRNAs, and functional annotations into final GFF3.'''
import subprocess
subprocess.run([
'agat_sp_merge_annotations.pl',
'--gff', braker_gff,
'--gff', trna_gff,
'-o', output_gff
], check=True)
subprocess.run([
'agat_sp_statistics.pl',
'--gff', output_gff,
'-o', output_gff.replace('.gff3', '_stats.txt')
], check=True)
print(f'Merged annotations written to {output_gff}')
Common Errors
| Symptom | Cause | Fix |
|---|
| Chimeric gene set; training corrupted genome-wide | Annotated a contaminated assembly | CheckM2/FCS-GX gate BEFORE annotation; decontaminate first |
| Genes split at recoded stops; low coding density, high hypothetical | Wrong genetic-code table | Set the table from GTDB-Tk taxonomy, not a guess |
| Real NLR/immune gene families deleted; suspiciously gene-poor | Over-masking with an uncurated repeat library | Filter the RepeatModeler library against a protein DB; confirm conserved families survive; soft-mask -xsmall |
| Accessory genome inflated ~10x in a pangenome | Frozen-DB annotation drift (Panaroo) | Re-annotate all assemblies with ONE pipeline + DB version from FASTA; existing GenBank annotations are unusable for pangenomics |
| BUSCO looks great but models are wrong | BUSCO-only quality claim (certifies ~1000 easy genes) | Proteome-mode BUSCO on delivered proteins; add mono-exonic fraction + length distribution + mRNA:gene ratio |
| Low gene count | Repeat masking too aggressive | Soft-mask (-xsmall) not hard-mask; curate the TE library first |
| Isoform-naive annotation (mRNA:gene = 1.00) | No RNA-seq evidence | Add RNA-seq (Iso-Seq for isoforms/UTRs) to BRAKER3 |
Complete Pipeline Script
#!/bin/bash
set -e
GENOME="assembly.fasta"
RNASEQ_BAM="rnaseq_sorted.bam"
PROTEINS="orthodb_proteins.fa"
BAKTA_DB="/path/to/bakta_db"
THREADS=8
ORGANISM_TYPE="${1:-eukaryotic}"
if [ "$ORGANISM_TYPE" == "prokaryotic" ]; then
echo "Running prokaryotic annotation with Bakta"
echo "Step 0: Contamination/completeness gate (CheckM2 before annotation is non-negotiable)"
checkm2 predict --input $GENOME --output-directory checkm2_out --threads $THREADS
bakta --db $BAKTA_DB --output bakta_out --prefix genome --translation-table 11 \
--locus-tag MYORG --threads $THREADS $GENOME
echo "Done. Results in bakta_out/"
else
echo "Running eukaryotic annotation pipeline"
echo "Step 0: Assembly QC (contiguity + completeness before committing to annotation)"
quast.py $GENOME -o quast_out --threads $THREADS
busco -i $GENOME -l "${LINEAGE:?set the DEEPEST applicable clade dataset, e.g. insecta_odb10 / embryophyta_odb10; eukaryota_odb10 is too shallow and inflates completeness}" -o busco_asm -m genome --cpu
BuildDatabase -name mygenome
RepeatModeler -database mygenome -threads -LTRStruct
RepeatMasker -lib mygenome-families.fa -pa -xsmall -gff - repeat_out
braker.pl --genome=repeat_out/$( ).masked \
--bam= --prot_seq= \
--softmasking --threads --gff3 --workingdir=braker_out
busco -i braker_out/braker.aa -l -o busco_check -m proteins --cpu
emapper.py -i braker_out/braker.aa --output eggnog_out --cpu -m diamond
tRNAscan-SE -E --thread -o trna_out.txt --gff trna.gff
cmscan --cpu --cut_ga --rfam --nohmmonly --tblout rfam.tbl -- 2 --clanin Rfam.clanin Rfam.cm
Related Skills
- genome-annotation/prokaryotic-annotation - Bakta and Prokka details
- genome-annotation/eukaryotic-gene-prediction - BRAKER3 and AUGUSTUS options
- genome-annotation/repeat-annotation - Soft-masking before gene prediction
- genome-annotation/functional-annotation - eggNOG-mapper and InterProScan
- genome-annotation/ncrna-annotation - Infernal/Rfam and tRNAscan-SE detail
- rna-structure/ncrna-search - Covariance-model search, gathering thresholds, and clan resolution
- genome-annotation/annotation-qc - BUSCO genome-vs-proteome, OMArk, CheckM2 gates
- genome-assembly/assembly-qc - Pre-annotation assembly quality checks
- genome-intervals/gtf-gff-handling - GFF3/GTF hierarchy traversal, AGAT sanitizing/validation, coordinate conversion, and seqid-consistency checks on the merged annotation
- workflows/genome-assembly-pipeline - Upstream: hands off the decontaminated, QC-passed FASTA (with its QV/BUSCO)
References
- Salzberg SL (2019) Next-generation genome annotation: we still struggle to get it right. Genome Biology 20:92. DOI 10.1186/s13059-019-1715-2. (error propagation.)
- Gabriel L, Bruna T, Hoff KJ, et al (2024) BRAKER3: fully automated genome annotation using RNA-seq and protein evidence with GeneMark-ETP, AUGUSTUS, and TSEBRA. Genome Research 34:769-777. DOI 10.1101/gr.278090.123. (high-confidence training-set mining.)
- Tonkin-Hill G, MacAlasdair N, Ruis C, et al (2020) Producing polished prokaryotic pangenomes with the Panaroo pipeline. Genome Biology 21:180. DOI 10.1186/s13059-020-02090-4. (annotation-drift accessory inflation.)
- Schwengers O, Jelonek L, Dieckmann MA, et al (2021) Bakta: rapid and standardized annotation of bacterial genomes via alignment-free sequence identification. Microbial Genomics 7:000685. DOI 10.1099/mgen.0.000685.
