【原】lncRNA组装流程的软件介绍之diamond

健明 2021-07-14

展开全文

咱们《生信技能树》的B站有一个lncRNA数据分析实战，缺乏配套笔记，所以我们安排了100个lncRNA组装案例文献分享，以及这个流程会用到的100个软件的实战笔记教程！

下面是100个lncRNA组装流程的软件的笔记教程

DIAMOND是一种高通量比对程序，可将DNA测序reads文件与蛋白质参考序列文件（如NCBI-nr）进行比较。它以C ++实现，旨在多核服务器上快速运行。

该程序明确地设计为，利用具有大内存容量和许多内核的现代计算机体系结构。那么为什么它那么快呢，因为它使用了种子和延伸方法。额外的算法成分是使用缩小的字母，间隔种子和双索引。算法简单了解一下就可以了，具体的算法的内容比较难懂就不深入讨论了。

其主要的特点包括(为什么该工具那么受欢迎):

在比较NCBI-nr数据库的短DNA reads时，DIAMOND比BLASTX 快4个数量级，并且在e-value < 0.001 的比对上保持相当的灵敏度水平。简单一句话就是又快又准。
可用于移位分析的长reads比对
资源要求低，适合在标准台式机或笔记本电脑上运行(入门门槛低，适合各类的玩家)
各种输出格式，包括BLAST对比格式，例如格式6的tabular分隔形式和格式5的XML格式 (很多下游分析，或者脚本都是基于BLAST的输出结果，diamond能够直接输出和Blast相同的格式不能不说是最大的优点之一)

一、软件安装

使用conda安装

conda install diamond

二、diamond的用法

安装完成以后，可以使用diamond --help来查看软件的帮助文档。

1. 软件用法：

常用参数：

 输入参数：
--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  该模式比默认模式要灵敏，一些碎的弱比对结果也会输出。通常建议用于长序列比对。而默认模式主要用于短读比对，如30~40个氨基酸长度，得分值>50的。
--more-sensitive 相比--sensitive,会输出更多的弱比对结果。
--seg (yes, no) 屏蔽query序列中的低复杂序列片段（如SSR序列），若query为蛋白序列，默认no, 若query为CDS序列，默认yes。
--max-target-seqs 输出query比对到database中目标序列的条数(默认值为25)。设置为0将输出所有的比对结果。一般设为10就够用了。
--top  与--max-target-seqs作用相同，控制输出结果数目，但top按最优百分比，如某query能比对100条目标序列 –top 20输出最优的20条。注意若设置该参数会导致 --max-target-seqs失效，故建议二者选一。
--evalue/-e  比对得分期望的E 值，默认0.001。
--min-score  设置最小得分值，注意若设置该参数会导致--evalue失效，建议二者选一。
--id  设置identity, 只输出>identity的比对结果。
--query-cover  设置query比对长度覆盖阈值，只输出高于该覆盖率的结果。
 
 内存和性能参数
--block-size/-b  程序运行每次处理的数据量，以Gb为单位。是控制程序内存和存储的重要参数，默认为2。内存大概用到12Gb,是其6倍。大家要结合服务器内存大小来设置该参数
--tmpdir  程序会产生很多隐藏的垃圾文件（.PS_*），默认输出到当前目录下，建议用临时文件夹存放，方便删除。
 
 输出格式参数
--outfmt  与BLAST类型提供多种输出格式 
               0  = BLAST pairwise
               5  = BLAST XML6
               6  = BLAST tabular
              100 = DIAMOND alignment archive (DAA)
              101 = SAM
这里详细介绍格式6，即tabular格式。与BLAST的m8格式完全一致，可能因为BLAST太常用了，所以整成一样，好方便下游的其他分析。但DIAMOND输出结果信息中，可以提供很多额外的信息，可以根据需求自由组合。我的建议是在m8格式的基础上添加qlen和qcovhsp两列信息，可在结果中直接查看query的覆盖率，有助于判断比对结果

三、软件运行命令

mkdir diamond
# 构建数据库索引
nohup diamond makedb --in ./nr -d nr  &
# 比对
nohup diamond blastx -e 1e-5 \
-d ~/database/blastDB/nr/diamond/nr \
-q ~/lncrna/test/filter4_by_pfam_exon.fa \
-f 6 -o ./dna_matches.txt &

命令参数解读：

-e 1e-5 # 比对得分期望的E 值为0.00005
-d ~/database/blastDB/nr/diamond/nr # nr数据库
-q ~/lncrna/test/filter4_by_pfam_exon.fa # 待比对fasta序列
-f 6 # 输出格式6
-o ./dna_matches.txt & # 输出文件

四、结果解读

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

DIAMOND提供输出信息项如下，可根据你的需求进行自由组合。

qseqid Query Seqid    # query的ID名
qlenQuery sequence length    # query序列长度
sseqid Subject Seqid        # 比对上的序列ID
sallseqid All subject Seqid(s)         # 所有比对上的序列ID
slen Subject sequence length          # 目标序列长度
qstart Start of alignment in query         # query起始位置
qend End of alignment in query            # query终止位置
sstart Start of alignment in subject       # 目标序列起始位置
send End of alignment in subject          # 目标序列终止位置
qseq Aligned part of query sequence           # query中部分比对上的序列
sseq Aligned part of subjectsequence   # 目标中部分比对上的序列    
evalue Expect value   # 比对E值
bitscore Bit score             # 比对得分值
score Raw score 
length Alignment length  # 比对长度
pident Percentage of identicalmatches      # 比对结果相同hits的比例
nident Number of identical matches    #  结果相同hits数量
mismatch Number of mismatches        # 错配数
positive Number of positive - scoringmatches   # 得分值大于0的比对片段数量
gapopen Number of gap openings        # 比对结果有几处gap
gaps Total number of gaps      # 所有gap处的氨基酸数总和
ppos Percentage of positive -scoring matches   # 得分值大于0的片段数比例
qframe Query frame        # query的移码
btop Blast tracebackoperations(BTOP)      # Blast中的回溯操作
stitle Subject Title      # 目标序列标题信息
salltitles All Subject Title(s)           # 所有目标序列的标题信息
qcovhsp Query Coverage Per HSP              #query比对长度占序列总长百分比（经实测：qcovhsp=(qend-qstart+1)/qlen，并非是 length/qlen ,但二者结果相差不大）
qtitle Query title        # query的标题信息

五、常见问题：

1. Diamond适合并行运行多个蛋白质fasta的比对吗？

建议不要同时运行多个DIAMOND的任务在同一台机器上，因为如果将更多的资源分配给单个任务，效率其实会更高。

2. Diamond比对好像不够准确，它没有找到我期待中的比对结果。

有许多原因导致无法找到期待中的比对结果。原因如下:

- 比对的hints可能没有满足最低的e-value的要求，即使他们是100%match的。
- 低复杂度掩蔽和组成偏差校正(Low complexity masking and compositional bias correction)可能导致hints被过滤。
- 由于算法的特点，diamond无法比对低于10AA蛋白质序列。另外，有一些高重复的input，会被过滤掉。