macでインフォマティクス

macでインフォマティクス

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

超高速でDNAとタンパク質のアライメントを行う AC-DIAMOND

 

 シーケンシング技術の急速な進歩により、微生物の大量シーケンシングデータを作成することが可能になった。このようなデータの解析では、コンティグやリードを大規模にタンパク質データベースに合わせることがしばしばある(例えば、土壌のような複雑なメタゲノミクスデータは何百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

https://github.com/Maihj/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モード、 Sensitive-2 modeを紹介する。

 

1、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)

 

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

 

  

2、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 30 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

ERR2173372.5 YP_001123680.1 100.0 83 0 0 2 250 1173 1255 3.6e-43 178.7

ERR2173372.5 YP_001123702.1 100.0 83 0 0 2 250 1173 1255 3.6e-43 178.7

ERR2173372.5 YP_001123505.1 100.0 83 0 0 2 250 1173 1255 3.6e-43 178.7

ERR2173372.5 YP_001123527.1 100.0 83 0 0 2 250 1173 1255 3.6e-43 178.7

ERR2173372.5 YP_001123768.1 100.0 83 0 0 2 250 1173 1255 3.6e-43 178.7

ERR2173372.5 YP_001123790.1 100.0 83 0 0 2 250 1173 1255 3.6e-43 178.7

ERR2173372.5 NP_051101.1 100.0 83 0 0 2 250 1174 1256 3.6e-43 178.7

ERR2173372.5 NP_051121.1 100.0 83 0 0 2 250 1174 1256 3.6e-43 178.7

ERR2173372.5 YP_006666335.1 100.0 83 0 0 2 250 1171 1253 3.6e-43 178.7

ERR2173372.5 YP_006666353.1 100.0 83 0 0 2 250 1171 1253 3.6e-43 178.7

ERR2173372.5 YP_007889985.1 100.0 83 0 0 2 250 1171 1253 3.6e-43 178.7

ERR2173372.5 YP_007890003.1 100.0 83 0 0 2 250 1171 1253 3.6e-43 178.7

ERR2173372.5 YP_009438097.1 98.8 83 1 0 2 250 1173 1255 6.1e-43 177.9

ERR2173372.5 YP_009438115.1 98.8 83 1 0 2 250 1173 1255 6.1e-43 177.9

ERR2173372.5 YP_009409508.1 98.8 83 1 0 2 250 1173 1255 1.0e-42 177.2

ERR2173372.5 YP_009409527.1 98.8 83 1 0 2 250 1173 1255 1.0e-42 177.2

ERR2173372.5 YP_009409594.1 98.8 83 1 0 2 250 1173 1255 1.0e-42 177.2

ERR2173372.6 NP_063977.1 52.5 80 32 3 27 248 62 141 1.6e-11 73.6

ERR2173372.6 YP_004222306.1 52.5 80 32 3 27 248 62 141 1.6e-11 73.6

ERR2173372.6 YP_004842111.1 52.5 80 32 3 27 248 62 141 1.6e-11 73.6

ERR2173372.9 XP_023644275.1 44.8 87 41 4 7 246 92 178 1.1e-12 77.4

ERR2173372.10 NP_051049.1 89.2 83 9 0 1 249 672 754 8.9e-34 147.5

 

 

200万のシーケンスリードを600万プロテインデータベース に対して検索すると、4分半で結果が出力された。初回だけ必要なindex作成時間を含めても、20分以内で結果を出せることになる。blastxで相同性検索を行うと、ラフに見積もってもおそらく24時間以上待たなければいけない。状況によっては数日かかるかもしれない(大規模データベースでは数百~数千倍の時間がかかる)。簡単なテストだったが、AC-DIAMOND v1が圧倒的な処理時間を達成していることがわかった。AC-DIAMOND はわざわざ"v1"と記載されており、今後さらなる高速化アルゴリズムやオプションが実装される可能性もある。

 

追記

論文表1と表2で、処理時間及び、blastxと比べての感度が議論されています。

https://academic.oup.com/bioinformatics/advance-article/doi/10.1093/bioinformatics/bty391/4996593

 

*1 timeコマンドで"real"の時間を測定。E5  2650 v3 20コア40スレッドの2wayサーバで実施した (DIAMONDはデフォルトではあるだけスレッドを使う)。

引用

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/