使用DIAMOND将全基因组蛋白序列比对到Nr数据库

1. DIAMOND简介

对全基因组的基因进行Nr注释是必不可少的一步。由于Nr数据库非常大,导致使用BLAST会消耗巨大的计算资源和时间。使用DIAMOND则能快500-20000倍,而获得和BLAST比较一致的结果。特别是对于长度为100-150bp,数量超过1M的核酸序列,DIAMOND的速度比BLASTX快20000倍;当e-value设置为1e-5时,DIAMOND能报告80~90%的BLASTX匹配结果;若设置sensitivie模式,DIAMOND速度是BLASTX的2500倍,DIAMOND能报告94%的BLASTX匹配结果。因此,对大批量的序列进行Nr注释或规模较大的蛋白数据库的比对,使用DIAMOND是优先选择的方式。

虽然DIAMOND软件的运行速度快,但是其结果中没有匹配的Subject序列的描述信息。若需要蛋白名称和物种的描述信息,需要自行编写软件从数据库的FASTA文件中提取并整合到DIAMOND结果中。

2. DIAMOND下载和安装

$ wget https://github.com/bbuchfink/diamond/releases/download/v0.9.24/diamond-linux64.tar.gz -P ~/software/
$ mkdir /opt/biosoft/diamond
$ tar zxf ~/software/diamond-linux64.tar.gz -C /opt/biosoft/diamond
$ echo 'PATH=$PATH:/opt/biosoft/diamond' >> ~/.bashrc
$ source ~/.bashrc

3. DIAMOND使用

DIAMOND软件的主命令是diamond,它的使用包含几个子命令。DIAMOND最常用的使用方法:

使用DIAMOND软件的子命令makedb将FASTA格式的蛋白序列创建成后缀为dmnd的数据库文件:
$ diamond makedb --in nr_eukaryon.fasta -d nr_eukaryon_20190405
    将Nr数据库(版本20190405)的真核子集创建成DIAMOND数据库,数据库包含37.4M条蛋白序列,服务器是24线程2.0GHz的CPU资源,耗时17min。

将Illumina测序数据比对到数据库:
$ diamond blastx --db nr_eukaryon_20190405 --query reads.fq.gz --out reads.tab

将全基因组蛋白序列比对到数据库:
$ diamond blastp --db nr_eukaryon_20190405 --query proteins.fasta --out nr.tab --outfmt 6 --sensitive --max-target-seqs 20 --evalue 1e-5 --id 30 --block-size 20.0 --tmpdir /dev/shm --index-chunks 1
    对18700个蛋白序列进行Nr注释,服务器160线程2.4GHz的CPU资源,耗时91min(不加--sensitive参数只需要16min)。

此外,DIAMOND还包含其它几个子命令:

blastx
    将核酸序列比对到蛋白序列数据库
view
    从DAA文件生成格式化的输出结果
getseq
    从DIAMOND数据库文件获取FASTA序列
version
    打印DIAMOND软件的版本号
dbinfo
    打印数据库的信息
help
    打印DIAMOND软件的详细帮助信息

DIAMOND的参数:

数据库参数:
--in <string>    default: STDIN
    输入FASTA格式的蛋白序列数据库文件。
--db | -d <string>
    设置数据库文件路径和前缀。创建数据库时,会生成一个后缀为.dmnd的数据库文件。比对时,则是输入相应的数据库文件。
--taxonmap <string>
--taxonnodes <string>
    这两个参数输入NCBI的物种分类信息。前者输入Nr数据库中蛋白编号和物种编号的对应关系信息( ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/prot.accession2taxid.gz);后者输入物种编号的层次信息(ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdmp.zip)。输入这两个文件有利于DIAMOND和Nr数据库进行比对时,可以选择和数据库中指定的物种大类的数据进行比对。

输入参数:
--db | -d <string>
    设置数据库文件路径和前缀。创建数据库时,会生成一个后缀为.dmnd的数据库文件。比对时,则是输入相应的数据库文件。
--query | -q <string>    default: STDIN
    输入需要注释的FASTA或FASTQ格式的序列文件。可以是带.gz后缀的压缩文件。
