- 追加された行はこのように表示されます。
- 削除された行は
このように表示されます。
{{category 系統解析}}
*MrBayesによるベイジアン系統推定
MrBayesによる系統樹のベイズ推定について簡単にまとめてみる。
*ベイズ法とは
私もよく分かってません(爆死)。最尤法とはまた別の推定方法ということで。ま、最尤法と同様、モデル選択は重要なことに変わりはないので気を付けましょう。東大の仲田さんが[MrBayesで用いられている方法の解説|http://www2.tba.t-com.ne.jp/nakada/takashi/bayes/idea.html]を書いて下さっていますので詳しく知りたい方はこちらをどうぞ。
*塩基置換モデルの選択
PAUP*を持っていない場合はMrAIC.pl + PHYMLという手があります。PAUP*をお持ちの方は[[複数遺伝子領域を用いた系統解析のための塩基置換モデル選択]]に掲載した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の項を参照のこと。[[複数遺伝子領域を用いた系統解析のための塩基置換モデル選択]]で紹介しているスクリプトは、データを読み込んでモデルの設定まで行うNEXUSコマンドファイルを出力しますので、これを1度Executeするだけですぐに解析できます。
MrBayesでは標準設定の場合、枝長を保存しますが、このMCMCの段階で
MrBayes > MCMC SaveBrLens=No
とすれば枝長保存を省くことができます。最尤法の初期系統樹決定に使うような場合は樹形さえ決められればよく枝長は不要なのでこのオプションを指定しておくと良いかもしれません。
Ver.3.1では、同時に2つ以上ののMCMCを走らせることで実行中にある程度MCMCの収束の程度が判断できるようになりました(「Average standard deviation of split frequencies (ASDSF)」として画面に情報が表示される。*.mcmcファイルにも出力される)。ただし、これは2つ以上のMCMCが独立したスレッドとして動作するという意味ではありませんので、SMPマシンにおける演算の高速化ができるわけではありません。SMPマシンやクラスタではMPI版のMrBayesを入れることでその能力を生かすことができます。また、それぞれで設定した世代数分のMCMCが走ります。世代数がMCMCのRun数分だけ分割されるわけではありません。
なお、MCMCの同時実行数は
MrBayes > MCMC NRuns=8
などとオプションを指定することで設定できます。前述のように標準は2つ同時実行です。2つ以上実行しないとASDSFの値が求められませんので必ず2以上にします。これは複数のMCMCが樹形空間内の同じ辺りに収束しているのかを表す値で、0.01未満かどうかを基準とします。どれかのMCMCが局所最適解に捕まっており、それ以外が実際の最適解に収束している場合には値が大きくなるワケですが、0.01を下回った後で再度増大するのは問題ありません。一度0.01を下回れば十分です。最適解近傍の第n最適解もそれなりに評価される方がより正当なクレード支持率を得られるわけですから収束してから再度別解にどちらかのMCMCが捕まるのは妥当と言えます。ただし、真の最適解へ収束する前に、なかなか脱出できない局所最適解へ収束してしまっている場合はこの限りではありません。そのため、ASDSFと同時に尤度の値も収束していることを確認すべきです。
MCMCが終了したら、まず*.mcmcを開き、StdDev(s)の項(ASDSFの値)を見て、何世代目で収束しているのかを探します。通常、各行の最後の値がそれです。0.01未満になっている世代を探します。ただ、データセットによってはそもそも収束しない場合もありますのでご注意下さい。何度もやり直せばたまたま局所最適解に捕まらずに収束することもありますが、その局所最適解も考慮すべき樹形である場合にはその結果は正当なものではありません(特に事後確率)。MCMCが走り終わると、続行するかどうかの質問が出ますので、その質問に答える前に、収束後の世代数が十分かどうかを確認し、もし足りなければ不足する世代数分だけ続けるように指定します。
5,000世代目(=100世代毎にサンプリングする標準設定なら51本目の系統樹)で収束しているなら、
MrBayes > SumT BurnIn=51
とすると、初めの51本(NRuns=2なら各51x2=102本)の系統樹は使わずに*.con(多数決合意樹にベイズ事後確率がマッピングされた樹形ファイル)などのファイルが出力される。筆者はMCMCが収束してから10,000本(デフォルトのNRuns=2 SampleFreq=100なら50万世代)の系統樹を用いて多数決合意樹と事後確率を得るようにしています。
MrBayes > SumP BurnIn=51
の結果もよく見ましょう(特に尤度プロットとPSRF)。
MrBayes > Log Start Filename=hogehoge
MrBayes > SumT
MrBayes > SumP
MrBayes > Log Stop
などとしてメッセージをログファイルに保存しておくとよい。
また、生成されるファイル群の中身をよく確認しておきましょう。特にASDSFの値も書き込まれる*.mcmcファイルは、各パラメータと系列交換試行のAcceptance rateが出力されていますので、これらの値が収束世代以降に極端に低い値にはなっていないことを確認します。Acceptance rateがあまりに低い場合は、温度間隔(Temp)を狭めたりサンプリング間隔(SampleFreq)を大きめにしたりして再解析するなどの対処をとります。
*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 SwapAdjacent=Yes NSwaps=19
としています。これで温度レンジも広がります。SwapAdjacentは温度序列的に隣接するMCMC系列間でのみ交換試行を行うようにするオプションです。NSwaps=19なのは、20系列のMCMCの隣り合う組み合わせの全て(=19通り)で、1世代毎に平均1回の状態交換試行が行われるようにするためです。NChainsをいくら増やしても交換試行数が少ないと温度の高い系列の状態がサンプリングするcold chainへなかなか伝わらないためNSwapsを増やす必要があるわけです。ただ、下記のバグによりSwapAdjacentとNSwapsは設定しても反映されませんので、ソースコードを修正してコンパイルしたバイナリを用意する必要があります。修正方法は本ページ下部のリンクにあるdiffファイルをご参照下さい。有効化済みWin32バイナリも配布しています。
温度序列的に隣接する系列間での状態交換成立確率を高めるために温度間隔を狭くしている場合がありますが、これは確かに確率は上がるものの、MCMC系列の数はそのままだと樹形探索能力とパラメータ値最適化能力の低下を招くため望ましくありません。要するに局所最適解に捕まりやすくなってしまうということです。また、系列数を増やした際に系列間の交換試行も増やさなければ、実際には交換成立の絶対数が増加しません(成立確率も上がりません)。むしろ、系列の全組み合わせ数当たりの交換成立数が少なくなり、その効果で最適化の進行が遅くなる(より低温の系列への状態伝播が遅くなる)ため、実際には収束していなくても収束しているように見えてしまいやすくなります。というわけで、温度間隔を狭くするなら必ず系列数も増やすべきで、系列数を増やしたら状態交換試行数も増やすべきだということです。状態交換試行は1世代に各系列間で平均1回を私は目安にしています。
*SumTが実行できないとき
UNIX・MacOS X上ではSumTがしばしばできなくなることがあります。ファイルを読み込み終わった途端に落ちるという症状です。これはNRunsおよびNTreesオプションを指定して実行することで改善することがあるようです。しかし、やっぱりどうしようもないことがあることを確認しています。このような場合、Windows版MrBayesで実行するとうまくいきます。Windowsマシンが無い場合は、枝長が得られませんがPAUP*に読み込んでConTreeするしか無さそうです。
別の問題として、Separate modelによる解析を行った場合にSumTできない問題があります。NTreesを指定すれば各領域で別々にSumTすることはできますが、全領域の合計もしくは平均枝長の系統樹のSumT結果を得るようなことはできません。そこで、一旦平均or合計枝長の系統樹を作成してからSumTするためのスクリプトを用意しました。[[同樹形・不等枝長の複数系統樹から平均枝長の系統樹を作成する]]をご参照下さい。
*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へと状態が伝わるようにすることができます。よって、デフォルト設定よりも高速に収束状態へと向かうものと思われます。SwapDescendAllオプションの値にはYesかNoを与えます。何度も言うようですが無保証ですので自己の責任において利用して下さい。ソースコードの差分はページ下部のdiffファイルに含まれています。
*MrBayesのインストール
Windows / MacOS用にはコンパイル済のバイナリがあるのでそちらを利用する。Windowsではmb.exeを標準でPATHの通っているC:\Windows\system32に置いておけばいいでしょう。PAUP*もwin-paup4b10-console.exeをpaup.exeにリネームして同じ場所に入れておきましょう。
**通常版のコンパイル・インストール(UNIX)
% tar xvzf mrbayes-3.1.2.tar.gz
% cd mrbayes-3.1.2
% make
% su
# cp ./mb インストール先
# exit
Makefileを適当に編集すること。最適化オプションはCPU最適化(-march=hogehoge)はともかく、あまり危険なオプションは付けない方が良いでしょう。私は以下のオプションを指定しています。
OPTFLAGS = -O3 -march=athlon-xp -pipe -fforce-addr -funroll-loops \
-frerun-loop-opt -falign-functions=4
**通常版のコンパイル・インストール(Windows)
Cygwin・MinGW環境では基本的にUNIX環境と同じです。USEREADLINEをnoにした上で、Cygwinでは-mno-cygwinをOPTFLAGSに加える程度です。以下はVisual C++の場合です。
まず、.tar.gzの展開できるアーカイバでソースの配布ファイルを展開します。Visual C++用のプロジェクトファイルが含まれていますのでそれを使われても構いませんが、私は無償配布されていたVisual C++ Toolkit 2003(現在は配布終了)を使いました。VCToolkit2003では、プロジェクトファイルが扱えないので、まずは[Platform SDK|http://www.microsoft.com/downloads/details.aspx?FamilyId=0BAF2B35-C656-4969-ACE8-E4C0C0716ADB&displaylang=en]を入れてMakefileを扱うためのnmakeを入手します。これにはMrBayesのWindows版コンパイルに必要なヘッダも含まれています。インストールしたら、環境変数PATH・INCLUDE・LIBにPlatform SDKとVCToolkit2003のBin・Include・Libをそれぞれ追加します(Toolkitのコマンドプロンプト起動に使われるバッチファイルを編集するのが楽)。そして、ソースを展開したディレクトリに以下の内容のファイルをメモ帳などのエディタで作成します。ファイル名はとりあえずmakemrbayesとしましょう。
CC = cl
CFLAGS = /O2 /Ox /TC /ML /DWIN_VERSION
all: mrbayes.exe
mrbayes.exe: mb.obj bayes.obj command.obj mbmath.obj \
mcmc.obj model.obj plot.obj sump.obj sumt.obj
link /OPT:REF /OUT:mrbayes.exe *.obj
mb.obj: mb.c
bayes.obj: bayes.c
command.obj: command.c
mbmath.obj: mbmath.c
mcmc.obj: mcmc.c
model.obj: model.c
plot.obj: plot.c
sump.obj: sump.c
sumt.obj: sumt.c
ファイルを用意したら、ソースのあるディレクトリをカレントにして
nmake -f makemrbayes
とするとコンパイルが実行されてmrbayes.exeができあがります。CFLAGSには、Pentium Pro以降なら/G6、Pentium 4・Athlon以降なら/G7などの最適化オプションを付けると高速化できます。さらにSSEを使うなら/arch:SSEを、SSE2なら/arch:SSE2を付けると良いようです。
**MPI並列化版をインストールする場合(UNIX)
まずはMPIの実装を用意します。UNIXでは主に以下の3つが考えられます。
-[LAM/MPI|http://www.lam-mpi.org/]
-[MPICH2|http://www-unix.mcs.anl.gov/mpi/mpich/]
-[Open MPI|http://www.open-mpi.org/]
MPICHはリファレンス実装で安定していますが遅いです。LAMが速度・安定性の両面で優れています。Open MPIはLAMの後継です。MrBayesはMPI-2を使っていませんので特に必要はありません。好きなように選んでいただいて構いませんが個人的にはLAMがお薦めです。
あとはMrBayesのMakefileでMPI=yesを指定してコンパイルすればよい。ICCなどのGCC以外のコンパイラでコンパイルしたい場合、LAMならばmake時に環境変数LAMHCCにコンパイラの実行ファイルを指定する。環境によってmpiccにPATHが通っていなかったりするので注意。
LAMの場合は事前に
lamboot -v
としてデーモンを起動しておく。演算終了後は
lamhalt
にてデーモン停止。この辺りの実行環境関連は実装によって異なるので各実装のマニュアルを参照。
実行環境が整ったら、MPI版MrBayesをmpirunから以下のように起動する。
% mpirun -np 使用するCPU数 mb
mpirunの起動オプションは
% mpirun --help
で表示される。
**MPI並列化版をインストールする(Windows)
現在のところ当方の環境では正常に動作していません。ご注意下さい。
MPIの実装には以下の2つのどちらかを使います。
-[MPICH2|http://www-unix.mcs.anl.gov/mpi/mpich/]
-[DeinoMPI|http://mpi.deino.net/]
DeinoMPIはMPICH2の改造版なので元は同じです。
インストールしたら、Visual C++ならPATH・INCLUDE・LIBにMPI実装のインストール先にあるBin・Include・Libを加えます。そして、通常版と同様に以下のようなmakemrbayesを用意します。
CC = cl
CFLAGS = /O2 /Ox /TC /MT /DMPI_ENABLED /DWIN_VERSION
all: mrbayes.exe
mrbayes.exe: mb.obj bayes.obj command.obj mbmath.obj \
mcmc.obj model.obj plot.obj sump.obj sumt.obj
link /OPT:REF /OUT:mrbayes.exe *.obj mpi.lib
mb.obj: mb.c
bayes.obj: bayes.c
command.obj: command.c
mbmath.obj: mbmath.c
mcmc.obj: mcmc.c
model.obj: model.c
plot.obj: plot.c
sump.obj: sump.c
sumt.obj: sumt.c
それ以外は通常版と同様にコンパイルします。
Cygwin・MinGWではMakefileの
CC = mpicc
を
CC = gcc
に書き換え、CFLAGSに-I/hogehoge/Includeを、LIBSに/hogehoge/Lib/mpi.libを加えるなどしてmakeします。
DeinoMPIでは管理ツールをスタートメニューから起動して、「Crudential Store」タブ→「enable create store options」にチェックを入れる→「Password」の設定(クラスタを組むわけではないならNo password)→「Encryption」の設定(クラスタを組むわけではないならno encryption)→「Location」の設定(とりあえずRegistry)→「Create Credential」の順でCredentialを作成。詳しくは[マニュアル|http://mpi.deino.net/manual.htm#_Toc147560663]を参照。
MPICH2では設定ツールをスタートメニューから起動するとアカウントとパスワードを訊かれるので適当に入力してRegistすると設定画面が開きます(以降はアカウントは訊かれません)。channelの設定などを[マニュアル|http://www-unix.mcs.anl.gov/mpi/mpich/downloads/mpich2-doc-windev.pdf](pdf)を見つつ決めて下さい(分からなければデフォルトで)。
MPI並列化版の実行ファイルは、MPI実装のインストール先にあるmpiexec.exe経由で起動します。クラスタのノードにするわけでなくマルチCPUを使うだけなら、MPI実装のサービスを起動した状態(インストールすると自動的に起動しています)で、
mpiexec -localonly プロセス数 実行ファイル
として起動します。
**本ページ添付のdiffファイルを利用する
何度も言うようですが内容は保証できませんので十分に注意して使って下さい。ソースコードの差分はページ下部のリンクから取得できます。Windows環境の方は各種CPUに最適化したコンパイル済バイナリがzipファイルに入っていますのでそれを利用できます(nmake用のmakefileも入っています)。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 < ../mrbayes-3.1.2.mod.diff
これで改造版をコンパイルできるようになりますのであとは上記の手順の通りです。ただし、MPI関連は特に気を遣っていませんのでエンバグしている可能性がかなりあると思います。
*参考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 法による系統解析」にて日本語で丁寧に説明されています。
[分子系統学演習 - データセットの作成から仮説検定まで|http://www.fifthdimension.jp/documents/molphytextbook/]にまとめ直しました。リンク先をご覧下さい。