Phylogenetics

Overview

After pangenome analysis, we can use the core gene alignment to construct a phylogenetic tree.

Phylogenetic analysis helps us understand the evolutionary relationship between bacterial isolates. Closely related samples usually cluster together in the tree, while more distant samples appear farther apart.

In this lesson, we will use:

  • snp-sites to extract variable sites from the core gene alignment
  • snp-sites -C to count constant sites
  • IQ-TREE to construct a phylogenetic tree

Learning objectives

By the end of this lesson, you should be able to:

  • activate the phylogeny environment
  • identify the core gene alignment file from Panaroo
  • extract variable SNP sites from the alignment
  • count constant sites
  • run IQ-TREE to build a phylogenetic tree
  • identify important tree output files
  • download the tree file for visualization

Input file

In the previous lesson, Panaroo produced a filtered core gene alignment:

pangenome/core_gene_alignment_filtered.aln

Login to server

ssh genomevm@172.28.28.12
cd tanzim

Check that the file exists:

ls pangenome

You should see:

core_gene_alignment_filtered.aln
NoteImportant

For this lesson, we will use:

pangenome/core_gene_alignment_filtered.aln

as the input alignment.

Activate the phylogeny environment

Activate the phylogeny environment:

conda activate phylogeny

Check that the required tools are available:

snp-sites --help
iqtree --help

Create tree output directory

Create a directory named tree:

mkdir -p tree

Check that the directory was created:

ls
grep ">" pangenome/core_gene_alignment_filtered.aln

Extract variable sites

The core gene alignment contains both variable and constant sites.

Variable sites are positions that differ among samples. These positions are useful for building a phylogenetic tree. Although you could use the alignment generated by Panaroo directly as input to IQ-TREE, this would be quite computationally heavy, because the core genome alignments tend to be quite big. Instead, what we can do is extract the variable sites from the alignment, such that we reduce our FASTA file to only include those positions that are variable across samples.

Here is a small example illustrating what we are doing. For example, take the following three sequences, where we see 3 variable sites (indicated with an arrow):

S1  C G T A G C T G G T
S2  C T T A G C A G G T
S3  C T T A G C A G A T
      ↑         ↑   ↑

For the purposes of phylogenetic tree construction, we only use the variable sites to look at the relationship between our sequences, so we can simplify our alignment by extract only the variable sites:

S1  G T G
S2  T A G
S3  T A A

This example is very small, but when you have a ~4Mb genome, this can make a big difference. To extract variable sites from an alignment we can use the SNP-sites software:

Use snp-sites to extract variable sites:

snp-sites \
  pangenome/core_gene_alignment_filtered.aln \
  -o pangenome/core_gene_alignment.snp.aln

Check SNP alignment

Check that the SNP alignment was created:

