Introduction

GTDB（Genome Taxonomy Database）是一個(gè)專注于基因組分類學(xué)的數(shù)據(jù)庫，旨在提供高質(zhì)量的細(xì)菌、古細(xì)菌、葉綠體和線粒體的分類信息。通過使用全基因組信息，GTDB提供了多級別的分類，包括域、界、綱、目、科、屬等級別。其系統(tǒng)發(fā)育分類系統(tǒng)使得研究人員能夠深入了解生物的演化和功能。GTDB主要以細(xì)菌中普遍存在的120個(gè)單拷貝蛋白質(zhì)（bac120）和古菌中53個(gè)（arc53，從R07-RS207開始）標(biāo)記蛋白為基礎(chǔ)，在對多分組類別消歧后，根據(jù)相對演化散度標(biāo)準(zhǔn)化和分級，得到基因組分類結(jié)果。

GTDB定期更新，與其他生物信息數(shù)據(jù)庫整合，為基因組學(xué)研究和微生物生態(tài)學(xué)等領(lǐng)域的研究提供了全面且可靠的分類信息。GTDB數(shù)據(jù)庫已經(jīng)于2018/2020連續(xù)發(fā)表兩篇Nature Biotechnology，2021年又發(fā)表在Nature Microbiology和Nucleic Acids
Research上，是目前最權(quán)威的微生物分類數(shù)據(jù)庫之一。
GTDB配套軟件在2019/2022發(fā)表了兩篇Bioinformatcs，時(shí)提供了物種分類、進(jìn)化樹構(gòu)建的全套工具。

官網(wǎng)：https://gtdb.ecogenomic.org/

GTDB數(shù)據(jù)庫

可以看到上面的官網(wǎng)主頁，目前最新版的GTDB數(shù)據(jù)庫包含了584,382個(gè)細(xì)菌基因組和12,477個(gè)古菌基因組（2024.04.24，Release 09-RS220），在屬、種分辨率水平上描述了190多個(gè)門。

Browsers

GTDB支持以Alphabetical
table和Tree兩種形式來探索或查詢整個(gè)GTDB數(shù)據(jù)庫，包括分類信息和基因組信息。

此外，可以在Taxon History輸入 GTDB 或 NCBI
分類單元，以可視化分配給該分類單元的基因組在版本之間如何變化，比如p__Proteobacteria從Release
214開始改名為p__Pseudomonadota，在不斷壯大的同時(shí)也有一部分genome被分配給了其他門。

Tools

GTDB提供了一些配套的工具：

• GTDB-Tk：下一章節(jié)詳細(xì)介紹
• Third-party tools：AnnoTree，IDTAXA，SingleM，sourmash
• Third-party GTDB
  databases：HUMANn2，Kraken2，Bracken2，Centrifuge，Kraken，DADA2等
• Third-party
  scripts：tax_from_gtdb.py，genome_updater，gtdb_to_taxdump
• FastANI calculator：可用于原核生物 NCBI
  基因組的ANI比較，一次最多可以進(jìn)行 1,000 次成對比較，例如下圖。

Downloads

該頁面提供了GTDB數(shù)據(jù)庫和配套軟件的下載鏈接。

Statistics

統(tǒng)計(jì)頁面提供了各個(gè)版本的GTDB數(shù)據(jù)庫的各類基礎(chǔ)信息：

1. Taxon overview：GTDB R220 包含 596,859 個(gè)基因組，組織成 113,104
   個(gè)物種簇。

2. Species overview：GTDB R220 由 584,382 個(gè)細(xì)菌和 12,477
   個(gè)古細(xì)菌基因組組成，分為 107,235 個(gè)細(xì)菌和 5,869 個(gè)古細(xì)菌物種簇。

3. Genome categories：GTDB
   分類群由分離基因組（isolate）、宏基因組組裝基因組 (MAG)
   和單擴(kuò)增基因組 (SAG) 組成。

4. GTDB species representatives：每個(gè) GTDB 物種簇都由一個(gè)基因組表示。

5. Quality of GTDB representative genomes：被選為 GTDB
   物種代表的基因組的質(zhì)量。使用 CheckM 估計(jì)基因組完整性和污染，并根據(jù)
   MIMAG
   基因組標(biāo)準(zhǔn)進(jìn)行著色。一般來說，代表性基因組僅限于滿足completeness - 5*contamination >50，存在一些例外情況，如下面所示。
   

6. Taxa with the largest number of species：給出了每個(gè)分類等級包含最多
   GTDB
   物種簇的分類單元，門水平是常見的Pseudomonadota，Bacillota，Bacteroidota，Actinomycetota等。

7. Taxa with the largest number of sequenced genomes，與6類似。

