An anvi'o workflow for microbial pangenomics

The latest version of anvi’o is v7.1. See the release notes.

If you find a mistake on this page or would you like to update something in it, please feel free to edit its source by clicking the edit button at the top-right corner (which you will see if you are logged in to GitHub) 😇

With the anvi’o pangenomic workflow you can,

• Identify gene clusters across your genomes, by providing anvi’o with two or more FASTA files,
• Combine metagenome-assembled genomes from your anvi’o projects genomes with genomes from other sources.
• Interactively visualize your pangenomes,
• Estimate relationships between your genomes based on gene clusters,
• Interactively (or programmatically) partition gene clusters into bins,
• Perform phylogenomic analyses on-the-fly given a set of gene clusters,
• Annotate your genes, and inspect amino acid alignments within your gene clusters,
• Extend your pangenome with contextual information about your genomes or gene clusters,
• Quantify geometric and biochemical homogeneity of your gene clusters,
• Perform functional enrichment analyses on groups of genomes in your pangenome,
• Compute and visualize average nucleotide identity scores between you genomes, and more.

You can use anvi’o for pangenomic workflow even if you haven’t done any metagenomic work with anvi’o. All you need is an anvi’o installation, and a FASTA file for each of your genomes.

Introduction

The anvi’o pangenomic workflow described here will walk you through the following steps:

• Generate an anvi’o genomes-storage-db using the program anvi-gen-genomes-storage.

• Generate an anvi’o pan-db using the program anvi-pan-genome (which requires a genomes storage)

• Display results in anvi’o interactive interface using the program anvi-display-pan.

You can then use the interactive interface to bin your gene clusters into collections, or use the program anvi-import-collection to import bins for your gene clusters, and finally you can use the program anvi-summarize to create a static HTML summary of your pangenome. Easy peasy!

The following sections will detail each step, and culminating in an example run will follow, but let’s first make sure you have all the required dependencies installed and your installation is good to go.

Dependencies

If your system is properly setup, this anvi-self-test command should run without any errors:

$ anvi-self-test --suite pangenomics

Generating an anvi’o genomes storage

A genomes storage is a special anvi’o database that stores information about genomes. A genomes storage can be generated only from external-genomes, only from internal-genomes, or a combination of both. Before we go any further, here are some definitions to clarify things:

• An external genome is anything you have in a FASTA file format (i.e., a genome you have downloaded from NCBI, or obtained through any other way). Which means, you will need to convert each of your FASTA file into an anvi’o contigs database first using the program anvi-gen-contigs-database. Please read the contigs-db artifact to make sure you populate your contigs database with most useful information (such as annotating your genes with functions, and so on).

• An internal genome is any genome bin you stored in an anvi’o collection (after binning and/or refining genomes from metagenomes in anvi’o). Working with internal genomes is quite straightforward since you already have an anvi’o contigs and an anvi’o profile database for them, but don’t worry if you are reading this tutorial and this does not yet make sense to you.

You can create a new anvi’o genomes storage using the program anvi-gen-genomes-storage, which will require you to provide descriptions of genomes to be included in this storage. File formats for external genome and internal genome descriptions differ slightly. For instance, this is an example --external-genomes file:

and this is an example file for --internal-genomes:

For names in the first column please use only letters, digits, and the underscore character.

Thanks to these two files, genomes described in anvi’o collections and genomes obtained via other sources can be combined and analyzed together seamlessly. The most essential need for the coherence within the genomes storage is to make sure each internal and external genome is generated identically with respect to how genes were called, how functions were assigned, etc. Anvi’o will check for most things, but it can’t stop you from doing mistakes. For instance, if the gene caller that identified open reading frames is not identical across all contigs databases, the genes described in genomes storage will not necessarily be comparable. If you are not sure about something, send us an e-mail, and we will be happy to try to clarify.