ls pangenome/*snp*

You should see:

pangenome/core_gene_alignment.snp.aln

Check the first few lines:

head -n 2 pangenome/core_gene_alignment.snp.aln

We can see this contains a sequence named “S3”, which corresponds to one of the samples we used in day 2. We can look at all the sequence names in the FASTA file:

We can see each input genome appears once

However, before we move on to that step, we need another piece of information: the number of constant sites in the initial alignment (sites that didn’t change). Phylogenetically, it makes a difference if we have 3 mutations in 10 sites (30% variable sites, as in our small example above) or 3 mutations in 1000 sites (0.3% mutations). The IQ-TREE software we will use for tree inference can accept as input 4 numbers, counting the number of A, C, G and T that were constant in the alignment. For our small example these would be, respectively: 1, 2, 2, 2.

Fortunately, the snp-sites command can also produce these numbers for us (you can check this in the help page by running snp-sites -h). This is how you would do this:

Count constant sites

When building a tree from a SNP-only alignment, we need to tell IQ-TREE how many constant sites were removed.

Use snp-sites -C to count constant sites:

snp-sites \
  -C \
  pangenome/core_gene_alignment_filtered.aln \
  -o pangenome/core_gene_alignment_countconstant

View the constant site counts:

cat pangenome/core_gene_alignment_countconstant

You should see something like:

716495,650538,717930,710604

As we said earlier, these numbers represent the number of A, C, G, T that were constant in our original alignment. We will use these numbers in the tree inference step detailed next So, copy the value printed to somewwhere safe. This value will be used in the -fconst option in IQ-TREE.

Tree inference: IQ-TREE

There are different methods for inferring phylogenetic trees from sequence alignments. Regardless of the method used, the objective is to construct a tree that represents the evolutionary relationships between different species or genetic sequences. Here, we will use the IQ-TREE software, which implements maximum likelihood methods of tree inference. Phylogenetic tree inference using maximum likelihood is done by identifying the tree that maximizes the likelihood of observing the given DNA sequences under a chosen evolutionary model.

In this process, DNA substitution models describe how DNA sequences change over time due to mutations. These models consider how frequently different bases (A, T, C, G) are replaced by each other. Another parameter these models can include is rate heterogeneity, which accounts for the fact that different DNA sites may evolve at different rates. Some sites might change rapidly, while others remain more stable.

Maximum likelihood aims to find the tree topology and branch lengths that make the observed DNA sequences most probable, given the chosen model. It does this by exploring various tree shapes and lengths to calculate the likelihood of the observed sequences. The tree with the highest likelihood is considered the best representation of the evolutionary relationships among the sequences. The process involves making educated guesses about the tree’s parameters, calculating the likelihood of the data under these guesses, and refining the parameters iteratively to find the optimal tree that best explains the observed genetic variations.

IQ-TREE offers various sequence evolution models, allowing researchers to match their analyses to different types of data and research questions. Conveniently, this software can identify the most fitting substituion model for a dataset (using a tool called ModelFinder), while considering the complexity of each model.

Construct the phylogenetic tree

We run IQ-TREE on the output from SNP-sites, i.e. using the variable sites extracted from the core genome alignment and replace FCONST_VALUE with the value from:

pangenome/core_gene_alignment_countconstant

Run iq-tree:

iqtree \
  -s pangenome/core_gene_alignment.snp.aln \
  -fconst FCONST_VALUE \
  -seed 12345 \
  -B 1000 \
  -m GTR+F \
  -nt 4 \
  -pre tree/iqtree
iqtree \
  -s pangenome/core_gene_alignment.snp.aln \
  -fconst 716495,650538,717930,710604 \
  -seed 12345 \
  -B 1000 \
  -m GTR+F \
  -nt 4 \
  -pre tree/iqtree

Understanding the IQ-TREE command

Option Meaning
iqtree Runs IQ-TREE
-s pangenome/core_gene_alignment.snp.aln Input SNP alignment
-fconst FCONST_VALUE Constant site counts
-seed 12345 Random seed for reproducibility
-B 1000 Ultrafast bootstrap replicates
-m GTR+F Substitution model
-nt 4 Use 4 CPU threads
-pre tree/iqtree Output prefix
NoteAbout FCONST_VALUE

FCONST_VALUE is not typed literally.

Replace it with the constant site counts produced by:

cat pangenome/core_gene_alignment_countconstant
WarningModel note

The model with +ASC does not always work in this setup, so this lesson uses:

GTR+F

Check IQ-TREE outputs

List the tree directory:

ls tree

You should see files such as:

iqtree.bionj   iqtree.contree  iqtree.log     iqtree.splits.nex
iqtree.ckp.gz  iqtree.iqtree   iqtree.mldist  iqtree.treefile

Important IQ-TREE output files

File Description
iqtree.treefile Main phylogenetic tree file
iqtree.contree Consensus tree file
iqtree.iqtree Detailed IQ-TREE report
iqtree.log Log file

The main output file for tree visualization is:

tree/iqtree.treefile

View the tree file

The tree file is stored in Newick format.

View it using:

cat tree/iqtree.treefile

A Newick tree looks like this:

((S1:0.01,S2:0.01):0.02,(S3:0.03,S4:0.03):0.02);
NoteExercise

Run IQ-TREE and identify the main tree output file.

Which file will you use for tree visualization?

The main tree output file is:

tree/iqtree.treefile

Download tree file from the server

If you are working on a remote server, download the tree directory to your local computer.

From your local desktop terminal, run:

scp -r genomevm@172.28.28.12:~/tanzim/tree .
scp -r genomevm@172.28.28.12:~/tanzim/`tree` .

Open iqtree.treefile using a text editor or upload it to an online tree viewer such as iTOL.

Use:

open iqtree.treefile

Use:

xdg-open iqtree.treefile

Check tree support values

If bootstrap analysis completed successfully, support values may be included in the tree.

You can inspect the IQ-TREE report:

cat tree/iqtree.iqtree

The report includes information about:

  • input alignment
  • substitution model
  • likelihood values
  • bootstrap settings
  • tree search summary

Compare tree with metadata

After constructing the tree, compare the clustering pattern with sample metadata.

Possible metadata variables include:

  • sample ID
  • location
  • collection date
  • source
  • sequence type
  • AMR profile
  • virulence factor profile
NoteQuestion

Do samples with similar MLST or AMR profiles cluster together in the tree?

Closely related samples may cluster together and may share similar MLST, AMR, or virulence profiles. However, clustering should be interpreted carefully together with metadata and quality-control results.

Final directory structure

At the end of this lesson, your directory should look like this:

.
├── pangenome
│   ├── core_gene_alignment_filtered.aln
│   ├── core_gene_alignment.snp.aln
│   └── core_gene_alignment_countconstant
└── tree
    ├── iqtree.treefile
    ├── iqtree.contree
    ├── iqtree.iqtree
    └── iqtree.log

Practical exercise

Complete the following tasks:

  1. Activate the phylogeny environment.
  2. Check that core_gene_alignment_filtered.aln exists.
  3. Extract variable sites using snp-sites.
  4. Count constant sites using snp-sites -C.
  5. Replace FCONST_VALUE in the IQ-TREE command.
  6. Run IQ-TREE.
  7. Identify the main tree file.
  8. Download the tree file to your local computer.
CautionExercise: phylogeny summary

Create a small table in your notes:

Item File or value
Input alignment pangenome/core_gene_alignment_filtered.aln
SNP alignment
Constant site value
IQ-TREE output folder
Main tree file

Fill in the table using the files created during this lesson.

Key points

NoteImportant
  • Phylogenetic analysis shows relationships among bacterial isolates.
  • The input is the filtered core gene alignment from Panaroo.
  • snp-sites extracts variable sites from the alignment.
  • Constant site counts are needed when building a tree from SNP-only alignments.
  • IQ-TREE constructs the phylogenetic tree.
  • The main tree file is tree/iqtree.treefile.

Next step

In the next lesson, we will visualize and interpret the phylogenetic tree.

Back to top