Trees should be placed in:

/srv/projects/gtdb/{release}/{domain}/pre_curation/bac120/{date}

The u_cluster_de_novo method of the GTDB Species Cluster Toolkit creates the files gtdb_reps_bac.lst and gtdb_reps_ar.lst which list the NCBI accession numbers for all GTDB bacterial and archaeal reference genomes.

Updates: starting with R232 WitChi is used to filter both the bac120 and ar53 MSA. Filtering is performed by first applying WitChi in touchdown mode, followed by removal of columns with large number of gaps or high/low AA conservation, and finally random subsampling to obtain a MSA of ~5,000 AA.

Two bacterial reference trees are inferred with FastTree as follows:

• Bacterial 120 (bac120) marker set:

  • MSA (./msa directory): gtdb -t 32 tree create --genome_batchfile gtdb_reps_bac.lst --guaranteed_batchfile gtdb_reps_bac.lst --marker_set_ids 1 --output ./msa --no_tree --individual --prefix gtdb_r<cur>_bac120
  • Witchi Chi2 filtering (./witchi directory): gtrim prune-chi2 -i ../msa/gtdb_r232_bac120_concatenated.faa -o . --max-residue-pruned 100000 --residues-to-prune-perc 1 --touchdown --threads-permutations 10 --threads-pruning 96
  • Filtering of gappy and high/low AA conserved columns (./filtered directory): genometreetk msa_filter ../msa/gtdb_r232_bac120_markers_info.tsv ../witchi-rust/gtdb_r232_bac120_concatenated_p100_s410_m100000_touchdown_pruned.tsv ../witchi-rust/gtdb_r232_bac120_concatenated_p100_s410_m100000_touchdown_pruned.fasta .
  • Tree inference (./snapshot directory): /usr/bin/time -v FastTreeMP -nosupport -wag -log gtdb_r232_bac120_fasttree.log ../filtered/filtered_msa.faa > gtdb_r232_bac120_fasttree.tree 2> gtdb_r232_bac120_fasttree.out
    • R232: required ~15 hrs and 87 GB (on Page)
  • Non-parametric bootstraps are calculated with the GenomeTreeTk. GTDB R214 required ~75 GB per bac120 bootstrap tree so it is possible to do 6 trees at a time on a 500 GB machine (e.g. Page). It is possible and recommended to calculate bootstrap trees across several machines. The code is smart enough to skip to the next replicate if intermediate files for a given replicate already exist (required genometreetk >=0.1.8). Bootstraps must be done with the WAG model and without the GAMMA model. For GenomeTreeTk v0.1.8 the command is as followed, but be warned older version of GenomeTreeTk are loaded on some machines:
    • Bootstrap trees (./bootstrap directory):
      • export OMP_NUM_THREADS=8
      • genometreetk bootstrap gtdb_r<cur>_bac120_fasttree.tree ../filtered/filtered_msa.faa . --mp -c 5
      • genometreetk bootstrap ../snapshot/gtdb_r232_bac120_fasttree.tree NONE . --boot_dir replicates -c 64

• Ribosomal protein set 2 (23 universal ribosomal proteins; modified from Rinke et al. 2013):

  • MSA: gtdb -t 32 tree create --genome_batchfile gtdb_reps_bac.lst --marker_set_ids 11 --output ./msa_rp2 --custom_msa_filters --quality_threshold 0 --completeness_threshold 0 --contamination_threshold 100 --cols_per_gene 10000 --min_perc_aa 40 --min_rep_perc_aa 40 --no_tree --individual --prefix gtdb_r<cur>_bac_rp2
  • Tree inference and non-parametric bootstraps are calculated in an analogous manner to the bac120 marker set.

Three archaeal reference trees are inferred with IQ-Tree using initial guide trees inferred with FastTree as follows:

• Archaeal 53 marker set (ar53; modified from Dombrowski et al. 2020):

  • MSA (./msa directory): gtdb -t 40 tree create --genome_batchfile gtdb_reps_ar.lst --guaranteed_batchfile gtdb_reps_ar.lst --marker_set_ids 19 --output msa --no_tree --prefix gtdb_r232_ar53_full --no_trim
  • Witchi Chi2 filtering (./witchi directory): witchi prune --file ./msa/gtdb_r<cur>_ar53_concatenated.faa --num_workers_chisq 4 --num_workers_permute 4 --pruning_algorithm quartic --top_n 135 --touchdown --max_residue_pruned 100000
  • Filtering of gappy and high/low AA conserved columns (./filtered directory): genometreetk msa_filter ./msa/gtdb_r232_ar53_full_markers_info.tsv ./witchi/gtdb_r232_ar53_full_concatenated_quartic_s135_touchdown_pruned.tsv ./witchi/gtdb_r232_ar53_full_concatenated_quartic_s135_touchdown_pruned.fasta filtered
  • Guide Tree inference (./snapshot directory): FastTreeMP -nosupport -wag -gamma -log gtdb_r<cur>_ar53_fasttree.log ../filtered/filtered_msa.faa > gtdb_r<cur>_ar53_fasttree.tree 2> gtdb_r<cur>_ar53_fasttree.out
  • Infer IQ-Tree (./snapshot directory): /usr/bin/time -v iqtree2 -nt 96 -s ../filtered/filtered_msa.faa -m LG+C10+F+G -ft gtdb_r<cur>_ar53_fasttree.tree -pre gtdb_r<cur>_ar53_iqtree
    • Chris recommended moving from C10 to at least C20
    • R232: required 303 GB and 127 hrs days (run on Howard)
  • Non-parametric bootstraps: ar53 IQ-Tree requires ~300 GB of memory for GTDB R232 which allows 1 tree to be inferred on a machine with 500 GB of memory. A simple script, run_iqtree_bootstraps.py, was used to infer the bootstrap trees. This script can currently be obtained from the previous GTDB release and needs to be customized for each release and reference tree. Bootstrap trees can be inferred across multiple servers and a simple canary file scheme is used to protect against processes on different servers working on the same bootstrap tree (not a guarantee, but highly unlikely). Ultimately, this script should be built into one of the toolkits. Once all bootstrap trees have been inferred, support values can be determined on the original tree with:
    • genometreetk bootstrap ../snapshot/gtdb_r232_ar53_iqtree.treefile NONE bootstraps/trees -c 50