8. Relative evolutionary
   divergence（RED）：提供了從門到屬的每個(gè)分類等級的類群的相對進(jìn)化分歧（RED）。
   RED 值提供了相對時(shí)間的操作近似值，其中現(xiàn)有分類群存在于當(dāng)前
   (RED=1)，最后一個(gè)共同祖先出現(xiàn)在過去的固定時(shí)間
   (RED=0)，內(nèi)部節(jié)點(diǎn)根據(jù)以下值在這些值之間線性插值譜系特定的進(jìn)化速率。

9. Comparison of GTDB and NCBI taxa：GTDB
   物種代表性基因組和所有具有分配的 NCBI 分類法的 GTDB 基因組的 GTDB 和
   NCBI
   分類學(xué)分配的比較。對于每個(gè)分類等級，如果一個(gè)分類單元在兩個(gè)分類標(biāo)準(zhǔn)中名稱相同，則該分類單元被分類為unchange；如果
   GTDB 分類標(biāo)準(zhǔn)提供了 NCBI
   分類標(biāo)準(zhǔn)中缺少的名稱信息，則該分類單元被分類為passive change；如果兩個(gè)分類標(biāo)準(zhǔn)之間的名稱不同，則該分類單元被分類為actively changed?？梢钥吹?，在genome水平上GTDB和NCBI還是有很多名稱不對應(yīng)的，可以考慮使用gtdb_to_taxdump進(jìn)行轉(zhuǎn)化。
   

10. Genomic statistics：GTDB 物種代表性基因組和 GTDB
    中所有基因組的關(guān)鍵基因組統(tǒng)計(jì)數(shù)據(jù)。

    
    

11. Nomenclatural types per
    rank：顯示了每個(gè)分類等級的占位符與拉丁化名稱的細(xì)分。

Forum

GTDB還專門搭建了一個(gè)用戶社區(qū)和討論區(qū)：

• Announcement：用于公布最新 GTDB
  新聞、即將舉行的會(huì)議和以分類學(xué)為重點(diǎn)的手稿的開放論壇。
• Taxonomy
  (Bacteria/Archaea)：這個(gè)分類專注于討論與GTDB分類法相關(guān)的內(nèi)容。用戶在這里討論分類法的具體細(xì)節(jié)、與其他數(shù)據(jù)庫（如NCBI）的比較、分類標(biāo)準(zhǔn)的更新等。
• GTDB-Tk：專注于GTDB-Tk工具包的使用和討論。用戶可以在這里提出關(guān)于GTDB-Tk的使用問題、分享使用經(jīng)驗(yàn)、討論工具包的功能和性能等。
• GTDB：與基因組分類數(shù)據(jù)庫 (GTDB) 相關(guān)的一般問題。
• Website feedback：討論該網(wǎng)站、其組織、其運(yùn)作方式以及如何改進(jìn)它。

Help

幫助頁面也提供了一些非常有用的信息。

• About：GTDB簡要介紹，管理團(tuán)隊(duì)以及發(fā)表論文

• FAQ：列出了一些常見問題，例如：

  1. 為什么門名稱的后綴更改為-ota？Whitman et al. (2018)
     建議像其他等級一樣規(guī)范門等級的后綴。
  2. 為什么一些family和更高級別的名字以字母后綴結(jié)尾？屬等級以上的分類單元名稱附加字母后綴表示屬于以下類別的類群：
     1. GTDB 參考樹中不是單系類群，但存在其他證據(jù)表明它們是單系類群；
     2. 在版本之間位置不穩(wěn)定的組。

• Methods：對GTDB數(shù)據(jù)庫構(gòu)建過程中的方法進(jìn)行具體描述，例如：

  1. Gene
     identification，使用Prodigal鑒定基因，HMMER鑒定標(biāo)記基因（來自Pfam和TIGRFAM）并序列比對。
  2. Multiple sequence alignment，細(xì)菌和古菌多重序列比對 (MSA) 分別由
     120 個(gè) (bac120) 或 53 個(gè) (arc53) 系統(tǒng)發(fā)育信息標(biāo)記串聯(lián)而成。
  3. Tree inference，細(xì)菌參考樹是在 WAG 模型下使用
     FastTree。古菌參考樹是根據(jù) PMSF 模型下的 IQ-Tree，并使用
     FastTree推斷初始樹以恢復(fù)單系譜系。
  4. Identifying 16S rRNA sequences，使用 nhmmer以及來自 RFAM
     數(shù)據(jù)庫的 16S rRNA 模型（RF00177）。
  5. Average nucleotide identity，skani（從版本 220
     開始），之前是FastANI。
  6. RED
     normalization：相對進(jìn)化分歧用于幫助描繪物種等級之上的分類群，RED的概念是相同等級的類群在過去應(yīng)該同時(shí)出現(xiàn)。RED
     提供了相對時(shí)間的操作近似值，其中現(xiàn)有分類群存在于當(dāng)前 (RED =
     1)，最后一個(gè)共同祖先出現(xiàn)在過去的固定時(shí)間 (RED =
     0)，內(nèi)部節(jié)點(diǎn)根據(jù)譜系在這些值之間線性插值-特定的進(jìn)化速率。

