Build a core-genome phylogeny

Overview

This tutorial describes step by step the reconstruction of a core genome phylogeny for a given clade of bacteria. The main steps are:

1. getting the list of available genomes from the NCBI website
2. download genomes from the NCBI website
3. identify orthologous proteins between the different genomes with orthofinder
4. Orthofinder will reconstruct the phylogeny of the species based on identified orthologous proteins

Setup

1. clone eutils repository

   1 # download repo from metagenlab
   2 git clone https://github.com/metagenlab/eutils.git
   3 # move to eutils directory
   4 cd eutils
   5 # add path to pythonpath and path
   6 echo "export PATH=:$PWD"':$PATH' >> ~/.bashrc
   7 echo "export PYTHONPATH=:$PWD"':$PYTHONPATH' >> ~/.bashrc
   8 # reload bashrc
   9 source ~/.bashrc
   

2. install miniconda
3. create a dedicated environment and install orthofinder

   1# create a new environment
   2conda create -n orthofinder
   3source activate orthofinder
   4# install orthofinder
   5conda install orthofinder
   6# install biopython (used to query the NCBI website and manipulate genbank/fasta files)
   7conda install biopython
   

Get table of genome assemblies

The list of available genome assemblies can be downloaded with the script get_ncbi_genome_assemblies.py. Various options are available to filter the data in various ways:

• include only complete genomes
• include only representative or reference genomes (as classified on the NCBI)
• exclude incompletely classified genomes
• exclude genomes derived from metagenomes

If your taxon of interest includes too many genomes, it is also possible to get only one representative genome per species (or per genus/family/order). If you do so, one genome is selected among all assemblies available for each species. The genome is not selected randomly: priority is given to 1) complete genomes 2) reference genomes 3) representative genomes.

usage: get_ncbi_genome_assemblies.py [-h] [-t TAXON_ID] [-c] [-rf] [-rp] [-rs]
                                     [-gb] [-md] [-e] [-d] [-ra RANK]
                                     [-o OUTPUT]

optional arguments:
  -h, --help            show this help message and exit
  -t TAXON_ID, --taxon_id TAXON_ID
                        taxon_idd
  -c, --complete        Filter complete genomes (default: False)
  -rf, --reference      Filter reference genomes (default: False)
  -rp, --representative
                        Filter representative genomes (default: False)
  -rs, --refseq         Filter refseq assemblies (default: False)
  -gb, --genbank        Filter genbank assemblies (default: False)
  -md, --exclude_metagenome_derived
                        Exclude metagenome derived assemblies (default: False)
  -e, --exclude_unclassified
                        Exclude species labelled as "unclassified xxx" or as
                        "xxxx sp. yyy" (e.g. Blautia sp. TF11-31AT)
                        (default=False)
  -d, --complete_list   Report complete list of assemblies (no taxonomy
                        clustering) (default: False)
  -ra RANK, --rank RANK
                        Keep one prepresentative/clade of specified taxonomic
                        <rank>
  -o OUTPUT, --output OUTPUT
                        Output name (default: assemblies_<taxid>.tab )

Example: get one assembly per species for the Lachnospiraceae family

• go on the NCBI taxonomy website. Search for Lachnospiraceae, then click to go to the main page of the family to find the unique taxonomic identifier of the family. It will look like this (the taxonomic ID is indicated with the red arrow):

• execute the script using the family taxon ID. Since the Lachnospiraceae family encompass > 100 different species with multiple sequenced genome for each species, we will restrict the search to refseq genome, exclude metagenome derived assemblies, exclude incompletely classified genomes and only keep one genome per species:

1# -t taxid
2# -rs refseq assemblies only
3# -ra <rank> (only keep one representative of this taxonomic rank)
4# -e exclude unclassified
5# -o output file
6get_ncbi_genome_assemblies.py -t 186803 -rs -ra species -e -o lachnospiraceae_assemblies.tab

The search will take some time because the script perform a search for each species to identify the highest quality genome (complete genome whenever possible). In this particular case (in october 2018), this search reported a total of 706 genomes assemblies classified in 478 different taxon identifiers. After filtration, the number of genomes was reduced to 106 (keeping only 1 genome/species and excluding metagenome dervied assemblies and incompletely classified genomes).

The resulting table contain the following colums:

1. target_rank ==> name of the target rank
2. rank_taxid ==> NCBI taxid of considered rank (the species in this case)
3. representative_species ==> the name of the selected representative genome
4. taxid_representative_species ==> the taxid of the selected representative genome
5. assembly_name ==> NCBI assembly name
6. RefSeq ==> refseq accession (if available)
7. Genbank ==> genbank accession
8. contig_count (number of contigs of the representative genome)
9. status ==> status of the selected genome (Complete, Scaffold, Contig,...)
10. repres ==> is the selected genome an NCBI reference or representative genome?
11. length ==> length of the assembly
12. clade
13. n_subtaxid_available
14. no_rank
15. superkingdom
16. phylum
17. class
18. order
19. family
20. genus

• the table can be manually edited to add/remove assemblies

• For instance, if you are interested in one particular genus (eg. Robinsoniella genus, NCBI taxid 588605), you can add all available assemblies for this specific genus. For this you need to use the "-d" parameter (Report complete list of assemblies).

1# -t taxid
2# -rs only refseq assemblies
3# -d get all assemblies available
4# -o output file
5get_ncbi_genome_assemblies.py -t 588605 -rs -d -o robinsoniella_assemblies.tab

• there is curently 4 assemblies available for this genus. You can add them to the initial table if not already present.

Download gbff files

• The 106 (+ eventually 3 Robinsoniella) assemblies can be downloaded based on their refseq accessions (e.g GCA_000424105.1).
• With the option -g, only the gbff file will be downloaded (sequence + annotations in genbank format)

1assembly_accession2genome.py -i GCA_000424105.1 -g

• a bash loop can be used to download all gbff from the lachnospiraceae_assemblies.tab table:

1for i in `sed -n '1!p' lachnospiraceae_assemblies.tab | cut -f 6`; do 
2    echo "Downloading $i ..."
3    assembly_accession2genome.py -i $i -g
4done

Prepare fasta files for orthofinder

Once downloaded, move all gbff files into a single directory and uncompress the files

1mkdir gbk_files
2mv GC*/*gz gbk_files
3cd gbk_files
4gunzip *gz

Now we need to extract predicted protein sequences from the gbff file for each genome

1for i in `ls *gbff`; do
2    # -i input gbff file
3    # -f put locus tags as headers in the fasta file
4    gbk2faa.py -i $i -f
5done

We can move all faa files into a separate folder

1mkdir -fp orthofinder/faa
2mv *faa orthofinder/faa
3cd orthofinder

Execute Orthofinder

1# -t number of threads
2orthofinder -t 6 -f faa/

Results

Orthofinder will take a while because it needs to BLAST all pairs of genomes (>24h in my single test with 6 threads). Once terminated, the phylogeny can be found in:

1cd faa/Results_Oct*/WorkingDirectory/Orthologues_*/SpeciesTree_rooted.txt

The tree can be visualized with figtree. Leaf labels will be the faa file names. They can be renamed using gbk files with the follwing script:

1phylo_convert_leaf_labels.py -i SpeciesTree_rooted.txt  -g <path/to/gbff/files/*gbff>

The resulting tree will be named "SpeciesTree_rooted_renamed.txt".