{{category 系統解析}} *MrBayesによるベイジアン系統推定 　MrBayesによる系統樹のベイズ推定について簡単にまとめてみる。 *ベイズ法とは 　私もよく分かってません(爆死)。最尤法とはまた別の推定方法ということで。ま、最尤法と同様、モデル選択は重要なことに変わりはないので気を付けましょう。東大の仲田さんが[MrBayesで用いられている方法の解説|http://www2.tba.t-com.ne.jp/nakada/takashi/bayes/idea.html]を書いて下さっていますので詳しく知りたい方はこちらをどうぞ。 *塩基置換モデルの選択 　PAUP*を持っていない場合はMrAIC.pl + PHYMLという手があります。PAUP*をお持ちの方は[[Site-Specific Rateなモデルも含めたモデル選択]]に掲載したPerlスクリプトでモデル選択を行うことが出来ます。アミノ酸配列には[ProtTest|http://darwin.uvigo.es/software/prottest.html]か[ModelGenerator|http://bioinf.may.ie/software/modelgenerator/]をどうぞ。 *解析の実行 　当たり前ですが、以下のものが必要です。 -MrBayes (以下ではVer.3.1を使用) 　以下、インストール済として説明します。 % mb MrBayes > Execute data.nexus MrBayes > Execute preference.nexus などとしてデータ(自分で用意する)と設定(モデル選択用ソフトが出力するものなど)を読み込んだ後、 MrBayes > MCMC とするとMCMCが標準設定(MCMC100世代毎にファイルへ出力)で開始される。オプションについては[MrBayesのコマンドリファレンス|http://mrbayes.csit.fsu.edu/commref_mb3.1.pdf]にあるmcmcの項を参照のこと。 　MrBayesでは標準設定の場合、枝長を保存するが、このmcmcの段階で MrBayes > MCMC SaveBrLens=No とすれば枝長保存を省くことができる。最尤法の初期系統樹決定に使うような場合は樹形さえ決められればよく枝長は不要なのでこのオプションを指定しておくと良いかもしれません。 　Ver.3.1では、同時に2つのMCMCを走らせることで実行中にある程度MCMCの収束の程度が判断できるようになった(「Average standard deviation of split frequencies」として画面に情報が表示される。*.mcmcファイルにも出力される)。ただし、これは2つ以上のMCMCが独立したスレッドとして動作するという意味ではありませんので、SMPマシンにおける演算の高速化ができるわけではありません。SMPマシンやクラスタではMPI版のMrBayesを入れることでその能力を生かすことができます。また、それぞれで設定した世代数分のMCMCが走ります。世代数がMCMC数分だけ分割されるわけではありません。 　なお、MCMCの同時実行数は MrBayes > MCMC NRuns=8 などとオプションを指定することで設定できます。前述のように標準は2つ同時実行です。2つ以上実行しないと「Average standard deviation of split frequencies」の値が求められませんので必ず2以上にします。これは複数のMCMCが樹形空間内の同じ辺りに収束しているのかを表す値で、0.01未満かどうかを基準とします。どれかのMCMCが局所最適解に捕まっており、それ以外が実際の最適解に収束している場合には値が大きくなるワケです。 　MCMCが終了したら、まずdata.nexus.mcmcを開き、StdDev(s)の項(「Average standard deviation of split frequencies」の値)を見て、何世代目で収束しているのかを探します。通常、各行の最後の値がそれです。それ以降0.01未満になっている世代を探します。ただ、データセットによってはそもそも収束しない場合もありますのでご注意下さい。何度もやり直せばたまたま局所最適解に捕まらずに収束することもありますが、その局所最適解も考慮すべき樹形である場合にはその結果は正当なものではありません(特に事後確率)。つまり、MCMCは複数回繰り返すか、NRunsの値を始めから大きな値にするべきです(私は8以上にします)。MCMCが走り終わると、続行するかどうかの質問が出ますので、その質問に答える前に、収束後の世代数が十分かどうかを確認し、もし足りなければ不足する世代数分だけ続けるように指定します。 　5,000世代目(=100世代毎にサンプリングする標準設定なら51本目の系統樹)で収束しているなら、 MrBayes > SumT BurnIn=51 とすると、初めの51本(NRuns=2なら51x2=102本)の系統樹は使わずにdata.nexus.con(多数決合意樹にベイズ事後確率がマッピングされた樹形ファイル)などのファイルが出力される。筆者は100世代毎にサンプリングする設定で、MCMCが収束してから10,000本(NRuns=2なら50万世代)の系統樹を用いて多数決合意樹と事後確率を得るようにしている。 MrBayes > SumP BurnIn=51 の結果もよく見る。また、data.nexus.*ファイル群を適当にテキストエディタで中身を見るべし。 MrBayes > Log Start Filename=hogehoge MrBayes > SumT MrBayes > SumP MrBayes > Log Stop などとしてメッセージをログファイルに保存しておくとよい。 *MCMCのパラメータについて 　MrBayesでは、温度(=乱数の乱れの大きさ)の異なる複数のMCMCを走らせ(パラレルテンパリング法・Metropolis-Coupled MCMC・Multi-Chain MCMC・MC3)、温度の低いcold chainから系統樹をサンプリングします。同じようなことをNRunsの説明でも書きましたが、前述の「複数のMCMC(=run)を走らせる」というのは「温度が少しずつ異なっている複数のMCMC(=chain)」をさらに複数走らせる、という意味です。この温度の異なるMCMC間では、一定のルールに従って状態交換が行われます。この温度の少しずつ異なるMCMCの数を決めるのがNChainsの値でありデフォルト値は4で、各MCMC間の温度間隔を決めるのがTempでデフォルト値は0.2です。つまり、デフォルトでは0.2度ずつ異なる4つのMCMC(温度レンジは0.8度)が2つ走っている(NRuns=2)わけです。状態交換を何世代毎に行うかはSwapFreqの値で決まっておりデフォルトでは1で、その際に各runで何回状態交換試行を行うかはNSwapsで決まっておりこれもデフォルトでは1です。要するに1世代毎に、1回の状態交換試行が行われます。しかし、パラレルテンパリングではなるべく温度間隔を狭くして沢山のMCMCを回し、頻繁に状態交換が行われるようにした方が良いということで、 MrBayes > MCMCP NRuns=8 NChains=20 Temp=0.05 NSwaps=190 としています。これで温度レンジも広がります。NSwaps=190なのは、20個のMCMCの全組合せ(=190通り)で、1世代毎に平均1回の状態交換試行が行われるようにするためです。NChainsをいくら増やしても交換試行数が少ないと温度の高いchainの状態がサンプリングするcold chainへなかなか伝わらないためNSwapsを増やす必要があるのです。ただ、下記のバグによりNSwapsは設定しても反映されませんので、ソースコードを修正してコンパイルしたバイナリを用意する必要があります。 　温度序列的に隣接するchain間での状態交換成立確率を高めるために温度間隔を狭くしている場合がありますが、これは確かに確率は上がるものの樹形探索能力とパラメータ値最適化能力の低下を招くため望ましくありません。要するに局所最適解に捕まりやすくなってしまうということです。また、chain数を増やした際にchain間の交換試行も増やさなければ、実際には交換成立の絶対数が増加しません(成立確率も上がりません)。むしろ、chainの全組み合わせ数当たりの交換成立数が少なくなり、その効果で最適化の進行が遅くなる(より低温のchainへの状態伝播が遅くなる)ため、実際には収束していなくても収束しているように見えてしまいやすくなります。というわけで、温度間隔を狭くするなら必ずchain数も増やすべきで、chain数を増やしたら状態交換試行数も増やすべきだということです。状態交換試行は1世代に各chain間で平均1回を目安にして下さい。ソースコードを見ると、MrBayesには温度序列で隣接するchain間でのみ状態交換試行を行うSwapAdjacentというオプションがありますが、実際には使うことができません(設定しても有効にならない)。これの修正方法は本ページ下部のリンクにあるdiffファイルをご参照下さい。有効化済みWin32バイナリも配布しています。 *Ver.3.1のバグについて 　最新リリース版のVer.3.1で、NRuns=1を指定せずにMCMC演算中にCtrl + Cで中止すると、*.tファイルにend;が複数出力されてしまい、SumTをしようとすると文法エラーになります(私の環境でしか起きないのかもしれませんが)。これには、*.tファイルをエディタで開いて、末尾のend;を1つ残して他は削除してしまえばよい。 　Windows版では上記の不具合は起きないことを確認しました。環境によって変わるようです。 　また、MCMCの設定パラメータのNSwapsの値を変えても有効にならないという現象もあります。Linux上でGCC 3.4.5およびICC 8.1でコンパイルしたバイナリ、配布されているコンパイル済みのWindows版バイナリのいずれにおいても確認しています。これは、MrBayesのソースコードのバグで、mcmc.cの33,547行目でfor構文を使ってNSwapsの値だけ繰り返し処理するはずが、構文ミスによりこの繰り返し処理が働かないためです。33,548行目から33,557行目までを囲うように{}を追加することでfor構文が有効になり、NSwapsの値が正常に反映されるようになります。修正用のdiffをこのページに添付してありますので下部のリンクから参照できます。また、修正済みソースを用いてMinGWのgccで各種CPUに最適化してコンパイルしたWindows用バイナリも置いておきましたのでご自由にどうぞ。ただし最適化をきつめにかけていますので計算結果は保証できません。予めご了承下さい。 *MCMCで最尤系統推定 　ご存じのように、MrBayesはMCMCを使って樹形とパラメータの最適化を行っています。それに対し、PAUP*などの最尤系統推定ソフトウェアはNewton-Raphson法などの数値計算法を用いて検討している樹形における各パラメータの最適化を行います。このNewton-Raphson法による最適化は非常に計算エネルギーのかかる反復改善法で、MCMCの方が遙かに計算は楽です(その代わり、良い初期値を与えられた場合のNewton-Raphson法ほどの最適化能力は無いでしょう)。また、樹形の変更はNearest Neighbor InterchangeやSubtree Pruning and Regrafting、Tree Bisection and Reconnectionを利用しますが、MrBayesのMCMCにはNNIとTBRが組み込まれています。ここで、MrBayesのMCMCは尤度を規準としていますので、最尤法におけるパラメータの最適化＋樹形探索のアルゴリズムとしてMCMCを使うことは理論上可能です。つまり、MrBayesを高速に広い範囲の樹形探索を行う最尤系統推定ソフトウェアとして用いることができるということです。しかも、PAUP*では扱えない複雑な分子進化モデルに基づいた尤度計算が可能ということもメリットです(PAUP*で扱えるのにMrBayesでは扱えないモデルもありますが)。 　MrBayesで最尤系統推定を行うことはさほど難しいことではなく、cold chainからのサンプリングを1世代に1回にし、良いheated chainの状態が高速にcold chainに伝わるように状態交換試行を増やし、収束したら最初から計算することを何度も繰り返すだけです(MCMCのStopRuleおよびStopValオプションを使うと便利です)。全chainから系統樹をサンプリングしてもいいのですが、MrBayesにはそのような設定は無いようです。また、ランダムに生成した系統樹を使う方が良いとは思いますが、ratchet法的に生成した初期系統樹を与えてやることでratchet法も可能です。そうして得られた系統樹群から、ベストな尤度を持つ系統樹を最尤系統樹として採用するわけです。 　ただしこの方法では得られた樹形の信頼性を調べるには別途ブートストラップ解析などを行う必要があります。Bayes法の利点は信頼性の値(事後確率)の計算まで含めての高速性にあるわけです。その信頼性値が信頼できるかどうかはさておいて。しかし最尤法でも、樹形としての上位いくらかをリストアップし座位毎の尤度を計算すれば、CONSELを用いてそこそこ高速に、近似的に不偏なブートストラップ確率(Approximately Unbiased Bootstrap Probability＝AUBP)を計算することができるようになりました。上位の系統樹をどれくらい使って比較すべきかという点が難しいですが、重要な対立仮説間の比較ができているかどうか、という観点からそれは決められるべきでしょう。「重要な対立仮説」がどれくらいの範囲かは場合によって異なるので一意に決めることは難しいと思います。 *MrBayesの系列交換試行を拡張する 　上記のように、ここで配布している改造版を使うことでMCMC/MCMCPのNSwapsオプションとSwapAdjacentオプションが利用できます。NSwapsは正の整数、SwapAdjacentはYesかNoを指定します。改造版ではさらに、SwapDescendAllオプションを追加し、これによって高温chainから降順に全隣接chain間の交換試行をSwapFreqに設定した世代数毎に行うようにできるようにしました。これにより、高温chainが尤度の良好な状態になった際に高速にcold chainへと状態が伝わるようにすることができます。よって、デフォルト設定よりも高速に収束状態へと向かうものと思われます。ただし事後確率も高くなりやすいはずなので最尤系統樹探索用の設定だと考えた方が良いかもしれません。ベイズ推定にはNSwapsを増やす以外はデフォルト設定の方が良いでしょう。SwapDescendAllオプションの値にはYesかNoを与えます。何度も言うようですが無保証ですので自己の責任において利用して下さい。ソースコードの差分はページ下部のdiffファイルに含まれています。 *番外編 **MrBayesのインストール 　Windows / MacOS用にはコンパイル済のバイナリがあるのでそちらを利用する。Windowsではmb.exeを標準でPATHの通っているC:\Windows\system32に置いておけばいいでしょう。PAUP*もwin-paup4b10-console.exeをpaup.exeにリネームして同じ場所に入れておきましょう。 　その他の環境ではソースコードをダウンロードしてインストールする。 % tar xvf mrbayes-3.1_src.tar % cd mrbayes-3.1 % make % su # cp mb インストール先 # exit 　Makefileを適当に編集すること。最適化オプションはCPU最適化(-march=hogehoge)はともかく、あまり危険なオプションは付けない方が良いでしょう。私は以下のオプションを指定しています。 OPTFLAGS = -O3 -march=athlon-xp -pipe -fforce-addr -funroll-loops \ -frerun-cse-after-loop -frerun-loop-opt -falign-functions=4 -m3dnow \ -mfpmath=sse -msse -mmmx ***MPI並列化版をインストールする場合 　まずはMPIの実装を用意します。とりあえずここでは[MPICH|http://www-unix.mcs.anl.gov/mpi/mpich/]ということで書きます。 % tar xvzf mpich.tar.gz % cd mpich-* % ./configure --prefix=インストール先 --with-arch=hogehoge --with-device=ch_shmem % make % su # make install # exit などとする(これは1台のSMPマシン用)。FreeBSDの場合はportsからmakeするだけ。ただしportsの標準ではクラスタ用にmakeされるので、1台のSMPマシンで動作させる場合は「--with-device=ch_shmem」オプションを付けてconfigureするようにMakefileを書き換える(「--prefix=/usr/local/mpich」の後ろに書き加えればいい)。 　クラスタ用にコンパイルするには % ./configure --prefix=インストール先 --with-arch=hogehoge --with-device=ch_p4 --with-comm=shared などとする。クラスタを構成する各ノードがシングルコアシングルプロセッサマシンの場合は「--with-comm=shared」は必要無い。インストール先/shared/machines.*の編集も忘れずに行うこと。その他各種オプションは % ./configure --help で確認できる。 　あとはMrBayesのMakefileでMPI=yesを指定してコンパイルすればよい。FreeBSDの場合、portsからインストールするとmpiccがPATHの通った場所にインストールされないためフルパス指定すること。 　MPI版MrBayesはmpirunから以下のように起動する。 % mpirun -np 使用するCPU数 mb 　mpirunの起動オプションは % mpirun --help で表示される。 ***本ページ添付のdiffファイルを利用する 　Windows環境用バイナリとソースコードの差分をページ下部のリンクから取得できます。Windows環境の方は各種CPUに最適化したコンパイル済バイナリがzipファイルに入っていますのでそれを利用できます。Visual C++かCygwin・MinGW環境を用意すればコンパイルも可能です。自前でビルドする方は以下のようにdiffを当ててコンパイルして下さい。UNIX互換環境が前提です。Windowsでも.tar.gzが展開できるアーカイバとpatchコマンドがあればCygwinやMinGWが無くてもいけます。 　ソースコードの配布ファイルとdiffファイルをダウンロードしたら、作業ディレクトリに置きます。ソースコードの配布ファイルは作業ディレクトリで展開します。 tar xzf mrbayes-3.1.2.tar.gz 　展開したらディレクトリが生成されるのでディレクトリを移動します。 cd mrbayes-3.1.2 　patchコマンドでdiffファイルを当てます。 patch < ../mrbayes312.diff 　これで改造版をコンパイルできるようになりますのでMakefileを適当に書き換えてから make とすると実行ファイルができあがります。 *参考URL -http://mrbayes.csit.fsu.edu/mb3.1_manual.pdf 　公式マニュアルのPDF。本格的に研究に利用する方は読みましょう。とりあえず少し試してみたいという方も末尾のチャート図は見ておきましょう。 -http://mrbayes.csit.fsu.edu/Help/ 　公式HTML Help。 -http://mrbayes.csit.fsu.edu/commref_mb3.1.pdf 　公式コマンドリファレンス。 -http://www2.tba.t-com.ne.jp/nakada/takashi/index.html 　東大の仲田さんのサイト「きまぐれ生物学」。「Bayes 法による系統解析」にて日本語で丁寧に説明されています。