• API：功能較多，但文檔不是很詳細(xì)，未來也許會(huì)用上

GTDB-Tk

GTDB-Tk 是一個(gè)軟件工具包，用于根據(jù)基因組數(shù)據(jù)庫分類法 (GTDB)
對細(xì)菌和古細(xì)菌基因組進(jìn)行客觀分類。
它旨在與最新進(jìn)展相結(jié)合，允許直接從環(huán)境樣本中獲得數(shù)百或數(shù)千個(gè)宏基因組組裝基因組
(MAG)，還可以應(yīng)用于分離和單細(xì)胞基因組。

地址：https://github.com/Ecogenomics/GTDBTk

官方文檔：https://ecogenomics.github.io/GTDBTk/

Install

首先要安裝軟件和下載對應(yīng)的Reference
data，最新的R220數(shù)據(jù)庫需要使用v2.4.0的GTDB-Tk（2024.04.24）。如果要和之前項(xiàng)目一致，可以使用R214版本數(shù)據(jù)庫對應(yīng)v2.3.2的GTDB-Tk。

1.conda

conda安裝比較簡單，但最好新建一個(gè)環(huán)境，可能會(huì)比較慢（mamba代替conda會(huì)快一點(diǎn)）：

# using conda
conda create -n gtdbtk-2.3.2 -c conda-forge -c bioconda gtdbtk=2.3.2
# 下載 GTDB-Tk 參考數(shù)據(jù)并為其命名
download-db.sh
# 激活GTDB-Tk conda環(huán)境
conda 激活 gtdbtk-2.3.2
# 將環(huán)境變量設(shè)置為包含GTDB-Tk參考數(shù)據(jù)的目錄
conda env config vars set GTDBTK_DATA_PATH="/path/to/release/package/";

2.pip

pip可以快速安裝gtdbtk包，但是它的第三方依賴項(xiàng)需要手動(dòng)安裝一下
（不是很復(fù)雜，我用的pip安裝，因?yàn)椴幌胍嗟腸onda環(huán)境，太冗余了）：

python -m pip install gtdbtk==2.3.2

第三方依賴項(xiàng)：

• Prodigal（原核基因識別和翻譯起始位點(diǎn)鑒定）
• HMMER（加速族HMM 搜索）
• pplacer（線性時(shí)間最大似然和貝葉斯系統(tǒng)發(fā)育將序列放置到固定參考樹上）
• skani (v2.4.0開始依賴，通過 skani
  的稀疏鏈接進(jìn)行快速、穩(wěn)健的宏基因組序列比較，比fastANI還快，下次有機(jī)會(huì)講一下這個(gè))
• FastTree（用于大對齊的近似最大似然樹）
• Mash（快速估計(jì)基因組和宏基因組距離）

這些都是很常用的軟件，沒有的話用conda安裝一下。

GTDB-Tk 需要將 GTDBTK_DATA_PATH 環(huán)境變量設(shè)置為包含未存檔的 GTDB-Tk
參考數(shù)據(jù)的目錄。

echo 'export GTDBTK_DATA_PATH=/path/to/release/package/' >> ~/.bashrc
source ~/.bashrc

3.docker

如果用過docker的話，這個(gè)也非常簡單：

#下載解壓的 GTDB-Tk 參考數(shù)據(jù)到您選擇的目錄，例如：/host/release_data/
mv /path/to/release/package/ /host/release_data/
#為 GTDB-Tk 創(chuàng)建本地 I/O 目錄
mkdir /host/gtdbtk_io/
# 運(yùn)行容器
docker run -v /host/gtdbtk_io:/data -v /host/release_data:/refdata ecogenomic/gtdbtk classify_wf --genome_dir /data/genomes --out_dir /data/output

Reference data

GTDB-Tk 需要約 110G
的外部數(shù)據(jù)，需要下載和解檔。直接wget很慢的話，也可以先下載到本地，再通過FTP軟件或scp上傳到服務(wù)器：

數(shù)據(jù)地址：https://data.gtdb.ecogenomic.org/releases/latest/auxillary_files/gtdbtk_package/

對于完整包(full package)：

wget https://data.ace.uq.edu.au/public/gtdb/data/releases/latest/auxillary_files/gtdbtk_package/full_package/gtdbtk_data.tar.gz
tar xvzf gtdbtk_data.tar.gz

對于拆分包(split package)：

1. 確保 GTDB-Tk 發(fā)布數(shù)據(jù)的所有部分都位于同一目錄中。
2. 打開終端或命令提示符。
3. 導(dǎo)航到包含 GTDB-Tk 發(fā)布數(shù)據(jù)部分的目錄。
4. 使用以下命令將所有部分連接到一個(gè)存檔中：
   cat gtdbtk_r220_data.tar.gz.part_* > gtdbtk_r220_data.tar.gz
