DIAMOND 基因序列快速比对工具使用及超算集群并行计算指南

4. 结果解读

比对结束后，可以查看输出文件中的比对结果。DIAMOND 的输出文件包含每个待比对序列的匹配结果，包括匹配的参考序列名、匹配位置、匹配得分等信息。

结果字段表示与原生 BLAST 结果一致。

5. 帮助说明

以上就是 DIAMOND 的基本使用方法，更详细的说明可以参考官方文档：https://github.com/bbuchfink/diamond/wiki。

# downloading the tool，下载工具
wget http://github.com/bbuchfink/diamond/releases/download/v2.1.8/diamond-linux64.tar.gz
tar xzf diamond-linux64.tar.gz
# creating a diamond-formatted database file 创建 diamond 数据库
./diamond makedb --in reference.fasta -d reference
# running a search in blastp mode 使用 blastp 模式比对序列
./diamond blastp -d reference -q queries.fasta -o matches.tsv
# running a search in blastx mode 使用 blastx 模式比对序列
./diamond blastx -d reference -q reads.fasta -o matches.tsv
# downloading and using a BLAST database update_blastdb.pl --decompress --blastdb_version 5 swissprot
./diamond prepdb -d swissprot
./diamond blastp -d swissprot -q queries.fasta -o matches.tsv
Some important points to consider:
Repeat masking is applied to the query and reference sequences by default. To disable it, use --masking 0.
默认情况下是允许重复结果，如果只输出最优结果就加上 --masking 0
DIAMOND is optimized for large input files of >1 million proteins. Naturally the tool can be used for smaller files as well, but the algorithm will not reach its full efficiency.
The program may use quite a lot of memory and also temporary disk space. Should the program fail due to running out of either one, you need to set a lower value for the block size parameter -b.
DIAMOND 是大文件效率更好，对于小文件建议添加 -b 的参数
The sensitivity can be adjusted using the options --fast, --mid-sensitive, --sensitive, --more-sensitive, --very-sensitive and --ultra-sensitive.
比对敏感性，越往后其结果越接近原生 blast 结果，但速度也越慢，一般使用--more-sensitive 比较适中，计算资源不够的就使用 fast。

全参数帮助文件

下面是 DIAMOND 的较为详细的帮助说明：

