macでインフォマティクス

macでインフォマティクス

NGS関連のインフォマティクス情報についてまとめています。

BLASTとコンパチブルで高速なホモロジー検索ツール Diamond

2019 1/20 help追加 、コマンド追記

2019 6/9 -コマンド例から-max-target-seqs削除

2019 7/19 追記

  

Diamondはindexのつけ方を工夫することでBLASTXの解析速度を加速できるツール。blastと同等の機能を持つが、論文ではblastより最大20000倍高速化できると主張されている。特にクエリー配列が非常に多い場合に高速とされる。2015年にnature methodsに論文が発表された。

  

マニュアル

manual

ppt

https://www.donarmstrong.com/ld/dmnd2015/diamond_presentation_2015.pdf

インストール

Github

brewで導入できる。

#anacondaを使っているならcondaで導入可能
conda install -c bioconda -y diamond

#brewでも導入可能
brew install diamond

 > diamond help

$ diamond help

diamond v0.9.13.114 | by Benjamin Buchfink <buchfink@gmail.com>

Licensed under the GNU AGPL <https://www.gnu.org/licenses/agpl.txt>

Check http://github.com/bbuchfink/diamond for updates.

 

Syntax: diamond COMMAND [OPTIONS]

 

Commands:

makedb Build DIAMOND database from a FASTA file

blastp Align amino acid query sequences against a protein reference database

blastx Align DNA query sequences against a protein reference database

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

 

General options:

--threads (-p)         number of CPU threads

--db (-d)              database file

--out (-o)             output file

--outfmt (-f)          output format

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

sseq means Aligned part of 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)

staxids means unique Subject Taxonomy ID(s), separated by a ';' (in numerical order)

stitle means Subject Title

salltitles means All Subject Title(s), separated by a '<>'

qcovhsp means Query Coverage Per HSP

qtitle means Query title

 

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

 

Makedb options:

--in                   input reference file in FASTA format

 

Aligner options:

--query (-q)           input query file

--strand               query strands to search (both/minus/plus)

--un                   file for unaligned queries

--unal                 report unaligned queries (0=no, 1=yes)

--max-target-seqs (-k) maximum number of target sequences to report alignments for

--top                  report alignments within this percentage range of top alignment score (overrides --max-target-seqs)

--overlap-culling      delete hits only if a higher scoring hit envelops the given percentage of the query range

--compress             compression for output files (0=none, 1=gzip)

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

--sensitive            enable sensitive mode (default: fast)

--more-sensitive       enable more sensitive mode (default: fast)

--block-size (-b)      sequence block size in billions of letters (default=2.0)

--index-chunks (-c)    number of chunks for index processing

--tmpdir (-t)          directory for temporary files

--gapopen              gap open penalty

--gapextend            gap extension penalty

--matrix               score matrix for protein alignment (default=BLOSUM62)

--custom-matrix        file containing custom scoring matrix

--lambda               lambda parameter for custom matrix

--K                    K parameter for custom matrix

--comp-based-stats     enable composition based statistics (0/1=default)

--masking              enable masking of low complexity 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

--taxonmap             protein accession to taxid mapping file

--taxonnodes           taxonomy nodes.dmp from NCBI

 

Advanced options:

--algo                 Seed search algorithm (0=double-indexed/1=query-indexed)

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

--window (-w)          window size for local hit search

--xdrop (-x)           xdrop for ungapped alignment

--ungapped-score       minimum alignment score to continue local extension

--hit-band             band for hit verification

--hit-score            minimum score to keep a tentative alignment

--gapped-xdrop (-X)    xdrop for gapped alignment in bits

--band                 band for dynamic programming computation

--shapes (-s)          number of seed shapes (0 = all available)

--shape-mask           seed shapes

--index-mode           index mode (0=4x12, 1=16x9)

--rank-ratio           include subjects within this ratio of last hit (stage 1)

--rank-ratio2          include subjects within this ratio of last hit (stage 2)

--max-hsps             maximum number of HSPs per subject sequence to save for each query

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

 

View options:

--daa (-a)             DIAMOND alignment archive (DAA) file

--forwardonly          only show alignments of forward strand

 

Getseq options:

--seq                  Sequence numbers to display.

 

 

 

ラン

はじめにデータベースとなるアミノ酸配列のindexファイルを作成する。

diamond makedb --in input.faa -d nr

 

 