5. 命令執(zhí)行完成后，將在同一目錄中擁有一個(gè)名為“gtdbtk_r220_data.tar.gz”的存檔文件。

Check

檢查數(shù)據(jù)庫和軟件，結(jié)果顯示OK和Done表示正常（檢查fastani這一步要很久，文件非常多。所以如果三方軟件都OK的話也可以直接運(yùn)行下一個(gè)示例測試）:

gtdbtk check_install

[2024-05-20 14:10:37] INFO: GTDB-Tk v2.3.2
[2024-05-20 14:10:37] INFO: gtdbtk check_install
[2024-05-20 14:10:37] INFO: Using GTDB-Tk reference data version r214: /data/home/jianglab/share/GTDB_r214/release214/
[2024-05-20 14:10:37] INFO: Running install verification
[2024-05-20 14:10:37] INFO: Checking that all third-party software are on the system path:
[2024-05-20 14:10:37] INFO:          |-- FastTree         OK
[2024-05-20 14:10:37] INFO:          |-- FastTreeMP       OK
[2024-05-20 14:10:37] INFO:          |-- fastANI          OK
[2024-05-20 14:10:37] INFO:          |-- guppy            OK
[2024-05-20 14:10:37] INFO:          |-- hmmalign         OK
[2024-05-20 14:10:37] INFO:          |-- hmmsearch        OK
[2024-05-20 14:10:37] INFO:          |-- mash             OK
[2024-05-20 14:10:37] INFO:          |-- pplacer          OK
[2024-05-20 14:10:37] INFO:          |-- prodigal         OK
[2024-05-20 14:10:37] INFO: Checking integrity of reference package: /data/home/jianglab/share/GTDB_r214/release214/
[2024-05-20 14:10:38] INFO:          |-- pplacer          OK
[2024-05-20 14:10:38] INFO:          |-- masks            OK
[2024-05-20 14:10:39] INFO:          |-- markers          OK
[2024-05-20 14:10:39] INFO:          |-- radii            OK
[2024-05-20 14:10:41] INFO:          |-- msa              OK
[2024-05-20 14:10:41] INFO:          |-- metadata         OK
[2024-05-20 14:10:41] INFO:          |-- taxonomy         OK
[2024-05-20 15:13:45] INFO:          |-- fastani          OK
[2024-05-20 15:13:45] INFO:          |-- mrca_red         OK
[2024-05-20 15:13:45] INFO: Done.

注意GTDB-Tk的硬件要求，關(guān)鍵是Memory一定要給足夠，不然會(huì)程序會(huì)中途終止：

測試流程，無報(bào)錯(cuò)確定程序的軟件和數(shù)據(jù)庫正常（此示例是3個(gè)古菌genome，建議分配60GB內(nèi)存來運(yùn)行，另外gtdbtk_test必須是不存在的或空文件夾）:

gtdbtk test --out_dir gtdbtk_test

[2024-05-20 14:34:51] INFO: GTDB-Tk v2.3.2
[2024-05-20 14:34:51] INFO: gtdbtk test --out_dir gtdbtk_test
[2024-05-20 14:34:51] INFO: Using GTDB-Tk reference data version r214: /data/home/jianglab/share/GTDB_r214/release214/
[2024-05-20 14:34:51] INFO: Command: gtdbtk classify_wf --skip_ani_screen --genome_dir gtdbtk_test/genomes --out_dir gtdbtk_test/output --cpus 1 -f
    <TEST OUTPUT>     [2024-05-20 14:34:52] INFO: gtdbtk classify_wf --skip_ani_screen --genome_dir gtdbtk_test/genomes --out_dir gtdbtk_test/    
    <TEST OUTPUT>     [2024-05-20 14:34:52] INFO: Using GTDB-Tk reference data version r214: /data/home/jianglab/share/GTDB_r214/release214/      
    <TEST OUTPUT>     [2024-05-20 14:42:06] INFO: Note that Tk classification mode is insufficient for publication of new taxonomic designatio    
    <TEST OUTPUT>     [2024-05-20 14:42:06] INFO: Done.     
[2024-05-20 14:42:09] INFO: Test has successfully finished.

我分配了4核共60GB內(nèi)存運(yùn)行gtdbtk test，約439s完成，消耗了~45GB內(nèi)存。

使用方法

先看看help文檔：

