Hatena::ブログ(Diary)

それなりにマジメなメモ このページをアンテナに追加 RSSフィード Twitter

2015-12-28

高速なタンパク質配列相同性検索ツール「DIAMOND」

22:05 | 高速なタンパク質配列相同性検索ツール「DIAMOND」を含むブックマーク 高速なタンパク質配列相同性検索ツール「DIAMOND」のブックマークコメント

この記事は今年読んだ一番好きな論文 Advent Calendar 2015の17日の記事です。

今回紹介する論文タンパク質用の配列相同性検索に関する以下の論文です。

Buchfink, B., Xie, C., & Huson, D. H. (2015). Fast and sensitive protein alignment using DIAMOND. Nature Methods, 12(1), 59–60. doi:10.1038/nmeth.3176

まったく同じ問題に大学時代から取り組んでいて、いろいろな意味で刺激を受けた論文で、良い機会なので紹介しようと思います。


配列相同性検索とは

配列相同性検索とは、DNAタンパク質配列クエリとしてデータベース内の類似する配列を検索する解析です。この検索を行ううえで、クエリ配列データベース内の配列の類似度をどのように計算するか?が重要になります。理想的には、すべてのクエリデータベース内の配列のペアに対して、Smith-Watermanアルゴリズム[1]を利用して類似度を計算し、クエリごとにスコアの高い順に類似するデータベース配列名を出力することが望ましいです。しかし、この方法ではデータ量が多い場合、非常に多くの計算時間が必要となります。

このため、Seed-And-Extend strategyと呼ばれる近似手法が利用されます。この手法は以下の2つのステップで配列相同性検索を実行します。

  1. クエリデータベース配列で短く、かつ、類似度が非常に高い領域(seed)を検索
  2. Seedの周りのみSmith-Watermanベースのスコア計算を実行

ステップ1で、あらかじめ当たりをつけることで大部分のスコア計算を実行せずに済みます。これにより、Smith-Watermanアルゴリズムをそのまま使う場合よりも、大幅な高速化が可能となります。


このSeed-And-Extend strategyを利用した有名なツールとしてBLAST[2]があります。タンパク質配列相同性検索の場合、BLASTではステップ1でクエリデータベース配列間で3文字の類似する部分文字列を2つ検索し、seedとしています。BLASTは現在でも多くの解析で用いられる一般的なツールであり、配列相同性検索を実行する場合、これが利用されます。


DIAMONDとは

近年、DNAシーケンサーの性能向上に伴う、解析データの大規模化に対応するために開発された、Seed-And-Extend strategyベースのタンパク質用の配列相同性検索ツールです。DIAMONDではBLASTとは違い、spaced seedとreduced alphabetを利用することで高い検索精度を維持しつつ、BLASTの2万倍の高速化を実現しています。


Spaced seed

高い検索精度が必要な場合、ある程度ミスマッチを許容したseedの検索が必要になってきます。しかし、ミスマッチを許容すると、途端に計算量が大きくなります。この問題を解決する方法の一つにspaced seedがあります。

Spaced seedは連続した文字列のマッチングではなく、飛ばし飛ばしの文字列の完全一致の検索を行います。以下にspaced seedを利用した文字列のマッチングの例を示します。


f:id:ang65:20151228212900p:image:w360

Spaced seedでは、1となっている部分のみ一致しているかどうか確認し、0のところは無視します。


詳細は省きますが、spaced seedによって連続した文字列の完全一致よりも精度の高い検索を行うことができます。この結果、seedの検索の際、連続した文字列検索の場合よりも長い文字列検索を用いても精度を維持できるようになります。1文字でも長い文字列の検索ができると、seedの検索で見つかるseedの数が平均1/|¥Sigma|^2(¥Sigma:アルファベットの集合)となります。このため、その次のステップの計算時間が大幅に削減できます。


実際、DNA配列の場合は非常にうまくいくことが知られています。ただし、タンパク質配列の場合、これまでspaced seedを利用してもあまり精度が上がらないため、使われることはほとんどありませんでした。しかし、DIAMONDでは複数のspaced seedのパターンを利用するmultiple spaced seedと次に詳しく説明するreduced alphabetを組み合わせることで高い検索精度を維持しつつ、高速化を実現しています。


