Metabin refinement and population genetics using human tongue metagenomes

Iva Veseli

Email Twitter Github LinkedIn ORCiD

Florian Trigodet

Email Twitter Github LinkedIn ORCiD

Real world datasets are complicated, and it can be difficult to accurately capture this in tutorials due to time and space constraints. The infamous Infant Gut Tutorial covers the basics of both binning and population genetics on an extremely simple dataset, which is great for learning but the process and results are not necessarily representative of what one would actually end up doing/seeing in a large-scale ‘omics study. Take binning, for example – if you have dozens of large co-assemblies in your dataset, you won’t have the time to do manual binning (and even if you wanted to, chances are that your dataset would be too large to visualize in the anvi’o interactive interface). And microbial communities containing hundreds of populations rarely show perfectly clean, easy-to-bin patterns of sequence composition and differential coverage. As for population genetics, our existing tutorial shows you how to do the thing, but leaves out some context for the question and important parameter decisions along the way.

With that in mind, we created this tutorial to show a more realistic example of how to work on both metagenomic binning and population genetics in anvi’o, using real-world data from the human oral microbiome as generated in the paper “Functional and genetic markers of niche partitioning among enigmatic members of the human oral microbiome” by Shaiber et al (2020). Our goal is to show you the following:

• How to run a read recruitment workflow from many samples to a co-assembly, to generate differential coverage data for binning as well as variant data for population genetics
• How to manually refine a ‘metabin’ that was automatically generated by CONCOCT and encompasses several microbial populations into a few smaller yet high-quality bins
• How to set some key parameters in a population genetics analysis that answers a (somewhat) realistic research question

Of course, we are still constrained a little bit by the time and computational resources allowed by a typical tutorial session. We will be working with a smallish co-assembly that is entirely visualizable in the anvi’o interactive interface (unlike many real-world datasets). And it is unlikely that we will actually run the read recruitment workflow in a tutorial setting (but we will give you all the tools to do so in case you want to try it at home).

Without further ado, let’s begin :)

Our dataset and questions

We are taking advantage of the publicly-available data from Shaiber et al 2020. This study described a time-series of tongue and plague metagenomes taken from several human individuals, most of which were paired in male-female couples. This sampling strategy is particularly effective for metagenomic binning for the following reasons. First, having multiple samples from a single individual enabled the authors to do co-assembly, which effectively increases the coverage and likelihood of catching low-abundance populations (at the cost of increasing the overall complexity, which can break assembly algorithms). Second, we know that similar populations likely occur across different individuals and it therefore makes sense to map multiple individuals to a given co-assembly to leverage more differential coverage signals – this was not the strategy used in the Shaiber et al 2020 paper, but it will work quite well for our purposes today.

Here are the questions we want to answer using this dataset:

1. What microbial populations can we recover high-quality metagenome-assembled genomes (MAGs) for by refining ‘metabins’?
2. Do the individuals in a couple share similar populations in their oral microbiomes? To put it another way – can the population variants specific to a given individual help us distinguish between the couples?

In order to answer these questions at a smaller scale, we will be using a subset of the data from the 2020 study - a co-assembly of 5 tongue samples from the time-series of a single individual (the sample named T-B-M in the original paper), plus 20 tongue metagenomes to recruit reads from. These include the 5 samples used to make the co-assembly and 5 each from three other individuals: T-A-F, T-A-M, and T-B-F.

For the purposes of this tutorial, we will use the co-assembly that was already generated by the study authors. We will start our binning journey from the results of the automatic binning tool CONCOCT (Alneberg et al 2014) that was already run on the co-assembly.

In fact, CONCOCT was run twice – once by specifying -c 10 to force the tool to create exactly 10 ‘metabins’ (each containing multiple microbial populations), and once without enforcing a number of bins (in which case the tool would try to group sequences from exactly one population in each bin). The ‘metabin’ strategy is extremely useful for partitioning your large metagenomic datasets into manageable chunks that can be manually refined, and is described elsewhere in a blog post by Tom and Meren.

The mapping workflow

If all we wanted to do was bin refinement, a read recruitment step wouldn’t really be necessary here. After all, the automatic binning results were already generated using existing read recruitment results whereby the 5 time-series samples were mapped against their co-assembly. So we could have just taken the existing, publicly-available contigs and profile databases from this Figshare and done our bin refinement on that.

But we want more than that. To answer question (2), we also want to see how the tongue-associated populations from the individual T-B-M differ from the tongue-associated populations in the three other individuals (which means we need mapping information from 15 additional samples). And if we are going to map additional samples to the co-assembly anyway, might as well use them to help guide our bin refinement process.