$ gtdbtk -h
## GTDB-Tk v2.4.0 幫助文檔

  Workflows:
    classify_wf -> 通過在 GTDB 參考樹中的位置對基因組進(jìn)行分類
                     (ani_screening -> identify -> align -> classify)
    de_novo_wf  -> 推斷新的樹并使用 GTDB 分類進(jìn)行修飾
                     (identify -> align -> infer -> root -> decorate)

  Methods:
    identify -> 識別基因組中的標(biāo)記基因
    align    -> 創(chuàng)建多序列比對
    classify -> 確定基因組的分類
    infer    -> 從多序列比對中推斷樹
    root     -> 使用外群對樹進(jìn)行根定
    decorate -> 用 GTDB 分類修飾樹

  Tools:
    infer_ranks        -> 利用RED建立內(nèi)部節(jié)點(diǎn)的分類等級
    ani_rep            -> 計(jì)算ANI到GTDB的代表性基因組
    trim_msa           -> 基于掩碼修剪未修剪的MSA文件
    export_msa         -> 導(dǎo)出未修剪的古生菌或細(xì)菌MSA文件
    remove_labels      -> 從Newick樹中刪除標(biāo)簽(引導(dǎo)值，節(jié)點(diǎn)標(biāo)簽)
    convert_to_itol    -> 將GTDB-Tk Newick樹轉(zhuǎn)換為iTOL樹
    convert_to_species -> 將GTDB基因組id轉(zhuǎn)換為GTDB物種名稱


  Testing:
    test          -> 用3個(gè)古細(xì)菌基因組驗(yàn)證classify_wf管道
    check_install -> 驗(yàn)證第三方程序和GTDB參考包

對于每一個(gè)模塊可以單獨(dú)查看help，例如最重要的classify_wf工作流：

$ gtdbtk classify_wf -h
usage: gtdbtk classify_wf (--genome_dir GENOME_DIR | --batchfile BATCHFILE) --out_dir OUT_DIR
                          (--skip_ani_screen | --mash_db MASH_DB) [--no_mash] [--mash_k MASH_K] [--mash_s MASH_S]
                          [--mash_v MASH_V] [--mash_max_distance MASH_MAX_DISTANCE] [-f] [-x EXTENSION]
                          [--min_perc_aa MIN_PERC_AA] [--prefix PREFIX] [--genes] [--cpus CPUS] [--pplacer_cpus PPLACER_CPUS]
                          [--force] [--scratch_dir SCRATCH_DIR] [--write_single_copy_genes] [--keep_intermediates]
                          [--min_af MIN_AF] [--tmpdir TMPDIR] [--debug] [-h]

mutually exclusive required arguments:
  --genome_dir GENOME_DIR 包含F(xiàn)ASTA格式基因組文件的目錄
  --batchfile BATCHFILE 描述基因組的文件路徑 - 2或3列的制表符分隔文件（FASTA文件、基因組ID、翻譯表[可選]）

required named arguments:
  --out_dir OUT_DIR     輸出文件的目錄

mutually exclusive required arguments:
  --skip_ani_screen     跳過ani篩選步驟，使用mash和skani對基因組進(jìn)行分類。（默認(rèn)：False）
  --mash_db MASH_DB     保存/讀取Mash參考草圖數(shù)據(jù)庫（.msh）的路徑

optional Mash arguments:
  --no_mash             跳過使用Mash進(jìn)行的基因組預(yù)過濾（默認(rèn)：False）
  --mash_k MASH_K       k-mer大小[1-32]（默認(rèn)：16）
  --mash_s MASH_S       最大非冗余哈希數(shù)（默認(rèn)：5000）
  --mash_v MASH_V       保留的最大p值[0-1]（默認(rèn)：1.0）
  --mash_max_distance MASH_MAX_DISTANCE 選擇潛在GTDB基因組作為用戶基因組代表的最大Mash距離。（默認(rèn)：0.15）

optional arguments:
  -f，--full_tree       在分類步驟中使用未分割的細(xì)菌樹；這是原始的GTDB-Tk方法（版本<2），
                        需要超過320 GB的RAM來加載參考樹（默認(rèn)：False）
  -x，--extension EXTENSION 要處理的文件擴(kuò)展名，gz = gzipped（默認(rèn)：fna）
  --min_perc_aa MIN_PERC_AA 
                        排除在MSA中AA百分比不足此值的基因組（包含界限）（默認(rèn)：10）
  --prefix PREFIX       所有輸出文件的前綴（默認(rèn)：gtdbtk）
  --genes               表示輸入文件包含預(yù)測的蛋白質(zhì)（作為氨基酸），跳過基因調(diào)用。
                        警告：此標(biāo)志將跳過ANI比較步驟（ani篩選和分類）。（默認(rèn)：False）
  --cpus CPUS           要使用的CPU數(shù)量（默認(rèn)：1）
  --pplacer_cpus PPLACER_CPUS pplacer放置過程中使用的CPU數(shù)量
  --force               如果單個(gè)基因組發(fā)生錯(cuò)誤，則繼續(xù)處理（默認(rèn)：False）
  --scratch_dir SCRATCH_DIR
                        通過寫入磁盤來減少pplacer內(nèi)存使用（較慢）
  --write_single_copy_genes
                        輸出未對齊的單拷貝標(biāo)記基因（默認(rèn)：False）
  --keep_intermediates  在最終目錄中保留中間文件（默認(rèn)：False）
  --min_af MIN_AF       將基因組分配到物種簇的最小對齊比例（默認(rèn)：0.5）
  --tmpdir TMPDIR       指定臨時(shí)文件的備用目錄（默認(rèn)：/tmp）
  --debug               創(chuàng)建用于調(diào)試目的的中間文件（默認(rèn)：False）