--taxonlist <string>
    输入NCBI分类编号,仅对数据库中的目标子集进行比对。可以输入多个使用逗号分隔的编号ID。例如,古菌(2157)、细菌(2)、病毒(10239)、真菌(4751)、动物(33208)和植物(3193)。使用该参数,必须要使用--taxonmap和--taxonnodes参数构建数据库。
--query-gencode <int>    default:1
    设置在进行BLASTX模式进行比对时所使用的遗传密码。
--strand <string>    default: both
    设置对query序列的那条链进行比对。可以设置的值为:boh、plus和minus。
--min-orf | -l <int>
    在使用BLASTX模式进行比对时,若序列的某个ORF低于此值,则忽略该ORF。默认设置下:若核酸序列长度低于30,则值为1;若核酸序列长度低于100,则值为20;若核酸序列长度不低于100,则值为40。若该值设为1,则表示在BLASTX模式时,识别所有长度的ORF并进行比对。

比对参数:
--sensitive
    添加该参数,则能得到更多比对结果。该模式适合比对较长的序列。默认模式主要适用于比对short reads序列(Illumina reads),搜寻比对长度为30~40aa且bit得分大于50的匹配结果。
--more-sensitive
    相比于sensitive,能得到更全的比对结果。
--frameshifit | -F <int>
    在进行BLASTX模式时,设置对移码匹配的罚分,且一般推荐设置该参数值为15左右。设置该参数后,有利于比对较长的有INDEL错误的reads序列。默认情况是是禁用了此功能。
--gapopen <int>    default: BLOSUM62矩阵的默认值为1
    设置打开Gap罚分。
--gapextend <int>    default: BLOSUM62矩阵的默认值为1
    设置延长Gap罚分。
--matrix <string>    default: BLOSUM62
    设置计分矩阵。可以设置的值有:BLOSUM45,BLOSUM50,BLOSUM62,BLOSUM80,BLOSUM90,PAM250,PAM70和PAM30。
--comp-based-stats <bool>    default: 1
    设置是否对比对得分进行校正。
--algo <bool>
    设置种子序列搜索的算法。0表示使用double-indexed算法,1表示query-indexed算法。第一种算法是软件的主要算法,但不适合query数据较少的情况,这时需要使用第二种算法。默认设置下,软件会根据query和数据库的数据量自动选择相应的算法。两种算法的结果基本一致。
--freq-sd <int>
    设置种子序列最大的匹配频率。若种子序列在query或数据库种出现的频率过高,则忽略之。在一般模式下,该参数值为50;在sensitive模式下,该参数值为100,在more-sensitive模式下,该参数值为200。


输出参数:
--out | -o <string>    default:STDOUT
    设置输出文件。
--outfmt | -f <int>    default: 6
    设置输出格式。支持的格式有:0,两两比对格式;5,XML格式;6,BLAST表格格式;100,DIAMOND匹配存档(DDA)格式,该格式可以使用diamond view命令转换成其它格式;101,SAM格式;102,分类鉴定结果,该模式下结果文件分三列,QueryID、NCBI物种分类ID、最佳匹配evaule;103,PAF格式。
--salltitles
    在DDA格式结果中包含all subject titles。
--sallseqid
    在DDA格式结果中包含all subject ids。由于id比title要短,相比于--salltitles,本参数输出的结果文件更小。
--compress <bool>    default: 0
    设置是否对暑促文件进行压缩。0表示不压缩,1表示进行gzip压缩。
--max-target-seqs | -k <int>    default: 25
    设置每个query比对结果的最大匹配数目。
--top <int>
    若设置该参数的值,则报告一定数目的top匹配结果。例如设置该参数值为10,则报告匹配得分为top10的结果,即报告所有得分和最高得分相差10%以内的匹配结果。该参数会取代--max-target-seqs参数设置的值。
--range-culling
    添加该参数能剔除匹配到query序列相同位置的hit结果。该参数适合于较长的query DNA序列。当DNA序列较长时,它可能跨越多个基因,而一般情况下仅报告25个最佳匹配结果,若这25个hits结果可能evaule值很好但得分较低,且一起不能覆盖query序列的50%,会导致阵势的比对结果较差。因此,添加该参数能剔除匹配到query序列相同位置的hit结果,以使最终的比对结果对query序列的覆盖率超过50%。