• Archaeal 122 (ar122) marker set:

  • MSA: gtdb -t 32 tree create --genome_batchfile gtdb_reps_ar.lst --guaranteed_batchfile gtdb_reps_ar.lst --marker_set_ids 2 --output ./msa --no_tree --individual --prefix gtdb_r<cur>_ar122
  • Tree inference with support values are computed in an analogous manner to the ar53 marker set.

• Ribosomal protein set 2 (23 universal ribosomal proteins; modified from Rinke et al. 2013):

  • MSA: gtdb -t 16 tree create --genome_batchfile gtdb_reps_ar.lst --marker_set_ids 11 --output ./msa --custom_msa_filters --quality_threshold 0 --completeness_threshold 0 --contamination_threshold 100 --cols_per_gene 10000 --min_perc_aa 40 --min_rep_perc_aa 40 --no_tree --individual --prefix gtdb_r<cur>_ar_rp2
  • Tree inference with support values are computed in an analogous manner to the ar53 marker set.

Note: These instructions are for producing the final curation trees. Typically, we also provide curators with initial snapshot trees which lack non-parametric bootstraps so they can do an initial inspection. The following instructions can equally be applied for decorating and scaling snapshot trees. The u_pmc_species_names method of the GTDB Species Cluster Toolkit produces initial bacterial and archaeal taxonomy files which should be used to root and decorating the snapshot trees.

Trees can be rooted using:

• Bacteria: genometreetk outgroup ./bootstraps/gtdb_r<cur>_bac120_fasttree.bootstrap.tree ./u_pmc_species_names_bac/final_taxonomy.tsv p__Spirochaetota gtdb_r<cur>_bac120.rooted.tree

• Archaea: genometreetk outgroup gtdb_r<cur>_ar53_iqtree.bootstrap.tree ./u_pmc_species_names_ar/final_taxonomy.tsv p__Undinarchaeota gtdb_r<cur>_ar53.rooted.tree

  • It is important to verify the selected outgroup phylum is monophyletic in the tree. A warning will be produced if this is not the case. Other recommended outgroup phyla are p__Thermotogota, p__Patescibacteria, or p__Chloroflexota for the bacterial trees and p__QMZS01 for the archaeal trees.

Dummy leaf nodes used to provide an internal node for species labels useful during manual curation can be added to reference trees and their corresponding taxonomy files using:

• gtdb_validation_tk add_dummy gtdb_r<cur>_bac120.rooted.tree --taxonomy_file ./u_pmc_species_names_bac/final_taxonomy.tsv
• This produces the output files gtdb_r<cur>_bac120.rooted.dummy_nodes.tree and final_taxonomy.dummy_nodes.tsv.
• Note: perhaps this functionality should be moved to genometreetk or gtdb_migration_tk

Trees can be decorated with PhyloRank >= v0.1.11 as follows:

• phylorank decorate gtdb_r<cur>_bac120.rooted.dummy_nodes.tree final_taxonomy.dummy_nodes.tsv gtdb_r<cur>_bac120.decorated.tree --skip_rd_refine --gtdb_metadata gtdb_r<cur>_metadata.updated_reps.<date>.tsv

PhyloRank is used to calculated RED values, the RED plot, and RED scaled tree:

• phylorank outliers gtdb_r<cur>_bac120.decorated.tree final_taxonomy.dummy_nodes.tsv phylorank_outliers --fmeasure_table gtdb_r<cur>_bac120.decorated.tree-table --highlight_polyphyly --dpi 300

Curomatic should be run to perform initial adjustments that reduce manual curation efforts:

• cut -f1,101,154 /srv/db/gtdb/metadata/<cur>/representatives/sp_cluster_update/gtdb_r<cur>_metadata.updated_reps.tsv > curatomatic_metadata.tsv
• python3 -m curatomatic --meta curatomatic_metadata.tsv ./phylorank_outliers/gtdb_r232_bac120.decorated.red_decorated.tree ./phylorank_outliers/gtdb_r232_bac120.decorated.dict curatomatic 2>&1 |tee curatomatic.out