diamond --help
diamond v2.0.11.149 (C) Max Planck Society for the Advancement of Science Documentation, support and updates available at http://www.diamondsearch.org Please cite: http://dx.doi.org/10.1038/s41592-021-01101-x Nature Methods (2021)
Syntax: diamond COMMAND [OPTIONS]
Commands:
makedb Build DIAMOND database from a FASTA file #以 fasta 文件创建 diamond 格式数据库
blastp Align amino acid query sequences against a protein reference database #功能与原生 blastp 功能一致
blastx Align DNA query sequences against a protein reference database #功能与原生 blastx 一致
view View DIAMOND alignment archive (DAA) formatted file
help Produce help message
version Display version information
getseq Retrieve sequences from a DIAMOND database file
dbinfo Print information about a DIAMOND database file
test Run regression tests
makeidx Make database index
General options:
--threads (-p) number of CPU threads #指定需要运行的线程数，可尽量大
--db (-d) database file #diamond makedb 产生的 diamond 可使用格式的数据库
--out (-o) output file #比对结果输出命名
--outfmt (-f) output format #outfmt，一般选 6 表格格式，与原生 blast 一致
0 = BLAST pairwise
5 = BLAST XML
6 = BLAST tabular
100 = DIAMOND alignment archive (DAA)
101 = SAM Value
6 may be followed by a space-separated list of these keywords:
qseqid means Query Seq - id
qlen means Query sequence length
sseqid means Subject Seq - id
sallseqid means All subject Seq - id(s), separated by a ';'
slen means Subject sequence length
qstart means Start of alignment in query
qend means End of alignment in query
sstart means Start of alignment in subject
send means End of alignment in subject
qseq means Aligned part of query sequence
qseq_translated means Aligned part of query sequence (translated)
full_qseq means Query sequence full
qseq_mate means Query sequence of the mate
sseq means Aligned part of subject sequence
full_sseq means Subject sequence
evalue means Expect value
bitscore means Bit score
score means Raw score
length means Alignment length
pident means Percentage of identical matches
nident means Number of identical matches
mismatch means Number of mismatches
positive means Number of positive - scoring matches
gapopen means Number of gap openings
gaps means Total number of gaps
ppos means Percentage of positive - scoring matches
qframe means Query frame
btop means Blast traceback operations(BTOP)
cigar means CIGAR string
staxids means unique Subject Taxonomy ID(s), separated by a ';' (in numerical order)
sscinames means unique Subject Scientific Name(s), separated by a ';'
sskingdoms means unique Subject Super Kingdom(s), separated by a ';'
skingdoms means unique Subject Kingdom(s), separated by a ';'
sphylums means unique Subject Phylum(s), separated by a ';'
stitle means Subject Title
salltitles means All Subject Title(s), separated by a '<>'
qcovhsp means Query Coverage Per HSP
scovhsp means Subject Coverage Per HSP
qtitle means Query title
qqual means Query quality values for the aligned part of the query
full_qqual means Query quality values
qstrand means Query strand
Default: qseqid sseqid pident length mismatch gapopen qstart qend sstart send evalue bitscore
--verbose (-v) verbose console output
--log enable debug log
--quiet disable console output
--header Write header lines to blast tabular format.
Makedb options:
--in input reference file in FASTA format
--taxonmap protein accession to taxid mapping file
--taxonnodes taxonomy nodes.dmp from NCBI
--taxonnames taxonomy names.dmp from NCBI
Aligner options:
--query (-q) input query file
--strand query strands to search (both/minus/plus)
--un file for unaligned queries
--al file or aligned queries
--unfmt format of unaligned query file (fasta/fastq)
--alfmt format of aligned query file (fasta/fastq)
--unal report unaligned queries (0=no, 1=yes)
--max-target-seqs (-k) maximum number of target sequences to report alignments for (default=25)
--top report alignments within this percentage range of top alignment score (overrides --max-target-seqs)
--max-hsps maximum number of HSPs per target sequence to report for each query (default=1)
--range-culling restrict hit culling to overlapping query ranges
--compress compression for output files (0=none, 1=gzip, zstd)
--evalue (-e) maximum e-value to report alignments (default=0.001)
--min-score minimum bit score to report alignments (overrides e-value setting)
--id minimum identity% to report an alignment
--query-cover minimum query cover% to report an alignment
--subject-cover minimum subject cover% to report an alignment
--fast enable fast mode
--mid-sensitive enable mid-sensitive mode
--sensitive enable sensitive mode)
--more-sensitive enable more sensitive mode
--very-sensitive enable very sensitive mode
--ultra-sensitive enable ultra sensitive mode
--iterate iterated search with increasing sensitivity
--global-ranking (-g) number of targets for global ranking
--block-size (-b) sequence block size in billions of letters (default=2.0)
--index-chunks (-c) number of chunks for index processing (default=4)
--tmpdir (-t) directory for temporary files
--parallel-tmpdir directory for temporary files used by multiprocessing
--gapopen gap open penalty
--gapextend gap extension penalty
--frameshift (-F) frame shift penalty (default=disabled)
--long-reads short for --range-culling
--top 10 -F 15
--matrix score matrix for protein alignment (default=BLOSUM62)
--custom-matrix file containing custom scoring matrix
--comp-based-stats composition based statistics mode (0-4)
--masking enable tantan masking of repeat regions (0/1=default)
--query-gencode genetic code to use to translate query (see user manual)
--salltitles include full subject titles in DAA file
--sallseqid include all subject ids in DAA file
--no-self-hits suppress reporting of identical self hits
--taxonlist restrict search to list of taxon ids (comma-separated)
--taxon-exclude exclude list of taxon ids (comma-separated)
--seqidlist filter the database by list of accessions
--skip-missing-seqids ignore accessions missing in the database
Advanced options:
--algo Seed search algorithm (0=double-indexed/1=query-indexed/ctg=contiguous-seed)
--bin number of query bins for seed search
--min-orf (-l) ignore translated sequences without an open reading frame of at least this length
--freq-sd number of standard deviations for ignoring frequent seeds
--id2 minimum number of identities for stage 1 hit
--xdrop (-x) xdrop for ungapped alignment
--gapped-filter-evalue E-value threshold for gapped filter (auto)
--band band for dynamic programming computation
--shapes (-s) number of seed shapes (default=all available)
--shape-mask seed shapes
--multiprocessing enable distributed-memory parallel processing
--mp-init initialize multiprocessing run
--mp-recover enable continuation of interrupted multiprocessing run
--mp-query-chunk process only a single query chunk as specified
--ext-chunk-size chunk size for adaptive ranking (default=auto)
--no-ranking disable ranking heuristic
--ext Extension mode (banded-fast/banded-slow/full)
--culling-overlap minimum range overlap with higher scoring hit to delete a hit (default=50%)
--taxon-k maximum number of targets to report per species
--range-cover percentage of query range to be covered for range culling (default=50%)
--dbsize effective database size (in letters)
--no-auto-append disable auto appending of DAA and DMND file extensions
--xml-blord-format Use gnl|BL_ORD_ID| style format in XML output
--stop-match-score Set the match score of stop codons against each other.
--tantan-minMaskProb minimum repeat probability for masking (default=0.9)
--file-buffer-size file buffer size in bytes (default=67108864)
--memory-limit (-M) Memory limit for extension stage in GB
--no-unlink Do not unlink temporary files.
--target-indexed Enable target-indexed mode
--ignore-warnings Ignore warnings
View options:
--daa (-a) DIAMOND alignment archive (DAA) file
--forwardonly only show alignments of forward strand
Getseq options:
--seq Sequence numbers to display.
Online documentation at http://www.diamondsearch.org