--range-cover <float>    default: 50.0
    和参数--range-culling配合使用,表示剔除hits以使所有的比对结果能对query序列的覆盖率超过此参数设定的百分比。
--long-reads
    对于long read DNA测序数据,推荐添加该参数,等同参数--range-culling --top 10 -F 15。
--evalue | -e <float>    default: 0.001
    设置比对的evalue阈值。
--min-score <int>
    设置比对的bit score阈值。设置该参数会替代--evalue参数。
--id <int>
    设置identity阈值。
--query-cover <int>
    设置对query序列的覆盖度阈值。该阈值是HSP的阈值。
--subject-cover <int>
    设置对subject序列的覆盖度阈值。该阈值是HSP的阈值。
--max-hsps <int>
    设置报告的HSPs数量。默认是报告所有的HSPs,若设置成1,则仅报告最优HSP。
--unal <bool>
    设置是否报告没有匹配的比对结果。默认设置下,--outfmt为1、5和101时,该参数值为1,其它输出方式该参数值为0。
--no-self-hits
    添加该参数,则不报告序列和自己的比对结果。

性能参数:
--threads | -p <int>    default: Max
    设置程序运行所使用的CPU线程数。默认是服务器可用的最大CPU线程数。
--block-size | -b <float>    default: 2.0
    设置每次处理多少G的序列字符数。该参数控制程序消耗的内存量,一般内存消耗该值的6倍。设置更大的值,会消耗更多的内存和临时磁盘空间,但能提高性能。默认设置下,程序一次处理2G个序列字符,消耗内存12G。
--tmpdir | -t <string>    default: directory of --out
    设置临时文件夹路径。推荐该文件所在剩余磁盘空间有100G及以上。若将该参数设置为/dev/shm,则会将临时文件存放在内存,会增加内存消耗和计算性能。
--index-chunks | -c <int>    default: 4
    将seed index分成指定的份数。推荐将该参数值设置成1,能增加计算性能和内存使用量。

DIAMOND软件运行过程中,若使用参数 –block-size 20.0 –tmpdir /dev/shm –index-chunks 1 能使用较大内存,不会对数据进行分块,程序运行较快。程序运行过程中可能完全不会生成过程文件。但在程序运行过程中给出一些进度信息。当shape到15的时候,则表示程序快运行完毕。通过查看shape值来确定程序运行的进度。

CPU threads: 160
Scoring parameters: (Matrix=BLOSUM62 Lambda=0.267 K=0.041 Penalties=11/1)
Temporary directory: /dev/shm
Opening the database…
Target sequences to report alignments for: 20
Opening the input file…
Opening the output file…
Loading query sequences…
Masking queries…
Building query seed set…
Algorithm: Double-indexed
Building query histograms…
Allocating buffers…
Loading reference sequences…
Building reference histograms…
Allocating buffers…
Initializing temporary storage…
Processing query chunk 0, reference chunk 0, shape 0, index chunk 0.
Building reference index…
Building query index…
Building seed filter…
Searching alignments…
Processing query chunk 0, reference chunk 0, shape 1, index chunk 0.
Building reference index…
Building query index…
Building seed filter…
Searching alignments…
...
Processing query chunk 0, reference chunk 0, shape 15, index chunk 0.
Building reference index…
Building query index…
Building seed filter…
Searching alignments…
Deallocating buffers…
Computing alignments…
Deallocating reference…
Loading reference sequences…
Deallocating buffers…
Deallocating queries…
Loading query sequences…
Closing the input file…
Closing the output file…
Closing the database file…
Deallocating taxonomy…
Total time = 659.515s
Reported 33888205 pairwise alignments, 33890340 HSPs.
122518 queries aligned.

DIAMOND默认设置下输出表格格式的结果。结果分12列,其结果信息和BLAST默认设置-outfmt 6输出的格式完全一致。

1. query序列ID
2. subject序列ID
3. Identity百分比
4. 匹配长度
5. 错配长度
6. 打开Gap的次数
7. query序列起始位点
8. query序列结束位点
9. subject序列起始位点
10. subject序列结束位点
11. E-vaule值
12. bitscore得分

发表评论

电子邮件地址不会被公开。 必填项已用*标注

此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据