2019 1/17 condaインストール追記
2019 1/29 追記
2020 7/23 誤字修正
2020 7/24 help更新、例修正
シーケンシング技術の急速な進歩により、微生物の大量シーケンシングデータを作成することが可能になった。このようなデータの解析では、コンティグやリードを大規模にタンパク質データベースに合わせることがしばしばある（例えば、土壌のような複雑なメタゲノミクスデータは何百Gbps、およびNCBI-nrのようなタンパク質データベースのアミノ酸20G以上を含む）。 BLASTX（Altschul et al、1990）は、その優れた感度からDNA-タンパク質アラインメントのゴールドスタンダードとみなされてきたが、大規模アラインメントには遅すぎる。たとえば、12コアのサーバーを使用して10Gbpの土壌コンティグをNCBI-nrにアライメントすることを検討すると、 BLASTXでは何年もかかる。過去10年間に、BLASTXをスピードアップするためのツールがいくつか発表されてきている。初期のツールはRAPSearch2（Zhao et al、2012）で、上記の例で言えば数週間しかかからない。 PAUDA（Huson and Xie、2014）やGHOSTZ（Suzuki et al、2015）のような他のツールは、スピードと感度のトレードオフで、アライメント時間をさらに短縮することができる。 DIAMOND（Buchfink et al、2015）（紹介）は、感度を犠牲にすることなく速度を大幅に向上させる最近のツールである。上記の例では、6日以上かかる。シーケンサーがデータを生成する時間よりも6日以上長いことに注意する。この論文では、DIAMONDと同じ感度を維持しながら、DIAMONDよりも6倍から7倍のスピードアップを持つ、より効率的なDNA-タンパク質アライメントツールAC-DIAMOND（バージョンv1）を紹介する。 上記の10Gbpの土壌コンティグをアライメントさせる例では、AC-DIAMOND を使うと20時間超で終わる。
AC-DIAMONDは、複雑なデータ構造とアルゴリズムを設計してアライメントシードの位置を特定し、シードからのアライメントを計算するためのダイナミックプログラミングを並列化することによって、DIAMONDをスピードアップしようとする。（一部略）。Seeding processの最初の改善点は、二重インデックスを圧縮して、より多くのシーケンスを処理し、参照の再ロードを減らすことであった（これはAC-DIAMOND（すなわちバージョンv0）の予備バージョンになった。最初にIWBBIO 2016で発表され、4倍のスピードアップを実現した）。より速いスピードアップを達成するために、クエリのindexingを避け、リファレンスのindexを2レベルのsuffix-array様の データ構造で再設計する。これにより、キャッシュ・メモリ内で非常に効率的にバイナリ検索を行うことができる。（以下 段落略）。まとめると、AC-DIAMONDは、128ビットのSIMD機能を備えた1つのCPUコアを使用して、最大16のDBテーブルを並列処理することで速度を大幅に向上させる。
インストール
cent os6でテストした。
本体 Github
#anacondaを使っているならcondaで導入可（linuxのみ）
conda install -c bioconda -y ac-diamond
git clone https://github.com/Maihj/AC-DIAMOND.git
cd AC-DIAMOND/src/
./install-boost
make
cd ../bin/
> ./ac-diamond
$ ac-diamond -h
Syntax:
ac-diamond COMMAND [OPTIONS]
Commands:
makedb Build AC-DIAMOND database from a FASTA file
align Align DNA query sequences against a protein reference database
view View AC-DIAMOND alignment archive (DAA) formatted file
General options:
-h [ --help ] produce help message
-p [ --threads ] arg (=0) number of cpu threads
-d [ --db ] arg database file
-a [ --daa ] arg AC-DIAMOND alignment archive (DAA) file
-v [ --verbose ] enable verbose out
--log enable debug log
Makedb options:
--in arg input reference file in FASTA format
-b [ --block-size ] arg reference sequence block size in billions of letters
(default=4)
--sensitive enable building index for sensitive mode
(default:fast)
Aligner options:
-z [ --query-block-size ] arg (=6) query sequence block size in billions of
letters (default=6)
-q [ --query ] arg input query file
-k [ --max-target-seqs ] arg (=25) maximum number of target sequences to
report alignments for
--top arg (=100) report alignments within this percentage
range of top alignment score (overrides
--max-target-seqs)
--compress arg (=0) compression for output files (0=none,
1=gzip)
-e [ --evalue ] arg (=0.001) maximum e-value to report alignments
--min-score arg (=0) minimum bit score to report alignments
(overrides e-value setting)
--id arg (=0) minimum identity% to report an alignment
--sensitive enable sensitive mode (default: fast)
-t [ --tmpdir ] arg (=/dev/shm) directory for temporary files
--gapopen arg (=-1) gap open penalty, -1=default (11 for
protein)
--gapextend arg (=-1) gap extension penalty, -1=default (1 for
protein)
--matrix arg (=blosum62) score matrix for protein alignment
--seg arg enable SEG masking of queries (yes/no)
Advanced options (0=auto):
-w [ --window ] arg (=0) window size for local hit search
--xdrop arg (=20) xdrop for ungapped alignment
-X [ --gapped-xdrop ] arg (=20) xdrop for gapped alignment in bits
--ungapped-score arg (=0) minimum raw alignment score to continue local
extension
--hit-band arg (=0) band for hit verification
--hit-score arg (=0) minimum score to keep a tentative alignment
--band arg (=0) band for dynamic programming computation
--index-mode arg (=0) index mode (1=10x1, 2=10x8)
--fetch-size arg (=4096) trace point fetch size
--single-domain Discard secondary domains within one target
sequence
View options:
-o [ --out ] arg output file
-f [ --outfmt ] arg (=tab) output format (tab/sam)
--forwardonly only show alignments of forward strand
ラン
3つのモード; fast mode、Sensitive-2 mode、Sensitive-1 modeがある。ここではfast mode、Sensitive-2 modeを紹介する。
１、fast mode
indexを作成。
#index作成
ac-diamond makedb --in input_pep.fa -d acd.fast.index -b 4
- -b reference sequence block size in billions of letters (default=4)
- --sensitive enable building index for sensitive mode (default:fast)
indexファイルサイズは場合によっては配列自身より巨大になる。作業ディスクの空き容量に注意する。Uniref90の場合
DNA配列をアライメントする。
ac-diamond align -d acd.fast.index -q read.fa -a acd.result -e 0.001 -z 6
- -d database file
- -q Path to query input file in FASTA or FASTQ format (may be gzip compressed).
- -a AC-DIAMOND alignment archive (DAA) file
- -e maximum e-value to report alignments
- -z query sequence block size in billions of letters (default=6)
- -p number of cpu threads (default=max)。
出力はblastの出力フォーマットになっていないので最後に変換する。
blastxコンパチブル（-outfmt 6）に変換。
ac-diamond view -a acd.result.daa -o acd.result.m8
- -a AC-DIAMOND alignment archive (DAA) file
- -o output
出力
> head -n 10 acd.result.m8
$ head -n 10 acd.result.m8
ERR2173372.2 AT2G08986.1 81.0 42 8 0 241 116 237 278 2.8e-13 72.4
ERR2173372.2 AT2G08986.1 76.2 42 10 0 241 116 481 522 2.3e-12 69.3
ERR2173372.2 AT2G08986.1 78.6 42 9 0 241 116 1020 1061 4.0e-12 68.6
ERR2173372.2 AT2G08986.1 76.2 42 10 0 241 116 1130 1171 6.8e-12 67.8
ERR2173372.2 AT2G08986.1 51.3 78 20 4 235 2 63 122 8.9e-12 67.4
ERR2173372.2 AT2G08986.1 76.2 42 10 0 241 116 914 955 2.0e-11 66.2
ERR2173372.2 AT2G08986.1 70.5 44 13 0 244 113 544 587 7.5e-11 64.3
ERR2173372.2 AT2G08986.1 82.8 29 5 0 206 120 114 142 1.3e-07 53.5
２、Sensitive-2 mode
#index
ac-diamond makedb --in ref.fa -d acd.sensitive.index --sensitive -b 4
#align
ac-diamond align -d acd.sensitive.index -q query.fa -a acd.result --sensitive -e 0.001 -z 6
#convert
ac-diamond view -a acd.result.daa -o acd.result.m8
他にindexをfastモードで行い、アライメントを--sensitiveをつけて実行するSensitive-1 modeがある。詳細はGithubで確認してください。
テスト１
NCBI のRefseqのプロテイン配列をデータベースに使って、fastモードの実行時間を調べてみる。
#refseq proteinのダウンロード（合計1430ファイル、43GB）
wget ftp://ftp.ncbi.nlm.nih.gov/refseq/release/complete/*.faa.gz
全部使うのはしんどいので、一部だけ使う（およそ3GB）。
gzip -dv complete.2*
cat complete.2* > refseq2000_protein.faa
index。
ac-diamond makedb --in refseq2000_protein.faa -d acd.fast.index -b 4
index作成はシングルスレッドなので、大きなデータベースだと多少時間がかかる。3GBのfaaで10分程度かかった。
シロイヌナズナのraw シーケンスデータ（リンク）（fasta）を200万リード（500MB）使って相同性検索。上限閾値は1e-10とした。
ac-diamond align -d acd.sensitive.index -q query.fa -a acd.result --sensitive -e 1e-10 -z 6
4分40秒で結果が出力された（*1）。出力のdaaファイルのファイルサイズは351MBだった。
blasstコンパチブルに 変換。
ac-diamond view -a acd.result.daa -o acd.result.m8
マルチスレッドなのですぐに終わる。10秒かからなかった。デフォルトでは1リードから複数ヒットするため、合計でリード数より多い1017万ヒットがあった。
> head -n 8 acd.result.m8
$ head -n 30 acd.result.m8
ERR2173372.5 YP_009409765.1 100.0 83 0 0 2 250 1171 1253 3.6e-43 178.7
ERR2173372.5 YP_009409784.1 100.0 83 0 0 2 250 1171 1253 3.6e-43 178.7
ERR2173372.5 YP_009409678.1 100.0 83 0 0 2 250 1175 1257 3.6e-43 178.7
ERR2173372.5 YP_009409698.1 100.0 83 0 0 2 250 1175 1257 3.6e-43 178.7
ERR2173372.5 YP_009409851.1 100.0 83 0 0 2 250 1162 1244 3.6e-43 178.7
ERR2173372.5 YP_009409868.1 100.0 83 0 0 2 250 1173 1255 3.6e-43 178.7
ERR2173372.5 YP_009438012.1 100.0 83 0 0 2 250 1173 1255 3.6e-43 178.7
ERR2173372.5 YP_009438030.1 100.0 83 0 0 2 250 1173 1255 3.6e-43 178.7
200万のシーケンスリードを600万プロテインデータベース に対して検索すると、４分半で結果が出力された。初回だけ必要なindex作成時間を含めても、20分以内で結果を出せることになる。
テスト2
uniprotのプロテオームをデータベースにしてメタゲノムのアセンブル配列（DNA配列）のベストヒットを探索。uniprotデータの準備はblobtoolsを参考にした（blobtools紹介）。（*2）。
#fastモードでindex作成
ac-diamond makedb --in uniref90.fasta.gz -d unbiref90 -b 4
#max target seq=1で実行
ac-diamond align -d unbiref90 -q metagenome_assmbly.fa -a result -e 1e-25 -z 6 --max-target-seqs 1
#blastxのoutfmt=6出力フォーマットに変換
ac-diamond view -a result.daa -o result
uniprot proteomeデータベースに対してcontig20000配列をクエリとしてベストヒット探索にかかった時間はわずか5分30秒ほどだった（*3）。diamondでは103分ほどだったので、20倍ほど早く終わったことになる。
blastxで相同性検索を行うと、ラフに見積もってもおそらく24時間以上待たなければいけない。状況によっては数日かかるかもしれない（大規模データベースでは数百~数千倍の時間がかかる）。簡単なテストだったが、AC-DIAMOND v1が圧倒的な処理時間を達成していることがわかった。AC-DIAMOND はわざわざ"v1"と記載されており、今後さらなる高速化アルゴリズムやオプションが実装される可能性もある。
引用
AC-DIAMOND v1: Accelerating large-scale DNA-protein alignment.
Mai H, Zhang Y, Li D, Leung HC, Luo R, Wong CK, Ting HF, Lam TW
Bioinformatics. 2018 May 16.
Biostars
Question: Download all RefSeq proteins from all organisms in one faa-file
https://www.biostars.org/p/130274/
*1 timeコマンドで"real"の時間を測定。E5 2650 v3 20コア40スレッドの2wayサーバで実施した （DIAMONDはデフォルトでは利用できる全CPUを使う（ワークステーション環境））。
*2
ランするにはある程度メモリが必要。メモリ180GBのサーバーではランできたが、メモリが128GBのサーバーではランできなかった。
*3
xeon E5 v3 2650x2、180GBメモリ、7200rpm 3TB HDDのサーバで実行。
追記
論文表１と表２で、処理時間及び、blastxと比べての感度が議論されています。
https://academic.oup.com/bioinformatics/advance-article/doi/10.1093/bioinformatics/bty391/4996593