新版本帮助文件

新版本帮助更简洁，不在一个层次的命令不显示出来，以免混淆。

diamond --help
diamond v2.1.8.162 (C) Max Planck Society for the Advancement of Science, Benjamin Buchfink, University of Tuebingen Documentation, support and updates available at http://www.diamondsearch.org Please cite: http://dx.doi.org/10.1038/s41592-021-01101-x Nature Methods (2021)
Syntax: diamond COMMAND [OPTIONS]
Commands:
makedb Build DIAMOND database from a FASTA file
prepdb Prepare BLAST database for use with Diamond
blastp Align amino acid query sequences against a protein reference database
blastx Align DNA query sequences against a protein reference database
cluster Cluster protein sequences
linclust Cluster protein sequences in linear time
realign Realign clustered sequences against their centroids
recluster Recompute clustering to fix errors
reassign Reassign clustered sequences to the closest centroid
view View DIAMOND alignment archive (DAA) formatted file
merge-daa Merge DAA files
help Produce help message
version Display version information
getseq Retrieve sequences from a DIAMOND database file
dbinfo Print information about a DIAMOND database file
test Run regression tests
makeidx Make database index
greedy-vertex-cover Compute greedy vertex cover
Possible [OPTIONS] for COMMAND can be seen with syntax: diamond COMMAND Online documentation at http://www.diamondsearch.org

要显示更具体的命令下的参数，直接增加功能命令回车即可显示，具体使用大家可在自己系统里面查看即可：

diamond makedb
diamond v2.1.8.162 (C) Max Planck Society for the Advancement of Science, Benjamin Buchfink, University of Tuebingen Documentation, support and updates available at http://www.diamondsearch.org Please cite: http://dx.doi.org/10.1038/s41592-021-01101-x Nature Methods (2021)
Options:
--threads number of CPU threads
--verbose verbose console output
--log enable debug log
--quiet disable console output
--tmpdir directory for temporary files
--db database file
--in input reference file in FASTA format/input DAA files for merge-daa
--taxonmap protein accession to taxid mapping file
--taxonnodes taxonomy nodes.dmp from NCBI
--taxonnames taxonomy names.dmp from NCBI
--file-buffer-size file buffer size in bytes (default=67108864)
--no-unlink Do not unlink temporary files.
--ignore-warnings Ignore warnings
--no-parse-seqids Print raw seqids without parsing
Error: Missing parameter: database file (--db/-d)

