利用GTDB-TK对细菌和古菌基因组进行物种分类
发表于:2021-12-10 | 分类: 生物信息
字数统计: 1.5k | 阅读时长: 7分钟 | 阅读量:

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

安装GTDB-Tk

  • 通过conda安装主程序

    1
    2
    3
    4
    5
    # 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。

    • 自动 (非常慢)

      1
      download-db.sh
    • 手动

      1
      2
      3
      4
      5
      6
      7
      8
      9
      10
      11
      12
      13
      14
      15
       	# 找到下载脚本的路径并打开,从中找到数据库的链接
      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命令完成整个工作流程。
:::

1
2
# 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]。

1
2
3
4
5
6
7
8
# 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)

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

+++success Results

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

1
ls /tmp/gtdbtk/identify/identify/intermediate_results/marker_genes/genome_a.fna/
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
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识别的标记。
:::

1
cat /tmp/gtdbtk/identify/gtdbtk.ar122.markers_summary.tsv
1
2
3
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。

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

+++success Results

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

根据domain的不同,ar122bac120的前缀文件将出现,其中包含用户基因组和GTDB基因组的MSA,或仅包含用户基因组(gtdbtk.ar122.MSA.fastagtdbtk.ar122.user_MSA.fasta)。

1
ls /tmp/gtdbtk/align
1
2
3
4
5
6
7
align/
gtdbtk.ar122.user_msa.fasta
identify/
gtdbtk.ar122.filtered.tsv
gtdbtk.log
gtdbtk.ar122.msa.fasta
gtdbtk.warnings.log

+++

Step 4 基因组分类 (classify)

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

1
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.tsvgtdbtk.ar122.classify.tree)。基因组的分类信息呈现于总结文件中。

1
ls /tmp/gtdbtk/classify
1
2
3
4
5
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了。

参考

上一篇:
计算蛋白质等电点并绘制全局pI图
下一篇:
宏病毒组分析流程1-VirSorter2