classify_wf 注釋物種

classify_wf全稱為 Classify workflow，即分類工作流程。

分類工作流程由四個(gè)步驟組成：ani_screen（v2.4.0開始加入）, identify, align, classify.

1. ani_screen 步驟將用戶基因組與由所有 GTDB 代表性基因組組成的 Mash 數(shù)據(jù)庫進(jìn)行比較，然后使用 skani 驗(yàn)證最佳的 mash 命中。使用 skani 分類的用戶基因組不會(huì)運(yùn)行管道的其余部分（識別、對齊、分類），并在摘要文件中報(bào)告。
2. identify步驟使用 Prodigal 調(diào)用基因，并使用 HMM 模型和 HMMER 包來識別用于系統(tǒng)發(fā)育推斷的 120 個(gè)細(xì)菌和 53 個(gè)古細(xì)菌標(biāo)記基因 (Parks et al., 2018)。通過將標(biāo)記基因與其各自的 HMM 模型進(jìn)行比對來獲得多重序列比對 (MSA)。
3. align步驟連接比對的標(biāo)記基因，并將連接的 MSA 過濾為大約 5,000 個(gè)氨基酸。
4. 最后，classify步驟使用 pplacer 查找 GTDB-Tk 參考樹中每個(gè)基因組的最大似然位置。 GTDB-Tk 根據(jù)每個(gè)基因組在參考樹中的位置、相對進(jìn)化分歧和/或與參考基因組的平均核苷酸同一性 (ANI) 對每個(gè)基因組進(jìn)行分類。

# 輸入基因組目錄
gtdbtk classify_wf --genome_dir genomes/ --out_dir classify_wf_out --cpus 16

可以根據(jù)下圖估計(jì)自己的任務(wù)運(yùn)行時(shí)間，評估資源分配：

output

GTDB-Tk 在classify_wf管道的各個(gè)階段輸出許多文件，具體格式和示例可以在官網(wǎng)查看：
https://ecogenomics.github.io/GTDBTk/files/index.html

結(jié)果目錄(classify_wf_out/output)中主要結(jié)果文件 (ar53為古菌相關(guān)文件，bac120為細(xì)菌相關(guān)文件)：

• classify/gtdbtk.ar53.classify.tree 樹文件，Newick格式，可以用于系統(tǒng)發(fā)育樹的繪制。
• classify/gtdbtk.ar53.summary.tsv 分類結(jié)果文件，主要結(jié)果，可以用于系統(tǒng)發(fā)育樹的注釋。
• identify/gtdbtk.ar53.markers_summary.tsv 使用標(biāo)記基因情況
• align/gtdbtk.ar53.msa.fasta.gz 多序列對齊文件(輸入+參考)
• align/gtdbtk.ar53.user_msa.fasta.gz 多序列對齊文件(輸入)

繪制系統(tǒng)發(fā)育樹可以參考我之前寫的R繪制優(yōu)美的進(jìn)化樹（基礎(chǔ)）和R繪制優(yōu)美的進(jìn)化樹（進(jìn)階），或者使用iPhylo網(wǎng)站來交互式繪圖：iPhylo 生成并繪制優(yōu)美的分類樹。

分步示例

獲取數(shù)據(jù)

此示例將使用以下兩個(gè)基因組：基因組A：GCF_003947435.1 [GTDB / NCBI]，基因組B：GCA_002011125.1 [GTDB / NCBI]

# 獲取基因組。
mkdir -p gtdbtk_test_steps/genomes
wget -q https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/003/947/435/GCF_003947435.1_ASM394743v1/GCF_003947435.1_ASM394743v1_genomic.fna.gz -O gtdbtk_test_steps/genomes/genome_a.fna.gz
wget -q https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/002/011/125/GCA_002011125.1_ASM201112v1/GCA_002011125.1_ASM201112v1_genomic.fna.gz -O gtdbtk_test_steps/genomes/genome_b.fna.gz

基因鑒定 (identify)

可以通過classify_wf單步完成，但是有時(shí)候分步流程更適合處理大數(shù)據(jù)（可以自行分配每一步的資源）：

ls -l gtdbtk_test_steps/genomes
# –extension，gz是基因組文件的后綴名
gtdbtk identify --genome_dir gtdbtk_test_steps/genomes --out_dir gtdbtk_test_steps/identify --extension gz --cpus 2

