読者です 読者をやめる 読者になる 読者になる

PSIPredをプロファイルして灯台の下を照らす困難なお仕事

PSIPred をビルドするよ!


生物学者でも研究室内で簡単に使用できるシーケンス検索アルゴリズムなんだぜというのが売りなのだが、こんなものが実際にビルドできる生物学者がいるとは思えない。そこで情報学科卒の出番がやってきた!configure -> make -> make install の流れは既知としようと思う。そいつぁすげぇや!


今回はプロファイルを取るというのが最終的な目標であります。プロファイルというのは、実行時にどの関数がどの位呼び出されてその実行にどの位の時間がかかりそれが全体の何%を占めるのか、という分析結果を見ることです。http://nenya.cis.ibaraki.ac.jp/TIPS/gprof.html とか見ながらやりました。プロファイルを知ったのは最近で、プログラミング作法を読み返していたら見つけた。

間隔をおいた木々の列が地平線を背にしてシルエットを浮かび上がらせていることもあれば、クローバーの荒地の上に熱い真昼がじっとたちこめていることもあり、霞む紺碧の空の遥かむこうに描き込まれたクロード・ロラン風の雲は、吸い込まれていくような無色の背景で積雲の部分だけを目立たせている。


Vladimir Vladimirovich Nabokov "Lolita" (若島正 訳)

というわけで、プロファイルを取ることで無駄な最適化を防ぐことができます。


