A primer on anvi'o

This tutorial is tailored for anvi’o v6 or later. You can learn the version of your installation by typing anvi-interactive -v. If you have an older version, some things will not work the way they should.

The goal of this tutorial is to explore some of the most fundamental aspects of anvi’o and its application to a real-world dataset. We organized it in multiple interconnected chapters, which all use the same dataset:

• Chapter I: Genome-resolved Metagenomics
• Chapter II: Automatic Binning
• Chapter III: Phylogenomics
• Chapter IV: Pangenomics
• Chapter V: Metabolism Prediction
• Chapter VI: Microbial Population Genetics
• Chapter VII: Genes and genomes across metagenomes
• Chapter VIII: From single-amino acid variants to protein structures

At the end of this tutorial, you should be able to

• Characterize high-quality genomes from metagenomes and manually curate them,
• Incorporate automatic binning results into your metagenomes,
• Perform phylogenomic analyses to learn about their approximate location in the tree of life,
• Determine core and accessory genes they contain compared to closely related genomes through pangenomics,
• Explore within population genomic diversity through single-nuleotide variants,
• Put a given metagenome-assembled, single-cell or isolate genome into the context of public metagenomes to investigate their core and accessory genes, and
• Investigate the genomic heterogeneity in the context of predicted protein structures.

All of these sections use the same publicly available metagenomic dataset generated by Sharon et al. (2013), which contains 11 feacal samples from a premature infant collected between Day 15 and Day 24 of her life in a time-series manner. After downloading raw sequencing reads, we co-assembled all samples and processed the resulting contigs longer than 1K using anvi’o following the anvi’o metagenomic workflow. The main reason we chose this dataset was its simplicity: more than 95% of metagenomic reads map back to the contigs, a great indication of the assembled contigs describe a very large fraction of this infant’s gut microbiome.

We hope you find the tutorial useful, and generously share your opinions or criticism should you have any.

Downloading the pre-packaged Infant Gut Dataset

To download and unpack the infant gut data-pack, copy-paste the following commands into your terminal:

curl -L https://figshare.com/ndownloader/files/45076909 -o INFANT-GUT-TUTORIAL.tar.gz
tar -zxvf INFANT-GUT-TUTORIAL.tar.gz && cd INFANT-GUT-TUTORIAL

This will download a 224 Mb compressed file, and unpack it, which will take an additional 502 Mb.

If you are using a newer version of anvi’o than was the one that was used to generate these databases (perhaps you are following the development branch), you may need to run anvi-migrate to get them up to date. If you are not sure whether you need this, do not worry - you could safely skip it and anvi’o would later remind you what exactly needs to be done.

anvi-migrate --migrate-safely *.db

If you were sent here somewhere from down below, now you can go back. If you have no idea what this means, ignore this notice, and continue reading. You’re OK :)

Chapter I: Genome-resolved Metagenomics

The purpose of this tutorial is to have a conversation about genome-resolved metagenomics (with a focus on manual binning) using the Infant Gut Dataset (IGD), which was generated, analyzed, and published by Sharon et al. (2013), and was re-analyzed in the anvi’o methods paper.

By the end of the tutorial, you will be able to:

• Familiarize with the interactive interface for binning (~25 minutes)
• Inspect contigs in the context of their metagenomic signal (~5 minutes)
• Characterize bins by manual binning (~20 minutes)
• Summarize manual binning results for downstream analyses (~20 minutes)
• Manually curate individual bins for quality control (~10 minutes)
• Import and visualize external binning results (~20 minutes)
• Combine manual and automatic binning (~10 minutes)

By the end of this chapter, you should (1) have a comprehensive understanding of genome-resolved metagenomics, and (2) be able to characterize and manually curate genomes from your own assembly outputs moving forwards.

A typical anvi’o genome-resolved metagenomic workflow starts with one or more BAM files and a FASTA file of your contigs. There are many ways to get your contigs and BAM files for your metagenomes and we have started implementing a tutorial that describes the workflows we often use. But in this tutorial we will start from a point in the workflow where you have used your BAM and FASTA files to generate anvi’o contigs and profile databases.

While this tutorial will take you through a simple analysis of a real dataset, there also is available a more comprehensive (but more abstract) tutorial on anvi’o metagenomic workflow.

Using the files in the data-pack directory, let’s take a first look at the merged profile database for the infant gut dataset metagenome. If you copy-paste this to your terminal:

anvi-interactive -p PROFILE.db -c CONTIGS.db

The anvi’o interactive interface should welcome you with this display (after you click “draw”):

When it is time to type other commands, you can close the window, go back to the terminal and press CTRL + C to kill the server.

Importing taxonomy for genes

Anvi’o can work with gene-level taxonomic annotations, but gene-level taxonomy is not useful for anything beyond occasional help with manual binning. Once gene-level taxonomy is added into the contigs database, anvi’o will determine the taxonomy of each contig based on the taxonomic affiliation of genes they describe, and display them in the interface whenever possible.

Centrifuge (code, paper) is one of the options to import taxonomic annotations into an anvi’o contigs database. Centrifuge files for the IGD are already in the directory additional-files/centrifuge-files.

If you import these files into the contigs database the following way,

anvi-import-taxonomy-for-genes -c CONTIGS.db \
                               -p centrifuge \
                               -i additional-files/centrifuge-files/centrifuge_report.tsv \
                               additional-files/centrifuge-files/centrifuge_hits.tsv

And run the interactive interface again,

anvi-interactive -p PROFILE.db -c CONTIGS.db

You will see an additional layer with taxonomy:

In the Layers tab find the Taxonomy layer, set its height to 200, then drag the layer in between DAY24 and hmms_Ribosomal_RNAs, and click Draw again. Then click Save State button, and overwrite the default state. This will make sure anvi’o remembers to make the height of that layer 200px the next time you run the interactive interface!

Inferring taxonomy for metagenomes

So at this point we don’t have any idea about what genomes do we have in this dataset, but anvi’o can make sense of the taxonomic make up of a given metagenome by characterizing taxonomic affiliations of single-copy core genes. The details of this anvi-run-scg-taxonomy is described here in greater detail.

You can take a very quick look at the taxonomic composition of the metagenome through the command line first:

anvi-estimate-scg-taxonomy -c CONTIGS.db \
                           --metagenome-mode

which should give us this output for the IGD:

╒════════════════════╤════════════════════╤══════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════╕
│                    │   percent_identity │ taxonomy                                                                                                                         │
╞════════════════════╪════════════════════╪══════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════╡
│ Ribosomal_S6_11655 │               98.9 │ Bacteria / Firmicutes / Bacilli / Staphylococcales / Staphylococcaceae / Staphylococcus / Staphylococcus epidermidis             │
├────────────────────┼────────────────────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Ribosomal_S6_12163 │               98.9 │ Bacteria / Firmicutes / Bacilli / Staphylococcales / Staphylococcaceae / Staphylococcus / Staphylococcus hominis                 │
├────────────────────┼────────────────────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Ribosomal_S6_15200 │               98.9 │ Bacteria / Firmicutes / Bacilli / Lactobacillales / Lactobacillaceae / Leuconostoc / Leuconostoc citreum                         │
├────────────────────┼────────────────────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Ribosomal_S6_25880 │               98.9 │ Bacteria / Firmicutes / Clostridia / Tissierellales / Helcococcaceae / Finegoldia / Finegoldia magna                             │
├────────────────────┼────────────────────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Ribosomal_S6_2915  │               97.9 │ Bacteria / Firmicutes / Bacilli / Lactobacillales / Enterococcaceae / Enterococcus / Enterococcus faecalis                       │
├────────────────────┼────────────────────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Ribosomal_S6_29818 │               98.9 │ Bacteria / Firmicutes / Bacilli / Lactobacillales / Streptococcaceae / Streptococcus / Streptococcus oralis                      │
├────────────────────┼────────────────────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Ribosomal_S6_30904 │               92.5 │ Bacteria / Firmicutes / Clostridia / Tissierellales / Helcococcaceae / Anaerococcus / Anaerococcus sp002359915                   │
├────────────────────┼────────────────────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Ribosomal_S6_4484  │               98.9 │ Bacteria / Firmicutes / Clostridia / Tissierellales / Peptoniphilaceae / Peptoniphilus / Peptoniphilus rhinitidis                │
├────────────────────┼────────────────────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Ribosomal_S6_6421  │                100 │ Bacteria / Actinobacteriota / Actinobacteria / Propionibacteriales / Propionibacteriaceae / Cutibacterium / Cutibacterium avidum │
├────────────────────┼────────────────────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┤
│ Ribosomal_S6_7660  │               98.9 │ Bacteria / Firmicutes / Bacilli / Staphylococcales / Staphylococcaceae / Staphylococcus /                                        │
╘════════════════════╧════════════════════╧══════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════════╛

Good, but could have been better. Why?

Pro tip: we have a profile database, what does it mean and how it could improve the information we see here?

Making use of our profile database the following way, will give us a little more information about our dataset:

anvi-estimate-scg-taxonomy -c CONTIGS.db \
                           -p PROFILE.db \
                           --metagenome-mode \
                           --compute-scg-coverages

which should give us the following output:

Taxa in metagenome "Infant Gut Contigs from Sharon et al."
===============================================
╒════════════════════╤════════════════════╤════════════════════════════════╤═══════════╤═══════════╤══════════╤═══════════╤═══════════╤══════════════╕
│                    │   percent_identity │ taxonomy                       │   DAY_15A │   DAY_15B │   DAY_16 │   DAY_17A │   DAY_17B │ ... 6 more   │
╞════════════════════╪════════════════════╪════════════════════════════════╪═══════════╪═══════════╪══════════╪═══════════╪═══════════╪══════════════╡
│ Ribosomal_S6_2915  │               97.9 │ (s) Enterococcus faecalis      │   372.512 │   699.853 │  663.241 │    186.34 │    1149.6 │ ... 6 more   │
├────────────────────┼────────────────────┼────────────────────────────────┼───────────┼───────────┼──────────┼───────────┼───────────┼──────────────┤
│ Ribosomal_S6_11655 │               98.9 │ (s) Staphylococcus epidermidis │   112.694 │   172.478 │  147.304 │   23.3901 │   140.769 │ ... 6 more   │
├────────────────────┼────────────────────┼────────────────────────────────┼───────────┼───────────┼──────────┼───────────┼───────────┼──────────────┤
│ Ribosomal_S6_4484  │               98.9 │ (s) Peptoniphilus lacydonensis │         0 │         0 │        0 │         0 │         0 │ ... 6 more   │
├────────────────────┼────────────────────┼────────────────────────────────┼───────────┼───────────┼──────────┼───────────┼───────────┼──────────────┤
│ Ribosomal_S6_7660  │               98.9 │ (g) Staphylococcus             │         0 │         0 │        0 │         0 │         0 │ ... 6 more   │
├────────────────────┼────────────────────┼────────────────────────────────┼───────────┼───────────┼──────────┼───────────┼───────────┼──────────────┤
│ Ribosomal_S6_6421  │                100 │ (s) Cutibacterium avidum       │   17.8935 │   5.87368 │  4.79909 │         0 │         0 │ ... 6 more   │
├────────────────────┼────────────────────┼────────────────────────────────┼───────────┼───────────┼──────────┼───────────┼───────────┼──────────────┤
│ Ribosomal_S6_12163 │               98.9 │ (s) Staphylococcus hominis     │   2.39322 │   22.5447 │  13.2806 │         0 │   9.94853 │ ... 6 more   │
├────────────────────┼────────────────────┼────────────────────────────────┼───────────┼───────────┼──────────┼───────────┼───────────┼──────────────┤
│ Ribosomal_S6_15200 │               98.9 │ (s) Leuconostoc citreum        │         0 │         0 │  1.86532 │         0 │    1.6936 │ ... 6 more   │
├────────────────────┼────────────────────┼────────────────────────────────┼───────────┼───────────┼──────────┼───────────┼───────────┼──────────────┤
│ Ribosomal_S6_29818 │               98.9 │ (s) Streptococcus oralis       │         0 │         0 │        0 │         0 │         0 │ ... 6 more   │
├────────────────────┼────────────────────┼────────────────────────────────┼───────────┼───────────┼──────────┼───────────┼───────────┼──────────────┤
│ Ribosomal_S6_25880 │               98.9 │ (s) Finegoldia magna           │         0 │         0 │        0 │         0 │         0 │ ... 6 more   │
├────────────────────┼────────────────────┼────────────────────────────────┼───────────┼───────────┼──────────┼───────────┼───────────┼──────────────┤
│ Ribosomal_S6_30904 │               98.1 │ (s) Anaerococcus nagyae        │  0.528958 │         0 │  0.14157 │         0 │         0 │ ... 6 more   │
╘════════════════════╧════════════════════╧════════════════════════════════╧═══════════╧═══════════╧══════════╧═══════════╧═══════════╧══════════════╛

These look like information that would have been useful to have in front of us in our interactive interface. Luckily, anvi’o can add these taxonomic insights into a given profile database, if you change the previous command just a bit:

anvi-estimate-scg-taxonomy -c CONTIGS.db \
                           -p PROFILE.db \
                           --metagenome-mode \
                           --compute-scg-coverages \
                           --update-profile-db-with-taxonomy

Tadaa. Now let’s take another look at our interactive interface and find the additional data for our layers:

anvi-interactive -c CONTIGS.db \
                 -p PROFILE.db

At this point we have an overall idea about the make up of this metagenome, but we don’t have any genomes from it. The following sections will cover some of multiple ways to do this.

Manual identification of genomes in the Infant Gut Dataset

Let’s first talk about the signal that leads to this specific organization of contigs we see in our display, and its biological implications.

Next, we can perform a round of manual binning. This should take about 10 minutes.

Here is an example of 16 bins we identified for comparison AFTER you performed your own binning:

Please save your bins as a collection. You can give your collection any name, but if you call it default, anvi’o will treat it differently.