wget https://ndownloader.figshare.com/files/28834476 -O Prochlorococcus_31_genomes.tar.gz
tar -zxvf Prochlorococcus_31_genomes.tar.gz
cd Prochlorococcus_31_genomes
anvi-migrate *.db --migrate-dbs-safely

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

Running a pangenome analysis

Once you have your genomes storage ready, you can use the program anvi-pan-genome to run the actual pangenomic analysis. This is the simplest form of this command:

$ anvi-pan-genome -g MY-GENOMES.db -n PROJECT_NAME

anvi-pan-genome -g PROCHLORO-GENOMES.db \
                --project-name "Prochlorococcus_Pan" \
                --output-dir PROCHLORO \
                --num-threads 6 \
                --minbit 0.5 \
                --mcl-inflation 10 \
                --use-ncbi-blast

anvi-import-misc-data layer-additional-data.txt \
                      -p PROCHLORO/Prochlorococcus_Pan-PAN.db \
                      --target-data-table layers

New layers additional data...
===============================================
Data key "clade" .............................: Predicted type: str
Data key "light" .............................: Predicted type: str

New data added to the db for your layers .....: clade, light.

When you run anvi-pan-genome, the program will,

• Use all genomes in the genomes storage. If you would like to focus on a subset, you can use the parameter --genome-names.

• Use only a single core by default. Depending on the number of genomes you are analyzing, this process can be very time consuming, hence you should consider increasing the number of threads to be used via the parameter --num-threads.

• Use DIAMOND (Buchnfink et al., 2015) in ‘fast’ mode by default (or you can ask DIAMOND to be ‘sensitive’ by using the flag --sensitive) to calculate the similarity of each amino acid sequence in every genome against every other amino acid sequence across all genomes (which clearly requires you to have DIAMOND installed). Alternatively you could use the flag --use-ncbi-blast to use NCBI’s blastp for amino acid sequence similarity search.

• Use every gene call, whether they are complete or not. Although this is not a big concern for complete genomes, metagenome-assembled genomes (MAGs) will have many incomplete gene calls at the end and at the beginning of contigs. Our experiments so far suggest that they do not cause major issues, but if you want to exclude them, you can use the --exclude-partial-gene-calls flag.

• Use the minbit heuristic that was originally implemented in ITEP (Benedict et al, 2014) to eliminate weak matches between two amino acid sequences. You see, the pangenomic workflow first identifies amino acid sequences that are somewhat similar by doing similarity searches, and then resolves gene clusters based on those similarities. In this scenario, weak similarities can connect gene clusters that should not be connected. Although the network partitioning algorithm can recover from these weak connections, it is always better to eliminate as much noise as possible at every step. So the minbit heuristic provides a mean to set a to eliminate weak matches between two amino acid sequences. We learned it from ITEP (Benedict MN et al, doi:10.1186/1471-2164-15-8), which is another comprehensive analysis workflow for pangenomes, and decided to use it because it makes a lot of sense. Briefly, If you have two amino acid sequences, A and B, the minbit is defined as BITSCORE(A, B) / MIN(BITSCORE(A, A), BITSCORE(B, B)). So the minbit score between two sequences goes to 1.0 if they are very similar over the entire length of the ‘shorter’ amino acid sequence, and goes to 0.0 if (1) they match over a very short stretch compared even to the length of the shorter amino acid sequence or (2) the match between sequence identity is low. The default minbit is 0.5, but you can change it using the parameter --minbit.

• Use the MCL algorithm (van Dongen and Abreu-Goodger, 2012) to identify clusters in amino acid sequence similarity search results. We use 2 as the MCL inflation parameter by default. This parameter defines the sensitivity of the algorithm during the identification of the gene clusters. More sensitivity means more clusters, but of course more clusters does not mean better inference of evolutionary relationships. More information on this parameter and it’s effect on cluster granularity is here http://micans.org/mcl/man/mclfaq.html#faq7.2, but clearly, we, the metagenomics people will need to talk much more about this. So far in the Meren Lab we have been using 2 if we are comparing many distantly related genomes (i.e., genomes classify into different families or farther), and 10 if we are comparing very closely related genomes (i.e., ‘strains’ of the same ‘species’ (based whatever definition of these terms you fancy)). You can change it using the parameter --mcl-inflation. Please experiment yourself, and consider reporting!