[2024-05-20 15:25:20] INFO: GTDB-Tk v2.3.2
[2024-05-20 15:25:20] INFO: gtdbtk identify --genome_dir gtdbtk_test_steps/genomes --out_dir gtdbtk_test_steps/identify --extension gz --cpus 2
[2024-05-20 15:25:20] INFO: Using GTDB-Tk reference data version r214: /data/home/jianglab/share/GTDB_r214/release214/
[2024-05-20 15:25:20] INFO: Identifying markers in 2 genomes with 2 threads.
[2024-05-20 15:25:20] TASK: Running Prodigal V2.6.3 to identify genes.
[2024-05-20 15:25:30] INFO: Completed 2 genomes in 9.31 seconds (4.65 seconds/genome).
[2024-05-20 15:25:30] TASK: Identifying TIGRFAM protein families.
[2024-05-20 15:25:38] INFO: Completed 2 genomes in 7.16 seconds (3.58 seconds/genome).
[2024-05-20 15:25:38] TASK: Identifying Pfam protein families.
[2024-05-20 15:25:38] INFO: Completed 2 genomes in 0.59 seconds (3.37 genomes/second).
[2024-05-20 15:25:38] INFO: Annotations done using HMMER 3.3.2 (Nov 2020).
[2024-05-20 15:25:38] TASK: Summarising identified marker genes.
[2024-05-20 15:25:38] INFO: Completed 2 genomes in 0.03 seconds (75.34 genomes/second).
[2024-05-20 15:25:38] INFO: Done.

耗時(shí)40s，內(nèi)存消耗131.76MB。

產(chǎn)生idetity/文件夾，關(guān)注摘要文件gtdbtk_test_steps/identify/identify/gtdbtk.ar53.markers_summary.tsv，其中詳細(xì)介紹了從古細(xì)菌 53 或細(xì)菌 120 標(biāo)記組中識別出的標(biāo)記。

還可以在每個(gè)基因組對應(yīng)的中間文件目錄下找到被調(diào)用的基因和標(biāo)記信息:

ls gtdbtk_test_steps/identify/identify/intermediate_results/marker_genes/genome_a.fna/

基因比對 (align)

比對步驟將對齊所有已識別的標(biāo)記，確定最可能的域并輸出串聯(lián)的MSA:

gtdbtk align --identify_dir gtdbtk_test_steps/identify --out_dir gtdbtk_test_steps/align --cpus 2

[2024-05-20 15:32:01] INFO: GTDB-Tk v2.3.2
[2024-05-20 15:32:01] INFO: gtdbtk align --identify_dir gtdbtk_test_steps/identify --out_dir gtdbtk_test_steps/align --cpus 2
[2024-05-20 15:32:01] INFO: Using GTDB-Tk reference data version r214: /data/home/jianglab/share/GTDB_r214/release214/
[2024-05-20 15:32:03] INFO: Aligning markers in 2 genomes with 2 CPUs.
[2024-05-20 15:32:04] INFO: Processing 2 genomes identified as archaeal.
[2024-05-20 15:32:06] INFO: Read concatenated alignment for 4,416 GTDB genomes.
[2024-05-20 15:32:06] TASK: Generating concatenated alignment for each marker.
[2024-05-20 15:32:06] INFO: Completed 2 genomes in 0.02 seconds (130.19 genomes/second).
[2024-05-20 15:32:06] TASK: Aligning 52 identified markers using hmmalign 3.3.2 (Nov 2020).
[2024-05-20 15:32:06] INFO: Completed 52 markers in 0.43 seconds (121.60 markers/second).
[2024-05-20 15:32:06] TASK: Masking columns of archaeal multiple sequence alignment using canonical mask.
[2024-05-20 15:32:10] INFO: Completed 4,418 sequences in 3.59 seconds (1,231.58 sequences/second).
[2024-05-20 15:32:10] INFO: Masked archaeal alignment from 13,540 to 10,135 AAs.
[2024-05-20 15:32:10] INFO: 0 archaeal user genomes have amino acids in <10.0% of columns in filtered MSA.
[2024-05-20 15:32:10] INFO: Creating concatenated alignment for 4,418 archaeal GTDB and user genomes.
[2024-05-20 15:32:12] INFO: Creating concatenated alignment for 2 archaeal user genomes.
[2024-05-20 15:32:12] INFO: Done.

耗時(shí)10s，內(nèi)存忽略不計(jì)（因?yàn)檫@里只有兩個(gè)genome，輸入的多了還是要考慮這步的內(nèi)存的）

產(chǎn)生align/文件夾，如果基因組識別出的標(biāo)記數(shù)量較少，它將被排除在這一步的分析之外。如果是這種情況，將會(huì)出現(xiàn)Warning。

根據(jù)基因組Domain的不同，將出現(xiàn) ar53 或 bac120 前綴文件。

• gtdbtk.ar53.msa.fasta.gz：其中包含用戶基因組和 GTDB基因組的 MSA。
• gtdbtk.ar53.user_msa.fasta.gz：僅包含用戶基因組

物種分類 (classify)

分類步驟會(huì)將基因組放入?yún)⒖紭渲?，然后確定其最可能的分類。