In the anvi’o lingo, a collection is something that describes one or more bins, each of which describe one or more contigs.

If you identified near-complete genomes, then congratulations, you have characterized genomic contents of microbial populations de novo.

Summarizing the binning results

Why do we do binning? Because we are interested in making sense of our metagenomes in the context of genomes we have recovered through binning. Understanding the distribution patterns of the genomes we have in a collection in a quantitative fashion, or getting back a table of function names found in each one of them, or even summarizing our bins as distinct FASTA files are critical for next steps of every binning analysis. After all, binning is a boring detail before you start doing your science.

To ensure that you have everything you need to continue working with the outcomes of your binning effort outside of anvi’o, we have a program called anvi-summarize. It is possible to summarize any collection stored in an anvi’o profile database through this program. The result is a static HTML page that can be viewed on any computer.

anvi-import-collection additional-files/collections/merens.txt \
                       --bins-info additional-files/collections/merens-info.txt \
                       -p PROFILE.db \
                       -c CONTIGS.db \
                       -C default

Let’s summarize the collection you have just created:

anvi-summarize -p PROFILE.db \
               -c CONTIGS.db \
               -C default \
               -o SUMMARY

Once the summary is finished, take a minute to look at its contents.

Renaming bins in your collection (from chaos to order)

As you can see from the summary file, at this point bin names are random, and we often find it useful to put some order on this front. This becomes an extremely useful strategy especially when the intention is to merge multiple binning efforts later. For this task we use the program anvi-rename-bins:

anvi-rename-bins -p PROFILE.db \
                 -c CONTIGS.db  \
                 --collection-to-read default  \
                 --collection-to-write MAGs  \
                 --call-MAGs \
                 --prefix IGD \
                 --report-file rename-bins-report.txt

With those settings, a new collection MAG will be created in which (1) bins with a completion >70% are identified as MAGs (stands for Metagenome-Assembled Genome), and (2) bins and MAGs are attached the prefix IGD and renamed based on the difference between completion and redundancy.

Now we can use the program anvi-summarize to summarize the new collection:

anvi-summarize -p PROFILE.db \
               -c CONTIGS.db \
               -C MAGs \
               -o SUMMARY_AFTER_RENAME

You can now visualize the results by double-clicking on the index.html file present in the newly created folder SUMMARY_MAGs.

Here are the MAGs we got:

Better.

Refining individual MAGs: the curation step

Why do we need to curate bins? Please see this paper to see how poorly refined bins can influence ecological and evolutionary insights we gain from MAGs. There is also a step-by-step workflow to describe steps of manual curation here.

Why not relying only on single-copy core genes to estimate purity of a genome bin? Let’s think about this altogether. But once we are done with this, please see the relevant section in this study.

To straighten the quality of the MAGs collection, it is possible to visualize individual bins and if needed, refine them. For this we use the program anvi-refine. For instance, if you were to be interested in refining one of the bins in our current collection, you could run this command:

anvi-refine -p PROFILE.db \
            -c CONTIGS.db \
            -C MAGs \
            -b IGD_MAG_00001

Now the interactive interface only displays contigs from a single bin. During this curation step, one can try different clustering strategies (i.e. by only relying on coverage, or only relying on sequence composition) to identify outliers and investigate carefully whether they may be contaminants. You can select everything, and remove those contigs you don’t want to keep in the bin before using the Bins panel to store your updated set of contigs in the database.

Here is an example of MAG we had to curate (we removed three contigs):

Storing the refined new bin in the database will modify the collection, but you will need to run anvi-summarize if you want the summary output to also be updated.

This is indeed a simple example. But in some cases refining a given MAG can take hours. But this is an extremely critical step of genome-resolved metagenomics studies, especially when important claims are made based on MAGs.

Chapter II: Automatic Binning

Even if you prefer manual binning over automatic binning for the sake of accuracy and control over your data, automatic binning is an unavoidable need due to performance limitations associated with manual binning. The actual purpose of this chapter is to talk about advantages and disadvantages of automatic binning, by comparing multiple binning approaches on the simplest real-world gut metagenome there is: the infant gut dataset. However, while we are going through this, you will also learn about how to incorporate automatic binning results into anvi’o.

The directory additional-files/external-binning-results contains a number of files that describe the binning of contigs in the IGD based on various automatic and manual approaches. These files include (1) outputs from some of the well-known binning algorithms (i.e., GROOPM.txt, MAXBIN.txt, METABAT.txt, BINSANITY_REFINE.txt, MYCC.txt, and CONCOCT.txt), (2) the original binning of this dataset (SHARON_et_al.txt), and the manual binning we performed in the anvi’o paper (MEREN_et_al.txt).

# GroopM v0.3.5 (followed the general workflow on their manual)
groopm parse groopm.db contigs.fa [list of bam files]
groopm core groopm.db -c 1000 -s 10 -b 1000000
groopm recruit groopm.db -c 500 -s 200
groopm extract groopm.db contigs.fa

# MetaBat v0.26.3 (used jgi_summarize_bam_contig_depths to get a depth file from BAM files).
metabat -i contigs.fa -a depth.txt -o bin

# MaxBin v2.1.1
run_MaxBin.pl -contig contigs.fa -out maxbin_IGM -abund_list [list of all coverage files in associated format]

# MyCC (ran via the docker image available here: https://sourceforge.net/projects/sb2nhri/files/MyCC/)
(used jgi_summarize_bam_contig_depths to get depth file from BAM files per the authors suggestion)
MyCC.py contigs.fa 4mer -a depth.txt

# BinSanity + refinement v0.2.1
Binsanity -f . -l contigs.fa -p -10 -c igm.coverage.lognorm

## After bin inspection using CheckM & Anvi'o, INFANT-GUT-ASSEMBLY-bin_18
## was refined using the following parameters
Binsanity-refine -f . -l -p -150 -c igm.coverage.lognorm

Now you have the background information about where these files are coming from. Moving on. But just before we continue to move on, let’s remove the default collection from our profile database, if there is one, to avoid any confusion.

You can use the program anvi-show-collections-and-bins to see all collections in your anvi’o profile database:

anvi-show-collections-and-bins -p PROFILE.db

And use this one to remove one that’s called default:

anvi-delete-collection -p PROFILE.db \
                       -C default

Importing an external binning result

You can create a collection by using the interactive interface (e.g., the default and MAGs collections you just created), or you can import external binning results into your profile database as a collection and see how that collection groups contigs. For instance, let’s import the CONCOCT collection:

anvi-import-collection additional-files/external-binning-results/CONCOCT.txt \
                       -c CONTIGS.db \
                       -p PROFILE.db \
                       -C CONCOCT \
                       --contigs-mode

you can immediately see what collections are available in a given profile database using the program anvi-show-collections-and-bins (which in this case should show us the collections CONCOCT, and MAGs):

anvi-show-collections-and-bins -p PROFILE.db

You can get a quick idea regarding the estimated completion of bins in a given collection:

anvi-estimate-genome-completeness -p PROFILE.db \
                                  -c CONTIGS.db \
                                  -C CONCOCT

You can also get an idea about their taxonomy:

anvi-estimate-scg-taxonomy -p PROFILE.db \
                           -c CONTIGS.db \
                           -C CONCOCT

OK. Let’s run the interactive interface again with the CONCOCT collection:

anvi-interactive -p PROFILE.db \
                 -c CONTIGS.db \
                 --collection-autoload CONCOCT

Alternatively you could load the interface without the --collection-autoload flag, and click Bins > Load bin collection > CONCOCT > Load to load the CONCOCT collection.

To turn off text annotation, go to Settings > Options > Bins Selection and then uncheck Show names. You will then see something like this:

So this is how you load and display an external collection. So far so good.

Comparing multiple binning approaches

Since we have all these results from different binning approaches, it clearly would have been interesting to compare them to each other (because benchmarking stuff is often very insightful). But how to do it? The simplest way to do it is to assume a ‘true organization of contigs’, and then investigate every other approach with respect to that.

Here we have multiple independent sources of information we could use. Including (1) the organization of contigs based on hierarchical clustering analysis, (2) per-contig taxonomy estimated from the gene-level taxonomic annotations by Centrifuge, and (3) results from the original publication from Sharon et al., in which authors did a very careful job to identify every genome in the dataset (even resolving the Staphylococcus pangenome, which is extremely hard for automatic binning approaches that work with a single co-assembly). So these are the things we can build upon for a modest comparison.

To include binning results in this framework, we could import each collection into the profile database the way we imported CONCOCT. But unfortunately at any given time there could only be one collection that can be displayed in the interface. Luckily there are other things we can do. For instance, as a workaround, we can merge all binning results into a single file, and use that file as an ‘additional data file’ to visualize them in the interactive interface.

Anvi’o has a script called anvi-script-merge-collections to merge multiple files from external binning results into a single merged file (don’t ask why):

