Orchestrates genome annotation from assembled contigs to functional annotation, forking prokaryotic (Bakta one-step, genetic-code table from GTDB-Tk) vs eukaryotic (RepeatMask -> BRAKER3 -> functional -> ncRNA), then eggNOG/InterProScan functional assignment and Infernal/tRNAscan ncRNA. Use when committing the pro-vs-eukaryotic path and the genetic-code table from taxonomy (never guessing), annotating ONLY a decontaminated QC-passed assembly (CheckM2 before prokaryotic annotation is non-negotiable), committing the evidence set (RNA-seq + protein drives BRAKER3 training), soft-masking with a curated repeat library before gene prediction, or pinning the tool + DB version for any pangenome comparison. Hands mechanism to the genome-annotation component skills; not a re-teach of any single step.
Instrucciones de origen · Vista previa de solo lectura
name
bio-workflows-genome-annotation-pipeline
description
Orchestrates genome annotation from assembled contigs to functional annotation, forking prokaryotic (Bakta one-step, genetic-code table from GTDB-Tk) vs eukaryotic (RepeatMask -> BRAKER3 -> functional -> ncRNA), then eggNOG/InterProScan functional assignment and Infernal/tRNAscan ncRNA. Use when committing the pro-vs-eukaryotic path and the genetic-code table from taxonomy (never guessing), annotating ONLY a decontaminated QC-passed assembly (CheckM2 before prokaryotic annotation is non-negotiable), committing the evidence set (RNA-seq + protein drives BRAKER3 training), soft-masking with a curated repeat library before gene prediction, or pinning the tool + DB version for any pangenome comparison. Hands mechanism to the genome-annotation component skills; not a re-teach of any single step.
[{"after_repeat_masking":"Repeat content within expected range for taxon"},{"after_gene_prediction":"Gene count plausible, BUSCO completeness >90%"},{"after_functional_annotation":">60% of genes with functional assignment"}]
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
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
# Set --translation-table from the GTDB-Tk classification, never a guess: table 11 for most# bacteria, but --translation-table 4 for Mycoplasma/Spiroplasma (UGA = Trp, not stop) --# annotating a Mycoplasma under table 11 splits every gene at its internal UGA codons.# Add --complete ONLY for finished replicons; omit it for draft contigs (the common input).
Prokaryotic QC Checkpoint
import subprocess
import json
defvalidate_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}
withopen(gff_file) as f:
for line in f:
if line.startswith('#'):
continue
fields = line.strip().split('\t')
iflen(fields) >= 3and fields[2] in feature_counts:
feature_counts[fields[2]] += 1
qc_pass = Trueifnot (expected_cds_range[0] <= feature_counts['CDS'] <= expected_cds_range[1]):
print(f'WARNING: CDS count {feature_counts["CDS"]} outside expected range {expected_cds_range}')
qc_pass = Falseif feature_counts['tRNA'] < 20:
print(f'WARNING: Only {feature_counts["tRNA"]} tRNAs detected (expect >= 20)')
qc_pass = Falseprint(f'Feature summary: {feature_counts}')
return qc_pass, feature_counts
Path B: Eukaryotic Annotation
Step 1: Repeat Masking
# Build the RepeatModeler database FIRST, then the species-specific library
BuildDatabase -name mygenome assembly.fasta
RepeatModeler -database mygenome -threads 8 -LTRStruct
# CURATE the de novo library against a protein DB before masking, or real multi-copy gene# families (NLR/R-genes, zinc-fingers) get masked and "discovered" as a gene-poor repertoire.# Then soft-mask with the curated library (RepeatMasker uses the bundled Dfam DB in addition).
RepeatMasker \
-lib mygenome-families.fa \
-pa 8 \
-xsmall \
-gff \
-dir repeat_out \
assembly.fasta
Repeat Masking QC Checkpoint
defcheck_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))
withopen(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
ifnot 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
# BRAKER3 combines GeneMark-ETP, AUGUSTUS, and TSEBRA# Uses both RNA-seq and protein evidence for best results
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
# If only RNA-seq evidence available
braker.pl \
--genome=repeat_out/assembly.fasta.masked \
--bam=rnaseq_sorted.bam \
--softmasking \
--threads 8 \
--species=my_species \
--gff3
# If only protein evidence available (use OrthoDB proteins)
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 completeness on predicted proteins. Use the DEEPEST applicable clade dataset# (e.g. insecta_odb10 / embryophyta_odb10), NOT the shallow eukaryota_odb10.# The diagnostic that matters: compare this proteome BUSCO to a genome-mode BUSCO on the# same assembly -- a large gap means the predictor missed present genes (see genome-annotation/annotation-qc).
busco \
-i braker_out/braker.aa \
-l <clade>_odb10 \
-o busco_annotation \
-m proteins \
--cpu 8
defcheck_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 = 0withopen(braker_gff) as f:
for line in f:
if line.startswith('#'):
continue
feature = line.strip().split('\t')[2] iflen(line.strip().split('\t')) >= 3else''if feature == 'gene':
gene_count += 1elif feature == 'exon':
exon_count += 1
mean_exons = exon_count / gene_count if gene_count > 0else0withopen(busco_summary) as f:
for line in f:
if line.strip().startswith('C:'):
completeness = float(line.strip().split('C:')[1].split('%')[0])
break
issues = []
ifnot (expected_genes_range[0] <= gene_count <= expected_genes_range[1]):
issues.append(f'Gene count {gene_count} outside expected range {expected_genes_range}')
if completeness < 90:
issues.append(f'BUSCO completeness {completeness:.1f}% < 90%')
if mean_exons < 2:
issues.append(f'Mean exons/gene {mean_exons:.1f} is low for eukaryote')
print(f'Genes: {gene_count}, Mean exons/gene: {mean_exons:.1f}, BUSCO: {completeness:.1f}%')
returnlen(issues) == 0, issues
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.
