A tutorial on pangenome graphs in anvi'o

The purpose of this reproducible tutorial is to describe step-by-step how gene-centric pangenome graphs are computed, interactively visualized, and summarized for downstream analyses in anvi’o.

Throughout this tutorial we will use 29 genomes that belong to Undatipelagibacter (a.k.a. SAR11 subclade Ia.3.VI), a recently named genus in the family Pelagibacteraceae.

The tutorial will start with raw FASTA files for each genome, and end with a pangenome graph and associated flat-text summaries of genomic regions identified in them for downstream analyses. We will demonstrate how to

• Re-orient FASTA files using anvi-reorient-genomes,
• Turn reoriented FASTA files into contigs-db files using anvi-gen-contigs-database,
• Annotate the resulting contigs-dbs with useful information using anvi-run-hmms, anvi-run-ncbi-cogs, anvi-run-kegg-kofams, and anvi-run-pfams,
• Compute a pangenome from these contigs-dbs using anvi-pan-genome,
• Compute a pangenome graph using the pan-db using anvi-pan-genome-graph,
• Visualize the resulting pan-graph-db using anvi-display-pan-graph, and,
• Summarize the pangenome graph into flat-text files for downstream analyses using anvi-summarize.

Why pangenome graphs?

If you are familiar with the anvi’o pangenomics workflow, you already know how to identify which genes are shared across a set of genomes and which are unique. But conventional pangenome analyses collapse genes into unordered clusters, discarding their chromosomal context. A pangenome graph preserves the syntenic relationships between genes: it tells you not only which genes are shared, but where they sit relative to each other on the chromosome. This makes it possible to precisely delineate structurally conserved backbone regions from variable regions, and to study how genomic neighborhoods differ across organisms. In short, a conventional pangenome answers “what genes do these genomes have?”, while a pangenome graph answers “how are these genomes organized?”.

This tutorial serves two purposes. The first one is to familiarize you to how pangenome graphs are implemented in anvi’o and walk you through each step for you to be able to apply it to your own datsaets. The second one is to generate the input files for a reproducible bioinformatics workflow for Henoch et al (2026) paper that officially introduces the gene-centric pangenome graph toolkit in anvi’o and its application to Undatipelagibacter. If you are interested in our downstream analyses of pangenome graphs and the kinds of insights they offer, please consider also reading our paper “Synteny-aware microbial pangenomes reveal blueprints of genomic variation”.

You can follow this tutorial with your own FASTA files, or using the FASTA files that are included below. If you have any questions, you can send an email, or find @alexh4180 and/or @meren.anvio at

Setting up the work environment

You will need two things to go through this tutorial:

• A working anvi’o installation
• FASTA files for your genomes

The next two subsections will help you ensure both of these needs are satisfied appropriately before you begin with the rest of the tutorial.

Anvi’o setup

The software tools for pangenome graphs in anvi’o are quite new and are not yet available with any of the major anvi’o releases. However, you can follow the instructions at https://anvio.org/install/ to install the development version of anvi’o (a.k.a., anvio-dev) on your computer.

If you are seeing a similar output when you run the following command on your terminal, you are golden:

anvi-self-test -v

Anvi'o .......................................: eunice (v9-dev)
Python .......................................: 3.10.15

Profile database .............................: 42
Contigs database .............................: 24
Pan database .................................: 21
Pangraph database ............................: 4
Genome data storage ..........................: 8
Structure database ...........................: 2
Metabolic modules database ...................: 4
tRNA-seq database ............................: 2
Genes database ...............................: 6
Auxiliary data storage .......................: 2
Workflow configurations ......................: 5

You also need to make sure that you have all the necessary annotation databases accessible on your system since some of the steps below will require them. If you have just installed anvi’o, or if you had it before but never ran any of the following programs, please run them in your terminal for anvi’o to set things up for you to avoid issues later:

anvi-setup-ncbi-cogs
anvi-setup-pfams
anvi-setup-kegg-data

If you are not sure, just run them anyway since if everything is already setup, anvi’o will tell you something along the lines of “things already seem to be set up – use the --reset flag if you want to download and rebuild everything from scratch”.

Data setup

Whether you will be using your own FASTA files, or those that this tutorial will provide, it is a good idea to start fresh in a work directory. Let’s first create a work directory for this tutorial:

mkdir ~/pangraph-tutorial/

cd ~/pangraph-tutorial/

We will do everything in this -currently empty- directory.

Download the FASTA files:

curl -L https://cloud.uol.de/public.php/dav/files/jSFTXG3cSQMBjYX -o FASTA-FILES-ORIGINAL.tar.gz

Unpack the archive and get rid of the tar.gz file:

tar -xzf FASTA-FILES-ORIGINAL.tar.gz && rm FASTA-FILES-ORIGINAL.tar.gz

Now we have all the FASTA files in the newly unpacked FASTA-FILES-ORIGINAL directory, and the rest of the tutorial will make use of this structure. Just to make sure we are on the same page at this point, this is how our current work directory looks like:

.
└── FASTA-FILES-ORIGINAL
    ├── HIMB122.fa
    ├── HIMB140.fa
    ├── HIMB1488.fa
    ├── HIMB1491.fa
    ├── HIMB1493.fa
    ├── HIMB1506.fa
    ├── HIMB1507.fa
    ├── HIMB1513.fa
    ├── HIMB1518.fa
    ├── HIMB1526.fa
    ├── HIMB1552.fa
    ├── HIMB1556.fa
    ├── HIMB1559.fa
    ├── HIMB1573.fa
    ├── HIMB1577.fa
    ├── HIMB1593.fa
    ├── HIMB1597.fa
    ├── HIMB1611.fa
    ├── HIMB1631.fa
    ├── HIMB1636.fa
    ├── HIMB1641.fa
    ├── HIMB1662.fa
    ├── HIMB1685.fa
    ├── HIMB1701.fa
    ├── HIMB1702.fa
    ├── HIMB1723.fa
    ├── HIMB1758.fa
    ├── HIMB1765.fa
    └── HIMB1770.fa

2 directories, 29 files

If you are intending to follow this tutorial with your own FASTA files and if you would like the rest of the tutorial to work for you seamlessly, please consider putting your FASTA files in a FASTA-FILES-ORIGINAL directory, and keep in mind the following critical suggestions:

• Give your FASTA files a one-word names without any spaces or funny characters (you can safely use _ character, and everything in this tutorial should still work).
• Make sure your FASTA files have the extension .fa in lowercase, and not .fasta, .fna, or any other extension.
• We assume that your FASTA files for this analysis correspond to relatively closely related organisms. Pangenomes can be computed for various levels of taxonomy, and one can go as far as an entire phylum for some branches of life. But pangenome graphs are much more conservative. We advise you to remain below genus level for your initial analyses. It is difficult to give precise numbers that will work equally well throughout the tree of life, but a minimum gANI above 90% would be a good starting expectation. In the case of Undatipelagibacter the average ANI was around 95%, and it went as low as 93%.
• This tutorial assumes that all your FASTA files represent complete genomes with a single chromosome. Working with pangenome graphs with fragmented genomes is also possible, but since preserved and explicit gene synteny is one of the critical features of this workflow, at this stage we are reluctant to encourage such applications. That said, please feel free to get in touch with us so we can share our best practices and two cents with you for a useful analysis of your genomes.

If you are not sure whether every single FASTA file you have has a single contig, you can always run this command and expect to see only number 1 next to each file:

for i in FASTA-FILES-ORIGINAL/*fa; do echo "$i: $(grep -c '>' $i)"; done

If you are here, and if each of the statements above make sense, we are ready to start the actual analysis.

Reorienting genomes

This purpose of this step is to make sure all sequences start from an evolutionarily meaningful position. Because bacterial chromosomes are circular, the position where a sequencing or assembly project “breaks” the circle into a linear FASTA sequence is arbitrary. If two genomes start at completely different positions on the chromosome, the graph algorithm could misinterpret genuinely conserved regions as rearrangements (even though we have means to avoid that, it is better to be safe than sorry). By reorienting all genomes so they start at the same landmark such as the origin of replication or any other highly conserved region across all genomes), we ensure that the linear representations of these circular chromosomes are directly comparable before gene calling takes place.

Depending on how genomes are finalized, the starting point of each sequence may differ from one another. While graph building step can take care of rotations of genomes, it is best practice to align genomes first, which we can do with the anvi’o program anvi-reorient-genomes. It is imperative that this step is implemented before the genes are called in FASTA files, which typically takes place when FASTA files are turned into anvi’o contigs-db files.

We first need to create a fasta-txt file, which is a two-column TAB-delimited file that can be passed to anvi-reorient-genomes to let the program know where to find all FASTA files of interest. You can generate this file for our original FASTA files by running the following command in our work directory:

anvi-script-gen-fasta-txt --input-dir FASTA-FILES-ORIGINAL/ \
                          -o FASTA-FILES-ORIGINAL/fasta.txt

The resulting fasta.txt will look like this:

Now we can run the program to reorient all these genomes (which should take about 15 to 20 seconds on a laptop computer):

anvi-reorient-genomes --fasta-txt FASTA-FILES-ORIGINAL/fasta.txt \
                      --use-dnaa-for-reference-orientation \
                      --threads 8 \
                      -o FASTA-FILES-FINAL

This alignment step (which produces a ton of visualizations for you to enjoy) will store all FASTA files under a new directory, FASTA-FILES-FINAL. Here is a fun experiment to test if things worked as expected: running the following command to see the first 50 nucleotides of each genome in the directory FASTA-FILES-ORIGINAL,

for f in FASTA-FILES-ORIGINAL/*.fa; do awk '!/^>/ {print substr($0, 1, 50); exit}' "$f"; done

results in the following output (where each line is the first 50 nucleotides of one genome):

AAAAATGTCATTTTGGATGACAAACAAAATAAAAATAAAATTTTCTCAAA
TATTTTGGGTGACTATACGATTACCCTGAGAGATTTTTGTGTTAAATTAT
AGTTGATGTAAGCGCTTTAAAACCTGGTAACACAATTACAGTTTTATGGA
AAGAATTAGTTTTTACCGAAGGAAACTCAAGAGCAATAGATGATCAAAAA
GTTTTTGTATTTAAAATATATCCCTCTTGTTGAATTATTTTATCTAAGTC
AAAATAAATTGATTTTTAATTATTTTTAGTTTTGTAAAAAAATTTTGATT
TAACAAATAAAGCAACTTATTTTAAAAATGATGAAATCATTTTTACAGAA
AAAACATTGAAATATTTGCTGAAGAAATAATTTATGAAAAAAATAAAGAA
CTTAGTCATGGGGGGATTTTAAATAATTGTGAAGGTTCTATTTCACTTTT
ATTTACCACCACCTGCCGTTTCAAGTAAGGCACCAAATAAAACGAACAGA
AGCTTGTTTAAAGCTTTCTTCAGATTCTTTTAATTTATTTTGTTCATGCA
AATCTGAAGAAAGCTTTAAACAAGCTATATCATTAAAGCCAAACTATGCT
TTATTATAAATAAATTTATCAGCCTCTATAATAATTCCTTCATTATTCTT
AAATGGAAATGTAGTAATAGATGATAAAGTTAAAGGATTAAAATTATTTG
AATTAGTATACAAATAAATACTATTAATTACTTCTTTTTAAAATGTATTC
ATTTAATAAAATTAAATATTCAATATTGCTAATTTTATTTACTTTTTATT
AAATAAATGAATCAAGAGGTAACAATTCTTATATTAACTCATAAAAGTAA
ATATTAGGATAATAAAAAAAAGAGTAAATTAAATAATTATTTTTATTTTT
CTCATTTCTTTTAAATTCAAATTCTTTTGCAGTTATTATTTTTTCATCAA
GCAAAGTATTTAAATTTAAAATTAAAATGGAAGGGTACAGGAATCAATAC
TTAAAGCCAAACTATGCTGAAGCTTACAGTAACTTGGGTAATATACTACG
CCAAACTATGCTGAAGCTTACAGTAACTTGGGTAATACATTACAAGAACA
TTAAAGGATATAAATATTCTTTAAATAAAGATATATTTCAACTTATTCCA
TATTTGTGACCCTGGTGTGACACATATTTTCTATTTGAAGTGGTAAGGGG
TTAATTCACTATATTAAACTATTAAAATTTTCTTCGTTTTAGTGATAAAA
TTAGAAAACCGATGCTCTATCCAGCTGAGCTAAGGGCGCAATATTTGAAA
GAAATAATTGATAAAAAAACAAACAATCAATTATTCACAAACAAGATAAC
GTAAGTAGCTTGCTCTGAAAAAATTTTATATTTATTTATCTGATCTTCTA
TTATTGAAGTTCAAAAAAATTCTTACAAAGAATTAACTGAATTCTATGCT

In contrast, running this command to see the first 50 nucleotides of each genome in the directory FASTA-FILES-FINAL,

for f in FASTA-FILES-FINAL/*.fa; do awk '!/^>/ {print substr($0, 1, 50); exit}' "$f"; done

results in the following output:

ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTAAAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGTTTCAGAGGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGACGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAACAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGTTTCAGAGGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA
ATGAATAAATCTTTTACAAATAAAAATATTAATAATTTTGCTTCAGAAGA

Demonstrating that anvi-reorient-genomes did its job: the near-identical starting sequences across all 29 genomes correspond to the DnaA region near the origin of replication, confirming that all chromosomes are now anchored at the same evolutionarily conserved position. We are ready to move on.

Generating contigs-db files

curl -L https://cloud.uol.de/public.php/dav/files/7bRpYznDNBedSRk -o CONTIGS-DBs.tar.gz

tar -xzf CONTIGS-DBs.tar.gz && rm CONTIGS-DBs.tar.gz

If you are here, it means you have your final FASTA files in the FASTA-FILES-FINAL directory. Next, is to make sure we have an anvi’o contigs-db file for each of the FASTA files.

Generating contigs-db is quite straightforward. For a single FASTA file, there are three things one typically needs to do: (1) fix FASTA deflines, (2) generate a contigs-db, and (3) annotate it with functions. To apply all these steps to all FASTA files in the FASTA-FILES-FINAL directory, here we can use a for loop in our terminal, since loops make everything much more fun and dynamic.

While this tutorial explains the following steps by running each tool one after another, our group typically generates contigs-db files and annotates them using the program anvi-run-workflow with the ‘contigs’ workflow. The anvi-run-workflow uses snakemake in the background, and it allows us to run a large number of tasks in parallel on our HPC environment, making life much easier. But this strategy makes the FASTA to contigs-db journey like a black box, so here in this tutorial we go through each step one by one to make sure what is really happening behind the scenes is more obvious.

Here is the for loop that will turn each FASTA file under FASTA-FILES-FINAL into a contigs-db and store it under a new CONTIGS-DBs directory:

If you are following this tutorial word-by-word, you must run the script below from within ~/pangraph-tutorial/.

# first, let's get rid of any CONTIGS-DBs directory we may have in our work
# directory that may be left behind from a previous attempt
rm -rf CONTIGS-DBs

# create a fresh CONTIGS-DBs directory so we store everything in it
mkdir CONTIGS-DBs

# let's create a variable that keeps all our FASTA files we are
# interested in, so we can 'loop through' them next:
fasta_files=( FASTA-FILES-FINAL/*.fa )

# now we are ready to go through each FASTA file, turn them into
# anvi'o contigs-db files, and annotate them with functions
# and HMMs. This is the part that will take some time and it is
# best to run it on an HPC
num_threads="8"
for fasta_file in "${fasta_files[@]}"
do
    # learn the genome name without the file extension
    genome=$(basename $fasta_file .fa)

    # simplify the headers. this script will simply ensure
    # that the FASTA deflines are anvi'o compatible. since
    # the original deflines are in the matching FASTA files
    # under FASTA-FILES-ORIGINAL, we are not losing any
    # information here
    anvi-script-reformat-fasta $fasta_file \
                               --simplify-names \
                               --prefix $genome \
                               --overwrite-input

    # this anvi'o program will generate a contigs database
    # for this file. note that it is reading from one directory
    # and writing its output into another
    anvi-gen-contigs-database -f FASTA-FILES-FINAL/$genome.fa \
                              -o CONTIGS-DBs/$genome.db

    # the following programs will annotate the newly
    # generated contigs-db file with (1) the default HMMs
    # ship with anvi'o, (2) COGs from NCBI, (3) KOfams from
    # KEGG, (4) Pfams from .. well Pfams, and will also call
    # rRNA and tRNA genes:
    anvi-run-hmms -c CONTIGS-DBs/$genome.db -T $num_threads
    anvi-run-ncbi-cogs -c CONTIGS-DBs/$genome.db -T $num_threads
    anvi-run-kegg-kofams -c CONTIGS-DBs/$genome.db -T $num_threads
    anvi-run-pfams -c CONTIGS-DBs/$genome.db -T $num_threads
    anvi-scan-trnas -c CONTIGS-DBs/$genome.db -T $num_threads
done

If you are here, regardless of whether you downloaded your contigs-db files, or generated them from scratch, you should see something like this when you look into the CONTIGS-DBs directory:

CONTIGS-DBs/
├── HIMB122.db
├── HIMB140.db
├── HIMB1488.db
├── HIMB1491.db
├── HIMB1493.db
├── HIMB1506.db
├── HIMB1507.db
├── HIMB1513.db
├── HIMB1518.db
├── HIMB1526.db
├── HIMB1552.db
├── HIMB1556.db
├── HIMB1559.db
├── HIMB1573.db
├── HIMB1577.db
├── HIMB1593.db
├── HIMB1597.db
├── HIMB1611.db
├── HIMB1631.db
├── HIMB1636.db
├── HIMB1641.db
├── HIMB1662.db
├── HIMB1685.db
├── HIMB1701.db
├── HIMB1702.db
├── HIMB1723.db
├── HIMB1758.db
├── HIMB1765.db
└── HIMB1770.db

1 directory, 29 files

And you should se an output like this if you were to run anvi-db-info on any of the contigs-db files in that directory:

anvi-db-info CONTIGS-DBs/HIMB122.db

DB Info (no touch)
===============================================
Database Path ................................: HIMB122.db
description ..................................: [Not found, but it's OK]
db_type ......................................: contigs (variant: unknown)
version ......................................: 24


DB Info (no touch also)
===============================================
project_name .................................: HIMB122
contigs_db_hash ..............................: hash09d8a1e8
split_length .................................: 20000
kmer_size ....................................: 4
num_contigs ..................................: 1
total_length .................................: 1453515
num_splits ...................................: 72
gene_level_taxonomy_source ...................: None
genes_are_called .............................: 1
external_gene_calls ..........................: 0
external_gene_amino_acid_seqs ................: 0
skip_predict_frame ...........................: 0
splits_consider_gene_calls ...................: 1
scg_taxonomy_was_run .........................: 0
scg_taxonomy_database_version ................: None
trna_taxonomy_was_run ........................: 0
trna_taxonomy_database_version ...............: None
reaction_network_ko_annotations_hash .........: None
reaction_network_kegg_database_release .......: None
reaction_network_modelseed_database_sha ......: None
reaction_network_consensus_threshold .........: None
reaction_network_discard_ties ................: None
creation_date ................................: 1778571941.70956
modules_db_hash ..............................: 68221bd12b30
gene_function_sources ........................: KEGG_Module,KEGG_Class,COG24_FUNCTION,KEGG_BRITE,COG24_CATEGORY,Pfam,COG24_PATHWAY,KOfam,Transfer_RNAs

* Please remember that it is never a good idea to change these values. But in some
  cases it may be absolutely necessary to update something here, and a
  programmer may ask you to run this program and do it. But even then, you
  should be extremely careful.


AVAILABLE GENE CALLERS
===============================================
* 'pyrodigal-gv' (1,509 gene calls)
* 'Transfer_RNAs' (30 gene calls)
* 'Ribosomal_RNA_16S' (1 gene calls)


AVAILABLE FUNCTIONAL ANNOTATION SOURCES
===============================================
* COG24_CATEGORY (1,278 annotations)
* COG24_FUNCTION (1,278 annotations)
* COG24_PATHWAY (513 annotations)
* KEGG_BRITE (1,042 annotations)
* KEGG_Class (329 annotations)
* KEGG_Module (329 annotations)
* KOfam (1,043 annotations)
* Pfam (2,775 annotations)
* Transfer_RNAs (30 annotations)


AVAILABLE HMM SOURCES
===============================================
* 'Archaea_76' (76 models with 32 hits)
* 'Bacteria_71' (71 models with 70 hits)
* 'Protista_83' (83 models with 3 hits)
* 'Ribosomal_RNA_16S' (3 models with 1 hit)
* 'Ribosomal_RNA_18S' (1 model with 0 hits)
* 'Transfer_RNAs' (68 models with 30 hits)

If all looks good on your end, we can move on to the next step.

Computing the pangenome

curl -L https://cloud.uol.de/public.php/dav/files/TN2bxBCbAS5DRDJ -o UNDATIPELAGIBACTER-GENOMES.db

curl -L https://cloud.uol.de/public.php/dav/files/ctRp8xRWwaPSnp5 -o UNDATIPELAGIBACTER-PAN.db

Now we have all the contigs-db files ready, we are ready to compute a conventional pangenome. Which is a necessary step prior to generating a pangenome graph. The topic of pangenomics is covered extensively in another tutorial, so here we will employ a pragmatic approach and use only what we need from it.

Generating an anvi’o pan-db is also straightforward: (1) put genomes of interest in a genomes-storage-db and (2) compute a pangenome from it.

We will use the program anvi-gen-genomes-storage for the first step, and for that, we first need to describe which genomes we are interested in using an external-genomes artifact. Since we are interested in all genomes now under CONTIGS-DBs, generating a file that describes all of them is as simple as running the following command:

anvi-script-gen-genomes-file --input-dir CONTIGS-DBs/ \
                             -o external-genomes.txt

This will create a two-column TAB-delimited file as explained here. Now we can complete the first step of computing a pangenome, and generate a genomes-storage-db for ourselves:

anvi-gen-genomes-storage -e external-genomes.txt \
                         -o UNDATIPELAGIBACTER-GENOMES.db

Once it is ready, we can generate the pangenome:

anvi-pan-genome -g UNDATIPELAGIBACTER-GENOMES.db \
                --project-name UNDATIPELAGIBACTER \
                --num-threads 8

And visualize it using anvi-display-pan:

anvi-display-pan -p UNDATIPELAGIBACTER-PAN.db \
                 -g UNDATIPELAGIBACTER-GENOMES.db

Which will give us the following interactive interface to play with this pangenome:

Undatipelagibacter pangenome with 29 complete genomes.

This interactive display and the pan-db itself offer a wide range of downstream analysis opportunities. But we will not go through them here since what to do with pangenomes is extensively covered elsewhere.

Instead, we will quickly move on to the primary purpose of this tutorial: turning a pangenome into a pangenome graph.

Computing the pangenome graph

curl -L https://cloud.uol.de/public.php/dav/files/8eZZYqNrAdXF4TA -o UNDATIPELAGIBACTER-PAN-GRAPH.db

As our paper also mentions and we reiterated at the beginning of this tutorial, pangenomics is indeed an extremely powerful tool to characterize the gene repertoire across a set of genomes and to quantify conserved and variable gene sets among closely related members of a taxon. The purpose of anvi-pan-genome-graph is to recover syntenic context of genes, and make them accessible for downstream exploratory or quantitative analyses, by turning pangenomes into pangenome graphs:

anvi-pan-genome-graph -e external-genomes.txt \
                      -g UNDATIPELAGIBACTER-GENOMES.db \
                      -p UNDATIPELAGIBACTER-PAN.db

The output messages from the program will walk you through a series of clearly labeled stages, each announced in the terminal with a green header, so you can follow the progress and make sense of what is going on. Without going too much in technical detail, here is a brief summary of what happens when you run the program, and what success (or failure) looks like:

Stage 1: Start-node selection. After displaying your settings, the program determines the first node in the graph. It will either use the gene you pointed to with the flag --start-gene, or it will automatically select the synteny gene cluster with the smallest average gene caller ID across the most genomes. The latter is a great choice that works well when you have reoriented your genomes so they all start at the same position like we did here in this tutorial.

Stage 2: SynGC contextualization. This is the most computationally intensive step. The algorithm splits conventional gene clusters calculated by anvi-pan-genome and described in pan-db into synteny gene clusters (SynGCs) by examining local genomic neighborhoods. At the end of this step, the program will print a simple summary table that will show the total number of genes across all genomes, number of gene clusters before contextualization, number of SynGCs after, and a breakdown of how many genes fell into core, rearrangement, accessory, singleton, duplication, and RNA categories.

When everything works well given your genomes (as in, the evolutionary or architectural distance between your genomes make them suitable for this kind of an analysis), the analysis should produce a SynGC count that is only modestly larger than the gene cluster count. If the SynGC count exceeds twice the gene cluster count, the program will stop and tell you that your data has too many repeats or multi-copy genes forcing over-splitting of the graph context. If things will fail, they will most likely fail at this stage. There is not a lot of ways to ameliorate this issue, but there are some tools in the program power users (or those who are really determined to figure things out) can try:

• Widen the context window by increasing --min-k, or by lowering --alpha and/or --n parameters.
• Account for inversions by using the --inversion-aware flag if chromosomal inversions are common across the set of genomes considered.
• Remove offending genomes from your analysis if there are clear offenders. Poorly assembled or poorly scaffolded genomes that contain chimeras or unrealistic inversions will ruin the graph for various reasons and have no place in a high-resolution analysis of genome architectures.
• Strip gene-level offenders by removing multi-copy genes with --max-num-multi-copy-genes and --max-num-multi-copy-genes-per-genome parameters.
• Press the big red YOLO button by declaring the --just-do-it flag, and bypass all the graph integrity checks and let the algorithm attempt to continue and see if the graph will magically converge (it likely won’t).

Stage 3: Graph connectivity check. Disconnected fragments (which are common with draft genomes assembled in many contigs) will be reported during this step with their sizes and the algorithm will automatically prune them to keep only the largest connected components. If the graph is still fragmented after pruning, the run will abort with a message about severe data inconsistencies. This is the second place where things can go horribly wrong if you have fragmented genomes that are fragmented in unlucky ways.

Stage 4: Edge directionality optimization. The program will optimize the directionality of edges throughout the graph to produce a directed acyclic graph that is suitable for downstream layout calculations, and will let you know how many edges had to be reversed. There may be zero reversals, which would be unusual but valid, or you may have many many reversals, which would indicate substantial amount of structural complexity in your pangenome. Which is what we are here to study, so all is well.

Stage 5: Genome distances and output. A successful run will end with a genome-distance matrix calculated from the graph topological properties, and a PAN-GRAPH.db file written to your disk with a final message summarizing the number of regions identified and whether they were classified as backbone or variable regions.

Here is the very end of a successful run in our case with the 29 Undatipelagibacter genomes:

(...)

COMPUTING GRAPH-BASED GENOME DISTANCES
===============================================
(...)
* d(HIMB1506,HIMB1491) = 0.303
* d(HIMB1506,HIMB140) = 0.31
* d(HIMB1506,HIMB1493) = 0.299
* d(HIMB1506,HIMB1488) = 0.284
* d(HIMB1506,HIMB122) = 0.311
* d(HIMB1506,HIMB1507) = 0.292
* d(HIMB1491,HIMB140) = 0.23
* d(HIMB1491,HIMB1493) = 0.225
* d(HIMB1491,HIMB1488) = 0.317
* d(HIMB1491,HIMB122) = 0.257
* d(HIMB1491,HIMB1507) = 0.252
* d(HIMB140,HIMB1493) = 0.244
* d(HIMB140,HIMB1488) = 0.313
* d(HIMB140,HIMB122) = 0.237
* d(HIMB140,HIMB1507) = 0.235
* d(HIMB1493,HIMB1488) = 0.299
* d(HIMB1493,HIMB122) = 0.253
* d(HIMB1493,HIMB1507) = 0.258
* d(HIMB1488,HIMB122) = 0.311
* d(HIMB1488,HIMB1507) = 0.286
* d(HIMB122,HIMB1507) = 0.237
* Done.

STORING PAN-GRAPH-DB
==============================================================
Pan Graph database ..........................................: A new database, ./UNDATIPELAGIBACTER-PAN-GRAPH.db, has been created.

✓ pan_genome_graph.py took 0:02:48.2

If everything went well, we have our pan-graph-db ready for visualization.

Interactive visualization of the graph

Interactive visualization of a pangenome graph is done through the program anvi-display-pan-graph. When you run it on the files we have just generated,

anvi-display-pan-graph -p UNDATIPELAGIBACTER-PAN-GRAPH.db \
                       -g UNDATIPELAGIBACTER-GENOMES.db

Your browser should pop-up with an interactive display that shows you something like this when you press the ‘Draw’ button:

Undatipelagibacter pangenome graph with 29 complete genomes.

Before diving into the interface, it helps to understand the visual grammar of the graph. Each node represents a synteny gene cluster (SynGC), which describes a group of homologous genes that occupy the same syntenic position across one or more genomes. Edges connect SynGCs that are adjacent on at least one chromosome, so the paths through the graph trace the actual gene order of individual genomes. Where all genomes share the same gene neighborhood, the graph collapses into a single linear track. When genomes diverge (through gene INDELs or rearrangements), the graph branches into parallel paths that eventually/hopefully reconverge. The coloring distinguishes backbone regions (shown in blue by default), where all or nearly all genomes follow the same path, from variable regions (shown in yellow), where genomes take different routes through the graph. And this is exactly the kind of information a conventional pangenome cannot offer: the pangenome can tell you that a gene cluster is shared by 6 of 29 genomes, but only the pangenome graph can tell you that those 6 genomes carry it in a specific syntenic context between two backbone regions, forming a coherent variable island.

There are a lot of things one can do with the interactive interface, and it is difficult to cover all functionality it offers. One of the ways to explore the interface is through clicking on region IDs that are automatically assigned by anvi’o. In the screenshot above you already can see some regions numbered. How many of them are shown is a function of the size of the region and zoom level. By zooming into the graph with your mouse, you can see more of these region IDs to appear:

The effect of zooming in and zooming out on the interface.

By clicking on these region IDs, or by putting them in bins entirely or manually by selecting nodes while clicking Alt, you can learn about the metabolic modules they are involved in, and the function annotation information for each SynGC. Here is the example output for region #32:

Region #32 in the Undatipelagibacter pangenome graph and the metabolic involvement of SynGCs encoded in this region.

The same window as above, scrolled down to show the functions per SynGC section of it.

You can also arbitrarily add regions into any number of bins,

Some backbone regions added to a bin.

And show metabolic involvement of these SynGCs or their functional annotations by clicking the nodes button on the Bins panel:

Some exploratory insights into a bin.

You can also left-click on individual nodes to learn more about them, and see the amino acid sequence alignments among the genes that contribute to them:

Clicking on nodes.

The interactive interface is still evolving, and will tremendously benefit from any input you may have to share with us.

Summarizing the graph

curl -L https://cloud.uol.de/public.php/dav/files/8snz92oDqJeARDK -o UNDATIPELAGIBACTER-PAN-GRAPH-SUMMARY.tar.gz

tar -xzf UNDATIPELAGIBACTER-PAN-GRAPH-SUMMARY.tar.gz && rm UNDATIPELAGIBACTER-PAN-GRAPH-SUMMARY.tar.gz

Browsing the graph and gaining metabolic and functional insights into backbone and variable regions is one of the most important activities one can do with pangenome graphs to wrap their minds around the nature of the genomes they are studying. But beyond these exploratory and largely qualitative insights, most research questions will also need text based summaries of this display so they can be fed into LLMs, or statistical analysis or visualization tools to put numbers and/or significance values on various claims.

To generate text summaries of a given pan-graph-db, one can use the program anvi-summarize. When you run it without any collection,

anvi-summarize -p UNDATIPELAGIBACTER-PAN-GRAPH.db \
               -g UNDATIPELAGIBACTER-GENOMES.db \
               -o UNDATIPELAGIBACTER-PAN-GRAPH-SUMMARY

This will generate an output directory with multiple output files.

REGIONS.txt

One of the most powerful insights pangenome graphs offer is the precise delineation of backbone and variable regions among a given set of genomes, what is shown in blue and yellow in default in the interactive display of the graph.

The REGIONS.txt file simply describe these regions by giving each independent region a unique identifier.

For instance, it is obvious from the display that one of the most variable regions in the Undatipelagibacter pangenome graph is the region #32. Using this file, one can use the quantitative estimates about the regions anvi’o have identified in the graph to reach the same conclusion:

Disclaimer, the following BASH one-liner is a copy-paste from ChatGPT answering the prompt: “Read the attached REGIONS.txt file, a documentation for which is available at https://anvio.org/help/main/artifacts/pan-graph-summary/, and give me a single BASH one-liner that will recover the most variable region ID in this pangenome graph”.

awk -F'\t' 'NR>1 && $16>max {max=$16; id=$1} END{print id}' UNDATIPELAGIBACTER-PAN-GRAPH-SUMMARY/REGIONS.txt

When we run it in the terminal, this is what I get:

32

Excellent. Or we could ask it to give us a BASH one-liner that finds the variable region with the lowest max_expansion and highest complexity to identify ‘tower’ like structures:

awk -F'\t' 'NR>1 && $2=="variable"' UNDATIPELAGIBACTER-PAN-GRAPH-SUMMARY/REGIONS.txt | sort -t$'\t' -k8,8n -k10,10rn | head -1 | awk '{print $1}'

When I run this in my terminal, this is what I get:

148

Which is indeed a very intersting region. One can go much further than this, of course using additional files, such as SYNGCs.txt.

SYNGCs.txt

Regions are made of SynGCs, and each SynGC is described in the SYNGCs.txt file, which looks like this:

Please note that there are many columns that may be hidded as a function of the resolution of your screen, but you can view them by sliding the display below towards right.

Combining this file with the previous one, one can do a lot through the command line tools, Python, or R.

Because the summary files combine both structural information (region assignments, variability scores, genome membership) and functional annotations (COG, KEGG, Pfam) in flat tables, they are also well-suited for natural-language querying through large language models.

Just to demonstrate how these files are linked together and how accessible the output is, let’s upload both files to Claude and ask the following question:

The two file attached, REGIONS.txt and SYNGCs.txt are produced by the summary of a pangenome graph generated by anvi’o. The column names and their meanings are explained in this anvi’o artifact you should read:

https://anvio.org/help/main/artifacts/pan-graph-summary/

Please find the most variable region in this pangenome graph, and describe the functions encoded in this region by solely relying on the data in these files and very briefly describe the likely role of this region in the lifestyle of this organism.

Here is Claude’s response:

The most variable region appears to be the region 32.

Out of 341 regions in this pangenome graph (170 variable, 171 backbone), region 32 stands out dramatically with a composite variability score of 0.994 (on a 0–1 scale) — well above the runner-up, region 90, at 0.861. Here are its key statistics:

Metric	Value
Composite variability score	0.994 (next highest: 0.861)
Synteny gene clusters	1,005 (721 found in only 1 of 29 genomes)
Dominant COG category	Cell wall/membrane/envelope biogenesis (53%)
Key functions	Sugar-nucleotide biosynthesis, glycosyltransferases, O-antigen ligases, LPS export, sialic acid-like sugar pathways

This region appears to encode genome-specific surface polysaccharide biosynthesis machinery.

You may recall that region #32 was the same region we visually explored in the interactive display above. The LLM arrived at the same conclusion independently from the flat-text summaries alone, correctly identifying the surface polysaccharide biosynthesis island that is also visually striking in the graph. Something that is visible, is also quantifiable.

Final words

We hope our tutorial was useful for you to gain some insights into how to run pangenome graphs on your own datasets.

Please feel free to let us know if you have any questions.

Questions? Concerns? Find us on

If you are interested, you can continue with the reproducible bioinformatics workflow behind the Henoch et al. paper that starts with the data products generated in this tutorial, and the downstream analyses that led to the observations we described in our work:

https://merenlab.org/data/undatipelagibacter-pangenome-graph

Thank you for your time and interest.