• Utilize every gene cluster, even if they occur in only one genome in your analyses. Of course, the importance of singletons or doubletons will depend on the number of genomes in your analysis, or the question you have in mind. However, if you would like to define a cut-off, you can use the parameter --min-occurrence, which is 1, by default. Increasing this cut-off will improve the clustering speed and make the visualization much more manageable, but again, this parameter should be considered in the context of each study.

• Use euclidean distance and ward linkage to organize gene clusters and genomes. You can change those using --distance and --linkage parameters.

• Try to utilize previous search results if there is already a directory. This way you can play with the --minbit, --mcl-inflation, or --min-occurrence parameters without having to re-do the amino acid sequence search to gain time. If you are having hard time making anvi’o use the existing output files, please see the explanation in this issue. If you are not interested in using the previous search results, you can remove the output directory, or use the --overwrite-output-destinations flag to redo the search from scratch.

You need another parameter? Well, of course you do! Let us know, and let’s have a discussion. We love parameters.

Once you are done, a new directory with your analysis results will appear. You can add or remove additional data items into your pan profile database using anvi’o additional data tables subsystem.

Displaying the pan genome

Once your analysis is done, you will use the program anvi-display-pan to display your results.

This is the simplest form of this command:

$ anvi-display-pan -p PROJECT-PAN.db -g PROJECT-PAN-GENOMES.db

The program anvi-display-pan is very similar to the program anvi-interactive, and the interface that will welcome you is nothing but the standard anvi’o interactive interface with slight adjustments for pangenomic analyses. Of course anvi-display-pan will allow you to set the IP address and port number to serve, add additional view data, additional layers, and/or additional trees, and more. Please familiarize yourself with it by running anvi-display-pan -h in your terminal.

Here is the pangenome for the 31 Prochlorococcus isolate genomes we have created in the previous sections of this tutorial:

anvi-display-pan -g PROCHLORO-GENOMES.db \
                 -p PROCHLORO/Prochlorococcus_Pan-PAN.db

Looks ugly. But that’s OK (for now). For instance, to improve things a little, you may like to organize your genomes based on gene clustering results by selecting the ‘gene_cluster frequencies’ tree from the Samples Tab > Sample Order menu:

This is what happens when you draw it again (note the tree that appears on the right):

Looks more meaningful .. but still ugly.

Well, this is exactly where you need to start using the interface more efficiently. For instance, this is the same pangenome after some changes using the additional settings items in the settings panel of the interactive interface:

If you downloaded the Prochlorococcus data pack, you will also find an anvi’o state file that you can import into your pan database using the program anvi-import-state:

anvi-import-state -p PROCHLORO/Prochlorococcus_Pan-PAN.db \
                  --state pan-state.json \
                  --name default

No excuses for bad looking pangenomes.

For an example, step-by-step instructions on how to improve a pangenome, see the polishing the pangenome section in this reproducible workflow.

Inspecting gene clusters

Every gene cluster in your analysis will contain one or more amino acid sequence that originate from one or more genomes. While there will likely be a ‘core’ section, in which all gene cluster will appear in every genome, it is also common to find gene clusters that contain more than one gene call from a single genome (i.e., all multi-copy genes in a given genome will end up in the same gene cluster). Sooner or later you will start getting curious about some of the gene clusters, and want to learn more about them. Luckily you can right-click on to any gene cluster, and you would see this menu (or maybe even more depending on when you are reading this article):

For instance, if you click ‘Inspect gene cluster’, you will see all the amino acid sequences from each genome that went into that gene cluster (with the same order and background colors of genomes as they are arranged in the main display):

You can learn more about the amino acid color coding algorithm here.