anvi-script-merge-collections -c CONTIGS.db \
                              -i additional-files/external-binning-results/*.txt \
                              -o collections.tsv

If you take a look at this file, you will realize that it has a very simple format:

head collections.tsv | column -t
contig                           BINSANITY_REFINE                CONCOCT  GROOPM     MAXBIN      METABAT               MYCC              SHARON_et_al
Day17a_QCcontig1000_split_00001  INFANT-GUT-ASSEMBLY-bin_19.fna  Bin_4    db_bin_11  maxbin.008  metabat_igm.unbinned  Cluster.5.fasta   Finegoldia_magna
Day17a_QCcontig1001_split_00001  INFANT-GUT-ASSEMBLY-bin_6.fna   Bin_7    db_bin_46  maxbin.006  metabat_igm.unbinned  Cluster.3.fasta   Staphylococcus_epidermidis_virus_014
Day17a_QCcontig1002_split_00001  INFANT-GUT-ASSEMBLY-bin_19.fna  Bin_4    db_bin_11  maxbin.007  metabat_igm.unbinned  Cluster.5.fasta   Finegoldia_magna
Day17a_QCcontig1003_split_00001  INFANT-GUT-ASSEMBLY-bin_14.fna  Bin_2    db_bin_1   maxbin.009  metabat_igm.7         Cluster.12.fasta
Day17a_QCcontig1004_split_00001  INFANT-GUT-ASSEMBLY-bin_16.fna  Bin_3    db_bin_8   maxbin.008  metabat_igm.10        Cluster.14.fasta
Day17a_QCcontig1005_split_00001  INFANT-GUT-ASSEMBLY-bin_13.fna  Bin_5    db_bin_47  maxbin.007  metabat_igm.unbinned  Cluster.3.fasta   Staphylococcus_epidermidis_viruses
Day17a_QCcontig1006_split_00001  INFANT-GUT-ASSEMBLY-bin_16.fna  Bin_3    db_bin_8   maxbin.008  metabat_igm.10        Cluster.14.fasta  Leuconostoc_citreum
Day17a_QCcontig1007_split_00001  INFANT-GUT-ASSEMBLY-bin_16.fna  Bin_3    db_bin_8   maxbin.008  metabat_igm.10        Cluster.14.fasta
Day17a_QCcontig1008_split_00001  INFANT-GUT-ASSEMBLY-bin_14.fna  Bin_2    db_bin_1   maxbin.009  metabat_igm.7         Cluster.8.fasta   Candida_albcans

Good. Now you can run the interactive interface to display all bins in all collections stored in collections.tsv as additional layers:

anvi-interactive -p PROFILE.db \
                 -c CONTIGS.db \
                 -A collections.tsv

anvi-import-misc-data collections.tsv \
                      -p PROFILE.db \
                      -t items

anvi-interactive -p PROFILE.db -c CONTIGS.db

anvi-delete-misc-data -p PROFILE.db \
                      -t items \
                      --just-do-it

Fine.

At this point you should be seeing a display similar to this (after setting the height of each additional layer to 200px):

The legends for each of the bin collections are available in the Legends tab of Settings. To visually emphasize relationships between bins, you can change the color of each bin manually by clicking on the colored boxes in the legends. Or, if you’re not a masochist, you can use the program anvi-import-state to import an anvi’o state we have created for you:

anvi-import-state --state additional-files/state-files/state-merged.json \
                  --name default \
                  -p PROFILE.db

and run the interactive interface again with the same command line,

anvi-interactive -p PROFILE.db \
                 -c CONTIGS.db \
                 -A collections.tsv

this time you should get this display:

So far so good?

Now we can discuss about different approaches of automatic binning.

Please note that the algorithms we have used here may have been improved since the time we did these analyses, therefore please don’t make any decisions about their performance or efficacy based on what you are seeing here.

Just a reminder, once you have the interactive interface in front of you, you can in fact investigate the taxonomy of contigs by BLASTing them against various NCBI collections using the right-click menu to have a second opinion about what do public databases think they are:

Mouse section moved under Settings > Data, after version 7.0.0.

We recently have added an option to quickly run them on BIGSI. Sometimes it takes a split second, sometimes (especially when you are in France) it takes minutes. So, no promises, but try it for sure! It is an excellent algorithm.

Manually curating automatic binning outputs

OK. Let’s assume, we didn’t see the interactive interface, and we have no idea about the dataset. We didn’t do any of the things we did up to this point. We just had profiled and merged the IGD, and we did binning of this dataset using MaxBin. Let’s start by importing MaxBin results into the profile database as a collection:

anvi-import-collection additional-files/external-binning-results/MAXBIN.txt \
                       -c CONTIGS.db \
                       -p PROFILE.db \
                       -C MAXBIN \
                       --contigs-mode

From here, there are two things we can do very quickly. First, we can create a summary of our new collection or we can take a quick look at the completion / redundancy estimates of bins described by this collection from the command line, using anvi-estimate-genome-completeness:

anvi-estimate-genome-completeness -p PROFILE.db \
                                  -c CONTIGS.db \
                                  -C MAXBIN

Alternatively, we can take a quick look at the binning results by initiating the interactive interface in collection mode:

anvi-interactive -p PROFILE.db \
                 -c CONTIGS.db \
                 -C MAXBIN

This command should give you a display similar to this:

All previous interactive displays were at the contig-level (each leaf in the center tree was a contig). However, this display is at the bin-level. Instead of contigs, this display shows us the distribution of bins MaxBin identified. We also have completion and redundancy estimates for each bin, which helps us make some early sense of what is going on.

Please read this post to learn more about completion and redundancy estimates: Assessing completion and contamination of metagenome-assembled genomes

It is clear that some bins are not as well-resolved as others. For instance, bins maxbin_007 and maxbin_008 have redundancy estimates of 22% and 91%, respectively, which suggests each of them describe multiple distinct populations. Well, clearly we would have preferred those bins to behave.

If you order bins based on their detection across metagenomes (by changing the ‘Items order’ to ‘detection’ from the menu in the Main tab), you can also see that bins maxbin_007 and maxbin_008 are right next to each other. This suggests that it may be a good idea to simply merge these bins first, and then use the program anvi-refine to avoid issues of over-splitting populations of interest. Let’s use the program anvi-merge-bins to merge them into a single bin first:

anvi-merge-bins -p PROFILE.db \
                --collection-name MAXBIN \
                --bin-names-list "maxbin_007, maxbin_008" \
                --new-bin-name maxbin_007_and_008

If you take a quick look at the collection again, you can see that the new bin maxbin_007_and_008 is the bin that needs immediate action based on single-copy core genes:

anvi-estimate-genome-completeness -p PROFILE.db \
                                  -c CONTIGS.db \
                                  -C MAXBIN

Bins in collection "MAXBIN"
===============================================
╒════════════════════╤══════════╤══════════════╤════════════════╤════════════════╤══════════════╤════════════════╕
│ bin name           │ domain   │   confidence │   % completion │   % redundancy │   num_splits │   total length │
╞════════════════════╪══════════╪══════════════╪════════════════╪════════════════╪══════════════╪════════════════╡
│ maxbin_001         │ BACTERIA │            1 │            100 │           4.23 │          148 │        2969341 │
├────────────────────┼──────────┼──────────────┼────────────────┼────────────────┼──────────────┼────────────────┤
│ maxbin_002         │ BACTERIA │          0.9 │          94.37 │              0 │           88 │        1801068 │
├────────────────────┼──────────┼──────────────┼────────────────┼────────────────┼──────────────┼────────────────┤
│ maxbin_003         │ BACTERIA │          0.7 │          81.69 │           2.82 │          144 │        2764617 │
├────────────────────┼──────────┼──────────────┼────────────────┼────────────────┼──────────────┼────────────────┤
│ maxbin_004         │ BACTERIA │          0.9 │          97.18 │           1.41 │          188 │        2571878 │
├────────────────────┼──────────┼──────────────┼────────────────┼────────────────┼──────────────┼────────────────┤
│ maxbin_005         │ BACTERIA │            1 │          98.59 │              0 │          151 │        2555414 │
├────────────────────┼──────────┼──────────────┼────────────────┼────────────────┼──────────────┼────────────────┤
│ maxbin_006         │ BACTERIA │          0.6 │          78.87 │           5.63 │          305 │        2901149 │
├────────────────────┼──────────┼──────────────┼────────────────┼────────────────┼──────────────┼────────────────┤
│ maxbin_009         │ EUKARYA  │          0.7 │          83.13 │           7.23 │         1379 │       13915165 │
├────────────────────┼──────────┼──────────────┼────────────────┼────────────────┼──────────────┼────────────────┤
│ maxbin_007_and_008 │ BACTERIA │            1 │          95.77 │         154.93 │         2381 │        6287535 │
╘════════════════════╧══════════╧══════════════╧════════════════╧════════════════╧══════════════╧════════════════╛

Fine. Let’s refine that bin:

anvi-refine -p PROFILE.db \
            -c CONTIGS.db \
            -C MAXBIN \
            -b maxbin_007_and_008

Which should give us this display, on which we see the distribution of contigs that were originally binned into maxbin_007 and maxbin_008 across samples. The hierarchical clustering picked up some trends, and you can see clusters one can identify quickly:

We can now make the following selections to split these two bins into six, and update our database by storing these refined bins from the Bins panel in the interface:

If you take another look after this, you can see that maxbin_007_and_008 is replaced with smaller bins, and the total redundancy in the collection is much lower:

anvi-estimate-genome-completeness -p PROFILE.db \
                                  -c CONTIGS.db \
                                  -C MAXBIN

Bins in collection "MAXBIN"
===============================================
╒══════════════════════╤══════════╤══════════════╤════════════════╤════════════════╤══════════════╤════════════════╕
│ bin name             │ domain   │   confidence │   % completion │   % redundancy │   num_splits │   total length │
╞══════════════════════╪══════════╪══════════════╪════════════════╪════════════════╪══════════════╪════════════════╡
│ maxbin_001           │ BACTERIA │            1 │            100 │           4.23 │          148 │        2969341 │
├──────────────────────┼──────────┼──────────────┼────────────────┼────────────────┼──────────────┼────────────────┤
│ maxbin_002           │ BACTERIA │          0.9 │          94.37 │              0 │           88 │        1801068 │
├──────────────────────┼──────────┼──────────────┼────────────────┼────────────────┼──────────────┼────────────────┤
│ maxbin_003           │ BACTERIA │          0.7 │          81.69 │           2.82 │          144 │        2764617 │
├──────────────────────┼──────────┼──────────────┼────────────────┼────────────────┼──────────────┼────────────────┤
│ maxbin_004           │ BACTERIA │          0.9 │          97.18 │           1.41 │          188 │        2571878 │
├──────────────────────┼──────────┼──────────────┼────────────────┼────────────────┼──────────────┼────────────────┤
│ maxbin_005           │ BACTERIA │            1 │          98.59 │              0 │          151 │        2555414 │
├──────────────────────┼──────────┼──────────────┼────────────────┼────────────────┼──────────────┼────────────────┤
│ maxbin_006           │ BACTERIA │          0.6 │          78.87 │           5.63 │          305 │        2901149 │
├──────────────────────┼──────────┼──────────────┼────────────────┼────────────────┼──────────────┼────────────────┤
│ maxbin_009           │ EUKARYA  │          0.7 │          83.13 │           7.23 │         1379 │       13915165 │
├──────────────────────┼──────────┼──────────────┼────────────────┼────────────────┼──────────────┼────────────────┤
│ maxbin_007_and_008_1 │ BACTERIA │          0.4 │          38.03 │           2.82 │          417 │         696607 │
├──────────────────────┼──────────┼──────────────┼────────────────┼────────────────┼──────────────┼────────────────┤
│ maxbin_007_and_008_2 │ BACTERIA │          0.6 │          47.89 │              0 │          369 │         753160 │
├──────────────────────┼──────────┼──────────────┼────────────────┼────────────────┼──────────────┼────────────────┤
│ maxbin_007_and_008_3 │ BLANK    │            1 │              0 │              0 │          206 │         290126 │
├──────────────────────┼──────────┼──────────────┼────────────────┼────────────────┼──────────────┼────────────────┤
│ maxbin_007_and_008_4 │ BACTERIA │          0.3 │          32.39 │           5.63 │          488 │         810476 │
├──────────────────────┼──────────┼──────────────┼────────────────┼────────────────┼──────────────┼────────────────┤
│ maxbin_007_and_008_5 │ BACTERIA │          0.3 │          38.03 │           2.82 │          247 │        2333489 │
├──────────────────────┼──────────┼──────────────┼────────────────┼────────────────┼──────────────┼────────────────┤
│ maxbin_007_and_008_6 │ BACTERIA │          0.8 │          71.83 │           2.82 │          556 │        1237568 │
╘══════════════════════╧══════════╧══════════════╧════════════════╧════════════════╧══════════════╧════════════════╛

The take home message here is that even when automatic binning approaches yield poorly identified bins, it is possible to improve the final results through a manual refinement step. Clearly these extra steps require a lot of expertise, intuition, attention, and decision making. And fortunately you are all familiar with each one of them because science.

Thank you for following the tutorial this far!

anvi-display-contigs-stats CONTIGS.db

anvi-cluster-with-concoct -p PROFILE.db \
                          -c CONTIGS.db \
                          --num-clusters 5 \
                          -C CONCOCT_C5

anvi-interactive -p PROFILE.db \
                 -c CONTIGS.db \
                 --title 'Infant Gut Time Series by Sharon et al. [w/ 5 CONCOCT clusters]' \
                 --collection-autoload CONCOCT_C5

anvi-refine -p PROFILE.db \
            -c CONTIGS.db \
            -C CONCOCT_C5 \
            -b Bin_1

Chapter III: Phylogenomics

This is more of a practical tutorial to do phylogenomic analyses on metagenome-assembled genomes described in anvi’o collections. For a more abstract tutorial on phylogenomics, please consider first reading ‘An anvi’o workflow for phylogenomics’.

To see a practical application of phylogenomics see this workflow.

If you haven’t followed the previous sections of the tutorial, you will need the anvi’o merged profile database and the anvi’o contigs database for the IGD available to you. Before you continue, please click here, do everything mentioned there, and come back right here to continue following the tutorial from the next line when you read the directive go back.

Please run the following command in the IGD dir, so you have everything you need. We will simply import our previously generated collection of bins in the IGD dataset as the default collection:

anvi-import-collection additional-files/collections/merens.txt \
                       --bins-info additional-files/collections/merens-info.txt \
                       -p PROFILE.db \
                       -c CONTIGS.db \
                       -C default

At this point, you have in your anvi’o profile database a collection with multiple bins:

anvi-show-collections-and-bins -p PROFILE.db

Collection: "default"
===============================================
Collection ID ................................: default
Number of bins ...............................: 13
Number of splits described ...................: 4,451
Bin names ....................................: Aneorococcus_sp, C_albicans, E_facealis, F_magna, L_citreum, P_acnes, P_avidum, P_rhinitidis, S_aureus, S_epidermidis, S_hominis, S_lugdunensis, Streptococcus

Putting genomes in a phylogenomic context is one of the common ways to compare them to each other. The common practice is to concatenate aligned sequences of single-copy core genes for each genome of interest, and generate a phylogenomic tree by analyzing the resulting alignment.

Let’s assume we want to run a phylogenomic analysis on all genome bins we have in the collection merens in the IGD (you may have your own collections somewhere, that is fine too).

In order to do the phylogenomic analysis, we will need a FASTA file of concatenated genes. And to get that FASTA file out of our anvi’o databases, we will primarily use the program anvi-get-sequences-for-hmm-hits.

Selecting genes from an HMM Profile

We first need to identify an HMM profile, and then select some gene names from this profile to play with.

Going back to the IGD, let’s start by looking at what HMM profiles are available to us:

anvi-get-sequences-for-hmm-hits -c CONTIGS.db \
                                -p PROFILE.db \
                                -o seqs-for-phylogenomics.fa \
                                --list-hmm-sources

* Bacteria_71 [type: singlecopy] [num genes: 71]
* Archaea_76 [type: singlecopy] [num genes: 76]
* Protista_83 [type: singlecopy] [num genes: 83]
* Ribosomal_RNAs [type: Ribosomal_RNAs] [num genes: 12]

As you know, you can use anvi-run-hmms program with custom made HMM profiles to add your own HMMs into the contigs database.

Alright. We have two. Let’s see what genes do we have in Bacteria_71:

anvi-get-sequences-for-hmm-hits -c CONTIGS.db \
                                -p PROFILE.db \
                                -o seqs-for-phylogenomics.fa \
                                --hmm-source Bacteria_71 \
                                --list-available-gene-names

* Bacteria_71 [type: singlecopy]: ADK, AICARFT_IMPCHas, ATP-synt, ATP-synt_A,
Chorismate_synt, EF_TS, Exonuc_VII_L, GrpE, Ham1p_like, IPPT, OSCP, PGK,
Pept_tRNA_hydro, RBFA, RNA_pol_L, RNA_pol_Rpb6, RRF, RecO_C, Ribonuclease_P,
Ribosom_S12_S23, Ribosomal_L1, Ribosomal_L13, Ribosomal_L14, Ribosomal_L16,
Ribosomal_L17, Ribosomal_L18p, Ribosomal_L19, Ribosomal_L2, Ribosomal_L20,
Ribosomal_L21p, Ribosomal_L22, Ribosomal_L23, Ribosomal_L27, Ribosomal_L27A,
Ribosomal_L28, Ribosomal_L29, Ribosomal_L3, Ribosomal_L32p, Ribosomal_L35p,
Ribosomal_L4, Ribosomal_L5, Ribosomal_L6, Ribosomal_L9_C, Ribosomal_S10,
Ribosomal_S11, Ribosomal_S13, Ribosomal_S15, Ribosomal_S16, Ribosomal_S17,
Ribosomal_S19, Ribosomal_S2, Ribosomal_S20p, Ribosomal_S3_C, Ribosomal_S6,
Ribosomal_S7, Ribosomal_S8, Ribosomal_S9, RsfS, RuvX, SecE, SecG, SecY, SmpB,
TsaE, UPF0054, YajC, eIF-1a, ribosomal_L24, tRNA-synt_1d, tRNA_m1G_MT,
Adenylsucc_synt

OK. A lot. Good for you, Bacteria_71.

Bacteria_71 is a collection that anvi’o developers curated by taking Mike Lee’s bacterial single-copy core gene collection first released in GToTree, which is an easy-to-use phylogenomics workflow.

For the sake of this simple example, let’s assume we want to use a bunch of ribosomal genes for our phylogenomic analysis: Ribosomal_L1, Ribosomal_L2, Ribosomal_L3, Ribosomal_L4, Ribosomal_L5, Ribosomal_L6.

OK. The following command will give us all these genes from all bins described in the collection default:

anvi-get-sequences-for-hmm-hits -c CONTIGS.db \
                                -p PROFILE.db \
                                -o seqs-for-phylogenomics.fa \
                                --hmm-source Bacteria_71 \
                                -C default \
                                --gene-names Ribosomal_L1,Ribosomal_L2,Ribosomal_L3,Ribosomal_L4,Ribosomal_L5,Ribosomal_L6

Init .........................................: 4451 splits in 13 bin(s)
Hits .........................................: 668 hits for 1 source(s)
Filtered hits ................................: 64 hits remain after filtering for 6 gene(s)
Mode .........................................: DNA sequences
Genes are concatenated .......................: False
Output .......................................: seqs-for-phylogenomics.fa

If you look at the resulting FASTA file, you will realize that it doesn’t look like an alignment (I trimmed the output, you will see the full output when you look at the file yourself):

less seqs-for-phylogenomics.fa
>Ribosomal_L1___Bacteria_71___9a3cc bin_id:E_facealis|source:Bacteria_71|e_value:8.4e-59|contig:Day17a_QCcontig1|gene_callers_id:116|start:106378|stop:107068|length:690
ATGGCTAAAAAGAGCAAAAAAATGCAAGAAGCATTGAAAAAAGTTGATGCTACAAAAGCTTACTCAGTTGAAGAAGCAGTAGCTTTAGCAAAAGATACAAACATCGCGAAATTTGACGCA
(...)
>Ribosomal_L6___Bacteria_71___05162 bin_id:S_epidermidis|source:Bacteria_71|e_value:1.1e-43|contig:Day17a_QCcontig7|gene_callers_id:2325|start:13243|stop:13780|length:537
ATGAGTCGTGTTGGTAAGAAAATTATTGACATTCCTAGTGACGTAACAGTAACTTTTGACGGAAGTCATGTCACTGTAAAAGGTCCAAAAGGTGAATTAGAAAGAACTTTAAATGAAAGA
(...)
>Ribosomal_L5___Bacteria_71___ec17f bin_id:S_epidermidis|source:Bacteria_71|e_value:1.8e-31|contig:Day17a_QCcontig7|gene_callers_id:2328|start:14441|stop:14981|length:540
TTGAACCGTTTAAAAGAAAAATTTAATACAGAAGTTACTGAAAACTTAGTGAAAAAATTCAATTATAGTTCAGTGATGGAAGTACCAAAAATTGAGAAAATCGTTGTGAATATGGGTGTA
(...)

To exit less mode, press q.

Every sequence for every HMM hit is for itself :/ Hmm.

Concatenated we stand, divided we fall.

Concatenating genes

Although you can see how and why the previous output could be very useful for many other purposes, it is kinda useless for a phylogenomic analysis since we need a single concatenated alignment of gene sequences per genome.

If you look at the help menu of anvi-get-sequences-for-hmm-hits, you will see that there is a flag, --concatenate-genes, to get your genes of interest to be concatenated. Let’s do that:

anvi-get-sequences-for-hmm-hits -c CONTIGS.db \
                                -p PROFILE.db \
                                -o seqs-for-phylogenomics.fa \
                                --hmm-source Bacteria_71 \
                                -C default \
                                --gene-names Ribosomal_L1,Ribosomal_L2,Ribosomal_L3,Ribosomal_L4,Ribosomal_L5,Ribosomal_L6 \
                                --concatenate-genes

Config Error: If you want your genes to be concatenated into a multi-alignment file, you must
              also ask for the best hit (using the `--return-best-hit`) flag to avoid issues
              if there are more than one hit for a gene in a given genome. Anvi'o could have
              set this flag on your behalf, but it just is not that kind of a platform :/

Well. That didn’t go well.

The reason why it didn’t go so well is because even in the most complete genomes, there may be multiple HMM hits for a given ‘single-copy gene’. Here is an evidence for that coming from gold standard genomes for skeptics. As a solution to this problem, anvi’o asks you to use the --return-best-hit flag, which will return the most significant HMM hit if there are more than one gene that matches to the HMM of a given gene in a given genome. Fine. Let’s do that, then:

anvi-get-sequences-for-hmm-hits -c CONTIGS.db \
                                -p PROFILE.db \
                                -o seqs-for-phylogenomics.fa \
                                --hmm-source Bacteria_71 \
                                -C default \
                                --gene-names Ribosomal_L1,Ribosomal_L2,Ribosomal_L3,Ribosomal_L4,Ribosomal_L5,Ribosomal_L6 \
                                --concatenate-genes \
                                --return-best-hit

If you take a look at the output file, you can see that we are getting somewhere.

But the output is in DNA alphabet, which may not be the best option for phylogenomic analyses, especially if the genomes you have are coming from distant clades (which happens to be the case for IGD). Fortunately, you can easily switch to AA alphabet with an additional flag --get-aa-sequences (and no, there is no end to anvi’o flags, and the earlier you start getting used to the idea of reading those help menus, the sooner you will master your anvi’o game).

anvi-get-sequences-for-hmm-hits -c CONTIGS.db \
                                -p PROFILE.db \
                                -o seqs-for-phylogenomics.fa \
                                --hmm-source Bacteria_71 \
                                -C default \
                                --gene-names Ribosomal_L1,Ribosomal_L2,Ribosomal_L3,Ribosomal_L4,Ribosomal_L5,Ribosomal_L6 \
                                --concatenate-genes \
                                --return-best-hit \
                                --get-aa-sequences

If you look at the resulting file again, you will see how everything looks just so lovely. Congratulations. You did it.

In v8, anvi’o forgets to tell you that two bins in the default collection, P_acnes and S_lugdunensis, have none of the 6 requested genes, and instead silently removes them from the analysis (bad anvi’o). This is fixed in the development version and will hopefully be updated in this tutorial at the time of the next stable release. But for now, just keep it in mind since the loss of those two bins will become relevant in the next section.

Computing the phylogenomic tree

Once you have your concatenated genes, which you now have them in seqs-for-phylogenomics.fa if you followed the previous section, it is time to perform the phylogenomic analysis.

There are multiple ways to do this. Here we will use the program anvi-gen-phylogenomic-tree, which accepts a FASTA file and uses one of the programs it knows about to compute the tree. Currently the only option is FastTree, which infers approximately-maximum-likelihood phylogenetic trees from FASTA files that look like yours. Send us your favorite, and we will happily consider expanding the collection of available tools for this analysis.

Computing a phylogenomic tree from our FASTA file is as simple as this:

anvi-gen-phylogenomic-tree -f seqs-for-phylogenomics.fa \
                           -o phylogenomic-tree.txt

If anvi’o complains that you don’t have FastTree, you need to do two things. First, feel ashamed. Second, get FastTree: click me. (If you are using a mac and have brew, quietly type brew install fasttree and pretend it never happened).

The resulting file phylogenomic-tree.txt is a proper newick tree. If you want, you can visualize it immediately with anvi-interactive in manual mode:

anvi-interactive --tree phylogenomic-tree.txt \
                 -p temp-profile.db \
                 --title "Pylogenomics of IGD Bins" \
                 --manual

Which should give you this (after clicking Draw):

You can replace the colors with the bin names by selecting Text from Main > Legends > Item Names > bin_name and re-clicking Draw:

We can do much more with this phylogenomic tree of our bins than visualizing it in manual mode.

For instance, we could use it immediately to organize our bins in our collection while showing their distribution across samples.

anvi-interactive -p PROFILE.db \
                 -c CONTIGS.db \
                 -C default \
                 --tree phylogenomic-tree.txt

You will likely get an error from this command due to the loss of 2 bins that are missing the ribosomal proteins we asked for (see warning box at the end of the previous section). Those 2 bins never made it into the seqs-for-phylogenomics.fa, and therefore weren’t included in the tree. We’ll hopefully fix this section of the tutorial soon, but in the meantime please keep reading for the solution.

If you got an error from the previous command saying something like “the 11 items in your tree file do not match to the 13 items the ‘default’ collection describes”, then you can fix the mismatch between the two by making a new collection containing only the 11 genomes that made it into the phylogenomic tree. Here is how you can do that by exporting the original ‘default’ collection, removing the 2 bins missing from the tree (here done using the inverse grep command, but other strategies are possible), and importing the reduced collection back into the profile database:

anvi-export-collection -p PROFILE.db -C default
grep -v P_acnes collection-default.txt | grep -v S_lugdunensis > collection-default-reduced.txt
anvi-import-collection -c CONTIGS.db -p PROFILE.db -C default_reduced collection-default-reduced.txt

After this, you can visualize the tree like so:

anvi-interactive -p PROFILE.db \
                 -c CONTIGS.db \
                 -C default_reduced \
                 --tree phylogenomic-tree.txt

Which would give you the following display after selecting the ‘Phylogenomic tree’ from the ‘orders’ combo box in the ‘Settings’ tab (note that the name of your collection might not match to the one in the screenshot).

The tree in the middle shows the phylogenomic organization of bins we identified in the IGD.

Now you know how to organize distantly related genomes using universally conserved genes.

Chapter IV: Pangenomics

Both phylogenomics and pangenomics are strategies under the umbrella of comparative genomics, and they are inherently very similar despite their key differences. In this chapter we will discuss pangenomics and use anvi’o to have a small pangenomic analysis using our famous E. faecalis bin we recovered from the infant gut dataset and a bunch of others from the interwebs.

You can find a comprehensive tutorial on the anvi’o pangenomic workflow here.

If you haven’t followed the previous sections of the tutorial, you will need the anvi’o merged profile database and the anvi’o contigs database for the IGD available to you. Before you continue, please click here, do everything mentioned there, and come back right here to continue following the tutorial from the next line when you read the directive go back.

Please run following command in the IGD dir. They will set the stage for us to take a look at the E. faecalis bin:

anvi-import-collection additional-files/collections/e-faecalis.txt \
                       --bins-info additional-files/collections/e-faecalis-info.txt \
                       -p PROFILE.db \
                       -c CONTIGS.db \
                       -C E_faecalis

Generating a genome storage

For this example I downloaded 6 E. faecalis, and 5 E. faecium genomes to analyze them together with our E. faecalis bin. For each of these 11 external genomes, I generated anvi’o contigs databases. You can find all of them in the additional files directory:

ls additional-files/pangenomics/external-genomes/*db
additional-files/pangenomics/external-genomes/Enterococcus_faecalis_6240.db
additional-files/pangenomics/external-genomes/Enterococcus_faecalis_6250.db
additional-files/pangenomics/external-genomes/Enterococcus_faecalis_6255.db
additional-files/pangenomics/external-genomes/Enterococcus_faecalis_6512.db
additional-files/pangenomics/external-genomes/Enterococcus_faecalis_6557.db
additional-files/pangenomics/external-genomes/Enterococcus_faecalis_6563.db
additional-files/pangenomics/external-genomes/Enterococcus_faecium_6589.db
additional-files/pangenomics/external-genomes/Enterococcus_faecium_6590.db
additional-files/pangenomics/external-genomes/Enterococcus_faecium_6601.db
additional-files/pangenomics/external-genomes/Enterococcus_faecium_6778.db
additional-files/pangenomics/external-genomes/Enterococcus_faecium_6798.db

The post Accessing and including NCBI genomes in ‘omics analyses in anvi’o explains how to download sets of genomes you are interested in from the NCBI and turn them into anvi’o contigs databases.

There also are two files in the additional-files/pangenomics directory to describe how to access to the external-genomes:

and the internal one:

It is this simple to combine MAGs and isolates.

So everything is ready for an analysis, and the first step in the pangenomic workflow is to generate an anvi’o genomes storage.

anvi-gen-genomes-storage -i additional-files/pangenomics/internal-genomes.txt \
                         -e additional-files/pangenomics/external-genomes.txt \
                         -o Enterococcus-GENOMES.db

Computing and visualizing the pangenome

Now we have the genomes storage, we can characterize the pangenome:

anvi-pan-genome -g Enterococcus-GENOMES.db \
                -n Enterococcus \
                -o PAN \
                --num-threads 10

Now the pangenome is ready to display. This is how you can run it:

anvi-display-pan -g Enterococcus-GENOMES.db \
                 -p PAN/Enterococcus-PAN.db \
                 --title "Enterococccus Pan"

To get this ugly looking display:

OK. I made the previous display a bit prettier for you. If you kill the server, and import the state file the following way, and re-run the server,

anvi-import-state -p PAN/Enterococcus-PAN.db \
                  --state additional-files/state-files/state-pan.json \
                  --name default

anvi-display-pan -g Enterococcus-GENOMES.db \
                 -p PAN/Enterococcus-PAN.db \
                 --title "Enterococccus Pan"

It will look much more reasonable:

Now not only can we see how our E. faecalis genome looks like compared to available genomes, we can also see that it is not missing or carrying a great number of proteins compared to other genomes. The clustering of genomes based on gene clusters indicate that it is most similar to the genome Enterococcus faecalis 6512, which, according to the 00_INFO_ABOUT_EXTERNAL_GENOMES.txt under additional-files/pangenomics/external-genomes directory, corresponds to the assembly ID ASM17257v2 if you were to be interested in exploring further.

Adding average nucleotide identity

The pangenome tells us about the similarities and dissimilarities between those genomes given the amino acid sequences of open reading frames we identified within each one of them. We could also compare them to each other by computing the average nucleotide identity between them. Anvi’o comes with a program to do that, which uses PyANI by Pritchard et al. to pairwise align DNA sequences between genomes to estimate similarities between them. See this tutorial for details.

We can compute average nucleotide identity among our genomes, and add them to the pan database the following way:

anvi-compute-genome-similarity -e additional-files/pangenomics/external-genomes.txt \
                               -i additional-files/pangenomics/internal-genomes.txt \
                               --program pyANI \
                               -o ANI \
                               -T 6 \
                               --pan-db PAN/Enterococcus-PAN.db

anvi-import-misc-data additional-files/pangenomics/ANI_percentage_identity.txt \
                      -p PAN/Enterococcus-PAN.db \
                      --target-data-table layers \
                      --target-data-group ANI_percentage_identity

If you visualize your pangenome again once you have your ANI results in your database,

anvi-display-pan -g Enterococcus-GENOMES.db \
                 -p PAN/Enterococcus-PAN.db \
                 --title "Enterococccus Pan"

And take a careful look at your ‘Layers’ tab, you can see visualize your pangenome with the new ANI results:

Organizing gene clusters based on a forced synteny

By default, anvi’o will organize gene clusters in a given pangenome based on their distribution patterns across genomes. Which, by definition, will break gene synteny, i.e., two gene clusters that are next to each other in the display will not necessarily carry genes that are next to each other in respective genomes. If you have a relatively complete genome, you can ask anvi’o to organize gene clusters based on the organization of genes in that genome.

Here, for instance, we can force the organization of gene clusters based on the synteny of genes in our E. faecalis genome:

And here we go:

Please take a look at this study and its reproducible workflow if you would like to see a useful practical application of this strategy by Julie Reveillaud.

Functional enrichment analyses

One of the things we often are interested in is this question: which functions are associated with a particular organization of genomes due to their phylogenomic or pangenomic characteristics.

For instance, genomes in this pangenome are organized into two distinct groups based on differential distribution of gene clusters. That is of course not surprising, since these genomes are classified into two distinct ‘species’ within the genus Enterococcus. So you can certainly imagine more appropriate or interesting examples where you may be wondering about functional enrichment across groups of genomes that do not have such clear distinctions at the level of taxonomy. The actual tutorial to make sense of functions in pangenomes is here. Here we will make a quick pass to demonstrate its relevance.

First, you need to define how you would like to group your genomes in the pangenome so anvi’o can find out which functions are characteristic of each group. Here is our TAB-delimited file to divide these genomes into two groups:

Which is also in your data-pack. Let’s import it into the pan database as an additional layer data:

anvi-import-misc-data -p PAN/Enterococcus-PAN.db \
                      --target-data-table layers \
                      additional-files/pangenomics/additional-layers-data.txt

Once it is imported, we can re-run the pangenome,

anvi-display-pan -g Enterococcus-GENOMES.db \
                 -p PAN/Enterococcus-PAN.db \
                 --title "Enterococccus Pan"

and see the new layer there that correspond to our clades:

Now we can use the program anvi-compute-functional-enrichment-in-pan to ask anvi’o to identify and report functions that are enriched in either of these clades along with the gene clusters they are associated with:

anvi-compute-functional-enrichment-in-pan -p PAN/Enterococcus-PAN.db \
                                          -g Enterococcus-GENOMES.db \
                                          --category-variable clade \
                                          --annotation-source COG20_FUNCTION \
                                          -o functional-enrichment.txt

Which would generate a new file, functional-enrichment.txt, which is an anvi’o artifact called functional-enrichment-txt, in our work directory that is just filled with stuff like this:

So it turns out killer toxin insensitivity protein is only encoded by E. faecium… I am not surprised. One should expect anything from these microbes :(

Binning gene clusters

There are multiple ways to identify gene clusters that match to a given set of criteria. If you like, you can use a combination of filters that are available through the interface:

The command line program anvi-get-sequences-for-gene-clusters can also give you access to these filters and more to get very precise reports. Another option is the good’ol interactive interface, and using the dendrogram it produces to organize gene clusters based on their distribution across genomes. From this display you can make manual selections of gene clusters. I already made some selections and stored them in a file for your convenience. If you import them the following way,

anvi-import-collection additional-files/pangenomics/pan-collection.txt \
                       --bins-info additional-files/pangenomics/pan-collection-info.txt \
                       -p PAN/Enterococcus-PAN.db \
                       -C default

and re-run the interactive interface,

anvi-display-pan -g Enterococcus-GENOMES.db \
                 -p PAN/Enterococcus-PAN.db \
                 --title "Enterococccus Pan"

you will see the following selections:

We used a collection to store bins of contigs in the first section (and that’s how we identified that E. faecalis population from the Sharon et al. metagenomes anyway), and now the same concept serves us as a way to store bins of gene clusters.

If you right-click on any of the gene clusters, you will see a menu,

which, among other things, will give you access to the inspection page:

Inspecting gene clusters can be a lot of fun and very depressing at the same time. But of course it is hard to get a comprehensive understanding by going through things one by one. But can we summarize the pangenome?

Yes we can!

Summarizing a pangenome

For instance, in my tentative selection above, there is a bin called CORE ALL, which describes all gene clusters that seems to be in all genomes in this analysis. You can in fact summarize the collection default to access all the information about each gene described in each gene cluster selected as CORE ALL.

You can summarize the pangenome using the collection we have the following way:

anvi-summarize -p PAN/Enterococcus-PAN.db \
               -g Enterococcus-GENOMES.db \
               -C default \
               -o PAN_SUMMARY

If you run the the following command on a Mac system (or simply open the index.html file in your browser),

open PAN_SUMMARY/index.html

You will see a web page with some information about your pangenome:

The most important part of this output is this one:

You can unzip this file,

gzip -d PAN_SUMMARY/Enterococcus_gene_clusters_summary.txt.gz

And play with it to see how it will solve all your problems. You can import it into R, or open it in EXCEL to have a quick look at its contents. But here is a quick look at the first 10 lines of this file that contains 35,175 gene entries:

I’m sure you need no help to know what to do with this file.

Chapter V: Metabolism Prediction

If you haven’t followed the previous sections of the tutorial, you will need the anvi’o contigs databases for the E. faecalis and E. faecium genomes used in the pangenomics chapter of the tutorial. Before you continue, please click here, do everything mentioned there, and come back right here to continue following the tutorial from the next line when you read the directive go back.

A more comprehensive standalone tutorial for metabolism prediction in anvi’o can be found here. It covers more features than this chapter does, so check it out if you want more. :)

Microbiologists care a lot about what microbes are doing, and for this reason we spend a lot of time looking at functional annotations in our ‘omics data. But what if we told you that you could go a step further than that, and look at functional annotations in the larger context of metabolism? This chapter is about how to leverage known, structured information on metabolic pathways to predict the metabolic capacity encoded by microbial genomes, MAGs, or metagenomes.

In the pangenomics chapter, we explored functions in two species of Enterococcus. To demonstrate the utility of metabolism reconstruction in anvi’o, we are now going to predict the metabolic potential of these organisms.

Some obligatory background on metabolism prediction

Metabolism prediction, also known as metabolic reconstruction, is the practice of guessing (estimating) what metabolic pathways an organism can use to build or break down molecules, based on the proteins it encodes in its genome. It typically involves integrating knowledge about metabolic pathways from well-curated, publically-available databases such as KEGG or MetaCyc.

In some cases, metabolic reconstruction is a modeling approach that uses fancy math to compute a metabolic network from a genome and refine it based on experimental data in a process called Flux Balance Analysis (FBA). You cannot currently do FBA in anvi’o, but you can generate metabolic models for use with your favorite FBA tool with the program anvi-get-metabolic-model-file.

Here in this tutorial, we will cover a type of metabolism estimation called pathway prediction, which involves collecting gene annotations and matching them to defined metabolic pathways to compute a ‘completeness’ score for each pathway. Statistics about each pathway and the genes that contribute to them are then summarized in a variety of output files that can be analyzed either by reading them directly or by using them as input to downstream programs/scripts.

As of v8, anvi’o can use either KEGG or user-defined metabolic pathways as the source of metabolism information for estimation purposes. Our examples below will show how to use KEGG data, which is the default. But you can find instructions on user-defined pathways elsewhere. In particular, we recommend this tutorial as a starting point.

Estimating metabolism in the Enterococcus genomes

We’ve already prepped the Infant Gut Dataset to be ready for the metabolism estimation commands below, but in case you are working on your own dataset, please note that before you can run metabolism estimation on a fresh contigs database, you would need to first set up KEGG data on your computer using anvi-setup-kegg-data and then annotate your database using anvi-run-kegg-kofams.

Let’s get to it. To start off, we want an overview picture of what metabolisms are encoded in the Enterococcus genomes. We could visualize this nicely if we had a matrix of the completeness scores for each metabolic module (from the KEGG MODULE database) in each bin. anvi-estimate-metabolism can provide output files in matrix format if we run it in “multi-mode”, which means we will need an external genomes file. Luckily, there is already one in the datapack - the same one that is described above in the pangenomics chapter - in the additional-files/pangenomics folder.

Here is the command to run metabolism estimation on each bin, and produce matrix-formatted output:

anvi-estimate-metabolism -e additional-files/pangenomics/external-genomes.txt \
                         -O Enterococcus \
                         --matrix-format

When this program runs, it will look at the KOfam annotations (for KEGG Orthologs, or KOs) within each genome, match them up to the KEGG module definitions to estimate the completeness of each module, and produce multiple output matrices. For simplicity, we’ll use only the matrix containing pathwise completeness scores (what is pathwise completeness, you ask? Check out the technical documentation here).

We can use anvi’o to visualize the pathwise completeness matrix as a heatmap. First, we generate a newick tree from the matrix with the program anvi-matrix-to-newick:

anvi-matrix-to-newick Enterococcus-module_pathwise_completeness-MATRIX.txt

And then we use the --dry-run flag to ask anvi-interactive to give us a new profile database, import our state in it, and THEN run anvi’o interactive interface to visualize the distribution of modules across the Enterococcus genomes:

# dry run to get the profile db:
anvi-interactive -d Enterococcus-module_pathwise_completeness-MATRIX.txt \
                 -p Enterococcus_metabolism_PROFILE.db \
                 --manual-mode \
                 --dry-run

# import the state file:
anvi-import-state -s additional-files/state-files/state-metabolism.json \
                  -p Enterococcus_metabolism_PROFILE.db \
                  -n default

# run for reals:
anvi-interactive --manual-mode \
                 -d Enterococcus-module_pathwise_completeness-MATRIX.txt \
                 -t Enterococcus-module_pathwise_completeness-MATRIX.txt.newick \
                 -p Enterococcus_metabolism_PROFILE.db \
                 --title "Enterococcus Metabolism Heatmap"

Which should give us this display:

Excellent. We can already see that the E. faecalis and E. faecium genomes form two distinct groups, with the distinguishing metabolic pathways on the rightmost side of the heatmap. But what exactly are those pathways? The module numbers, which are IDs from the KEGG MODULE database, aren’t very informative. We can fix that by adding additional layers of text data describing each metabolic module. Let’s first generate a useful miscellaneous data file to import into our profile database:

# learn where the MODULES.db is:
export ANVIO_MODULES_DB=`python -c "import anvio; import os; print(os.path.join(os.path.dirname(anvio.__file__), 'data/misc/KEGG/MODULES.db'))"`

# start an empty file:
echo -e "module\tclass\tcategory\tsubcategory\tname" > modules_info.txt

# get module classes:
sqlite3 $ANVIO_MODULES_DB "select module, data_value from modules where data_name='CLASS'" | \
    sed 's/; /|/g' | \
    tr '|' '\t' >> module_class.txt

# get module names:
sqlite3 $ANVIO_MODULES_DB "select module,data_value from modules where data_name='NAME'" | \
    tr '|' '\t' > module_names.txt

# join everything
paste module_class.txt <(cut -f 2 module_names.txt ) >> modules_info.txt

# empty the trash bin:
rm module_names.txt module_class.txt

Throughout these steps we essentially have used the program sqlite3 and some terminal magic for a bit of text manipulation to generate a tab-delimited file from the anvi’o modules database that looks like this:

We can now import this file into the established profile database using the following program:

anvi-import-misc-data additional-files/metabolism/modules_info.txt \
                      -p Enterococcus_metabolism_PROFILE.db \
                      -t items

To painlessly make your interactive visualizations look like the ones in the screenshots, you can import the following state file into your profile database:

anvi-import-state -s additional-files/metabolism/metabolism_state.json \
                  -p Enterococcus_metabolism_PROFILE.db \
                  -n default

All this command does is save some settings for the interactive interface to the profile database. The settings will be automatically loaded the next time you run anvi-interactive so that you don’t have to manually change the settings yourself to get the same figure that you see in the screenshots.

Now we can run the interface command yet again to visualize the heatmap with labels. Here is the same command from above, copied here for your convenience:

anvi-interactive --manual-mode \
                 -d Enterococcus-module_pathwise_completeness-MATRIX.txt \
                 -t Enterococcus-module_pathwise_completeness-MATRIX.txt.newick \
                 -p Enterococcus_metabolism_PROFILE.db \
                 --title "Enterococcus Metabolism Heatmap"

And here is a (rotated) screenshot of the leftmost part, where we can see which pathways distinguish the two species:

You can see that E. faecalis genomes are on the left side and E. faecium ones are on the right. Based on the estimated pathwise completeness scores, it looks like the two species differ in select energy metabolism, amino acid metabolism, and vitamin/cofactor metabolism pathways. Yet many of those are faint bands indicating a low completeness score, so perhaps these genomes share a select few KOs that contribute to several of the species-specific pathways.

If you open the ‘Settings’ on the left side of the interface and choose the ‘Data’ panel, then you can hover your mouse over any part of the heatmap and see the underlying data. For example, one interesting observation is that all the E. faecalis genomes have near-complete pathways for Menaquinone Biosynthesis, Molybdenum Cofactor Biosynthesis, and Phylloquinone Biosynthesis (the completeness scores for these pathways are 88%, 80%, and 70%, respectively, in each E. faecalis genome), while these pathways are incomplete in all the E. faecium genomes (the corresponding completeness scores are 11%, 0%, and 0% in each E. faecium) genome. It would not be surprising if the same KOs from these pathways are present in each genome.

So now that we know what to look for, let’s get some more detailed metabolism estimation output. We’ll run the metabolism estimation again, but this time we will get long-format output - specifically ‘modules’ mode output, which will print information about each module in each genome, and ‘hits’ mode output, which will print information about each KEGG Ortholog that was annotated in each genome. ‘modules’ mode is the default, so if that was all we wanted, we wouldn’t need to specify it on the command line, but since we are also asking for ‘hits’ mode here we pass both of them, in a comma-separated list, to the --output-modes parameter. (Side note: you can find more details about the possible outputs of anvi-estimate-metabolism here).

anvi-estimate-metabolism -e additional-files/pangenomics/external-genomes.txt \
                         -O Enterococcus_metabolism \
                         --output-modes modules,hits

This produces two output files: Enterococcus_metabolism_modules.txt, and Enterococcus_metabolism_hits.txt. Here is a sample from the top of the modules file:

As you can see, each row contains one metabolic module in one genome, and includes for that module:

• some details about its name, categorization, and definition (all of which are sourced from KEGG, in this case)
• followed by information specific to this genome:
  • completeness score (two types, depending on how you interpret the module definition)
  • whether or not each score is above the completeness score threshold
  • how many enzymes (and which ones) were annotated that are unique to the current module
  • which enzymes (here, KEGG Orthologs) were annotated that contribute to this module
  • the gene caller IDs of those genes
  • some caveats to be aware of (in the ‘warnings’ column)

Let’s take a look at the Menaquinone Biosynthesis pathway that was differentially present between the two species. This is module M00116, so we can search for that ID number in the first column. Since we are specifically interested in which genes from this pathway were present in each genome, we will filter the output so that we only see the module number, genome name (2nd column), module definition (8th column), pathwise completeness score (11th column), list of annotated KOs (16th), and gene caller ids (17th) (scroll right to see the latter columns). Here is the command to do that search and filter from your terminal:

grep '^M00116' Enterococcus_metabolism_modules.txt | cut -f 1,2,8,11,16,17

Turns out everyone has got K02548, the second-to-last step in the module definition. What’s more, all the E. faecalis genomes have the same set of 8 KOs to make this module 88% complete (and they all have two copies of K02548). It also looks like the several of these KOs that are not shared with E. faecium are always right next to each other in the genome, because their gene caller ids are off-by-one from each other (this means that, for instance, these enzymes seem to be encoded in an operon with the following order: K19222, K01661, K01911, K02552, K02551, K08680, and K02549). All genomes are missing the enzyme for the last step in the pathway (K03183).

Let’s see what these KOs are by checking out the second file, Enterococcus_metabolism_hits.txt. That file looks like this:

In this file, each row is an annotation in one genome. The columns are hopefully self-explanatory from the headers. :) Most of the enzyme-related information will be the same across every genome, with the exception of the gene caller ID and the contig that the gene is located on. Let’s look for each of the KOs in module M00116. For simplicity, we’ll only look for the annotations in one of the E. faecalis genomes. Here is a loop to do the search for each KO number in this file:

for ko in K02552 K02551 K08680 K02549 K01911 K01661 K19222 K02548 K03183; do \
  grep $ko Enterococcus_metabolism_hits.txt | grep Enterococcus_faecalis_6255; \
done

then we see the following:

You can see that most of these enzymes are actually useful for two different metabolic pathways: M00116 and M00932. In fact, out of all of these KOs, only K02548 (which is also present in the E. faecium genomes) is unique to the menaquinone biosynthesis pathway. M00932 is the phylloquinone biosynthesis pathway that we also noticed earlier when looking at the heatmap.

So, are these E. faecalis populations capable of producing menaquinone? It appears likely, considering the 88% completeness score of the menaquinone biosynthesis module. However, we cannot know for sure given the absence of the enzyme required for the last step of the pathway. It is entirely possible that we failed to annotate this enzyme for some reason, and perhaps a more thorough search for K03183, would turn up something. Or perhaps menaquinone-producing capabilities have already been confirmed in these organisms, and a literature search would turn up evidence for this. However, KEGG’s own reference pathway for M00116 in E. faecalis indicates that this pathway is incomplete (their reference genome is actually missing two enzymes – K19222 and K03183. It is likely that we were able to annotate K19222 because anvi-run-kegg-kofams has a heuristic for finding missing annotations, which you can read about in the documentation).

But, as often happens in science, every scientist must decide for themself at what point they feel comfortable accepting results and at what point they need to go deeper to confirm them. And it is of course important to keep in mind that while genomically-encoded metabolic pathways are the first requirement for an organism to be able to perform some metabolism, even that does not necessarily mean the organism is doing so in all environments or stages of life. Perhaps anvi-estimate-metabolism is thus most reliable as a hypothesis-generating tool used to inform wet-lab experiments or to guide literature searches. It is up to you, and your science. :)

Metabolism Enrichment

Just like we looked at functional enrichment in the pangenomics chapter, we can look for enriched metabolic modules in the two Enterococcus species. In fact, we can use a program that will run the same statistical tests - except that instead of functional annotations, it will use modules which have a completeness score above the threshold that the user sets (by default, 0.75). All we need is the ‘modules’ mode output, and a file indicating which group each genome belongs to. We already have the modules file from before. Here is the groups file (which is also in the datapack):

And here is the command to run the enrichment script on modules:

anvi-compute-metabolic-enrichment -M Enterococcus_metabolism_modules.txt \
                                  -G additional-files/metabolism/entero_groups.txt \
                                  -o Enterococcus_enriched_modules.txt

We get from this a file called Enterococcus_enriched_modules.txt, in which the modules are organized from highest to lowest enrichment score. So if we look at the top 10 or so rows in the file, we will see the metabolic pathways that are most enriched in either group:

Well, it seems like there are only a few modules that could be enriched in either group, considering that the p-value and adjusted q-value jump to high, non-significant values (ie, q < 0.01) starting with the Thiamine salvage pathway. But of the 6-7 modules with high enrichment scores and low p/q-values, it seems that the E. faecium genomes are enriched with a couple of sugar degradation pathways while the E. faecalis genomes are enriched with several biosynthesis pathways - including the Menaquinone and Molybdenum cofactor biosynthesis pathways that we found earlier when looking at the heatmap. In fact, if you look back at the zoomed heatmap you can find a few more of these enriched modules and visualize the difference in completeness score between the two clades. The enrichment script is a good way to quantify these differences and assign a significance value to them.

If you are interested in learning more details about this enrichment analysis and its output, the tutorial for it is here.

Chapter VI: Microbial Population Genetics

Here we will profile the single-nucleotide variations (SNVs) in the E. faecalis bin found in Sharon et al.’s Infant Gut Dataset (IGD).

This is more of a practical tutorial for hands on experience to recover and make sense of SNVs. For a more theoretical one on the same topic, please consider first reading the tutorial Analyzing sequence variants with anvi’o.

If you haven’t followed the previous sections of the tutorial, you will need the anvi’o merged profile database and the anvi’o contigs database for the IGD available to you. Before you continue, please click here, do everything mentioned there and come back right here to continue following the tutorial from the next line when you read the directive go back.

First of all, if you haven’t already, run this command to load the collection containing the E. faecalis bin (no harm done running it twice):

anvi-import-collection additional-files/collections/merens.txt \
                       --bins-info additional-files/collections/merens-info.txt \
                       -p PROFILE.db \
                       -c CONTIGS.db \
                       -C default

Please run the following commands in the IGD dir. They will set the stage for us to take a look at the E. faecalis bin:

# importing taxonomy for gene calls
anvi-import-taxonomy-for-genes -c CONTIGS.db \
                               -i additional-files/centrifuge-files/centrifuge_report.tsv \
                                  additional-files/centrifuge-files/centrifuge_hits.tsv \
                               -p centrifuge

# importing the state file so things look pretty
anvi-import-state --state additional-files/state-files/state-merged.json \
                  --name default \
                  -p PROFILE.db

OK. Let’s first remind ourselves where this bin is.

You remember the E. faecalis bin from previous sections of this tutorial:

The red selection in the most outer layer represents the E. faecalis bin, which is very abundant in every sample, and it is stored in the collection merens under the bin name E_facealis (yes, the name has a typo, we know, it is all Tom’s fault). In this section we will focus on this bin, and in fact we don’t really need the rest of this metagenome.

Since we have identified the genome bin of focus, we can set aside the rest of the metagenome, and ‘split’ that bin out of this profile database to create a self-contained anvi’o profile that only describes this single bin. This is also a great way to share your research with others. Here is how we can get that bin out of this metagenome:

anvi-split -p PROFILE.db \
           -c CONTIGS.db \
           -C default \
           -b E_facealis \
           -o MAGs

Once it is done, this is what you will see in the resulting directory:

ls MAGs/
E_facealis
ls MAGs/E_facealis/
AUXILIARY-DATA.h5 CONTIGS.db        CONTIGS.h5        PROFILE.db

Look familiar? This generated a new anvi’o profile database and anvi’o contigs database that pertains only to the E_facealis bin. While the program anvi-split does its magic, it also adds into the resulting profile database a collection (called DEFAULT), with a single bin in it (called ALL_SPLITS) that contains all the contigs in this MAG. This way we can access to its contents from all anvi’o programs. Here is an example:

anvi-estimate-genome-completeness -p MAGs/E_facealis/PROFILE.db \
                                  -c MAGs/E_facealis/CONTIGS.db \
                                  -C DEFAULT

╒════════════╤══════════╤══════════════╤════════════════╤════════════════╤══════════════╤════════════════╕
│ bin name   │ domain   │   confidence │   % completion │   % redundancy │   num_splits │   total length │
╞════════════╪══════════╪══════════════╪════════════════╪════════════════╪══════════════╪════════════════╡
│ ALL_SPLITS │ BACTERIA │            1 │            100 │           4.23 │          140 │        2865861 │
╘════════════╧══════════╧══════════════╧════════════════╧════════════════╧══════════════╧════════════════╛

OK. If you were to visualize this bin,

anvi-interactive -p MAGs/E_facealis/PROFILE.db \
                 -c MAGs/E_facealis/CONTIGS.db

this is what you should get:

Expert question. Do you see that one split that has a different taxonomy? Who does not think it is not contamination? Why?

wget https://ndownloader.figshare.com/files/8248433 \
   -O TARA_ANW_MAG_00006.tar.gz

tar -zxvf TARA_ANW_MAG_00006.tar.gz

# we had compressed the auxiliary data to save space,
# you need to unzip it:
gzip -d TARA_ANW_MAG_00006/AUXILIARY-DATA.h5.gz

# if you are not using anvi'o v2.3.0, you need to migrate
# anvi'o databases to your version:
ANVIO_SAMPLES_DB=SKIP anvi-migrate TARA_ANW_MAG_00006/*.db --migrate-safely


anvi-interactive -p TARA_ANW_MAG_00006/PROFILE.db \
                 -c TARA_ANW_MAG_00006/CONTIGS.db

Alright.

Here is a different representation of the coverage of this bin across samples (which is a screenshot from our paper):

See, it really is abundant (every dot here is the coverage of a nucleotide position that was reported as a variable position).

OK. Clearly, we have no way of knowing the extent of variation within this bin through this perspective. But anvi-gen-variability-profile is exactly for that purpose, and that’s what we will do here using two different methods to visualize its report (using R, and using anvi’o).

Let’s first generate the SNV profile output file, details of which were described extensively here. Here it is:

anvi-gen-variability-profile -c MAGs/E_facealis/CONTIGS.db \
                             -p MAGs/E_facealis/PROFILE.db \
                             -C DEFAULT \
                             -b ALL_SPLITS \
                             --samples-of-interest additional-files/samples.txt \
                             --min-coverage-in-each-sample 20 \
                             --min-occurrence 3 \
                             --include-split-names \
                             --quince-mode \
                             -o E-faecalis-SNVs.txt

This command simply requests anvi-gen-variability-profile to select all E. faecalis nucleotide positions that were identified as ‘variable’ in at least one sample, and (1) covered by more than 20X in all 8 samples of interest (2) and display a minimum scattering power of 3 (minimum scattering power is a very simple attribute of a nucleotide position, and described here). Instead of --min-scatter 3 you could use --min-occurrence 3, however, in that case the program would have returned every SNV position that was reported in more than 3 samples, which would have included nucleotide positions that were variable in every sample. Finally we also used the --samples-of-interest parameter with the following file:

cat additional-files/samples.txt
DAY_15B
DAY_16
DAY_17B
DAY_18
DAY_19
DAY_22A
DAY_23
DAY_24

which contains only the sample names with larger number of reads for every duplicate day to simplify things.

OK. Our beloved anvi-gen-variability-profile will kindly store results from 466 nucleotide positions in the output file, E-faecalis-SNVs.txt. Because we used the flag --quince-mode, there will be a total of 3,728 entries in the file (=466 x 8).

At this point we now know that the E. faecalis is not a monoclonal bin, and does maintain some level of heterogeneity. What we don’t know is whether it is just noise, or if there is any signal that may tell us something about the nature of this variability.

Now we will visualize the information two different ways.

Visualizing SNV profiles using R

First, we will use R to recapitulate what we did in our paper. Of course, please feel absolutely free to look at the Figure 3 behind that link if you are fine with RUINING THE SURPRISE :(

This section is here due to historical reasons. This is not our preferred way of visualizing SNV data.

For this step we will need this R script Meren had written before. You can download it the following way:

wget https://raw.githubusercontent.com/meren/anvio-methods-paper-analyses/a57b0cee07e9dd6fc59892114f2ad5bb9df78215/SHARON_et_al/VARIABILITY_REPORTS/02_GEN_FIGURE_SUMMARY.R \
      -O visualize-SNVs.R
chmod +x visualize-SNVs.R

Here is an extra step, and why we need it: E-faecalis-SNVs.txt contains filtered SNVs, but we also want to know about each sample’s variation density, the number of nucleotide positions reported as a variable position for each 1,000 nts. The R script can already do it, but we need another file that reports ALL SNVs for every sample, without any filters to get the raw counts. Let’s first create that the following way:

anvi-gen-variability-profile -c MAGs/E_facealis/CONTIGS.db \
                             -p MAGs/E_facealis/PROFILE.db \
                             -C DEFAULT \
                             -b ALL_SPLITS \
                             --samples-of-interest additional-files/samples.txt \
                             -o E-faecalis-SNV-density.txt

Now we can use the R script to visualize this information the following way,

This R script will require some libraries to be installed. You can install all of them by typing R in your terminal, and then entering this command: install.packages(c('ggplot2', 'reshape2', 'reshape', 'gridExtra', 'grid', 'plyr', 'gtools')). After everything installs, exit R by typing q(). My R version that worked with this script was v3.2.2.

./visualize-SNVs.R E-faecalis-SNVs.txt E-faecalis-SNV-density.txt 150 2870000

where 150 is the number of random positions we want to show in the heatmap, and 2780000 is the genome size to calculate the SNV density per sample. If you have all necessary R libraries installed, you should see a new PDF file in the directory that looks like this:

With some rearrangement and polishing using Inkscape, you can see how this output fit into our Figure 3:

The figure shows that the variation density changes quite dramatically from one day to another, despite the rather stable coverage. Plus, the emergence of this pattern is not random: the heatmap shows that the nucleotide positions that show variation re-occur, and competing base identities remains the same.

Investigating what causes this, is of course when things start to get exciting. However, we will not go there. Instead, we would like to leave you with this thought: by using patterns of variability, we can start characterizing changes in microbial population structures across environments, and generate better-informed hypotheses to investigate mechanisms that drive these shifts.

Visualizing SNV profiles using anvi’o

R visualization is useful, but the heatmap in that figure can’t effectively visualize more than a couple hundred positions. That’s why there is a random subsampling step. But we can use the anvi’o interactive interface to display up to 25,000 nucleotide positions easily.

For this, we have a little program called anvi-script-snvs-to-interactive which reads an anvi’o variability profile, and creates an anvi’o interactive interface-compatible output from it.

So when you run this script the following way,

anvi-script-snvs-to-interactive E-faecalis-SNVs.txt -o e_faecalis_snvs

it will do its magic, and create an output directory with material that can directly be used with anvi-interactive with the --manual flag.

A little note for people who are interested in programming: Feel free to take a look at the relevant line of the source code of this script to see how easy it is to generate an anvi’o-compatible visualizable output from any TAB-delimited matrix file.

If you run the interactive interface on these results the following way,

anvi-interactive --profile e_faecalis_snvs/profile.db \
                 --tree e_faecalis_snvs/tree.txt \
                 --view-data e_faecalis_snvs/view.txt \
                 --title "SNV Profile for the E. faecalis bin" \
                 --manual

You will get this view:

This view can definitely be improved. I prepared a state file to match colors of competing nucleotides to the R results. If you import that state file and run the interactive interface the following way,

anvi-import-state -p e_faecalis_snvs/profile.db \
                  --state additional-files/state-files/state-snvs.json \
                  --name default

anvi-interactive -d e_faecalis_snvs/view.txt \
                 -t e_faecalis_snvs/tree.txt \
                 -p e_faecalis_snvs/profile.db \
                 --title "SNV Profile for the E. faecalis bin" \
                 --manual

This time you will get this display:

As we’ve seen before, occurrence of SNVs follow a bi-daily fashion. Not that it needs any further convincing, but just to show off here, if you were to click Main > Layers > Order by > Default(tree) > Draw, you can see that even days and odd days nicely separate from each other:

For instance while I was on this page, I selected two SNV positions that showed different patterns:

Position 80 (the selection towards the right) shows variability only in odd days, position 686 shows variability only in even days. This will probably at some point will become a right-click menu function, but even now we can do some tricks to explore the context of these SNVs. Because that is what anvi’o is all about. Exploring, sometimes to disturbing depths. Well, here are two tiny AWK one-liner to get some specific information about these positions from our files. Just to have a bit more insights into these SNVs, here is one nucleotide position that is varible only in odd days:

awk '{if(NR == 1 || $2 == 142) print $2 "\t" $6 "\t" $14 "\t" $22 "\t" $17 "\t" $18 "\t" $19 "\t" $21 "\t" $23 "\t" $24}' E-faecalis-SNVs.txt | anvi-script-as-markdown

And here is another one that is only variable in even days:

awk '{if(NR == 1 || $2 == 186) print $2 "\t" $6 "\t" $14 "\t" $22 "\t" $17 "\t" $18 "\t" $19 "\t" $21 "\t" $23 "\t" $24}' E-faecalis-SNVs.txt | anvi-script-as-markdown

Good, everything checks out. Now since we know the split names and positions in splits, we can in fact see where they actually are using the interactive interface to visualize the merged profile database again, and look at the wider context using the ‘inspect’ option. Which I have already done for you:

There are many directions you can go once you have the gene caller IDs associated with a question you have. Just take a look at this post and see some of the hidden talents of anvi’o: Musings over a Nitrospira genome that can do complete nitrification.

Here I will stop, but still we have a lot to talk about!

Visualizing SNV profiles as a network

You can ALSO generate an XML description of the SNV profiles you have generated using anvi-gen-variability-profile program, using the program anvi-gen-variability-network:

anvi-gen-variability-network -i E-faecalis-SNVs.txt \
                             -o E-faecalis-SNVs.gexf

You can use Gephi to play with the resulting file to visualize or to analyze the network properties of nucleotide positions further.

Here is a screenshot from Gephi for SNV profiles in the E. faecalis genome bin:

Networks. Everyone likes networks.

Measuring distances between metagenomes with F_ST

The fixation index, F_ST, is a measure of differentiation revealed by SNVs between two populations, and is a variation of the most widely used approach to describe population structures, F-statistic. Anvi’o includes a program, anvi-gen-fixation-index-matrix, which implements F_ST in accordance with Schloissnig et al.’s 2013 work to permit multi-allelic variant positions.

As an alternative to previous strategies, we can use F_ST to estimate distances between metagenomes using the SNV data following the footsteps of Tom Delmont and Evan Kiefl et al’s previous work by providing this program with the variability profile anvi’o generates:

anvi-gen-fixation-index-matrix --variability-profile E-faecalis-SNVs.txt \
                               --output-file FST_E_facealis.txt

The resulting file is a distance matrix, in which 0 indicates maximum similarity, and 1 indicates maximum distance between two metagenomes as a function of, in our case, SNVs they have in the context of a given genome. We can quickly turn that distance matrix into a hierarchical clustering dendrogram,

anvi-matrix-to-newick FST_E_facealis.txt \
                      --output-file FST_E_facealis.newick

And visualize the output in anvi’o,

anvi-interactive -t FST_E_facealis.newick \
                 -p FST_E_facealis.db \
                 --manual

And voila:

Even though F_ST is quite a sane way to investigate associations between metagenomes, it is important to look at the entirety of the SNV profile by first looking at it using the second method.

OK.

You just read about four different ways to visualize the single-nucleotide variability profiles with anvi’o. We are done here.

Chapter VII: Genes and genomes across metagenomes

In the previous sections of this tutorial we recovered an E. faecalis genome from an infant gut metagenome.

Increasing availability of shotgun metagenomes offer quite powerful ways to characterize MAGs and isolate genomes across different systems. Harnessing the power of available metagenomes, in fact we could ask the question “how does our E. faecalis MAG look like in other infant gut metagenomes”. To provide a preliminary answer to this question, and how similar questions can be answered using in anvi’o, here we will use another Banfield group study by Brooks et al.

In your data-pack you have a merged anvi’o profile that describes the read recruitment statistics of our E. faecalis MAG across 20 gut metagenomes from the Brooks et al. study.

{
    "samples_txt": "samples.txt",
    "fasta_txt": "fasta.txt",
    "references_mode": true,
    "import_percent_of_reads_mapped": {
        "run": false
    },
    "bowtie": {
        "additional_params": "--no-unal",
        "threads": 4
    },
    "anvi_profile": {
        "--overwrite-output-destinations": true,
        "--min-contig-length": 0,
        "--profile-SCVs": true,
        "--skip-SNV-profiling": false,
        "threads": 4
    },
    "anvi_script_reformat_fasta": {
        "run": true
    },
    "anvi_run_ncbi_cogs": {
        "run": true,
        "threads": 15
    },
    "output_dirs": {
        "FASTA_DIR": "02_FASTA",
        "LOGS_DIR": "00_LOGS",
        "QC_DIR": "01_QC",
        "PROFILE_DIR": "05_ANVIO_PROFILE",
        "MERGE_DIR": "06_MERGED",
        "MAPPING_DIR": "04_MAPPING",
        "CONTIGS_DIR": "03_CONTIGS"
    },
    "anvi_run_hmms": {
        "run": false
    },
    "iu_filter_quality_minoche": {
        "run": false
    },
    "anvi_merge": {
        "--skip-concoct-binning": true
    }
}

anvi-run-workflow -w metagenomics \
                  -c config.json

You can visualize the merged profile in your data-pack by typing the following command,

anvi-interactive -p additional-files/e_faeealis_across_hmp/PROFILE.db \
                 -c additional-files/e_faeealis_across_hmp/CONTIGS.db

This output simply puts our contigs in the E. faecalis genome we recovered from Sharon et al.’s infant gut in the context of 20 other infant gut metagenomes.

The gene mode: studying distribution patterns at the gene-level

Furthermore, one could want to do it per gene basis for a more higly resolved understanding of distribution patterns. Luckily, if you have a collection and a bin in your profile database, anvi’o allows that.

If you don’t have a collection, anvi-script-add-default-collection will help you add one.

The ‘gene mode’ can be invoked using the following command:

anvi-interactive -p additional-files/e_faeealis_across_hmp/PROFILE.db \
                 -c additional-files/e_faeealis_across_hmp/CONTIGS.db \
                 -C DEFAULT \
                 -b EVERYTHING \
                 --gene-mode

which should produce the following output, where the coverage of each gene is independently visualized in the context of other metagenomes:

If you right click any of the genes, you can see that this menu has more options,

Including the inspection of a given gene in its context with upstream and downstream coverage information:

Chapter VIII: From single-amino acid variants to protein structures

This section of the tutorial uses the data described in the above section Putting the E. faecalis genome in the context of HMP gut metagenomes. You already have this data, but if you’re wondering “where did this come from?”, you should read that section and then come back.

This section briefly showcases 2 separate, yet integrated features. The first is binding site prediction prediction with InteracDome, and the second is protein structure visualization with anvi’o structure. This is a showcase of these features, not an in-depth tutorial on either topic. For more comprehensive material on InteracDome, check out this blog. For more detailed info and in depth tutorial on anvi’o structure and what it can do for you, please refer to this post.

Predicting binding sites with InteracDome

This is considered an experimental feature at this point in time. It has not been sufficiently tested, yet we are very excited to experiment with it and make it available to you to experiment with. Please report any bugs you may find.

In the previous chapters of this tutorial we used single-nucleotide variants (SNVs) to explore sequence heterogeneity of an E. faecalis population. This demonstrated the ease with which SNVs can be explored and characterized, however offered zero insight into their potential impact on fitness. This is an extremely difficult task, yet as a step towards this direction, the first section of this chapter is dedicated to exploring how we can predict binding sites using InteracDome via the program anvi-run-interacdome, and thereby attribute more meaning and context to observed distributions of SNVs.

For this demo, we once again will use the E. faecalis genome as a medium to estimate per-residue binding frequencies. We start (and end) the process with the following command:

In order for the command below to work, you will need to have at some point ran anvi-setup-interacdome, which downloads all necessary files. This will require an internet connection, ~250MB of free space, and about 10 minutes. So if the command didn’t work, first type anvi-setup-interacdome and let anvi’o loose.

anvi-run-interacdome -c additional-files/e_faeealis_across_hmp/CONTIGS.db -T 2

As you may be getting used to, -c specifies the contigs database you would like to predict binding sites for, and -T means to run the command with 2 threads (feel free to use more if you have the hardware).

In brief, running this command searches each gene in your contigs database against a subset of the Pfam database that have been annotated with per-residue “binding frequencies”. These binding frequencies are values between 0 and 1 which indicate how likely it is a residue is involved in interacting with a ligand. The search is ran with HMMER and the binding frequencies are from the InteracDome.

If you want to understand the inner workings as to how the InteracDome dataset was created, you should read this paper which introduces the InteracDome resource.

After running the command, you will see a bunch of output. Let’s take a look at the output, one chunk at a time. Note that if you are unsatisfied with any of my descriptions of these outputs (Evan speaking here), check out my implementation diary of this anvi’o feature for way more info. Here is the first chunk of output:

CITATION
===============================================
Anvi'o will use 'InteracDome' by Kobren and Singh (DOI: 10.1093/nar/gky1224) to
attribute binding frequencies. If you publish your findings, please do not
forget to properly credit their work.

InteracDome is the resource published by Kobren and Singh, that has annotated per-residue binding frequencies to Pfams. Without them, this would not be possible, so don’t forget to cite. Next:

INITIALIZATION
===============================================
Interacdome dataset used .....................: representable
Minimum hit fraction .........................: 0.5
Pfam database directory ......................: /Users/evan/Software/anvio/anvio/data/misc/Interacdome
Pfam database version ........................: 31.0 (2017-02)
Contigs DB ...................................: Initialized: additional-files/e_faeealis_across_hmp/CONTIGS.db (v. 18)
HMM profiles loaded from .....................: /Users/evan/Software/anvio/anvio/data/misc/Interacdome/Pfam-A.hmm
Number of HMM profiles loaded ................: 2375 profiles
num genes that HMM will be run on ............: 2,754
Output .......................................: /var/folders/58/mpjnklbs5ql_y2rsgn0cwwnh0000gn/T/tmpz7zrm92p/AA_gene_sequences.fa

This simply provides some information about where the Pfams are located, which contigs database is being used, as well as the amino acid sequences exported from the contigs database. It also lets us know that we will be attempting to compute binding frequencies for 2754 E. faecalis genes since that’s how many genes in the contigs database had amino acid sequences. To compute these binding frequencies, we will use an InteracDome database of 2375 Pfams that have annotated binding frequencies. This number comes from the fact that we are using the ‘representable’ subset of Pfams, which “correspond[s] to domain-ligand interactions that had nonredundant instances across three or more distinct PDB structures. [Kobren and Singh] recommend using this collection to learn more about domain binding properties”. If you want to include only highly confident interactions that have been cross-validated, supply the flag --interacdome-dataset confident. Next up:

HMM Profiling for InteracDome
===============================================
Reference ....................................: unknown
Kind .........................................: unknown
Alphabet .....................................: AA
Context ......................................: DOMAIN
Domain .......................................: N/A
HMM model path ...............................: /Users/evan/Software/anvio/anvio/data/misc/Interacdome/Pfam-A.hmm
Number of genes in HMM model .................: 2,375
Noise cutoff term(s) .........................: --cut_ga
Number of CPUs will be used for search .......: 2
HMMer program used for search ................: hmmsearch
Temporary work dir ...........................: /var/folders/58/mpjnklbs5ql_y2rsgn0cwwnh0000gn/T/tmpeqd9lraw
Log file for thread 0 ........................: /var/folders/58/mpjnklbs5ql_y2rsgn0cwwnh0000gn/T/tmpeqd9lraw/AA_gene_sequences.fa.0_log
Log file for thread 1 ........................: /var/folders/58/mpjnklbs5ql_y2rsgn0cwwnh0000gn/T/tmpeqd9lraw/AA_gene_sequences.fa.1_log

Rather boringly, this output just lets us know how we are running HMMER. Next:

HMMER results
===============================================
Loaded HMMER results from ....................: /var/folders/58/mpjnklbs5ql_y2rsgn0cwwnh0000gn/T/tmpeqd9lraw/hmm.standard
num total domain hits ........................: 2,898
num unique genes .............................: 1,521
num unique HMMs ..............................: 956
filtered hits (no InteracDome data) ..........: 0
filtered hits (disagrees with conserved match states) : 1
filtered hits (hit is too partial hit) .......: 261

These are the results from running HMMER. 2898 domain hits were found for 1521 genes. That means that about half the E. faecalis genome had at least one domain hit, and that the average number of domain hits for these genes was about 2. There were 956 Pfams (HMMs) that contributed domain hits, which is about 40% of the InteracDome “representable” dataset. The remaining lines indicate some filtering info. Most notably, 261 hits were removed because they corresponded to domain hits that had Pfams that aligned less than half of their length to the E. faecalis gene. You can change this fraction threshold with the parameter --min-hit-fraction, which by default is 0.5.

Next up:

InteracDome Results
===============================================
* (Pre-filter) 1488 genes with at least 1 binding freq
* (Pre-filter) 728542 total position-ligand combos with binding freqs
    - position-ligand combos with binding freqs in range [0.0,0.1): 458055
    - position-ligand combos with binding freqs in range [0.1,0.2): 96063
    - position-ligand combos with binding freqs in range [0.2,0.3): 49816
    - position-ligand combos with binding freqs in range [0.3,0.4): 35944
    - position-ligand combos with binding freqs in range [0.4,0.5): 22678
    - position-ligand combos with binding freqs in range [0.5,0.6): 18016
    - position-ligand combos with binding freqs in range [0.6,0.7): 16228
    - position-ligand combos with binding freqs in range [0.7,0.8): 10809
    - position-ligand combos with binding freqs in range [0.8,0.9): 7704
    - position-ligand combos with binding freqs in range [0.9,1.0]: 13229

Filter by minimum binding frequency ..........: True
Minimum binding frequency ....................: 0.2

* (Post-filter) 1488 genes with at least 1 binding freq
* (Post-filter) 174424 total position-ligand combos with binding freqs
    - position-ligand combos with binding freqs in range [0.0,0.1): 0
    - position-ligand combos with binding freqs in range [0.1,0.2): 0
    - position-ligand combos with binding freqs in range [0.2,0.3): 49816
    - position-ligand combos with binding freqs in range [0.3,0.4): 35944
    - position-ligand combos with binding freqs in range [0.4,0.5): 22678
    - position-ligand combos with binding freqs in range [0.5,0.6): 18016
    - position-ligand combos with binding freqs in range [0.6,0.7): 16228
    - position-ligand combos with binding freqs in range [0.7,0.8): 10809
    - position-ligand combos with binding freqs in range [0.8,0.9): 7704
    - position-ligand combos with binding freqs in range [0.9,1.0]: 13229

The following shows a text-based histogram of the binding frequency scores attributed to the E. faecalis genes (Reminder: the InteracDome dataset attributes to each residue a score between 0 and 1, where 0 means unlikely and 1 means likely to be involved in binding a ligand). The first section shows before filtering low-frequency hits, and the second shows after filtering. As we can see, most binding frequencies attributed fall under 0.2, and by default there is a filter cutoff of 0.2. You can change this with --min-binding-frequency.

Finally, there is:

* Binding frequencies for 49195 unique positions successfully stored in
additional-files/e_faeealis_across_hmp/CONTIGS.db
* Domain hit summaries stored in INTERACDOME-domain_hits.txt
* Match state contributors stored in INTERACDOME-match_state_contributors.txt

After all filtering, the results are in: 49195 residues in the E. faecalis genome have been implicated in ligand binding. The positions, ligand types, and binding frequencies are stored in the contigs database, and a slew of additional information is stored in the text outputs INTERACDOME-domain_hits.txt and INTERACDOME-match_state_contributors.txt.

With the addition of this information, there is seemingly an uncountable number of directions to go, and the truth is that we ourselves have not even scratched the surface. We therefore can’t predict how this may be useful to you. I, Evan, have some ideas of my own, but I haven’t even had time to explore them. That’s how fresh these cookies are. That said, in the next section we will transition to talking about anvi’o structure, a tool to visualize protein structures, and as you may have guessed, there will be a section in which we visualize the predicted binding sites we calculated in this section.

Visualizing protein structures, SAAVs, SCVs, and binding sites with anvi’o structure

So far we have been referring to sequence variants as SNVs, yet when analyzing sequence variation in the context of proteins, it makes a lot more sense to characterize sequence variation with single codon variants (SCVs) or single amino acid variants (SAAVS) rather than SNVs. If you are unfamiliar with these concepts, or want to learn more, click me. Point is, we will be working with SAAVs and SCVs moving forward.

In this section we will introduce anvi’o structure, an integrated component of anvi’o that enables one to predict protein structures for genes in your contigs database and visualize them in an interactive interface. The motivation was seeded by wanting more direct insight into the potential effect that sequence variants had on fitness, by visualizing them directly on the structure, since it is the physical properties of these structures that define function.

The structure database

To analyze and visualize where nucleotide-level variation occurs in the three-dimensional structure of an encoded protein, you need a protein structure. For this, anvi’o structure has a dedicated structure database that stores all structural information requested by the user. Users can create their own structure database by providing a contigs database, and a list of gene IDs in the contigs database they want to predict structures for.

Very briefly, structures are predicted using MODELLER, which attempts to find homologous proteins with experimentally-solved structures. If any exist, they are used as templates to model your proteins of interest.

You already have a pre-made structure database in the data-pack, and you can make sure it exists with the following command:

ls additional-files/e_faeealis_across_hmp/STRUCTURE.db

anvi-gen-structure-database -c additional-files/e_faeealis_across_hmp/CONTIGS.db \
                            --max-number-templates 1 \
                            --gene-caller-ids 58,620,2105,2371,2602 \
                            -o additional-files/e_faeealis_across_hmp/STRUCTURE_NEW.db

Visualizing proteins with anvi-display-structure

With a structure database in hand, you are ready to visualize SAAVs and SCVs on protein structures.

anvi-display-structure -s additional-files/e_faeealis_across_hmp/STRUCTURE.db \
                       -p additional-files/e_faeealis_across_hmp/PROFILE.db \
                       -c additional-files/e_faeealis_across_hmp/CONTIGS.db \
                       --min-departure-from-consensus 0

After running the above command, you should see the following view.

This is the predicted structure for gene 58. The left panel indicates that it is a Hydroxyethylthiazole kinase of the sugar kinase family, and the template protein used to model the structure (PDB code 3DZV) was crystallized from an E. faecalis isolate by the Joint Center for Structural Genomics. The template shares 96.5% sequence identity with our gene, and so we are very confident in the accuracy of the predicted structure. On the display, each red sphere indicates the position of a SAAV that occurred in at least one of the metagenomes. You can learn more about each SAAV by hovering the mouse above it. Try and find the SAAV occurring at the 87th residue. To help, you can click the Views tab at the top, find the Backbone Color section, and then select the Rule to be residueindex. The SAAV at position 87 is orange in color.

Under ‘Reference info’ you can see that the MAG had an asparagine at this position, however under ‘Variant frequencies’, 88% of reads were threonine within the infants sampled. If you inspect the amino acid frequencies of more of these SAAVs, you will notice that most of them have insignificant levels of variation. These may be true biological variation in the environment, but they are at such low frequencies they are indistinguishable from sequencing error. To avoid fooling ourselves, we can filter out these reads by going to the Filters tab, and conservatively setting the minimum departure from reference value to 0.05:

(I turned off the translucent surface by going to Views|Structure|Surface|Show and unchecking the box)

To get a better idea of the characteristics of each SAAV, you can set the color and size for each according to various metrics of interest in the Perspectives tab. In the animation below, we color each SAAV according to its entropy to quantify the degree of variation at each site (click here and search for ‘entropy’ for an explanation):

(I turned off rainbow coloring by going to Views|Structure|Backbone Color|Rule and selecting Static)

The 4 white-ish spheres are sites in which the metagenomes disagree with the MAG, but where there is little to no variation within the environment. In contrast, the 3 red-ish spheres indicate sites in which the environment has a mixed population of 2 or more amino acids. Verify this for yourself by looking at the variant frequencies and comparing their identities to the reference residue.

SAAVs and SCVs

So far we’ve focused exclusively on SAAVs, however we can also study variation with single codon variants (SCVs). You can visualize SCVs by clicking CDN (read: codon) in the Main tab:

In this view, the first thing you will notice is that there are many more SCVs than SAAVs. That’s because unlike SAAVs, SCVs also report on synonymous variation. For example, this residue is completely synonymous, but impressively houses 4 of the 6 arginine codons:

If you inspect enough of these genes, you’ll start to notice that SCVs at solvent inaccessible residues (relative solvent accessibility = 0) are almost always entirely synonymous (synonymity = 1). Do you have any explanation for why?

Case study: sequence variation within a kinase binding pocket

Since gene 58 is a kinase, it may not surprise you it has a binding pocket for ADP. Remember when we ran anvi-run-interacdome and predicted ligand sites? Let’s visualize where this ADP pocket is. First, let’s navigate to the Views tab, turn on the surface, and increase the opacity to 1. Then color the surface dynamically according to ADP binding frequency values. Your settings should mimic this:

Here is what you should be seeing:

Rather gorgeously, we can see where the ADP is predicted to bind. Since binding pockets are meticulously constrained to maintain optimal binding kinetics, they are under extreme purifying selection. To see whether are there are any SAAVs in the binding pocket, I turned off the surface and changed the settings so that the backbone color represents ADP binding frequency, where purple represents higher frequencies:

Here I’ve zoomed into the binding pocket, where you can see there are no SAAVs in “site”.

However, when we switch to a view of SCVs, the pocket lights up like a Christmas tree.

Visually it’s clear that only variants that are fully synonymous are permissible in the binding pocket. This is a fine example of how non-synonymous sequence evolution is usually more highly constrained than synonymous sequence evolution due to fitness-decreasing alterations in the biophysical properties of the encoded protein. Since these synonymous SCVs do not influence the binding site, they have no fitness effect on the protein’s ability to function. It’s impossible to know without conducting fitness experiments, but most likely they are neutral mutations, and have entered the population through drift and/or draft.

Grouping metagenomes

You now know the basics. Let’s talk about how to get the most out of this interface. Sometimes it is desirable to partition metagenomes into groups based on the experimental design. For example, these 20 metagenomes are from 7 infants, and it would be useful to group which metagenomes belong to which infant. Unless the user provides custom categories for groupings, only two are by default available: merged and samples. Currently you are using the merged category, which is why you see all variant data merged onto a single protein view. Go to the Main tab and select samples from within the Category dropdown menu. You should see this:

Now you can visualize each metagenome’s variants individually. Let’s import some metadata so we can create our own custom groupings. Return to your terminal and import the additional data table into the profile database:

anvi-import-misc-data additional-files/e_faeealis_across_hmp/additional_layers_table.txt \
                      -p additional-files/e_faeealis_across_hmp/PROFILE.db \
                      -t layers

This is just a TAB-delimited table that is compatible with anvi’o misc data table inputs:

Afterwards, we can rerun the interactive display:

anvi-display-structure -s additional-files/e_faeealis_across_hmp/STRUCTURE.db \
                       -p additional-files/e_faeealis_across_hmp/PROFILE.db \
                       -c additional-files/e_faeealis_across_hmp/CONTIGS.db \
                       --min-departure-from-consensus 0

Now, there are many more grouping categories available to us.

cohort_infant categorizes metagenomes into the infants from which they were sampled. Select it to reveal one protein view for each of the 7 infants. Now the display can inform us about subtle differences in SAAVs between each infant.

We hope this tutorial give you enough motivation to explore your own data with this framework! Feel free to click on the other 4 genes in the Main tab and further explore, or better yet, do it with your own metagenomes.

Let’s face it: it is a lot easier to live in a world of A, C, T, and G. But in this restricted space we can’t explore the biochemical consequences of our sequences, and we therefore can’t understand why DNA sequences are the way they are. Anvi’o structure has been designed to help make this undertaking more accessible, but it is only in its infancy. If you have suggestions, please let us know.

Final words

Do you want more examples? Do you have questions? Please don’t hesitate to get in touch with us!