This is probably not a step that you want to run on a laptop. We did it on our high-performance computing cluster (HPC), using the snakemake workflow for metagenomics implemented in anvi-run-workflow. We will show you how we set this workflow up, but there is no need for you to run it yourself - you will be able to download the output of the workflow that is required for the rest of the tutorial in the next section.

# keep only accession and sample name
cut -d ',' -f 1,45 Oral_Microbiome_SraRunTable.csv | tr ',' '\t' > oral_samples.txt
# search for sample names of interest and extract to a new file
for n in T_B_M T_B_F T_A_F T_A_M; do \
  grep $n oral_samples.txt >> oral_samps_to_download.txt ; \
done
# extract only the accession for the samples of interest
cut -f 1 oral_samps_to_download.txt > SRA_accession_list.txt

anvi-run-workflow -w sra_download --get-default-config download.json

anvi-run-workflow -w sra_download -c download.json

So, how did we run the mapping workflow?

First, we got our reference metagenome co-assembly for individual T-B-M by downloading the corresponding contigs and profile database from this Figshare link. Then, we extracted the FASTA file of the contigs, as well as both existing CONCOCT collections that were already stored in the profile database:

# extract le data
tar -xvf T-B-M.tar.gz

# update the dbs to match your current anvi'o version
anvi-migrate --migrate-safe T-B-M/*.db

# get our reference FASTA file
anvi-export-contigs -c T-B-M/CONTIGS.db -o T_B_M-contigs.fa

# export the metabin collection
anvi-export-collection -p T-B-M/PROFILE.db -C CONCOCT_c10 -O CONCOCT_c10

# export the standard CONCOCT collection
anvi-export-collection -p T-B-M/PROFILE.db -C CONCOCT -O CONCOCT

We will use the FASTA file of the co-assembly contigs as the reference in our mapping workflow. The mapping workflow will create a new contigs database out of this FASTA file, and associate the mapping data to that new contigs database via a new (merged) profile database containing the mapping info from all 20 of our samples. Because we’ll get a fresh profile database that won’t have any of the existing collection information in it, we will need to import the collection data (using the collection files we just exported) so that we can do our bin refinement from there. Luckily, because we exported our reference FASTA directly from the contigs database, the split names in the exported collection files will match to the names in the new contigs database. But we are getting ahead of ourselves. First we actually have to run the workflow.

On our HPC, we navigated to the folder where we downloaded the 20 metagenome samples that we want to map against this co-assembly. We created a samples-txt file that describes the paths to each of the samples (you can see the format here). Then we generated a fasta-txt file with the path to our reference co-assembly FASTA, and we got a default config file for the metagenomics workflow:

echo -e "name\tpath\nT_B_M\tT_B_M-contigs.fa" > fasta.txt
anvi-run-workflow -w metagenomics --get-default-config map.json

After modifying the config file a bit to turn on ‘reference mode’, turn off unnecessary rules (like functional annotation),, and increase the number of threads for each rule as much as our cluster could handle, we did a dry run (using the -n dry run flag and -q quiet flag passed directly to snakemake using -A or --additional-params) to see what the workflow would actually run:

anvi-run-workflow -w metagenomics -c map.json -A -n -q

And it showed us the following list of jobs:

job                                  count
---------------------------------  -------
annotate_contigs_database                1
anvi_gen_contigs_database                1
anvi_init_bam                           20
anvi_merge                               1
anvi_profile                            20
anvi_run_hmms                            1
anvi_run_scg_taxonomy                    1
bowtie                                  20
bowtie_build                             1
gen_qc_report                            1
gzip_fastqs                             40
import_percent_of_reads_mapped          20
iu_filter_quality_minoche               20
iu_gen_configs                           1
metagenomics_workflow_target_rule        1
samtools_view                           20
total                                  169

Everything looked okay in the dry run. We have 20 samples, so it makes sense that we are doing all the sample-specific jobs 20 times (the only exception, gzip_fastqs, works individually on R1 and R2 files, so it runs twice per samples).

We very confidently started the workflow using the command:

anvi-run-workflow -w metagenomics -c map.json

And after it finished successfully, we took the resulting contigs database and merged profile database, and we imported the original collections of CONCOCT bins that we extracted previously:

anvi-import-collection -p PROFILE.db -C CONCOCT_c10 CONCOCT_c10.txt -c T_B_M-contigs.db
anvi-import-collection -c T_B_M-contigs.db -p PROFILE.db -C CONCOCT CONCOCT.txt

These are the databases that were put into the datapack for this tutorial, which you can download in the next section.

Refining CONCOCT metabins

First, let’s download the tutorial datapack. Just so you know, the archived databack is about half a gigabyte in size, and once unpacked it will take up ~1.6Gb of space on your computer. If you are okay with that, here are the download commands:

wget https://figshare.com/ndownloader/files/52819088 -O BINNING_POPGEN_TUTORIAL.tar.gz
tar -xvf BINNING_POPGEN_TUTORIAL.tar.gz && cd BINNING_POPGEN_TUTORIAL/

You should now be inside the datapack directory on your terminal. The directory should contain 3 databases – the contigs database containing the co-assembly of the 5 time-series tongue metagenomes from one individual (T-B-M), the profile database containing the mapping data from 20 tongue metagenome samples (from the same T-B-M individual as well as three others), and an auxiliary database with lots of extra info related to the mapping (which we will use later).

An initial look at the co-assembly

Let’s take a look at what we have in the database.

anvi-display-contigs-stats T_B_M-contigs.db

A little interactive window should pop up in your browser, filled with useful tables of statistics about the co-assembly. At the top you will see a bar chart of the number of annotations to each single-copy core gene (SCGs; which can be annotated by running anvi-run-hmms) in the assembly contigs. Here is the chart for the Bacteria_71 set of bacterial SCGs:

A bar chart of hits to bacterial SCGs in the T-B-M co-assembly

Since we expect that each microbial population should have one copy of each of these SCGs, it looks like we have at least 20 bacterial populations in this co-assmbly (and at most 35). In fact, the most common number of hits to a given bacterial SCG is 25, which is exactly the number of bacterial genomes that anvi’o predicted this assembly to have:

Predicted number of populations in the T-B-M co-assembly, based on SCG counts

If you want to learn more about how anvi’o estimates the number of genomes, you can click the little [?] link in the interactive interface, or just click on this conveniently-placed link to the same documentation.

Anyway, we expect to find at least 25 bacterial populations in this co-assembly (there will probably be even more than that because many of the populations are likely to be incomplete). We’re going to do our best to bin at least some of these populations.

We will be refining metabins that were automatically binned by CONCOCT. Let’s check which collections we have available for this co-assembly. Collections of bins are stored in the profile database, so we pass that database as a parameter to anvi-show-collections-and-bins:

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

You should see two collections. One of them is just called CONCOCT and it contains over 50 bins:

Collection: "CONCOCT"
===============================================
Collection ID ................................: CONCOCT
Number of bins ...............................: 58
Number of splits described ...................: 12,499
Number of contigs described ..................: 12,173

These are the binning results from a standard run of CONCOCT, whereby the program did its best to split each individual microbial population into its own bin. As you can see, the number of bins is much larger than the expected number of microbial populations based on single-copy core genes. It is always possible (and perhaps even likely) that we are missing some SCG annotations due to incomplete coverage of less abundant populations. However, it is also possible (and perhaps even very likely) that CONCOCT fragmented some of these population genomes across multiple bins. Which is one reason that we like to use the metabin approach described in this blog post - if you tell CONCOCT to make just a few bins (much less than the number of expected populations), then it is less likely to split sequences from the same population across different bins.

Happily, the second collection in the profile database, CONCOCT_c10, is the metabin collection. Shaiber et al. created this collection by running CONCOCT with the -c 10 parameter, to request 10 bins:

Collection: "CONCOCT_c10"
===============================================
Collection ID ................................: CONCOCT_c10
Number of bins ...............................: 10
Number of splits described ...................: 12,499
Number of contigs described ..................: 12,173
Bin names ....................................: Bin_1, Bin_10, Bin_2, Bin_3, Bin_4, Bin_5, Bin_6, Bin_7, Bin_8, Bin_9

This is the collection we are going to work on in just a moment.

First, let’s just take a look at the co-assembly and its mapping data.

anvi-interactive -c T_B_M-contigs.db -p PROFILE.db \
    --title "T-B-M co-assembly"

The T-B-M co-assembly in the interactive interface.

The profile database comes with pre-loaded settings to make the visualization a bit more digestible. The sample layers are colored according to which individual they come from - in particular, the 5 T-B-M samples that were used to make the co-assembly are in dark green. We can see coverage patterns in the samples from the other individuals that were mapped to this co-assembly – the other individuals’ tongue microbiomes include similar populations to those in T-B-M’s – which is good news for us because it means we can use the mapping data from the other samples to help guide our bin refinement.

If you check the ‘Items order’ label at the top, you will see that the contigs (technically, their splits) are organized according to their shared sequence composition (tetranucleotide frequency) and differential coverage patterns across all 20 samples. The inner dendrogram displays this organization, and we often use this dendrogram to select which contigs to put into a bin – because we expect that sequences originating from the same genome should have roughly the same composition and follow similar differential coverage patterns.

Just out of curiousity, let’s take a look at how the standard CONCOCT bins group the contigs from the co-assembly. You can either load the bin collection CONCOCT from the ‘Bins’ panel in the interface, or close and re-open the interface with --collection-autoload in the command:

anvi-interactive -c T_B_M-contigs.db -p PROFILE.db \
    --title "T-B-M co-assembly (CONCOCT collection)" \
    --collection-autoload CONCOCT

The regular CONCOCT bins

You can see a lot of mixing between the different bins. There are some sections of the figure where a clade on the inner dendrogram has been sorted into a single bin (such as Bin 2, Bin 10, and Bin 11), but for the most part, the bins are quite fragmented. This is exactly what we wanted to avoid by using the metabin technique.

What do the metabins look like? Again, you can load the collection from the interface, or re-start the interface with a new collection name to load:

anvi-interactive -c T_B_M-contigs.db -p PROFILE.db \
    --title "T-B-M co-assembly (CONCOCT_c10 collection)" \
    --collection-autoload CONCOCT_c10

The 10 CONCOCT metabins

It is a little bit better. Now multiple bacterial populations have been purposefully combined into just a few metabins, and we can see a bit less mixing; ie, the binning follows the dendrogram structure better. And in general, the mixing is not as chaotic – many clades of the inner dendrogram are sorted into just two bins. We should be able to work with this output to refine a metabin into some high-quality MAGs.

Why so much mixing? Is CONCOCT doing a really bad job? Maybe. After all, automatic binning of complex metagenomes is an extremely difficult task and algorithmic heuristics won’t always give us an answer we expect. But maybe not. It is important to remember that CONCOCT was given only a fraction of the information shown here to base its binning decisions on. It only used differential coverage signal from the 5 dark green samples from T-B-M – the samples used to make the co-assembly. Here, we are showing mapping information from 15 extra samples from 3 different individuals, and while it seems clear that those individuals harbor similar populations, there is likely still some variation like large-scale insertions and deletions or hyper-variable regions that could confuse the differential coverage signals. So give CONCOCT a break – it is trying its best to live up to our unrealistic expectations. :P

Refining Bin 3

We want to refine a metabin and split it into its component populations. But which metabin should we work on? If you take a look in the ‘Bins’ panel on the interactive interface, you should see the overall completeness and redundancy estimates for each one:

Completion and redundancy of the metabins

As expected, they have a lot of redundancy because each bin contains several microbial populations. The only potential exception is Bin 2, which clearly stands out as one of the only bins that cohesively groups several consecutive clades in the inner dendrogram and seems to represent a single population based on the coverage patterns (though the size of the bin is a bit large for most microbial genomes). So we would like to focus on refining one of the other bins.

Somewhat arbitrarily, we selected the second-largest bin, Bin 3, to refine for this tutorial. Based on its redundancy score, we should be able to obtain at least 4-5 bins out of this metabin (potentially more if the populations are incomplete).

Here is the command to refine Bin 3. It will open the interactive interface again, but this time only showing the contigs (technically, their splits) belonging to Bin 3.

anvi-refine -c T_B_M-contigs.db -p PROFILE.db -C CONCOCT_c10 -b Bin_3

Here is what Bin 3 looks like in the interface:

Refining metabin ‘Bin 3’ in the interface’s ‘refine mode’.

Now we can use the dendrogram organization, the differential coverage patterns, and the SCG-based completion/redundancy/taxonomy estimates to identify the subsets of this metabin that represent individual population genomes. Give it a try for yourself, and once you are ready, click the Show/Hide box to see if your results match ours.

Make sure to hit the button called ‘Store refined bins in database’ in the ‘Bins’ panel to save your work once you are done refining. Please note that this will overwrite your existing larger collection (in our case, this is the CONCOCT_c10 collection), replacing the original metabin (Bin_3) with the set of smaller bins you curated.

Congratulations! You just refined a metabin. :)

If you want to see how your bin looks in the context of the larger co-assembly, you can do that by going back to the regular interactive interface, using the following command to automatically open the (now partially-refined) CONCOCT_c10 collection.

anvi-interactive -c T_B_M-contigs.db -p PROFILE.db \
    --title "T-B-M co-assembly (CONCOCT_c10 collection)" \
    --collection-autoload CONCOCT_c10

anvi-import-collection -p PROFILE.db -C refined_c10 refined.txt -c T_B_M-contigs.db
anvi-import-collection -p PROFILE.db -C refined_Bin_3 refined_Bin_3.txt -c T_B_M-contigs.db

anvi-interactive -c T_B_M-contigs.db -p PROFILE.db \
    --title "T-B-M co-assembly (CONCOCT_c10 collection)" \
    --collection-autoload refined_c10

anvi-interactive -c T_B_M-contigs.db -p PROFILE.db \
    --title "T-B-M co-assembly (refined bins collection)" \
    --collection-autoload refined_Bin_3

As you can see, visualizing binning results in different contexts and different organizations can be, paradoxically, both confusing and helpful. We could spend ages trying to improve our binning results, and often only see a marginal improvement. Bin refinement is truly an art rather than an exact science. For the sake of this tutorial, we are done here.

Population genetics

Now that we have some manually-refined bins, we can answer our second question: Do the individuals in a couple share similar populations in their oral microbiomes? Now it becomes relevant for you to know that the four individuals represented in our datapack were not chosen randomly – they actually represent two couples. T-B-M is the partner of T-B-F (their sample layers were dark and light green in our interactive interface figures above), and T-A-F is the partner of T-A-M (their samples were the light and dark pink layers). So the question that we want to answer is: are the microbes from T-B-M more similar to the microbes from T-B-F than to either of the communities from T-A-F or T-A-M? We will answer this question using population genetics, whereby we will compare the single-nucleotide variant patterns of a microbial population across different samples to see if the variation is similar between the couples’ samples.

In theory, we could run our comparison using all of the sequence variants in the entire co-assembly, ignoring the boundaries of population genomes and utilizing all the data we have available. But in our case today, we are interested in visualizing the output with the anvi’o interactive interface, and unfortunately that places a limit on the number of variants we can handle (to around 10,000 variant positions). So to reduce the number of variants to a tractable amount, we’ll be focusing on comparing one microbial population across the different samples.

Choosing a population (and samples) to focus on

Which microbial population should we compare? In order to robustly call sequence variants, anvi-profile requires a nucleotide position to have a minimum coverage of 10x (by default). If coverage is less than that minimum threshold at a given position, it won’t identify any single-nucleotide variants (SNVs) at that position. This means that we need a population with enough mean coverage across our samples to be sure that we will have enough variant information for our comparison.

Let’s take a look at the coverage of our refined bins across our samples. If you haven’t already, you can import our refined collection into the contigs database of the datapack:

anvi-import-collection -p PROFILE.db -C refined_Bin_3 refined_Bin_3.txt -c T_B_M-contigs.db

These are the 7 bins we refined from the metabin ‘Bin 3’ in the previous chapter. To get a table of coverage for our bins across samples, we will use the program anvi-summarize:

anvi-summarize -p PROFILE.db -c T_B_M-contigs.db -C refined_Bin_3

Once this command finishes, you should see a new directory of output containing a bunch of summary tables describing our collection of refined bins. Within that directory is a subfolder called bins_across_samples:

ls SUMMARY/bins_across_samples/

In the subfolder’s content list, you should see files with coverage in the name. The mean coverage table (mean_coverage.txt) describes the average coverage across all sequences in a bin within each sample, and the Q2Q3 version (mean_coverage_Q2Q3.txt) does the same thing while ignoring positions with outlier coverage. std_coverage.txt is a table of the standard deviation of the coverage within a bin. We’ll take a look at the mean_coverage_Q2Q3.txt file to get a sense of which microbial populations are consistently well-covered across most samples.

Currently that table has bins in the rows and samples in the columns, which makes it quite wide and unwieldy to look at in the terminal. Let’s quickly transpose that matrix and print it out.

anvi-script-transpose-matrix SUMMARY/bins_across_samples/mean_coverage_Q2Q3.txt -o Q2Q3_transposed.txt
cat Q2Q3_transposed.txt

You should see a table like the following:

Remembering that we want to have a lot of nucleotide positions with at least 10x coverage in order to have decent SNV profiles, what we are looking for is a microbial population that has fairly high overall coverage in as many samples as possible. From this table, it seems clear that none of the microbial populations we binned are going to fit perfectly, since the overall coverages are rather low. But the microbial population with the overall highest coverage is Prevotella jejuni, so that bin will be the target of our population genetics analysis.

import pandas as pd
df = pd.read_csv("Q2Q3_transposed.txt", sep="\t", index_col=0)
df.mean()
df.median()

You might remember from our binning results that the Prevotella jejuni bin was rather small and incomplete (it didn’t even have any SCGs annotated in it). So in reality, we’ll only be using a subset of the P. jejuni population genome. It is not ideal, but for our particular question it’s more important that we have high enough coverage for proper SNV calling than to have a complete genome.

To that end, we should probably not work with any sample in which the Prevotella jejuni bin has low coverage. We can parse the coverage table to keep only samples in which P. jejuni has a non-outlier mean coverage of at least, let’s say 15x. A mean coverage of 15x doesn’t mean that all nucleotide positions in the bin will have >10x coverage (our minimum threshold for SNV calling), but at least it is likely that a majority of positions will cross that threshold.

To get a list of samples in which this bin has at least 15x mean coverage, we can go back to our coverage table (the table directly in the summary output, SUMMARY/bins_across_samples/mean_coverage_Q2Q3.txt). Below, I export a list automatically using python, but you can feel free to use another way to do it:

import pandas as pd
df = pd.read_csv("SUMMARY/bins_across_samples/mean_coverage_Q2Q3.txt", sep="\t", index_col=0)
df[df > 15].count(axis=1)
df[df > 15].loc['Prevotella_jejuni',].dropna()
samps = df[df > 15].loc['Prevotella_jejuni',].dropna().index.to_list()
with open('samples_of_interest.txt', 'w') as f:
     for s in samps:
             f.write(f"{s}\n")

(A small reminder that you can leave the Python interpreter in the command line by running exit()).

The output we should get is a text file, called samples_of_interest.txt, that contains the 11 samples in which the P. jejuni bin has >15x mean (Q2Q3) coverage. This is what it should look like:

T_A_F_01
T_A_F_05
T_A_M_01
T_A_M_02
T_A_M_03
T_A_M_04
T_A_M_05
T_B_M_01
T_B_M_02
T_B_M_03
T_B_M_05

As you can see, we only have 3 individuals left in our analysis. T-B-F dropped out because P. jejuni didn’t have enough coverage in any of their samples. However, we still have one couple (T-A-F and T-A-M), so we can still take a stab at seeing if these two individuals harbor similar P. jejuni populations.

Let’s move on to the population genetics comparison.

The population genetics part

Now that we know which bin to compare across which samples, we can ask anvi’o to profile the variation for us. Here is the command:

anvi-gen-variability-profile -c T_B_M-contigs.db \
               -p PROFILE.db \
               -C refined_Bin_3 \
               -b Prevotella_jejuni \
               --samples-of-interest samples_of_interest.txt \
               --min-occurrence 5 \
               --min-coverage-in-each-sample 20 \
               --include-split-names \
               --quince-mode \
               -o Prevotella-SNVs.txt

The first few inputs should hopefully be self-explanatory – we provide the co-assembly, the read recruitment information (which includes SNVs that were called by the anvi-profile program during the mapping workflow), our collection of refined bins, the name of the bin that we want to compare, and the list of samples that we selected in the previous section.

The --min-occurrence 5 argument requires that a given SNV position has to be variable in at least 5 samples (out of the 11 samples that we specified in samples_of_interest.txt) in order to be reported. This threshold ensures that we get rid of spurious variants only present in a small fraction of samples that won’t be very useful in our comparison anyway. The --min-coverage-in-each-sample 20 is yet another filter, requiring a variant position to have at least 20x coverage in each sample to be reported. The higher coverage requirement means that when a given position is variable in some samples but not in other samples, we can be sure we aren’t accidentally missing variation in the latter samples just due to low coverage. It also means that when a position is variable, we can be sure it is a true variant and not a random sequencing or alignment error (although the 10x min coverage threshold in anvi-profile should already take care of this to some extent). Both of these filters reduce the number of potential variable positions we will look at, and their goal is to reduce noise.

The other flags control the output of the program. --include-split-names asks for split names to be included in the output variability table, so that we could inspect splits with interesting SNV positions if we wanted to. --quince-mode ensures that each position that is variable in at least 1 sample (by default, or in our case, at least 5 samples due to our use of --min-occurrence 5) will have a corresponding entry from every single sample in the output variability table, even if the position is not variable in some of those samples. As the output file will report variable nucleotide positions on a per-sample basis, you can imagine that ‘Quince mode’ basically completes the matrix of variability across samples – in our case, every position in the co-assembly will have 11 corresponding entries in the output file (one entry for the position’s variation in each of our 11 samples of interest). Finally, -o specifies the name of our output file.

If you run the anvi-gen-variability-profile command, you will see how the filters are working from the output on the terminal. We’ve copied some of the relevant output lines here:

Variability data ..............................................: 9,514,685 entries in 12,499 splits across 11 samples
Entries after split name(s) of interest filter ................: 258,587 (9,256,098 were removed) [320 genes remaining]
Entries after minimum occurrence of positions filter ..........: 224,760 (33,827 were removed) [315 genes remaining]
Entries after quince mode filter ..............................: 281,259 (56,499 were added) [315 genes remaining]
Entries after minimum coverage across all samples filter ......: 109,527 (171,732 were removed) [222 genes remaining]
Num entries reported ..........................................: 109,527
Num NT positions reported .....................................: 9,957

In other words, while there were ~9.5 million SNV positions in the entire T-B-M co-assembly, there were only 258,587 in our Prevotella jejuni bin. When we kept only nucleotide positions that were variable (ie, SNVs at that position) in at least 5 samples, most of them (224,760) remained, which means there were a lot of SNVs that were consistently found across many samples. Filling out the matrix to include an entry for each position in each sample of interest (quince mode) obviously adds a bunch of entries, but many SNV positions drop out due to having <20x coverage in at least one sample.

So we ultimately end up with ~109k entries corresponding to about 10,000 variable nucleotide positions. That’s quite a lot of variation for such a little bin. This population must be quite variable across the different samples (which could include variation over time in one individual and/or variation between the P. jejuni populations living in different hosts).

Let’s take a look at the output. Here are the first 15 rows of the Prevotella-SNVs.txt output file:

The columns in this output file have been explained elsewhere, so we won’t repeat ourselves here. You can at least notice that we have 11 entries for the same nucleotide position (position 210 in the table), which is due to the quince mode filter. This position is a T in the co-assembly contig sequence, but most of the samples have at least a portion of reads mapping to the position that have a C nucleotide instead – in fact, C is the consensus nucleotide in samples like T_A_F_01, T_A_M_03, T_A_M_04, and all of the T_B_M samples. The changing consensus base in some samples suggests that the variation at this particular position changed over time in the individuals T-A-F and T-A-M, but as far as we know (from the samples we were able to include), it mostly stayed stable with a majority C base in the individual T-B-M.

Of course, manually examining all variable positions is not a tractable strategy to compare our samples. Ideally we would like to visualize the results to answer our target question.

Visualizing the variable positions

The SNV interactive interface

First, let’s do it using the anvi’o interactive interface. We can do this by converting our variability table into a view-data file, which just means a tab-delimited text file that the interactive interface can read from in --manual mode, and in this particular case will make the variable nucleotide positions into ‘items’ and decorating information about its variation in different samples into the ‘layers’.

Luckily, there is already a program that will consolidate and reformat our variability table for us, then create a profile database with a nice organization of the variable positions:

anvi-script-snvs-to-interactive Prevotella-SNVs.txt -o Prevotella-snvs

Note that this script by default makes the departure_from_consensus column the value shown in the interactive interface for each position. This helps you distinguish between positions that are extremely uniform within a given sample and those that are extremely variable. If you want to visualize the departure_from_reference values instead (to identify at which positions the short reads do not agree with the reference assembly), you can use the parameter --display-dep-from-reference.

We can then visualize the variable positions by running the following command:

anvi-interactive --manual -p Prevotella-snvs/profile.db \
                 --tree Prevotella-snvs/tree.txt \
                 --view-data Prevotella-snvs/view.txt \
                 --title "P. jejuni variable positions"

There is a technical limit to the number of variable positions that can be shown in the anvi’o interactive interface in this way (similar to how there is a limit for the number of splits when visualizing an assembly). We are below the limit in our case, but if you end up passing that limit when using your own data, you’ll either need to filter out more positions to get below the limit, find an alternative way to visualize the data, or analyze the variability table without visualizing it.

The initial SNV interactive interface.

The main layers of this interface show the departure_from_consensus value for each position in each sample that we analyzed. We also get layers to show the competing nucleotides, the corresponding gene caller ID (if the position is within a gene call), and its relative position within a codon (again, if it is within a protein-coding gene) for each variable position.

But currently the interface isn’t really pretty. Let’s make it match the colors from our co-assembly so that we can more easily distinguish between samples from different individuals (and different couples).

To do this, we’ll export the visualization settings from the co-assembly’s profile database, and import those settings into this profile database:

anvi-export-state -p PROFILE.db -s default -o default.json
anvi-import-state -p Prevotella-snvs/profile.db -n default -s default.json

Then we can simply re-run the interactive command:

anvi-interactive --manual -p Prevotella-snvs/profile.db \
                 --tree Prevotella-snvs/tree.txt \
                 --view-data Prevotella-snvs/view.txt \
                 --title "P. jejuni variable positions"

A much prettier version of the SNV interactive interface.

Much better. Now it is extra clear that the different individuals have distinct versions of the P. jejuni population. In particular, there is not as much variation in the dark green T-B-M samples, which makes sense because this is the population that we have binned from the T-B-M co-assembly. The positions that do vary in this individual generally seem to agree with each other (low departure from consensus) at time points 1 and 2, meaning that there isn’t much subpopulation variation in this individual at that time. But this seems to have shifted after time point 3, at which point extensive subpopulation variation of P. jejuni appeared.

We can also see that there is high departure from consensus in the samples from the two other (partnered) individuals, and more importantly, that these variable positions clearly match across those two individuals with a similar amount of departure from consensus. So we know that they are different from the consensus at these positions, but we don’t yet know if their specific nucleotides match between T-A-F samples and T-A-M samples. So we don’t quite have have an answer to our question just yet.

Fixation index to compare samples

To confirm whether the P. jejuni populations are similar between T-A-F and T-A-M, we can compute a matrix of fixation index instead, and automatically cluster this matrix to group similar samples together.

Let’s do it. All we need is to pass the variability table to anvi-gen-fixation-index-matrix, cluster the output matrix to get a dendrogram in Newick format, and open the dendrogram in the interface:

anvi-gen-fixation-index-matrix --variability-profile Prevotella-SNVs.txt \
                               --output-file FST_Prevotella.txt
anvi-matrix-to-newick FST_Prevotella.txt -o FST_Prevotella.nwk
anvi-interactive -t FST_Prevotella.nwk \
                 -p FST_Prevotella.db \
                 --manual

After a bit of fussing with the interface settings to convert the sample names into text (instead of the default colors), and binning the samples from individual T-B-M, you could see a tree that looks like the following:

The fixation index tree.

So, now we have our answer :) The samples from couple ‘A’ do truly share similar P. jejuni populations, which are distinct from the P. jejuni living in our individual from couple ‘B’.

A network visualization

Another way of displaying the same data would be to compute a network, where variable positions are small nodes, and each sample is a big node, and the big nodes are connected via edges to small nodes if the corresponding position is variable in the corresponding sample. Visualizing such a network is another strategy to answer the question of ‘which samples come together based on having similar P. jejuni populations’, like we did with the fixation index tree, but this time you’d also be able to see the SNVs that bring the samples together. And – you guessed it – we have a helper program in anvi’o to get you the proper network file.

Here is the command to convert the variability table into a graph file that is compatible with the network visualization software Gephi:

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

If you open that file in the Gephi software, you can choose the ‘Force Atlas 2’ option in the ‘Layout’ panel on the left side, and hit the ‘Run’ button so that it computes an optimal arrangement of the nodes and edges in the network.

The Force Atlast 2 layout option for visualizing the network.

We also recommend turning on the options ‘LinLog mode’ and ‘Prevent Overlap’ to make the network even easier to read. Once the layout algorithm has run for a bit and hopefully reached some sort of consensus arrangement, you can hit ‘Stop’ in the ‘Layout’ panel, and switch to the ‘Preview’ tab to view the final visualization. You can turn on the ‘Show Labels’ option and turn off the ‘Show Edges’ option in the settings on the left side to make it look less like a crazy spiderweb and more like a somewhat organized collection of named dots.

The network visualization should hopefully look similar to this one:

Our Gephi visualization of the P. jejuni SNV network.

As you can see, once again the samples from the ‘A’ couple loosely group together, and are clearly distinct from the samples of our T-B-M individual.

Conclusion

In our small yet realistically-complex tongue microbiome datapack, we’ve been able to achieve our two goals. We refined 7 bins from one CONCOCT metabin, including population genomes of expected oral microbes like TM7 and Prevotella. And, we used one of those refined bins to prove that humans in long-term relationships with each other share similar oral microbiota – the reasons for this are likely obvious to you, but now you can see the proof with your own eyes. ;)

We hope this little tutorial helped to demonstrate the reality of metabin refinement and the utility of population genetics to answer a realistic biological question. If you have any questions or suggestions, let us know in the comments.