It is not only fun but also very educational to go through gene clusters one by one. Fine. But what do you do if you want to make sense of large selections?

As you already know, the anvi’o interactive interface allows you to make selections from the tree. So you can select groups of gene clusters into bins (don’t mind the numbers on the left panel, there clearly is a bug, and will be fixed in your version):

You can create multiple bins with multiple selections, and even give them meaningful names if you fancy:

While selecting gene clusters manually using the dendrogram is an option, it is also possible to identify them using the search interface that allows you to define very specific search criteria:

You can highlight these selections in the interface, or you can add them to a collection to summarize them later.

In addition, you can search gene clusters also based on functions:

Similarly, you can add these gene clusters into collections with whatever name you like, and summarize those collections later.

Advanced access to gene clusters is also possible through the command line through the program anvi-get-sequences-for-gene-clusters. For more information see this issue or this vignette.

Inferring the homogeneity of gene clusters

This functionality is available since anvi’o v5.2 thanks to the efforts of Mahmoud Yousef.

Gene clusters are good, but not all gene clusters are created equal. By simply inspecting the alignments within just a few of gene clusters, you can witness differing levels of disagreements between amino acid sequences across different genomes.

Concept of homogeneity

A gene cluster may contain amino acid sequences from different genomes that are almost identical, which would be a highly homogeneous gene cluster. Another gene cluster may contain highly divergent amino acid sequences from different genomes, which would then be a highly non-homogeneous one, and so on.

One could infer the nature of sequence homogeneity within a gene cluster by focusing on two primary attributes of sequence alignments: functional homogeneity (i.e., how conserved aligned amino acid residues across genes), and geometric homogeneity (i.e., how does the gap / residue distribution look like within a gene cluster regardless of the identity of amino acids). While it is rather straightforward to have an idea about the homogeneity of gene clusters through the manual inspection of the aligned sequences within them, it has not been possible to quantify this information automatically. But anvi’o now has you covered.

Indeed, understanding the within gene cluster homogeneity could yield detailed ecological or evolutionary insights regarding the forces that shape the genomic context across closely related taxa, or help scrutinize gene clusters further for downstream analyses manually or programmatically. The purpose of this section is to show you how we solved this problem in anvi’o and to demonstrate its efficacy.

Functional and geometric homogeneity estimates in anvi’o

Anvi’o pangenomes contains two layers that summarize homogeneity indices for each gene cluster (and an additional one that combines the two).

Here is an example in the context of our Prochlorococcus pangenome (see the outermost two additional layers):

Homogeneity indices are calculated based on the alignment. If the alignment for some reason fails for a given set of genes in a gene cluster, the homogeneity indices for it will show -1.

Geometric homogeneity index

The geometric homogeneity index indicates the degree of geometric configuration between the genes of a gene cluster. Specifically, we look for the gap/residue patterns of the gene cluster that have been determined by the alignment. The alignment process by definition aligns similar sections of genes to each other, and denotes missing residues (or different structural configurations of gene contexts) through gap characters. If gap/residue distribution patterns are mostly uniform throughout the gene cluster, then this gene cluster will have a high geometric homogeneity, and the maximum value of 1.0 indicates that are no gaps in the alignment.

Anvi’o computes the geometric homogeneity index by combining analysis of the gene cluster content in two levels: the site-level analyses (i.e., vertically aligned positions) and the gene-level analyses (i.e., horizontally aligned positions because they are in the same gene). We convert the information in a gene cluster into a binary matrix, where gaps and residues are simply represented by 1s and 0s, and we utilize the logical operator xor to identify and enumerate all comparisons that have a different pattern. In site-level homogeneity checks, the algorithm sweeps from left to right and identifies these differences across aligned columns. In gene-level homogeneity checks, the algorithm sweeps from one gene to the next and identifies all differences in a gene where the pattern is not the same, effectively checking for the spread of gaps across the gene cluster. It is possible to skip the gene-level geometric homogeneity caluclations with the flag --quick-homogeneity. While this will not be as accurate or as comprehensive as the default approach, it may save you some time, depending on the number of genomes with which you are working.