PSIPred は、その筋では有名なアルゴリズムなのです。が、何せニッチな需要なだけに日本語情報がほとんど無い。PSIPred の動きは割と楽なんだけど、PSI-BLAST の方は論文も結構難しい。いつかエントリー書けたらいいなぁ、というか、書ける程度には詳しくならないとまずいなぁ(ホジホジ


さて、 PSIPred は内部で PSI-BLAST というツールを使っています。なので、まず PSI-BLAST をビルドせねばならないです。


情報のソースは ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.2.24/user_manual.pdf

【PSI-BLAST】1. ソースを取ってくる

ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/[version]/ncbi-blast-[version]+-src.tar.gz


間違っても blast 以下を全部ごそっととか取らないように。特に blast/db は巨大なので注意。一番新しい奴を落とせば良いような気がする。 blast+/LATEST というリンクが最新版に張られているはず。 zip もあるので zip が好きな人はそっちでも大丈夫。落としたら解凍。

【PSI-BLAST】2. ./configure

[src-dir]/c++/ の中に configure があるので configure します。 configure すると makefile ができます。が、その前にまず ./configure --help をします。というのも、今回はプロファイルを取るという重要な使命があるため、ビルドオプションを選ばないといけないからです。コンパイル時に -pg オプションを付けないといけません。


./configure --help を見るといろいろオプションがあるんですが、今回大事なのは --with-profiling オプションです。これを付けるとプロファイリング用にビルドされます。ソースをチラ見すると、中では -pg を付けているだけなので、最適化オプションを付けてはいけません。つまり、 --without-debug を付けてはいけないし、付ける場合は --without-optimization も付けないといけません。

$ ./configure --with-profiling

ってな感じになるかしら。

【PSI-BLAST】3. make, and make install

configure で致命的なエラーが出ているなら解決しましょう。たぶんライブラリが足りないとか、そんなエラーかなのではないかと。gcc のバージョンを上げないといけない可能性もある。


ビルドには結構時間かかるので注意。僕の PC では2、30分かかった気がする。

【PSI-BLAST】4. データベースの構築

データベースは ftp://ftp.ncbi.nlm.nih.gov/blast/db にあります。


PSI-BLAST はデータベースから与えられたクエリ配列に似た配列を見つけるアルゴリズムです。なので検索先のデータベースが必要です。が、ローカルに全部入れるのは趣味人以外はあまり意味が無いと思う。今回はアルゴリズムのボトルネックを確認するのが目的なので、何か手頃なデータベースが一つあればOK。中身は良く分からないので、サイズで選ぶと ftp://ftp.ncbi.nlm.nih.gov/blast/db/pataa.tar.gz あたりが手頃な気がするのでダウンロード。


どこでもいいので、落として解凍しておく。複数のデータベースを落とす場合はデータベースごとにフォルダに分けておくと吉。


以上で PSI-BLAST のビルドは終わりです。


PSI-BLAST 自体の使い方はよく知りません。たぶんすごくいろいろオプションがあります。 -help としてみたけど、やっぱりいろいろあった。でも論文を読む限り、基本はデフォルト値で使うのが良いんじゃないかなぁと思う。データベース、繰り返し回数、 E-Value の指定で良さそう。


続いてPSIPredのビルド作業。

【PSIPred】1. ソースを落とす

場所は http://bioinfadmin.cs.ucl.ac.uk/downloads/psipred/

【PSIPred】2. makefileを編集する

-pg オプションを付けないといけません。CFLAGS = -pg と書き加えます。

【PSIPred】3. make, and make install

../bin にインストールされます。


これで PSIPred もビルドできた!


ちなみに、 PSI-BLAST を含む BLAST パッケージは最近まとめ直しが行われて BLAST+ というパッケージに編成を変えたみたいです。と BLAST+ フォルダの README か何かに書いてありました。で、プログラム名なども少し変わったので、 PSIPred も BLAST+ 向けのが用意されています。それが BLAST+ ディレクトリに入っている実行ファイルです。実行ファイルと言っても、実行権限が付いた tcsh スクリプトです。実際に PSIPred を使う時は、この実行ファイルを編集してそれを実行するのが手っ取り早いと思います。

# The name of the BLAST+ data bank
set dbname = pataa/pataa

この部分を先にダウンロードして展開したデータベースを指定します。なお、この実行スクリプトはカレントディレクトリが BLAST+ であることを前提にして書かれているので注意しましょう。


ちなみに fedora のログインシェルは bash なので、僕の場合は tcsh を追加でインストールしないといけませんでした。似た境遇の人は適当に

# yum install tcsh

とかすれば tcsh が入るので、シェル入れてから実行してください。


実行は

$ ./runpsipredplus example.fasta

結果

で、プロファイルの結果ですが、実は PSIPred 自体は十分に早いということが分かりました。というのも、PSIPred を実行してみると PSI-BLAST の部分で時間がかかり、実際に PSIPred が走っている時間はほとんど1秒に満たない。ということは、次なる問題は PSI-BLAST のどの部分が最も時間を食っているのかということです。


そこで、先ほどの実行ファイル runpsipredplus の中身を読むと、実際に PSI-BLAST を呼び出している部分があります。

$ncbidir/psiblast -db $dbname
-query $tmproot.fasta
-inclusion_ethresh 0.001
-out_pssm $tmproot.chk
-num_iterations 3
-num_alignments 0 >& $tmproot.blast

つまり、これの $dbname と $tmproot を適度に置き換えて端末で実際に実行すればプロファイルが取れるというわけです。検索クエリの例として example.fasta が用意されているので、テストクエリはそれを使えばOKと思われます。実際は、

$ psiblast -db pataa
-query example.fasta
-inclusion_ethresh 0.001
-out_pssm example.chk
-num_iterations 3
-num_alignments 0 >& example.blast

こんな感じになるかと。

$ gprof gsiblast gmon.out > gmon.log

とすると、gmon.logに結果がリダイレクトされます。


で、プロファイルの結果が


f:id:makisyu:20101102235457p:image


つまり、というかやはり、ダブルヒットの発見とギャップ付きアライメントの作成で全実行時間の70%を占めているようです。他にもデータベースからシーケンスを取ってくる部分で20%ほど喰っているのでそこも気になるところです。後者はソースを見てみないとよく分からないけど、前者は「うん。まぁそうだよね。」という結果ですな。まぁそれでも確かめることができてよかった。


ちなみに、プロファイルの可視化には gprof2dot というのを使いました。これ便利ー。思ったより綺麗に出力してくれて素敵。重い部分は色が変わる辺りとか泣ける。d:id:ousttrue:20091017 とかが参考になりました。


さて、このアルゴリズム。いかに高速化すべきか。という研究を院で日々しているのです。そもそも PSIPred が BLAST を3回も回すからダメとか、ダブルヒットサーチが遅すぎるとか、あの GetSequence の20%は何なのかとか、ニューラルネットってもっと研究されてるし改善できるんじゃねとか、そんな感じです。


たまにはコンピュータネタで書いてみた。次は普通の話題で書きたいなぁ…。