blastxでホモロジー検索を行う。inputは塩基配列

diamond blastx -d nr -q query.fna -o matches.m8 
  • --threads   Number of CPU threads. By default, the program will auto-detect and use all available virtual cores on the machine.

 出力はタブ区切り形式である。

 user$ head matches.m8 

gi|451813329|ref|NC_020286.1|:3569362-3569561,1-772 gi|451813330|ref|YP_007449782.1| 100.0 323 0 0 1 969 1 323 1.8e-179 622.1

gi|451813329|ref|NC_020286.1|:3569362-3569561,1-772 gi|451813441|ref|YP_007449893.1| 33.0 233 142 4 82 747 34 263 8.8e-25 108.2

 Diamondの検出閾値はblastのdefaultの検出閾値よりずっと低いため、stringencyはblastより高くなっている。また、defaultのパラメータはショートリード向けの設定のため、クエリ配列が長い場合、--sensitiveや--more-sensitiveをつけることが推奨されている。

 

マニュアルに書いてあるが、タブ出力するには--outfmt 6をつける。さらに以下のような指定を行うことで、出力項目を好きに設定できる。

  • qseqid Query Seq - id
  • qlen Query sequence length
  • sseqid Subject Seq - id
  • sallseqid All subject Seq - id(s), separated by a ’;’ slen Subject sequence length
  • qstart Start of alignment in query
  • qend End of alignment in query
  • sstart Start of alignment in subject
  • send End of alignment in subject
  • qseq Aligned part of query sequence
  • sseq Aligned part of subject sequence
  • full sseq Full subject sequence
  • evalue Expect value
  • bitscore Bit score
  • score Raw score
  • length Alignment length
  • pident Percentage of identical matches
  • nident Number of identical matches
  • mismatch Number of mismatches
  • positive Number of positive - scoring matches
  • gapopen Number of gap openings
  • gaps Total number of gaps
  • ppos Percentage of positive - scoring matches
  • qframe Query frame
  • btop Blast traceback operations(BTOP)
  • staxids Unique Subject Taxonomy ID(s), separated by a ’;’ (in numerical order). This field requires setting the --taxonmap parameter for makedb.
  • salltitles All Subject Title(s), separated by a ’<>’
  • qcovhsp Query Coverage Per HSP
  • qtitle Query title

デフォルトでは

"qseqid sseqid pident length mismatch

gapopen qstart qend sstart send evalue bitscore"が出力されるようになっている。

例えばdiamondでblastxサーチを行う。tabularでqseqid、sseqid、evalueのみ出力する。

diamond blastx --query input.fa \
--db uniprot_ref_proteomes.diamond.dmnd \
--outfmt 6 qseqid sseqid evalue \
--sensitive \
> blast.out

 

追記

E-value ≤1E-100、amino acid identity ≥ 80%、minimum length (amino acid) ≥ 300 a.a。

diamond blastx --query input.fa \
--db uniprot_ref_proteomes.diamond.dmnd \
--outfmt 6 qseqid sseqid evalue \
--min-orf 300 --id 80 --evalue 1e-100
> blast.out

 

引用

Fast and sensitive protein alignment using DIAMONDFast and sensitive protein alignment using DIAMOND
Benjamin Buchfink, Chao Xie & Daniel H

Nature Methods 12, 59–60 (2015) doi:10.1038/nmeth.3176

PDF

https://lemosbioinfo.files.wordpress.com/2016/11/nmeth-3176.pdf

 

追記

Question about max-target-seqs option #29

Question about max-target-seqs option · Issue #29 · bbuchfink/diamond · GitHub

 

KMCのホームページで、diamondとKMCの連携について提案があります。詳細はKMCのHPからスクリプトをダウンロードして確認してください。

http://sun.aei.polsl.pl/REFRESH/index.php?page=projects&project=kmc&subpage=download

 

追記

 論文図1で色々なシーケンサー由来のリードを使ってタンパク質と相同性検索した時の処理時間と感度が比較されています。データによっては2万倍以上高速化しています。

https://lemosbioinfo.files.wordpress.com/2016/11/nmeth-3176.pdf

 

追記

ローカルデータベースのダウンロードについては、こちらが参考になります。

https://www.uppmax.uu.se/resources/databases/diamond-protein-alignment-databases/

 

追記

AC-DIAMOND