Functional homogeneity index

In contrast, the functional homogeneity index considers aligned residues (by ignoring gaps), and attempts to quantify differences across residues in a site by considering the biochemical properties of differing residues (which could affect the functional conservation of genes at the protein-level). To do this, we divided amino acids into seven different “conserved groups” based on their biochemical properties. These groups are: nonpolar, aromatic, polar uncharged, both polar and nonpolar characteristics (these amino acids also lie in groups besides this one), acidic, basic, and mostly nonpolar (but contain some polar characteristics) (for more information please see this article). Then, the algorithm goes through the entire gene cluster and assigns a similarity score between 0 and 3 to every pair of amino acids at the same position across genes based on how close the biochemical properties of the amino acid residues are to each other. The sum of all assigned similarity scores is indicative of the functional homogeneity index of the gene cluster, and will reach to its maximum value of 1.0 if all residues are identical.

Both indices are on a scale of 0 to 1, where 1 is completely homogeneous and 0 is completely heterogeneous. If the algorithm is interrupted by a runtime error (due to unexpected issues such as not all genes being the same length for whatever reason, etc.), it will default to error values of -1 for the indices. So if you see a -1 in your summary outputs, it means we failed to make sense of the alignments in that gene cluster for some reason :/

In reality, the complex processes that influence protein folding and the intricate chemical interactions between amino acid residues should remind us that these assessments of similarity are only mere numerical suggestions, and do not necessarily reflect accurate biochemical insights, which is not the goal of these homogeneity indices.

Using homogeneity indices in pangenomes

Given all of that, let’s look at our pangenome again. What are they good for? Well, there is a lot one could do with these homogeneity indices, and it would be unfair to expect this tutorial to cover all of them. Nevertheless, here is an attempt to highlight some aspects of it.

To get a quick idea about some of the least homogeneous gene clusters, one could order gene clusetrs based on increasing or decreasing homogeneity. Let’s say we want to order it by increasing functional homogeneity (as you likely know already, this can be done through the “Items order” combo box in the main settings panel):

If you open the “Mouse” panel on the right side of your display, you can actually hover over the layer and see the exact value of the index for all gene clusters. The gene cluster shown underneath the cursor in the figure above has a relatively low functional homogeneity estimate, but a higher geometric homogeneity. We can inspect the gene cluster to take a look for ourselves:

You can see why the relatively higher geometric homogeneity score makes sense. Three of these genes have the same gap/residue pattern, and the gaps at the end of the other gene throw things off slightly, making the geometric score close to 0.75. On the other hand, the color coding of the aligned amino acids also gives us a hint regarding the lack of functional homogeneity across them. We can do the same for the geometric homogeneity index. Try it yourself: order the pangenome based on the geometric homogeneity index, and inspect a gene cluster with a relatively low score.

So, what can we do for more in depth analyses of our gene clusters?

Anvi’o offers quite a powerful way to filter gene clusters both through the command line program anvi-get-sequences-for-gene-clusters, and through the interface for interactive exploration of the data:

Exploratory analyses

Let’s assume you wish to find a gene cluster that represents a single-copy core gene with very high discrepancy between its geometric homogeneity versus its functional homogeneity. As in, you want something that is highly conserved across all genomes, it is under structural constraints that keeps its alignment homogeneous, but it has lots of room to diversify in a way that impacts its functional heterogeneity. You want a lot. But could anvi’o deliver?

Well, for this very specific set of constraints you could first order all gene clusters based on decreasing geometric homogeneity index, then enter the following values to set up a filter before applying it and highlighting the matching gene clusters:

The first gene cluster in counter-clockwise order is the one that best matches to the listed criteria. Upon closer inspection,