6. 结果过滤

1. 根据比对的得分进行过滤

Diamond Blastp 会产生比对得分，你可以根据得分来筛选结果。通过设置一个阈值，只保留得分高于阈值的比对结果。例如：

diamond blastp -d database.fasta -q query.fasta -o output.m8 --min-score 100
### --min-score 100 表示只保留得分高于等于 100 的比对结果。

2. 根据期望的价值（E 值）进行过滤

E 值表示在随机情况下，期望观察到给定比对得分的次数。可以根据 E 值来过滤结果，只保留期望值低于设定阈值的比对结果。例如：

diamond blastp -d database.fasta -q query.fasta -o output.m8 --evalue 1e-5
### --evalue 1e-5 表示只保留 E 值低于或等于 1e-5 的比对结果。

3. 根据相似性阈值过滤

diamond blastp -d database.fasta -q query.fasta -o output.m8 --id 97
### --id 97 表示只保留相似性大于等于 97% 的比对结果。

4. 取唯一结果

### 使用 AWK 命令根据第一个列（query ID）或其他标识符来提取唯一结果
sort -k1,1 -u output.m8 > unique_output.m8

5. Python 脚本处理输出最优结果

# 打开 Diamond Blastp 输出文件
with open('output.m8', 'r') as file:
    best_hit = {}
    # 逐行读取文件
    for line in file:
        fields = line.strip().split('\t')
        # 根据制表符分割字段
        query_id, subject_id, percent_identity, alignment_length, e_value, bit_score = fields[:6]
        # 如果查询 ID 不在 best_hit 中或当前行比最佳结果更好，则更新最优结果
        if query_id not in best_hit or float(bit_score) > float(best_hit[query_id]['bit_score']):
            best_hit[query_id] = {
                'subject_id': subject_id,
                'percent_identity': float(percent_identity),
                'alignment_length': int(alignment_length),
                'e_value': float(e_value),
                'bit_score': float(bit_score)
            }
    # 输出最优结果
    for query_id, hit_info in best_hit.items():
        print(f"Query ID: {query_id}")
        print(f"Subject ID: {hit_info['subject_id']}")
        print(f"Percent Identity: {hit_info['percent_identity']}")
        print(f"Alignment Length: {hit_info['alignment_length']}")
        print(f"E-value: {hit_info['e_value']}")
        print(f"Bit Score: {hit_info['bit_score']}")
        print("-------------")

脚本说明：

• 脚本读取 Diamond Blastp 输出文件，并遍历每一行。
• 对于每个查询 ID，它保留具有最高比对分数（或其他选择的标准）的比对结果。
• 最终输出每个查询 ID 对应的最优匹配结果的相关信息。

6. 使用 R 提取最有结果

# 读取 Diamond Blastp 输出文件
data <- read.table("output.m8", header = FALSE, sep = "\t")
# 命名列名
colnames(data) <- c("query_id", "subject_id", "percent_identity", "alignment_length", "e_value", "bit_score")
# 根据查询 ID 获取最优结果
library(dplyr)
best_hits <- data %>% group_by(query_id) %>% slice(which.max(bit_score))
# 根据最高比对分数选择最优结果，可以根据其他标准替换 bit_score
# 显示最优结果
print(best_hits)

脚本说明：

• 脚本使用 read.table() 函数读取 Diamond Blastp 输出文件。
• 使用 dplyr 包中的 group_by() 和 slice() 函数按照查询 ID 分组，并选择每个查询 ID 的最高比对分数，以获取最优结果。
• 最终输出包含最优匹配结果的数据框。

7. 参考文献

Benjamin Buchfink, Chao Xie, and Daniel H. Huson. Fast and sensitive protein alignment using diamond. Nature methods, 12(1):59–60, Jan 2015.