Shiny从入门到入定——16-Escaping the graph

GTDB-Tk 是一个软件工具包，用于根据基因组数据库分类 GTDBGenome Database Taxonomy GTDB 为细菌和古菌基因组分配客观分类学分类assigning objective taxonomic classifications。可以用于宏基因组组装基因组 (MAG)、单菌基因组和单细胞基因组。

安装GTDB-Tk

• 通过conda安装主程序

  # latest version
  conda create -n gtdbtk -c conda-forge -c bioconda gtdbtk
  
  # specific version (replace 2.3.0 with the version you wish to install, recommended)
  conda create -n gtdbtk -c conda-forge -c bioconda gtdbtk=2.3.0

• 下载GTDB-Tk数据库 (~47 G)

  最新版本的数据库为R214，需要 >=2.1.0 版本的GTDB-Tk。

  • 自动 (非常慢)

    download-db.sh

  • 手动

     	# 找到下载脚本的路径并打开，从中找到数据库的链接
    which download-db.sh
      
    # 进入数据库默认路径
    cd ~/miniconda3/envs/gtdbtk/share/gtdbtk-2.3.0/db/
    	  
       # 用wget下载 (非常慢)，推荐在Windows下用其他软件下载，再上传服务器
    nohup wget -c https://data.gtdb.ecogenomic.org/releases/release214/214.0/auxillary_files/gtdbtk_r214_data.tar.gz --no-check-certificate &
    
    tar zxvf gtdbtk_r214_data.tar.gz
    
    mv release214/* ./
    
    # 也可自定义数据库的存放位置
    编辑`~/miniconda3/envs/gtdbtk/etc/conda/activate.d/gtdbtk.sh`，将`export GTDBTK_DATA_PATH=`后面的内容改成数据库所在的路径。

运行软件

:::info
注意：pplacer requires ~215 GB of RAM to fully load the bacterial tree into memory. 因此要么在服务器上运行，要么设置足够大的swap分区方可在PC中运行。
:::

classify_wf

:::default
可以通过classify_wf命令完成整个工作流程。
:::

# v2.3.0增加了两个参数（二选一）：--skip_ani_screen，--mash_db，前者无参数值，后者指定MASH参考数据库 (.msh)的输出或读取路径
gtdbtk classify_wf --genome_dir <my_genomes> --out_dir <output_dir> --extension fna --cpus 20 --force --skip_ani_screen

+++primary 参数解析

• --genome_dir
  directory containing genome files in FASTA format

• --batchfile
  path to file describing genomes - tab separated in 2 or 3 columns (FASTA file, genome ID, translation table [optional])

• --out_dir
  directory to output files

• --skip_ani_screen Skip the ani_screening step to classify genomes using mash and FastANI (default: False)

• --mash_db MASH_DB path to save/read (if exists) the Mash reference sketch database (.msh)

• -x, --extension
  extension of files to process, gz = gzipped (Default: fna)

• --min_perc_aa
  exclude genomes that do not have at least this percentage of AA in the MSA (inclusive bound) (Default: 10)

• --prefix
  prefix for all output files (Default: gtdbtk)

• --cpus
  number of CPUs to use (Default: 1)

• --pplacer_cpus
  use pplacer_cpus during placement (default: cpus)

• --force
  continue processing if an error occurs on a single genome

• --scratch_dir
  Reduce pplacer memory usage by writing to disk (slower).

• --min_af
  minimum alignment fraction to consider closest genome (Default: 0.65)

• --tmpdir
  specify alternative directory for temporary files （Default: /tmp

• --debug
  create intermediate files for debugging purposes
  +++

:::default
但下面我们演示分步运行。在处理大型管道时，单独运行这些步骤有时会很有用。
:::

Step 1 准备输入文件

以下两个基因组作为示例文件：* Genome A: GCF_003947435.1 [GTDB / NCBI] * Genome B: GCA_002011125.1 [GTDB / NCBI]。

# Create the directory.
mkdir -p /tmp/gtdbtk && cd /tmp/gtdbtk

# Obtain the genomes.
mkdir -p /tmp/gtdbtk/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 /tmp/gtdbtk/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 /tmp/gtdbtk/genomes/genome_b.fna.gz

Step 2 Gene calling (identify)

gtdbtk identify --genome_dir /tmp/gtdbtk/genomes --out_dir /tmp/gtdbtk/identify --extension gz --cpus 20

+++success Results

获得的基因和标记marker信息可以在每个基因组反应中间文件目录genomes respeective intermediate files directory下找到，如下所示。

ls /tmp/gtdbtk/identify/identify/intermediate_results/marker_genes/genome_a.fna/

genome_a.fna_pfam_tophit.tsv
genome_a.fna_protein.gff.sha256
genome_a.fna_pfam_tophit.tsv.sha256
genome_a.fna_tigrfam.out
genome_a.fna_pfam.tsv
genome_a.fna_tigrfam.out.sha256
genome_a.fna_pfam.tsv.sha256
genome_a.fna_tigrfam_tophit.tsv
genome_a.fna_protein.faa
genome_a.fna_tigrfam_tophit.tsv.sha256
genome_a.fna_protein.faa.sha256
genome_a.fna_tigrfam.tsv
genome_a.fna_protein.fna
genome_a.fna_tigrfam.tsv.sha256
genome_a.fna_protein.fna.sha256
prodigal_translation_table.tsv
genome_a.fna_protein.gff
prodigal_translation_table.tsv.sha256

:::default
大部分情况下只需要读取摘要文件，其中详细说明了从archaeal 122或bacterial 120标记集marker set识别的标记。
:::

cat /tmp/gtdbtk/identify/gtdbtk.ar122.markers_summary.tsv

name    number_unique_genes     number_multiple_genes   number_multiple_unique_genes    number_missing_genes    list_unique_genes       list_multiple_genes     list_multiple_unique_genes      list_missing_genes
genome_a.fna    109     3       0       10      PF00368.13,PF00410.14,PF00466.15,PF00687.16,PF00827.12,PF00900.15,PF01000.21,PF01015.13,PF01090.14,PF01092.14,PF01157.13,PF01191.14,PF01194.12,PF01198.14,PF01200.13,PF01269.12,PF01280.15,PF01282.14,PF01655.13,PF01798.13,PF01864.12,PF01868.11,PF01984.15,PF01990.12,PF02006.11,PF02978.14,PF03874.11,PF04019.7,PF04104.9,PF04919.7,PF07541.7,PF13656.1,PF13685.1,TIGR00021,TIGR00037,TIGR00042,TIGR00111,TIGR00134,TIGR00240,TIGR00264,TIGR00270,TIGR00279,TIGR00283,TIGR00291,TIGR00293,TIGR00307,TIGR00308,TIGR00323,TIGR00324,TIGR00335,TIGR00336,TIGR00337,TIGR00389,TIGR00392,TIGR00398,TIGR00405,TIGR00408,TIGR00422,TIGR00425,TIGR00432,TIGR00442,TIGR00448,TIGR00456,TIGR00463,TIGR00468,TIGR00471,TIGR00490,TIGR00491,TIGR00501,TIGR00521,TIGR00549,TIGR00670,TIGR00729,TIGR00936,TIGR00982,TIGR01008,TIGR01012,TIGR01018,TIGR01020,TIGR01025,TIGR01028,TIGR01038,TIGR01046,TIGR01052,TIGR01060,TIGR01077,TIGR01080,TIGR01213,TIGR01309,TIGR01952,TIGR02076,TIGR02153,TIGR02236,TIGR02258,TIGR02390,TIGR02651,TIGR03626,TIGR03627,TIGR03628,TIGR03636,TIGR03653,TIGR03665,TIGR03671,TIGR03672,TIGR03674,TIGR03677,TIGR03680,TIGR03683,TIGR03684    PF01496.14,TIGR00458,TIGR00658          PF01866.12,TIGR00064,TIGR00373,TIGR00522,TIGR02338,TIGR02389,TIGR03629,TIGR03670,TIGR03673,TIGR03722
genome_b.fna    118     2       0       2       PF00368.13,PF00410.14,PF00466.15,PF00687.16,PF00827.12,PF00900.15,PF01000.21,PF01015.13,PF01090.14,PF01092.14,PF01157.13,PF01191.14,PF01194.12,PF01198.14,PF01200.13,PF01269.12,PF01280.15,PF01282.14,PF01655.13,PF01798.13,PF01864.12,PF01866.12,PF01868.11,PF01984.15,PF01990.12,PF02006.11,PF02978.14,PF03874.11,PF04019.7,PF04104.9,PF04919.7,PF07541.7,PF13656.1,TIGR00021,TIGR00037,TIGR00042,TIGR00064,TIGR00111,TIGR00134,TIGR00240,TIGR00264,TIGR00270,TIGR00279,TIGR00283,TIGR00291,TIGR00293,TIGR00307,TIGR00308,TIGR00323,TIGR00324,TIGR00335,TIGR00336,TIGR00337,TIGR00373,TIGR00389,TIGR00392,TIGR00398,TIGR00405,TIGR00408,TIGR00422,TIGR00425,TIGR00432,TIGR00442,TIGR00448,TIGR00456,TIGR00458,TIGR00463,TIGR00468,TIGR00471,TIGR00490,TIGR00491,TIGR00501,TIGR00521,TIGR00522,TIGR00549,TIGR00658,TIGR00670,TIGR00729,TIGR00936,TIGR00982,TIGR01008,TIGR01012,TIGR01018,TIGR01020,TIGR01025,TIGR01028,TIGR01038,TIGR01046,TIGR01052,TIGR01060,TIGR01077,TIGR01080,TIGR01309,TIGR01952,TIGR02076,TIGR02153,TIGR02236,TIGR02258,TIGR02338,TIGR02389,TIGR02390,TIGR02651,TIGR03626,TIGR03628,TIGR03629,TIGR03636,TIGR03653,TIGR03665,TIGR03670,TIGR03671,TIGR03672,TIGR03673,TIGR03674,TIGR03677,TIGR03680,TIGR03683,TIGR03684,TIGR03722 PF01496.14,PF13685.1            TIGR01213,TIGR03627

+++

Step 3 Aligning genomes (align)

该步骤将对齐所有已识别的标记，确定最可能的域domain 并输出串联的MSA。

gtdbtk align --identify_dir /tmp/gtdbtk/identify --out_dir /tmp/gtdbtk/align --cpus 20

+++success Results

:::warning
要注意输出，如果一个基因组识别出的标记数量较少，那么在这一步，它将被排除在分析之外。如果出现这种情况，将出现警告。
:::

根据domain的不同，ar122或bac120的前缀文件将出现，其中包含用户基因组和GTDB基因组的MSA，或仅包含用户基因组（gtdbtk.ar122.MSA.fasta和gtdbtk.ar122.user_MSA.fasta）。

ls /tmp/gtdbtk/align

align/
gtdbtk.ar122.user_msa.fasta
identify/
gtdbtk.ar122.filtered.tsv
gtdbtk.log
gtdbtk.ar122.msa.fasta
gtdbtk.warnings.log

+++

Step 4 基因组分类 (classify)

将用户的基因组放置于参考树上，然后决定其最为可能的分类。

gtdbtk classify --genome_dir /tmp/gtdbtk/genomes --align_dir /tmp/gtdbtk/align --out_dir /tmp/gtdbtk/classify -x gz --cpus 20

+++success Results

两个主要的输出文件包括总结文件summary file和包含基因组的参考树 (gtdbtk.ar122.summary.tsv和gtdbtk.ar122.classify.tree)。基因组的分类信息呈现于总结文件中。

ls /tmp/gtdbtk/classify

classify/
gtdbtk.ar122.summary.tsv
gtdbtk.warnings.log
gtdbtk.ar122.classify.tree
gtdbtk.log

+++

错误信息处理

• OpenBLAS blas_thread_init: pthread_create failed for thread 109 of 128: Resource temporarily unavailable

  • 解决方案：少调用一些线程就ok了。

参考

• GTBD-Tk Documentation