One could see that there is tremendous amount of functional variability among aligned residues across genomes, despite the geometric homogeneity of the sequences in the gene cluster (only a portion of the alignment is shown here):

Those of you who read our study on the topic, would perhaps not be surprised to learn that the COG functional annotation of this particular gene cluster resolves to an enzyme that is related to cell wall/membrane/envelop biogenesis. It is always good when things check out.

Scrutinizing phylogenomics

Here is another example usage of homogeneity indices. We often use single-copy core gene clusters for phylogenomic analyses to estimate evolutionary relationships between our genomes. Identifying single-copy core gene clusters is easy either through advanced filters, or through manual binning of gene clusters:

But single-copy core gene clusters will contain many poorly aligned gene clusters that you may not want to use for a phylogenomic analysis so as to minimize the influence of noise that stems from bioinformatics decisions regarding where to place gaps. On the other hand, there will be many gene clusters that are near-identical in this collection, which would be rather useless to infer phylogenomic relationships. Luckily, you could use homogeneity indices and advanced search options to identify those that are geometrically perfect, but functionally diverse:

You could easily append those that match to these criteria to a separate bin:

And perform a quick-and dirty analysis on-the-fly to see how they would organize your genomes directly on the interface:

Or you could get the FASTA file with aligned and concatenated amino acid sequences for these gene clusters to do a more appropriate phylogenomic analysis:

$ anvi-get-sequences-for-gene-clusters -p PROCHLORO/Prochlorococcus_Pan-PAN.db \
                                       -g PROCHLORO-GENOMES.db \
                                       -C default -b Better_core \
                                       --concatenate-gene-clusters \
                                       -o better_core.fa

Needless to say, estimates for homogeneity indices per gene cluster will also appear in your summary files from anvi-summarize to satisfy the statistician within you.

Splitting the pangenome

In some cases one might want to split a given pangenome into multiple independent pangenomes, such as one that contains only core gene clusters, or one that contains only singletons, etc. This advanced functionality is most useful when one wishes to do very high-resolution phylogenomic analyses by scrutinizing core genes in the pangenome with an orchestrated use of functional and geometric homogeneity indices (which can also be done without splitting the pangenome, but we find that it sometimes it makes things easier).

For this we have anvi-split, a multi-talented program to split things, which also works with collections and bins in pan databases. Using this program you can focus on any group of gene clusters that are defined in a bin in a given pangenome, and split them into an independent and stand-alone pangenome.

If you are feeling lost, you will likely find the visual description of this functionality much more clear than the technical one. Following steps will demonstrate it using the Prochlorococcus pangenome. Let’s assume in the Prochlorococcus pangenome you happened to have three bins, Core, HL Core, and Singletons, all stored in the collection named default:

These bins can be split from this pangenome into their own little pan databases using anvi-split this way:

anvi-split -p Prochlorococcus-PAN.db \
           -g Prochlorococcus-GENOMES.db \
           -C default \
           -o SPLIT_PANs

The resulting directory would contain the following files:

$ ls -lR SPLIT_PANs
total 0
drwxr-xr-x  3 meren  staff  96 Apr 26 16:37 Core
drwxr-xr-x  3 meren  staff  96 Apr 26 16:37 Hl_Core
drwxr-xr-x  3 meren  staff  96 Apr 26 16:37 Singletons

SPLIT_PANs/Core:
total 3240
-rw-r--r--  1 meren  staff  1658880 Apr 26 16:37 PAN.db

SPLIT_PANs/Hl_Core:
total 1560
-rw-r--r--  1 meren  staff  798720 Apr 26 16:37 PAN.db

SPLIT_PANs/Singletons:
total 4672
-rw-r--r--  1 meren  staff  1384448 Apr 26 16:37 PAN.db

As you can see these are individual pangenomes that can indeed be visualized with anvi-display-pan. For instance running the following command,

anvi-display-pan -p SPLIT_PANs/Core/PAN.db \
                 -g Prochlorococcus-GENOMES.db

