Generate GTDB metadata - Ecogenomics/GTDBNCBI GitHub Wiki
Metadata for NCBI genomes needs to be generated after syncing with RefSeq and GenBank. Ideally, metadata would only be generated when required, but in some cases it is calculated de novo for all genomes for simplicity. Some data should be calculated de novo for all genomes on a regular bases. For example, 16S classification when there is a new SILVA/LTP database.
Note: commands should be run from
/srv/db/gtdb/metadata/<release>/processing
NCBI gene calling is conservative. This can negatively impact the CheckM estimates and results in valid marker genes being missed. As such, genes are called de novo using Prodigal to ensure they are called in a consistent manner:
- Predict translation tables using gTranslate and compare to NCBI translation table data:
gtdb_migration_tk trans_table -g ../../../genomes/ncbi/<release>/genome_dirs.tsv -t ../ncbi/taxonomy/standardised_taxonomy/ncbi_r237_standardized.tsv -o translation_tables -c 96 -l logs/trans_table-<server>.loggTranslate is run on batches of genomes (default = 10,000) so it can be run across multiple servers and to allow restarting if required. A comparison with NCBI translation table data is performed and conflicting gTranslate predictions written to ncbi_tt_conflict.tsv. These should be reviewed to ensure gTranslate is correct. If gTranslate appears to be incorrect, a corrections file named needs to be manually produced which is a two column TSV file with the headers genome_id and translation_table. This file should be named tt_corrections.tsv and placed in the translation_tables output directory. This file can be provided in the prodigal step and will override the gTranslate prediction. Hints on how to resolve conflicts is given on the Resolving gTranslate and NCBI conflicts page.
Requirements: took ~23 hrs for the ~1.3 million GTDB R237 genomes when run using 96 CPUs on Chifley, Howard, Holt, Hurley, and Page (i.e. effectively using 480 CPUs)
- Call genes using Prodigal:
Prodigal version: v2.6.3
gtdb_migration_tk prodigal -g ../../../genomes/ncbi/<release>/genome_dirs.tsv --trans_table translation_tables/gtranslate.translation_table_summary.tsv.gz --tt_override ./translation_tables/tt_corrections.tsv -o prodigal -c 96 -l logs/prodigal-<server>.logProdigal is run on batches of genomes (default = 10,000) so it can be run across multiple servers and to allow restarting if required. Final Prodigal results are stored in a prodigal directory under the directory paths specified in genome_dirs.tsv.
Note: Genomes without a predicted translation table from the previous step are skipped (see ./translation_tables/gtranslate_no_prediction.tsv).
Requirements: took ~10 hrs for the ~300,000 GTDB R237 genomes without prior Prodigal results when run using 96 CPUs on Chifley, Howard, Holt, Hurley, and Page (i.e. effectively using 480 CPUs)
hmmer version: v3.4
Note: it is important to use HMMER v3.4 to obtain the same results as previous GTDB versions
GTDB marker genes are identified using a reduced set of Pfam and TIGRfam HMMs (DB Lite) that have been evaluated to give near identical results to using the full set of HMMs. Using DB Lite is recommended and is the default of GTDB:
DB Lite (recommended)
gtdb_migration_tk hmmsearch -g ../../../genomes/ncbi/<release>/genome_dirs.tsv -r ../../../genomes/ncbi/<release>/report.log --db pfam --hmm_db_path ../../../marker_genes/hmms_extended_pfam33.1_tigr15 -o hmmsearch -c 96 -l logs/hmmsearch_pfam-<server>.log
gtdb_migration_tk hmmsearch -g ../../../genomes/ncbi/<release>/genome_dirs.tsv -r ../../../genomes/ncbi/<release>/report.log --db tigrfam --hmm_db_path ../../../marker_genes/hmms_extended_pfam33.1_tigr15/tigrfam.hmm -o hmmsearch -c 96 -l logs/hmmsearch_tigrfam-<server>.loghmmsearch is run on batches of genomes (default = 10,000) so it can be run across multiple servers and to allow restarting if required. Final results are stored under the prodigal directory of each genome.
Requirements: took ~9 hrs for the ~300,000 GTDB R237 genomes without prior Pfam results when run using 96 CPUs on Chifley, Howard, Holt, Hurley, and Page (i.e. effectively using 480 CPUs) Requirements: took ~15 hrs for the ~300,000 GTDB R237 genomes without prior TIGRfam results when run using 96 CPUs on Chifley, Howard, Holt, Hurley, and Page (i.e. effectively using 480 CPUs)
Full DB (for testing/evaluation purposes only)
gtdb_migration_tk hmmsearch -g ../../../genomes/ncbi/<release>/genome_dirs.tsv -r ../../../genomes/ncbi/<release>/report.log --db pfam --hmm_db_path /srv/db/pfam/33.1 -o hmmsearch -c 96 -l logs/hmmsearch_pfam.log
gtdb_migration_tk hmmsearch -g ../../../genomes/ncbi/<release>/genome_dirs.tsv -r ../../../genomes/ncbi/<release>/report.log --db tigrfam --hmm_db_path /srv/db/tigrfam/15.0/TIGRFAMs_15.0_HMM -o hmmsearch -c 96 -l logs/hmmsearch_tigrfam.logNote: if required, HMMER tophit files can be regenerated without recalling genes using the
gtdb_migration_tk top_hitcommand.
Calculate genome (e.g., GC) and protein (e.g., gene count) attributes:
gtdb_migration_tk genomic_metadata -g ../../../genomes/ncbi/<release>/genome_dirs.tsv -o genomic_metadata -c 96 -l logs/genomic_metadata.logQuestion for Pierre: I modified this command to take the previous release domain report and the NCBI taxonomy file for the current release. These are used to determine the domain for each genome. This works, but doesn't seem right. Presumably we should use our domain marker sets to determine the correct domain for each genome, but this doesn't happen until later in this workflow. I haven't look into if the workflow can be modified to get this data earlier, but does seem like we need to know the domain of each genome at this point. It is used here and for the 5S/15S/23S rRNA step.
Identification of transfer RNAs (tRNAs) is done using trnascan:
module load trnascan-se (optional)
gtdb_migration_tk trnascan -g ../../../genomes/ncbi/<release>/genome_dirs.tsv -d ../../<prev_release>/representatives/sp_cluster_update/gtdb_<prev_release>_domain_report.tsv -t ../ncbi/taxonomy/standardised_taxonomy/ncbi_<release>_standardized.tsv -o trnascan -c 96 -l logs/trnascan-<server>.logtRNAscan-SE is run on batches of genomes (default = 10,000) so it can be run across multiple servers and to allow restarting if required. Final results are stored in a trna directory under the directory paths specified in genome_dirs.tsv.
Requirements: took ~20 hrs for the ~300,000 GTDB R237 genomes without prior tRNA results when run using 96 CPUs on Chifley, Howard, Holt, Hurley, and Page (i.e. effectively using 480 CPUs)
Identification and classification of 5S, 16S, and 23S rRNA genes can be done following the steps at: Generate Ribosomal Metadata
Genome quality is estimated using CheckM v1 and v2.
- Estimate genome quality using CheckM v1:
gtdb_migration_tk checkm -g ../../../genomes/ncbi/<release>/genome_dirs.tsv --report ../../../genomes/ncbi/<release>/report.log -o checkm -c 64 -l logs/checkm-<server>.logCheckM v1 is run on batches of genomes (default = 10,000) so it can be run across multiple servers and to allow restarting if required.
Requirements: took ~3 days for the ~300,000 new/updated GTDB R237 genomes when run using 96 CPUs on Chifley, Howard, Holt, Hurley, and Page (i.e. effectively using 480 CPUs)
- Estimate genome quality using CheckM v2:
gtdb_migration_tk checkm2 -g ../../../genomes/ncbi/<release>/genome_dirs.tsv --report ../../../genomes/ncbi/<release>/report.log -o checkm2 -c 96 -l logs/checkm2-<server>.logCheckM v2 is run on batches of genomes (default = 50,000) so it can be run across multiple servers and to allow restarting if required.
Requirements: took ~21 hours for the ~300,000 new/updated GTDB R237 genomes when run using 96 CPUs on Chifley, Howard, Holt, Hurley, and Page (i.e. effectively using 480 CPUs)
Next steps of the update are under Update database metadata.