Reduced alphabet

タンパク質は一般的には20種類のアミノ酸で構成されています。配列相同性検索の場合、この20種類ごとに別々の文字を割り当てて、タンパク質文字列で表現します。しかし、近い性質をもつアミノ酸がいくつかあります。このため、近い性質のアミノ酸は同じ文字で表現し、解析するということが近年良く行われます。DIAMONDにおいても、seedの検索の段階では、20種類のアミノ酸を11種類に減らした配列を変換し、変換後の配列seedの検索を実行します。これにより、検索を行いやすい完全一致の文字列検索でも、ミスマッチを許した検索を実行した場合と同じように高い検索精度を保つことができます。


Seedの検索

DIAMONDのseedの検索では先ほど紹介した通り、multiple spaced seedとreduced alphabetを利用します。このうち、multiple spaced seedを利用する場合、以下の2つの問題点があります。


  1. Spaced seedのパターン数分、検索回数が増加
  2. Spaced seedのパターン数分、検索インデックスのメモリ量が増加

まず、DIAMONDでは検索回数の増加に対して、検索の際になるべくメモリのキャッシュヒット率を上げることで、seedの検索を高速化して、問題を解決しています。このキャッシュヒット率の向上に関係データベースなどで用いられるsort merge joinと同じような方法を利用します。具体的には以下のように行います。

  1. spaced seedに基づいてデータベース内の配列の部分文字列を生成
  2. できた部分文字列を辞書式順序に並び替え
  3. 上記と同様に複数のクエリ配列の部分文字列を生成し、並び替え
  4. データベース配列クエリ配列の部分文字列配列を上から順にマッチング

大量なクエリを一度に処理する場合は、この方法により高速に複数の部分文字列の検索を行うことができます。


そして二つ目の問題点のメモリ量の増加に関してですが、これに対しては「メモリが足りなければデータベースクエリを分割して処理する」だけです。DIAMONDの場合、データベースのサイズ×8バイト分のメモリがインデックスに必要になります。このため、なんらかの対処が必要になると思っていたのですが、DIAMONDでは分割するだけでインデックスなどはそのままメモリにのっけるという男らしい方法を取ります。

(このメモリを気にして、今までmultiple spaced seedを使わなかっただけに正直驚きました。)


実験結果

DIAMONDの評価ではBLASTと比較して評価 (Figure 1, Supplementary Table 1)を行っています。KEGGデータベースとIlluminaのreadを利用した場合、DIAMOND-sensitive (16パターンのspaced seedを利用)の場合、BLASTと92%のアラインメントが一致し、2000倍の高速化を達成しています。また、DIAMOND-fast(4パターンのspaced seedを利用)の場合、BLASTと69%のアラインメントが一致し、約20,000倍の高速化を達成しています。


感想

正直、論文を読んだだけでは良くわからない部分が多々あります。また、DIAMONDのソースコードを読むと、論文にかかれてないフィルタがいくつか入っているため、本当にspaced seedでよくなっているのか若干疑わしいです。

また、精度に関しても、この評価でいいのか?と疑問もあります。(いろいろなクエリで試すとDIAMOND-sensitiveでも精度が非常に落ちる場合があります。)

ただ、タンパク質配列でまじめにspaced seedを使って精度と速度を出しているすごいツールなので、今後どう発展していくのか楽しみです。私のツール(GHOSTX, GHOSTZ)も負けないように改良を続けていこうと思います。




Reference

1. Smith, T. F., & Waterman, M. S. (1981). Identification of common molecular subsequences. Journal of Molecular Biology, 147(1), 195–197. doi:10.1016/0022-2836(81)90087-5

2. Altschul, S. (1997). Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Research, 25(17), 3389–3402. doi:10.1093/nar/25.17.3389

2015-11-28

SC15のBioinformaticsセッションの論文紹介

12:18 | SC15のBioinformaticsセッションの論文紹介を含むブックマーク SC15のBioinformaticsセッションの論文紹介のブックマークコメント