Would give you this display:

Or runnig this one instead,

anvi-display-pan -p SPLIT_PANs/Singletons/PAN.db \
                 -g Prochlorococcus-GENOMES.db

Would give you this one:

And that’s that.

Quantifying functional enrichment in a pangenome

Once we have our pangenome, one of the critical things that we usually want to do is look at the functions associated with our gene clusters, especially those that are accessory to subsets of genomes in our pangenome. Anvi’o can help you identify functions that are enriched for some of the clades or sub clades that are included in your pangenome via the program This is done with our new and improved program anvi-compute-functional-enrichment-in-pan.

This program utilizes information you provide by splitting your genoems into multiple grups, and finds functions that are characteristic to any of those groups, i.e. functions that are enriched in genomes that belong to one group, but predominantly absent in genomes outside of it.

To use this feature you must have at least one categorical additional layer information (which can easily be done via anvi-import-misc-data), and at least one functional annotation source for your genomes storage (which will automatically be the case if every contigs database that were used when you run anvi-gen-genomes-storage was annotated with the same functional source). If you are unsure whether there are any functions available to you in your genomes storage, you can run anvi-db-info on it and thake a look at the output.

Let’s use the Prochlorococcus example to learn what we can do with this.

First we will compare the low-light vs. the high-light genomes in order to see if there are any functions that are unique to either group using the ‘light’ categorical additional layer data (if this doesn’t make sense to you please go back to one of the pangenome figures above and see the layer additional data ‘light’):

anvi-compute-functional-enrichment-in-pan -p PROCHLORO/Prochlorococcus_Pan-PAN.db \
                                          -g PROCHLORO-GENOMES.db \
                                          --category light \
                                          --annotation-source COG14_FUNCTION \
                                          -o enriched-functions.txt

Here is the structure of the output file enriched-functions.txt (there are more columns, scroll towards right to see them):

The following describes each column in this output:

• COG14_FUNCTION this column has the name of the specific function for which enrichment was calculated. In this example we chose to use COG14_FUNCTION for functional annotation, and hence the column title is COG14_FUNCTION. You can specify whichever functional annotation source you have in your PAN database using the --annotation-source, and then the analysis would be done according to that annotation source. Even if you have multiple functional annotation sources in your genome storage, only one source could be used for a single run of this program. If you wish, you could run it multiple times and each time use a different annotation source. If you don’t remember which annotation sources are available in your genomes storage, you can use the flag --list-annotation-sources.

• enrichment_score is a score to measure how much is this function unique to the genomes that belong to a specific group vs. all other genomes in your pangenome. This score was developed by Amy Willis. For more details about how we generate this score, see the note from Amy below.

• unadjusted_p_value is the p-value for the enrichment test (unadjusted for multiple tests).

• adjusted_q_value is the adjusted q-value to control for False Detection Rate (FDR) due to multiple tests (this is necessary since we are running the enrichment test each function separately).

• associated_groups is the list of groups (or labels) in your categorical data that are associated with the function. Notice that if the enrichment score is low, then this is meaningless (if the function is not enriched, then it is not really associated with any specific group/s). See Alon’s explanation above to understand how these are computed.

• function_accession is the function accession number.

• gene_clusters_ids are the gene clusters that were associated with this function. Notice that each gene cluster would be associated with a single function, but a function could be associated with multiple gene clusters.

• p_LL, p_HL for each group (in the case of this example there are two groups: LL and LH) there will be a column with the portion of the group members in which we detected the funtion.

• N_LL, N_HL for each group these columns specify the total number of genomes in the group.

Our example here includes only two categories (LL and HL), but you can have as many different categories as you want. Just remember that if some of your groups have very few genomes in them, then the statistical test will not be very reliable. The minimal number of genomes in a group for the test to be reliable depends on a number of factors, but we recommend proceeding with great caution if any of your groups have fewer than 8 genomes.

