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

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 1.7.0 with the version you wish to install, recommended)
    conda create -n gtdbtk -c conda-forge -c bioconda gtdbtk=1.7.0
  • 下载 GTDB-Tk 数据库 (~47 G)

    最新版本的数据库为 R202,需要 >=1.5.0 版本的 GTDB-Tk。

    • 自动 (非常慢)

      download-db.sh
    • 手动

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

运行软件

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

classify_wf

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

gtdbtk classify_wf --genome_dir <my_genomes> --out_dir <output_dir> --extension fna --cpus 20 --force
参数解析
  • --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

  • -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

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

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
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

大部分情况下只需要读取摘要文件,其中详细说明了从 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
Results

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

根据 domain 的不同, ar122bac120 的前缀文件将出现,其中包含用户基因组和 GTDB 基因组的 MSA,或仅包含用户基因组( gtdbtk.ar122.MSA.fastagtdbtk.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
Results

两个主要的输出文件包括总结文件summary file和包含基因组的参考树 ( gtdbtk.ar122.summary.tsvgtdbtk.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 了。

参考

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