gtdbtk classify --genome_dir gtdbtk_test_steps/genomes --align_dir gtdbtk_test_steps/align --out_dir gtdbtk_test_steps/classify -x gz --cpus 2 --skip_ani_screen

[2024-05-20 15:42:48] INFO: GTDB-Tk v2.3.2
[2024-05-20 15:42:48] INFO: gtdbtk classify --genome_dir gtdbtk_test_steps/genomes --align_dir gtdbtk_test_steps/align --out_dir gtdbtk_test_steps/classify -x gz --cpus 2 --skip_ani_screen
[2024-05-20 15:42:48] INFO: Using GTDB-Tk reference data version r214: /data/home/jianglab/share/GTDB_r214/release214/
[2024-05-20 15:42:49] TASK: Placing 2 archaeal genomes into reference tree with pplacer using 2 CPUs (be patient).
[2024-05-20 15:42:49] INFO: pplacer version: v1.1.alpha19-0-g807f6f3
[2024-05-20 15:48:46] INFO: Calculating RED values based on reference tree.
[2024-05-20 15:48:48] TASK: Traversing tree to determine classification method.
[2024-05-20 15:48:48] INFO: Completed 2 genomes in 0.00 seconds (8,621.39 genomes/second).
[2024-05-20 15:48:48] TASK: Calculating average nucleotide identity using FastANI (v1.33).
[2024-05-20 15:48:52] INFO: Completed 6 comparisons in 2.69 seconds (2.23 comparisons/second).
[2024-05-20 15:48:52] INFO: 2 genome(s) have been classified using FastANI and pplacer.
[2024-05-20 15:48:52] INFO: Note that Tk classification mode is insufficient for publication of new taxonomic designations. New designations should be based on one or more de novo trees, an example of which can be produced by Tk in de novo mode.
[2024-05-20 15:48:52] INFO: Done.

耗時(shí)368s，內(nèi)存消耗45GB（因?yàn)檫@里只有兩個(gè)古菌genome，如果還有細(xì)菌的話應(yīng)該分配100GB左右）

產(chǎn)生align/文件夾，摘要文件gtdbtk.ar53.summary.tsv和包含這些基因組的參考樹gtdbtk.ar53.classify.tree。

注意事項(xiàng)

• 為什么 GTDB-Tk 和 NCBI 或 Silva 分類名稱之間的命名系統(tǒng)存在差異？

GTDB-Tk 使用 GTDB 分類法 (https://gtdb.ecogenomic.org/)。此分類法與NCBI 和 Silva 類似，但不完全相同。
在許多情況下，GTDB分類法更嚴(yán)格地遵循等級后綴的命名規(guī)則，比如GTDB有 Nitrospirota 而沒有Nitrospirae。

• 可以將細(xì)菌樹和古細(xì)菌樹合并成一棵樹嗎？

細(xì)菌樹和古細(xì)菌樹是從不同的標(biāo)記基因推斷出來的。目前，這些樹的正確生根仍然是一個(gè)開放的研究領(lǐng)域。
GTDB-Tk 不提供合并樹的工具，但可以通過操作 Newick文件來人為地組合它們。

• GTDB-Tk 如何使用平均核苷酸同一性（ANI）驗(yàn)證物種分配？

GTDB-Tk 使用 skani（在 v2.3.2 之前一直使用 fastANI）來估計(jì)基因組之間的ANI。
僅當(dāng)基因組之間的 ANI 在物種 ANI 限制半徑內(nèi)（通常為 95%）且比對分?jǐn)?shù)(AF)>=0.5時(shí)，查詢基因組才被分類為與參考基因組屬于同一物種。
在某些情況下，查詢基因組的系統(tǒng)發(fā)育位置可能不支持物種分配。

• 互斥選項(xiàng) --mash_db 和 --skip_ani_screen 之間有什么區(qū)別？

從 GTDB-Tk v2.2+ 開始，classify_wf和分類函數(shù)需要額外的參數(shù)才能運(yùn)行：--mash_db 或 --skip_ani_screen。

使用新版本的 GTDB-Tk v2.4.0，分類管道（classify_wf 和 classify）的第一階段是將所有用戶基因組與所有參考基因組進(jìn)行比較，并在可能的情況下根據(jù)ANI 匹配對它們進(jìn)行注釋。
使用 --mash_db 選項(xiàng)將向 GTDB-Tk 指示 ANI篩選所需的草圖 Mash 數(shù)據(jù)庫的路徑。
如果沒有可用的數(shù)據(jù)庫（即這是第一次運(yùn)行分類），則 --mash_db選項(xiàng)將繪制一個(gè)可用于后續(xù)調(diào)用的新 Mash 數(shù)據(jù)庫。 --skip_ani_screen 選項(xiàng)將跳過預(yù)篩選步驟并對所有基因組進(jìn)行分類，類似于之前版本的 GTDB-Tk。