Now let’s search for one of the top functions in the table “Ser/Thr protein kinase RdoA involved in Cpx stress response, MazF antagonist”, which is enriched for the members of the LL group, and we can see in the table that it matches four gene clusters.

In fact, if we look carefully, then we find that this function matches a gene cluster that is unique for each one of the four low light clades, and is a single-copy core gene for the low light genomes in our pangenome. Cool!

Let’s look at another enriched function: Exonuclease VII, large subunit. When we search this function, we should be careful since the name of the function contains a comma, and the function search option in the interactive interface treats the comma as if it is separating multiple functions that are to be searched at the same time. Hence I just searched for Exonuclease VII, and here are the results:

We can see that the search matched hits for both Exonuclease VII, large and small subunits, with a total of 22 hits.

The large subunit matches a single gene cluster which is in the CORE LL, and the small subunit matches a gene cluster in each one of the clade specific cores (similar to what we saw above for the Ser/Thr protein kinase). Both of these genes are also part of the single copy core unique to low light members and absent from high light members.

Computing the average nucleotide identity for genomes (and other genome similarity metrics too!)

Anvi’o also contains a program, anvi-compute-genome-similarity, which uses various similarity metrics such as PyANI to compute average nucleotide identity across your genomes, and sourmash to compute mash distance across your genomes. It expects any combination of external genome files, internal genome files, or a fasta text file that points to the paths of FASTA files (each FASTA is assumed to be 1 genome). In addition, anvi-compute-genome-similarity optionally accepts a pan-db to add all results into it as additional layer data.

Here is an example with our Prochlorococcus pangenome:

anvi-compute-genome-similarity --external-genomes external-genomes.txt \
                               --program pyANI \
                               --output-dir ANI \
                               --num-threads 6 \
                               --pan-db PROCHLORO/Prochlorococcus_Pan-PAN.db

Once it is complete, we can visualize the pan genome using anvi-display-pan again to see what is new:

anvi-display-pan -g PROCHLORO-GENOMES.db \
                 -p PROCHLORO/Prochlorococcus_Pan-PAN.db

When you first look at it you will not see anything out of ordinary. But if you go to the ‘Layers’ tab, you will see the following addition:

If you click the checkbox for, say, ANI_percentage_identity, you will see a new set of entries will be added into the list of layer data entries. Then, if you click the button ‘Draw’ or ‘Redraw layer data’, you should see the ANI added to your display:

Yes. Magic.

You may need to change the min value from the interface for a better representation of ANI across your genomes in your own pangenome.

Summarizing an anvi’o pan genome

When you store your selections of gene clusters as a collection, anvi’o will allow you to summarize these results.

Even if you want to simply summarize everything in the pan genome without making any selections in the interface, you will still need a collection in the pan database. But luckily, you can use the program anvi-script-add-default-collection to add a default collection that contains every gene cluster.

The summary step gives you two important things: a static HTML web page you can compress and share with your colleagues or add into your publication as a supplementary data file, and a comprehensive TAB-delimited file in the output directory that describes every gene cluster.

You can summarize a collection using the program anvi-summarize, and a generic form of this command looks like this:

$ anvi-summarize -p PROJECT-PAN.db \
                 -g PROJECT-PAN-GENOMES.db \
                 -C COLLECTION_NAME \
                 -o PROJECT-SUMMARY

If you open the index.html file in the summary directory, you will see an output with some essential information about the analysis:

As well as the TAB-delimited file for gene clusters:

The structure of this file will look like this, and will give you an opportunity to investigate your gene clusters in a more detailed manner (you may need to scroll towards right to see more of the table):

This file will link each gene from each genome with every selection you’ve made and named through the interface or through the program anvi-import-collection, and will also give you access to the amino acid sequence and function of each gene.

Creating a pangenome with functions

It is also possible to define gene clusters based on gene functions rather than gene sequences. This may be a useful way to compare genomes that are too distant for high-resolution analyses with pangenomics. For such applications, please see the program anvi-display-functions