Case 2
Perfect match
Same start, same end, same strand and the same length as the aligned protein.
AMAP keepsBakta's call, named after the protein: Rv0002, dnaN.
A Nextflow pipeline for bacterial genomes
Bakta calls the genes. AMAP checks every one.
AMAP annotates assembled genomes with Bakta, aligns a trusted reference proteome onto them with miniprot, and settles every disagreement between the two: late starts, genes split in two, calls that are not genes, genes Bakta missed. Built for the Mycobacterium tuberculosis complex, where it changes about 30% of the genes, and for any species you have a reference for.
Reads assembled genomes. Writes, for each one, a corrected GFF3, GenBank and proteome with the tables behind every decision, and one report for the whole run.
M. tuberculosis CDC1551, its first 12.4 kb, against H37Rv. Point at a gene.
What AMAP keeps
Bakta's call
H37Rv protein, aligned
what AMAP keeps
same coordinates
coordinates differ
4calls agree with the reference
4genes get the reference's coordinates
3calls are dropped
Start here
From nothing to a corrected genome
Every tool AMAP calls lives in one container image, so Nextflow and a container engine are all you install.
-
Install Nextflow and a container engine
Nextflow needs Java 17 or newer, and conda brings both. Then Docker on a laptop, or Singularity or Apptainer on a cluster.
$ conda create -n nextflow-env -c conda-forge -c bioconda nextflow $ conda activate nextflow-env $ docker --version
On an Apple Silicon Mac the image runs under emulation: fine for the test below, slow for a real genome. Installation, in full.
-
Check the plumbing
The test profile annotates the H37Rv genome that ships inside the image, so it needs no data of your own. If it finishes without an error, the whole chain works.
$ nextflow run PathoGenOmics-Lab/AMAP -profile docker,testThe first run pulls the image: about 3.4 GB compressed, 12 GB on disk, once.
-
Run your genomes
A sample sheet with two columns and a real tab between them, a check that takes seconds, then the run. Pick what you have and copy the command.
samples.tsvsampleId fasta CDC1551 /data/cdc1551.fna Erdman /data/erdman.fna
Check it first$ nextflow run PathoGenOmics-Lab/AMAP \ -profile docker \ --tsv samples.tsv --validate_onlyWhere will it run?Which genomes?Anything else?First time? The tutorial runs one real genome from an empty folder, with the numbers to check yours against. Every option is in Parameters.
What AMAP decides
Each call, held up against the reference
Every Bakta call is compared with where the reference proteins align, on the same strand of the same contig. That puts it in one of four cases, and a reference gene with no call at all is a fifth. These are real genes of CDC1551: point at one for its coordinates.
Case 3
Imperfect match
The same coordinates, but the alignment's length differs: H37Rv's fhaA has 6 residues that CDC1551's lacks.
AMAP keepsBakta's call, named Rv0020c.
Case 4
Conflicting overlap
Bakta's dnaA starts 33 bp after the aligned protein does. Where the two disagree, the reference wins by default.
Bakta calls two short genes where H37Rv's Rv0007 aligns end to end. The reference is longer than either call, so it wins.
Bakta's call is 3,303 bp longer than the alignment, which stops before the end of its protein. A call longer by more than 550 bp wins over an alignment that is not complete.
AMAP keepsthe reference's coordinates, 1 to 1,524.one gene, 9,914 to 10,828.Bakta's call, 2,267,753 to 2,272,573, named Rv2024c.
Case 1
Potentially novel
Nothing from H37Rv aligns to these four calls: a stretch CDC1551 has and H37Rv lacks. None is a call Bakta is known to make wrongly, and none matches the H37Rv genome at 80% identity.
Nothing aligns here either, but the call is covered end to end by one that Bakta makes on H37Rv itself and that the curated annotation does not have.
AMAP keepsall four, under Bakta's locus tags.nothing: the call is dropped.
Missing
A reference gene Bakta missed
H37Rv has a 71-residue protein, Rv0031, between two of Bakta's calls, and Bakta called nothing there. No gene accounts for the alignment.
AMAP keepsthe gene, put back at the alignment: 33,526 to 33,738.
Then come the harder ones: a protein placed on its paralog's locus, proteins too small for miniprot, genes a contig end cuts, frameshifts and lost stops, and the real copies of a gene. How it works, step by step.
Does it work?
Checked against genomes whose answer is known
The accuracy benchmark runs the correction on seven genomes and scores it against an annotation AMAP did not make: the reference's own genes, carried onto the other genomes by a whole-genome alignment. Every gene a genome carries unchanged is found at exactly its coordinates, under its name. Of the genes that changed, which have no exact answer, 99.8 to 100% still sit where the alignment puts them, under their name.
Unchanged genes, at exactly their coordinatesChanged genes, found under their name
M. tuberculosis H37Rvagainst itself
H37Rv, another assemblyCP110619, against H37Rv
M. tuberculosis CDC1551against H37Rv
M. bovis AF2122/97against H37Rv
CDC1551 in 40 contigsa draft, against H37Rv
E. coli K-12 MG1655against itself
E. coli O157:H7 Sakaiagainst K-12
Scores of tests/benchmark, which anyone can rerun in the container; the truth comes from the references and minimap2, never from AMAP.
The pipeline
Where each step happens
Stations are the processes of main.nf; lines are what travels between them. Pick a station, or flip an option to see what a run leaves out.
Run with
Scroll sideways to follow the lines.
What you get
Files you can open, and the reasons behind them
Everything lands under results/, one folder per step. These answer most questions; Outputs has every file.
AMAP/{sample}/*-AMAP.gff3 · .gbk
The corrected annotation. Pseudogenes and genes a contig end cut are marked the way NCBI marks them.
AMAP/{sample}/*-AMAP.clean.faa
The proteins without pseudogenes or internal stops: what orthology, phylogenomics and dN/dS tools expect.
AMAP/{sample}/*-AMAPmastertable.tsv
One row per gene: coordinates, Bakta, eggNOG and, for the tuberculosis complex, Mycobrowser.
CDS_CASES/{sample}/
Every decision: a table per case, the reference genes Bakta missed, and the counts of each.
final_report.html
The whole run in one page that opens offline and travels by email.
pipeline_info/
The exact tool versions and parameters of the run, for your methods section.
Lines from real runs on CDC1551; the annotation comes from the genome cut into 40 contigs.
contig_4 AMAP CDS 164980 165406 . - 0 ID=L4_CDC1551d_Rv0354c;…;Dbxref=AMAP_ids:Rv0354c;pseudo=true;Note=frameshift or lost stop codon: its frame runs on 41 residues past the end of the reference protein contig_9 AMAP CDS 4160 4591 . - 0 ID=L4_CDC1551d_Rv1239c;…;Dbxref=AMAP_ids:Rv1239c;partial=true;end_range=4591,. contig_10 AMAP CDS 101478 102146 . + 0 ID=L4_CDC1551d_Rv1239c_2;…;Dbxref=AMAP_ids:Rv1239c;partial=true;end_range=102146,.
pseudo=true a pseudogene, with the reason in
Notepartial=true cut by the end of its contig, with the open end in end_range or start_range CDS complement(164980..165406)
/locus_tag="L4_CDC1551d_Rv0354c"
/db_xref="AMAP_ids:Rv0354c"
/pseudo=""
/note="frameshift or lost stop codon: its frame runs on 41
residues past the end of the reference protein"
CDS 101478..>102146
/locus_tag="L4_CDC1551d_Rv1239c_2"
/db_xref="AMAP_ids:Rv1239c"
/translation="VFPGFDALPEVLRPVARPQPPNAHPVAQPPAQALVDCGVYVCGQR…"
/pseudo no
/translation, and the reason in /note>102146 the 3' end runs off the contig (< marks a cut 5' end)$ column -t -s$'\t' results/CDS_CASES/summary/*-summary.tsv case2 2813 coordinates matched the reference exactly case3 24 coordinates matched, alignment length differed case4 1100 overlapped the reference but disagreed, resolved case1_initial 141 Bakta called it, the reference has nothing there rv_missing 97 the reference has a gene there, Bakta missed it total_cds 4077 what Bakta called before correction case1_removed 130 "novel" genes that were not really novel case1_kept 11 kept as genuinely absent from H37Rv
The tutorial's run on the whole CDC1551 genome. How to read it.
The report
The whole run on one page
This one is the example report: twelve made-up genomes across four lineages, built around a real CDC1551 run. Synthetic data, so read the layout, not the numbers. Click to explore.
The report is made for a large screen. This one is the example, on synthetic data.
Open it full screen ↗If something looks off
The questions people ask first
The sample sheet is rejected
It has to be separated by a real tab: spaces will not do, and this is the most common reason. Check with cat -t samples.tsv, where real tabs show up as ^I. Keep sample IDs to letters, digits and underscores.
Nextflow will not start: "unsupported class file version"
That is Java 8. Nextflow needs Java 17 or newer: conda install -c conda-forge openjdk=17.
A task dies with exit code 137
On macOS and Windows, Docker runs inside a virtual machine with a memory limit of its own, often 8 GB. Raise it in Docker Desktop under Settings, Resources, Memory, and keep --max_memory at or below it. What each step needs.
The first run seems stuck
The first run pulls a 12 GB image and builds the eggNOG database, which dominates the wall clock; both happen once. A long silence during the database build is not a failure. Point --eggnog_db somewhere you keep, and -resume picks a stopped run up where it left off.
Case 2 is at or near zero
Every gene disagreeing with the reference by a few bases means something systematic is wrong with the reference, not with 4000 genes. The usual cause is a reference annotation from another assembly than the genome given with it. What a healthy summary looks like.
Does it work for anything but tuberculosis?
Yes. Give it a trusted annotation and the genome it describes with --ref_annotation and --ref_genome, and turn off the two steps that are for the tuberculosis complex only. Using your own reference.
A task failed. Where do I look?
Nextflow prints the work directory of the failed task. Its .command.err is almost always the answer, and .command.sh is the exact script that ran. Troubleshooting, by error.
Open source, on GitHub
Read the code, follow the changes, ask a question or report a bug.