11/15-11/20にスパコンのランキングが発表されることで有名なSCがあった。今年もbioinformatics関係のセッションがあり、3つ論文があったのでせっかくなので読んでみたので簡単な論文紹介と感想の記事。

今回、発表のあった論文は以下の3つ。

  • Evangelos Georganas, Aydın Buluç, Jarrod Chapman, Steven Hofmeyr, Chaitanya Aluru, Rob Egan, Leonid Oliker, Daniel Rokhsar, and Katherine Yelick. 2015. HipMer: an extreme-scale de novo genome assembler. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC '15). ACM, New York, NY, USA, , Article 14 , 11 pages. DOI=http://dx.doi.org/10.1145/2807591.2807664
  • Patrick Flick, Chirag Jain, Tony Pan, and Srinivas Aluru. 2015. A parallel connectivity algorithm for de Bruijn graphs in metagenomic applications. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC '15). ACM, New York, NY, USA, , Article 15 , 11 pages. DOI=http://dx.doi.org/10.1145/2807591.2807619
  • Patrick Flick and Srinivas Aluru. 2015. Parallel distributed memory construction of suffix and longest common prefix arrays. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (SC '15). ACM, New York, NY, USA, , Article 16 , 10 pages. DOI=http://dx.doi.org/10.1145/2807591.2807609

HipMer: an extreme-scale de novo genome assemble

3本の中で一番bioinformaticsっぽい論文。HipMerというMeraculous[1]をマルチノードに拡張したツールの話。

昨年のSC14でも同じようにMeraculousの分散並列版の話をしていたが、contigの生成のところに新しい最適化を入れたのとscaffoldingの部分も分散並列させたという部分が新しい。

これによってMeraculousの全部の処理を分散環境で実行できるようになった。

実験として、ヒトのデータ(101 bpのread、2.9 billion本)のアセンブルを12コアのマシン1280ノード(合計15360コア)使って8.4分で処理させることができたとのこと。


スケーリング自体は悪いものの他のマルチノードアセンブラであるRayやABySSと比べても速度が10倍以上でているらしいのでうまくタスク分散できている模様。

Scaffoldingの部分の並列化はデータ並列である程度並列化がでるようなのであまり面白くなかったが、contig生成の部分で新しく追加された最適化である頻度の高いk-merの判定(Section 3.1)とde Bruijn Graphの探索時の通信回数削減方法(Section 3.2)が面白かった。

A parallel connectivity algorithm for de Bruijn graphs in metagenomic applications

AssemblerにGraph500でおなじみBreadth First Search (BFS)を適応するという話だったのでde Bruijn Graphの探索部分に使うのか?と思ったが違った。Metagenomeのデータだとhigh species level(?)でde Bruijn graphsがつながっていないことが多いらいし、ということが[2]で示されているとか。このため、最初にde Bruijn graphsの連結を調べて、グラフを分割してしまえばあとは個別にアセンブルできる。([2]を読んでないので私自身どれくらい分割できるかとかわかっていない)

この最初の分割をBFSでやってしまうというのがこの論文で紹介されていること。なのでメインの話はBFSを如何に高速に処理するかの部分で、BFSを使ってどのようにde Bruijn graphsを分割するかがおまけ程度に書かれている。

BFSの応用にこんなのがあったのかという発見があったので、その部分は面白かった。

Parallel distributed memory construction of suffix and longest common prefix arrays

3つの論文でバイオっぽい単語が一番出てこなかった論文。タイトルの通り、suffix arrayとlcpの構築を分散環境で高速に行う話。

分散環境においてのsuffix arrayの構築はすごい通信が多くなりそうで並列度でないのでは?と思ったがそうでもないようで128コアから1024コアの部分でみるとstrong scalingで0.5でているのでそこそこでている。

大本のアルゴリズムprefix doublingで並列用に少し手を加えたものを利用している。

並列化方法はall-to-all、all-reduce,、scan、exscanなどの並列化するときによく用いられている処理を駆使して行っている。

一番の重要な部分はsortのときの通信回数を減らすためにメモリの局所性を意識して実装されているところ。これによって速度がかなり上がっている模様。

個人的に驚いたのは1024コアでヒトゲノムsuffix arrayがたった5秒でできてしまうところ。

コア数が多い気もするがこれだけ速ければsuffix arrayの構築で時間がかかるから使うのをあきらめてた部分とかに使えないかな?


全体の感想

昨年よりもバイオ感がかなり減っていたが、面白い部分がいくつかあった。

ただ、問題設定がde novo assemblerに偏ってたのでもう少しいろいろbioinformaticsHPCを活用しなきゃならない問題がないのかな?と思った。



Reference

  1. Meraculous: De Novo Genome Assembly with Short Paired-End Reads Chapman JA, Ho I, Sunkara S, Luo S, Schroth GP, et al. (2011) Meraculous: De Novo Genome Assembly with Short Paired-End Reads. PLoS ONE 6(8): e23501. doi: 10.1371/journal.pone.0023501
  2. A. C. Howe, J. K. Jansson, S. A. Malfatti, S. G. Tringe, J. M. Tiedje, and C. T. Brown. Tackling soil diversity with the assembly of large, complex metagenomes. Proceedings of the National Academy of Sciences, 111(13):4904--4909, 2014.

2014-03-31

GTC2014の参加レポート!

00:18 | GTC2014の参加レポート!を含むブックマーク GTC2014の参加レポート!のブックマークコメント

GTCとは?

NVIDIAが主催しているGPUを使った技術についてのカンファレンスです。GPUといえば画像処理を思い浮かべるかもしれませんが、General-purpose computing on graphics processing units (GPGPU)と呼ばれるような画像処理だけでなくGPUを使ったものなら何でもありカンファレンスです。

日本でもここ最近は毎年5月くらいにGTC Japanとして開催されていますが、今回は2014年3月24日〜3月27日まで開催された本家のGTCのレポートです。

GTCの公式サイトは以下のところです。

http://www.gputechconf.com/page/home.html

また、以前のGTCのトークは以下のところで見れます。

http://www.gputechconf.com/gtcnew/on-demand-gtc.php


今年の注目トピック

今年はなんと言っても機械学習に関するセッションです。数年前からDeep Neural Network (DNN) が話題になってますが、このDNNの学習でGPUがよく使われます。

この関係で今回はKeynote機械学習について言及された他、機械学習に関する多くのトークがりました。

DNNについて詳しく知りたい人はPFIブログのDNNの記事を見るとわかりやすいのでお勧めです。

http://research.preferred.jp/2012/11/deep-learning/

DNNを使って画像認識をはじめ、音声認識などにも活用した話がありました。

DNNでTheano使うって話良く聞くけど音声認識ならkaldiを使うっててもあるらしいです。

また、私が見ていた限りではDNN以外にもRandom Forestの構築にGPUを使った話がいくつかあったようです。

個人的な感想

自分がバイオインフォな研究を行っている関係でそれっぽいものを見てましたが今年はバイオインフォっぽい話が去年よりも多かったように思います。

中でも個人的に気になったのは以下のもの。

GPU Accelerated Genomics Data Compression (BingQiang Wang)

この業界で超有名な中国BGIの方が発表してたfastqファイルの圧縮の話。いくつかの圧縮手法GPU版を実装して、それらを組み合わせて圧縮率、圧縮の計算時間、解凍の計算時間について発表してました。特に気になったのがクオリティスコアの部分でMarkov transform + Huffmanで40〜50%くらいの圧縮率を出してました。クオリティスコアってほぼランダムな配列って聞いてたんですが、Markov transform効くんですね。多少は前の文字と相関があるのかな?

BWT Indexing: Big Data from Next Generation Sequencing and GPU (Jeanno Cheung)

FM-indexをGPUで構築する話。BWAなど多くのDNAのmapperでFM-indexが使われますが、これ速くしても?って思ったら2010年くらいにFM-indexを使ったde novo assemblerが発表されてて、それを意識して行われた研究の話らしいです。このassemblerがとっても気になるので後で調べる予定。


Introducing NVBIO: High Performance Primitives for Computational Genomics (Jonathan Cohen, Nuno Subtil)

NVIDIAが作っているDNAなどの配列解析をGPUで手軽に行うためのライブラリで、アラインメントなど基本的なことのいくつかはすでにライブラリの中に入ってて簡単な関数呼び出しでGPUで計算することができます。去年発表で気になってたBowtie2のGPU版もこのライブラリで実装されたものらしいです。

このライブラリの場所は以下のところ。

https://github.com/NVlabs/nvbio

最後に

GPUがどんどん進化してきたり、新しい手法で前まではあまりGPU向きでないと思われてきたものもGPUで計算できるようになってきたりしているので今後もGPUに注目ですね。

2013-04-04

Xeon Phiのプログラミングについての情報(2) Offloadによるプログラミング

14:17 | Xeon Phiのプログラミングについての情報(2) Offloadによるプログラミングを含むブックマーク Xeon Phiのプログラミングについての情報(2) Offloadによるプログラミングのブックマークコメント

Xeon Phiを使ったプログラムの実行としては以下の方法がある。

今回はOffloadを使用したプログラムを作ったのでその解説を簡単にする。

Offloadを使用したプログラミングについて

Offloadを使用したプログラミングを簡単に説明するとGPUプログラミングのような使い方になる。

Intelが公開している資料等は以下のところにある。

実際のプログラミングOpenMPの様にディレクティブを使ってXeon Phiへのデータ転送やプログラミングの指示を行う。

基本的な実行の流れは以下の通り。

  1. CPU側でXeon Phiで使用するデータの準備
  2. CPUからXeon Phiにデータの転送
  3. Xeon Phiで計算
  4. Xeon PhiからCPUに結果を転送
  5. CPUXeon Phiの結果を用いて計算

今回Xeon Phiを使ったプログラムの練習とテストのために行列積のプログラムを書いた。コードは以下のところ。

Offloadによる記述の説明

Xeon Phiで実行する部分の指示

練習で書いたサンプルプログラムの中の./src/matrix_multiplication.cppのMultiplyMatricesMicImple()という関数の中で使っている。

まず以下のようなディレクティブによってXeon Phiで実行する部分を記述

#pragma offload target(mic) 
{
// Xeon Phiで実行される部分
}

Xeon Phiが複数あり、実行するXeon Phiを指定する必要がある場合は以下の様にする

#pragma offload target(mic:0) 
{
// Xeon Phiで実行される部分
}

上の例ではmic0というXeon Phiが使われる。この0で指定した部分には変数も使える。サンプルでは試しに変数で指定している。

Xeon Phi、CPU間のデータ転送

データを転送する場合は以下のようにして行う。

#pragma offload target(mic) \
  in(matrix_a:length(column_a*row_a) alloc_if(1) free_if(1)) \
  in(matrix_b:length(column_b*row_b) alloc_if(1) free_if(1)) \
  out(matrix_c:length(column_a*row_b) alloc_if(1) free_if(1))
{
// Xeon Phiで実行される部分
}

matrix_aはCPU側のポインタで、転送する要素数をlength(n)で指定、今は行列なのでcolumn_a*row_a個ある。その次のalloc_if(1) free_if(1)の部分ではディレクティブの開始時にメモリを確保、終了時に開放するように指定。指定しない場合と同じだが、今回は明示的にやっている。

メモリを使いまわすために確保、開放しない場合は1の部分を0にする。

そしてinはディレクティブの開始時にCPUからXeon Phiにメモリのコピーを、outは終了時にXeon PhiからCPUのメモリに転送を行うときに指定。

上の例ではmatrix_aとmatrix_bのデータを最初にCPUからXeon Phiに転送し、計算が終わり次第、matrix_cの内容をXeon PhiからCPUに転送していることになる。

他にもいろいろできるのでこれについては以下のチュートリアル資料等を参考に。

http://software.intel.com/sites/default/files/Beginning%20Intel%20Xeon%20Phi%20Coprocessor%20Workshop%20Offload%20Compiling%20Part%201.pdf

(ALLOC、FREEなどがdefineされてるらしいが、使うとコンパイル時にエラーになってできなかった・・・。だれか教えて)

Xeon Phiで処理するコード

ディレクティブ以下のブロック内については普通にOpenMPなどが使える。その部分についての説明は省略。

(今回のサンプルは何も考えずに普通にOpenMPで並列化したのでもっと頑張れよ!といわれるかもしれないがそのうち頑張ります。)


Offloadで書いたプログラムを実行する際の環境変数周り

OpenMPスレッド数を指定する際に環境変数のOMP_NUM_THREASにスレッド数をセットするのが一般的だが、何もしないとCPUXeon Phiで同じ環境変数を見にいってしまう。このため、Xeon Phiだけ指定したい場合はMIC_ENV_PREFIXを使う。

export MIC_ENV_PREFIX=M_

としておくと、M_から始まる環境変数Xeon Phiで使われる。

なのでXeon PhiでOpenMPスレッド数を指定する場合は以下のようにする。

export MIC_ENV_PREFIX=M_

export M_OMP_NUM_THREAS=128

また、CPUXeon Phiのデータのやり取りやXeon Phiの実行の様子を出力するには以下のように環境変数をセットする。

export OFFLOAD_REPORT=1 #2や3でより詳細なものが出せる

実験

2つの正方行列のサイズを1024*1024,2048*2048,4096*4096にしてその積を計算するときの計算時間を測定。

(CPUBLASなどのライブラリを使っていませんのでその分も考慮してください)

実験環境
行列サイズCPU 1コア (sec.)Xeon Phi 236スレッド (sec.)高速化率
1024*10242.11.41.5
2048*204817.81.710.5
4096*4096728.43.5208.1

最大約200倍、CPUのソケットで比較すると理想的には約26倍くらいは出る計算。なので割りと速い気も・・・GPUと比べたらどうかは試してみないとよくわからないが。

頑張ってチューニングしたらどうなるかはきっとHPC屋さんがわんさか?やってるのでそちらを参照。

まとめ

どういうふうに書けばいいかわかれば割りとすぐに使える感じ。そこそこ速くなりそうなので今回は満足。

今後研究に使うかは謎。

他にもいろいろな制御ができるので公式のサンプルコードやチュートリアルを参考に。

情報のリンクのまとめは以下のところにある。

http://d.hatena.ne.jp/ang65/20130402/1364873045

2013-04-02

Xeon Phiのプログラミングについての情報(1) 資料のまとめ

12:24 | Xeon Phiのプログラミングについての情報(1) 資料のまとめを含むブックマーク Xeon Phiのプログラミングについての情報(1) 資料のまとめのブックマークコメント

研究室でXeon Phi コプロセッサー 5110Pを買ったので試しに使おうとしたときに見た資料等の場所についてのメモです。

練習用にコードも書いて、それの解説は以下のところにあります。

http://d.hatena.ne.jp/ang65/20130404/1365052622

公式の資料

ここにインストールやらチュートリアルやらがいろいろあります。

特によく見るものとして以下のもの

micの設定制御のツールの説明

サンプルプログラム

サンプルプログラムについては公式の資料の場所にも少しありますが、インストール時に入る以下のところにあるもののほうが参考になります。

  • /opt/intel/composer_xe_2013/Samples/

この中にXeon Phi以外のサンプルもありますが、Xeon Phiのものは以下のところにあります。

C++
  • /opt/intel/composer_xe_2013/Samples/en_US/C++/mic_samples/
  • /opt/intel/composer_xe_2013/Samples/ja_JP/C++/mic_samples/
Fortran
  • /opt/intel/composer_xe_2013/Samples/en_US/Fortran/mic_samples/
  • /opt/intel/composer_xe_2013/Samples/ja_JP/Fortran/mic_samples/

その他の資料

公式以外の資料で参考にしたものは以下のものです。

日本語でインストール手順などの説明があるところ
コンパイルや実行について書かれているところ
自作のサンプル

今回のテストで作ったサンプルはGitHubの以下のところにあげてあります。