| トップへ | 生物の起源 | 雑記 | 新刊紹介 | 地質年代表 | 生物分類表 |
| Genera Cellulatum | 系統解析 | リンク集 | プロフィール | 掲示板 | 系統解析トップへ |
目次 |
この解説は MrBayes 3.1（最新版は 3.1.2）を使い始めるにあたって，マニュアル・ 参考書などを整理しながら書いています。 従って，細かいところで今後修正する箇所も出てくる可能性があることをご了承下さい。 特に世代数や burnin の値については他の研究者らがどのような値を用いるのか， 今後出てくる MrBayes 3.1 を用いた系統解析の論文なども参照したいところなので， このプログラムを実際に使ってみる場合にはご注意下さい。
MrBayes 3.1.2 には，当然ながらここで紹介した以外にも様々な機能や設定がありますので， 必要な場合には直接参考書や説明書をお読み下さい。
ここでは Windows XP，Vista の環境で系統解析を行うことを前提とした方法を書いていますが， Mac OS の場合でも OS に依存する操作の違いはありますが，大きな違いはないと思われます。 ただし，MrModeltest の操作については別のプログラムを用いますので，ソフトの説明文を参照してください。
MrBayes の使用法については，MrBayes 3.1 用のマニュアルも公開されており， また教科書（Hall, 2004; MrBayes 3.0b4 に準拠）もあるので，これらを参照しました。 余裕のある方は，原典に当たって見ることをお勧めします。
なお、MrBayes 3.1.1、MrBayes 3.1.2 に先立った公開された MrBayes 3.1 については、 プログラムの動作が従来のものに比べて非常に（4 倍程度）遅いため、 この欠陥が修正された MrBayes 3.1.2 を用いることを強くお勧めします。 なお、MrBayes 3.1.2 においては MrBayes 3.1.1 にて報告されたいくつかのバグが修正されています。 しかし、本解説で扱った範囲においては違いはありません。
以下のプログラムをダウンロードしてインストールする。（PAUP* を除き、全て無料でダウンロードできる）
作業用のフォルダを作成する（後々便利なので）。ここでは "Gene1" フォルダとする。
ClustalX などでアラインメントを作成し、NEXUS 形式で保存する。（ここでは Gene1.nxs とする。 ファイル名を入力することがあるので、シンプルな名称の方が便利である。またハイフン "-" などがファイル名に入っていると、 MrBayes がファイルを読み込めないため注意）
注：OTU（operational taxonomic unit，操作上の分類単位）の名称は、全て半角で英数字とアンダーバー "_" しか使わないようにすること。特にハイフン "-" は ClustalX では使えるが、MrBayes で使えないので要注意。
メモ帳などのテキストエディタを用いて、ファイルのトップにあるデータ定義の部分を MrBayes
で読み込める形式に編集する。
ClustalX で作成した NEXUS ファイルの出だしは以下のようになっている。
#NEXUS
BEGIN DATA;
dimensions ntax=32 nchar=1392;
format missing=?
symbols="ABCDEFGHIKLMNPQRSTUVWXYZ"
interleave datatype=DNA gap= -;
matrix
Chloromonas_rosae_UTEX1337 ATTTGGGATCCGCACTTTGGCCAACCAGCAGTTGAAGCATTTACACGTGG
Chloromonas_serbinowi ATCTGGGATCCTCATTTTGGTCAACCTGCTGTAGAAGCGTTTACACGTGG
（以下略）
この内 format の部分（青字）は MrBayes の形式に合わないため、編集する必要がある。
具体的にはsymbols=
の部分が不要である点と、
interleave
をinterleave=yes
と変更する必要がある点が重要である。
これらを修正すると、以下のようになる。
#NEXUS
BEGIN DATA;
dimensions ntax=32 nchar=1392;
format missing=? interleave=yes datatype=DNA gap= -;
matrix
Chloromonas_rosae_UTEX1337 ATTTGGGATCCGCACTTTGGCCAACCAGCAGTTGAAGCATTTACACGTGG
Chloromonas_serbinowi ATCTGGGATCCTCATTTTGGTCAACCTGCTGTAGAAGCGTTTACACGTGG
（以下略）
#NEXUS
BEGIN DATA;
dimensions ntax=32 nchar=376;
format missing=?
symbols="ABCDEFGHIKLMNPQRSTUVWXYZ"
interleave datatype=PROTEIN gap= -;
matrix
Cg_elongatum AGFKAGVKDYRLTYYTPDYVVRDTDVLAAFRMTPQPGVPPEECGAAVAAE
Cg_euchlorum AGFKAGVKDYRLTYYTPDYVVRDTDVLAAFRMTPQPGVPPEECGAAVAAE
（以下略）
この内 format の部分（青字）は MrBayes の形式に合わないため、編集する必要がある。
具体的にはsymbols=
の部分が不要である点と、
interleave
をinterleave=yes
と変更する必要がある点が重要である。
これらを修正すると、以下のようになる。
#NEXUS
BEGIN DATA;
dimensions ntax=32 nchar=376;
format missing=? interleave=yes datatype=PROTEIN gap= -;
matrix
Cg_elongatum AGFKAGVKDYRLTYYTPDYVVRDTDVLAAFRMTPQPGVPPEECGAAVAAE
Cg_euchlorum AGFKAGVKDYRLTYYTPDYVVRDTDVLAAFRMTPQPGVPPEECGAAVAAE
（以下略）
こうして編集したファイルを別名で保存する（Gene1.bay）。
（MrModeltest と PAUP* が必要。これを用いない場合については後述）
MrModeltest というプログラムを用いて，進化モデルを推定する。
まず，MrModeltest のフォルダから、Gene1 フォルダに MrModelblock と mrmodeltest2.exe をコピーする （このフォルダ上で解析を進めるのに便利なため。後で削除する）。
PAUP* を開いて、Gene1.bay を読み込む。
MrModelblock を読み込む。 これにより同じフォルダに mrmodel.scores と mrmodelfit.log の二つのファイルが作られる。 前者が解析結果で、後者が解析のログである。
コマンドプロンプト（DOS window）を通じて mrmodeltest2.exe を実行するため， 以下のバッチファイルを用意する（メモ帳などテキストエディタを用いて， 下記の命令文を同じフォルダに ".bat" という拡張子をつけて保存）。ここでは mrmod.bat とする （mrmod.bat http://www2.tba.t-com.ne.jp/nakada/takashi/bayes/dl/mrmod.bat からダウンロードできる）。
mrmodeltest2 < mrmodel.scores > MrMod.out
このバッチファイルをダブルクリック。すると MrModeltest が mrmodel.scores を解釈して、 結果を MrMod.out に出力する。
注：mrmod.bat が mrmodel.scores および mrmodeltest2.exe と同じフォルダにないと機能しないので注意。
MrMod.out をメモ帳などで開くと、hLRT（階層的尤度比検定）もしくは、AIC（赤池情報量規準）に基いて選択されたモデルと、 MrBayes 用の構文（MrBayesBlock）がそれぞれ出ている。この hLRT（または AIC）の構文の部分に注目。
[!最初の[*******]の部分は、構文の情報に関するメモで、削除しても良い。この中に選ばれたモデルと、規準が記されている。 BEGIN 〜 END; がプログラムの構文で、********の部分はモデルによって変わってくる。
MrBayes settings for the best-fit model (******) selected by hLRT in MrModeltest 2.2
]
BEGIN MRBAYES;
Prset ****************************;
Lset ********************;
END;
というわけで以上の構文を、上記のデータファイル（Gene1.bay）の末尾にコピーし、上書き保存する。
注：hLRT と AIC の結果のいずれを用いるべきか，特に好みがなければ，よりシンプルなモデルを選択する傾向の強い hLRT の結果に従えば問題はないだろう。
PAUP* が利用できない場合、モデルの検証を行うには手間がかかる。
モデルの選択は解析結果に大きく影響するため、検証を省くことには問題があるが、系統樹を描くことは出来る。
ここでは恣意的に GTR+I+G というモデルを用いた解析方法を紹介しておく。
データファイル（Gene1.bay）を開き、末尾に以下の構文をコピーし，上書き保存する。
BEGIN MRBAYES;
Prset statefreqpr=dirichlet(1,1,1,1);
Lset nst=6 rates=invgamma;
END;
解析する塩基配列がタンパク質をコードしている場合， 第 1，第 2，第 3 コドンごとに進化速度を区別するモデルを用いるのが良いと考えられる。
よって，ここでは General Time Reversible（GTR）モデルと，site-specific rates model（SS）を併せたモデルを紹介する。
（気になる場合は，「タンパク質をコードしていない塩基配列の場合」と同様のモデル推定も行い，
この結果の尤度や事後確率と比較しても良いかもしれない。MrModeltest と同様の尤度比検定や赤池情報量規準の比較を行うことも考えられるが，
ここではそこまで踏み込まない。なお，SS も含めたモデル選択を行うプログラム（modelselect）が，
Life is fifthdimension で公開されている）。
以下の構文を，上記のデータファイル（Gene1.bay）の末尾にコピーし，上書き保存する。
BEGIN MRBAYES;
charset 1st_pos = 1-.\3;
charset 2nd_pos = 2-.\3;
charset 3rd_pos = 3-.\3;
partition by_codon = 3:1st_pos,2nd_pos,3rd_pos;
set partition = by_codon;
lset nst=6;
prset ratepr=variable;
END;
編集する場合の参考に、簡単な解説を加えたものも示しておく。
BEGIN MRBAYES;
charset 1st_pos = 1-.\3;[ここから，]
charset 2nd_pos = 2-.\3;
charset 3rd_pos = 3-.\3;
partition by_codon = 3:1st_pos,2nd_pos,3rd_pos;
set partition = by_codon;[ここまでがコドンの仕分けの定義]
lset nst=6;[nst=6 で GTR モデルを選択]
prset ratepr=variable;[ratepr=variable で上記の仕分けに基づき， site-specific rates モデルを選択]
END;
アミノ酸置換のモデルが複数存在するので，その中から選択する。 ここでは WAG というモデルを用いた例を紹介する（他に BLOSUM や JTT などのモデルも必要に応じて用いること）。
以下の構文を，上記のデータファイル（Gene1.bay）の末尾にコピーし，上書き保存する。
BEGIN MRBAYES;
lset rates=gamma;
prset aamodelpr=fixed(wag);
END;
編集する場合の参考に、簡単な解説を加えたものも示しておく。
BEGIN MRBAYES;
lset rates=gamma;[rates=gamma で gamma 分布を指定]
prset aamodelpr=fixed(wag);[aamodelpr=fixed(wag) で WAG モデルを選択。 この () 内を変更すれば他のモデルも利用可能なので，詳しくは MrBayes の説明書を参照。 例えば BLOSUM モデルであれば "blosum"、JTT モデルであれば "jones" を () 内に入れる]
END;
MrBayes はコマンド入力によって動くが，事前にコマンドを MrBayes Block の形で保存しておくのが良い。 実際の解析は長時間に渡るため，時間のあたりをつける目的で少ない世代数で初めに解析を走らせ， その後，許容時間に応じて（例えば 1 晩，1 週間など）世代数を決定するのが良いだろう。 以下にそのような目的の MrBayes Block を示す。
MCMC（マルコフ連鎖モンテカルロ法；これにより、ベイズの事後確率を求める）
のコマンドを記した以下の構文を必要に応じて編集する（mcmcA.bay
http://www2.tba.t-com.ne.jp/nakada/takashi/bayes/mcmcA.bay からダウンロードできる）。
begin mrbayes;
log start filename=Gene1.log replace;
mcmcp ngen=10000 printfreq=1000 samplefreq=100
nchains=4 savebrlens=yes filename=Gene1;
mcmc;
log stop;
end;
編集する場合の参考に、簡単な解説を加えたものも示しておく。
begin mrbayes;
log start filename=Gene1.log replace;[filename= でログファイルの名称の定義]
mcmcp ngen=10000 printfreq=1000 samplefreq=100[ngen= で世代数の指定、 samplefreq= で何世代ごとにデータをサンプリングするかを指定]
nchains=4 savebrlens=yes filename=Gene1;[nchains= でマルコフ連鎖の数の指定。 filename= で出力ファイルの名称の定義]
mcmc;
log stop;
end;
こうして編集したファイルを mcmc.bay などとして保存する。
mrbayes-3.1.2 のフォルダから、Gene1 フォルダに mrbayes.exe をコピーする （このフォルダ上で解析を進めるのに便利なため。後で削除する）。
mrbayes.exe をダブルクリックして起動させる。
execute Gene1.bay としてデータファイルを読み込む （execute は簡単に exe でも良い）。
exe mcmc.bay としてマルコフ連鎖を開始する。この設定では 1,000 世代ごとに経過と Average standard deviation of split frequencies（以下 ASDSF）が表示される（後者は連鎖が収束しているかどうかの指標となる）。 解析時間は世代数にほぼ比例するので，これで終了時間が推定できる。
10,000 世代の試行が終了すると、解析を続行するか否かを尋ねられる。
ここで，経過時間と収束の程度から追加の世代数を考える。
解析を終了して信頼できる樹形を描くためには、いくつかの条件が満たされる必要があることに注意したい。
連鎖が収束している状態でサンプリングを行っていること。
サンプリング数が十分量得られていること。
連鎖の収束は、上述の ASDSF の値を見ることによって確認することができる。 サンプリングしているデータが収束している場合には、ASDSF の値は 0 に近づくことが予想される。 実際には 0.01 未満に達していることを指標にするのが適当ではないかと思っている。
次にサンプル数は、少なすぎれば事後確率の値が不安定になることから、最低でも収束に達してから 1,000 （この場合、世代数に直すと 100,000）程度は必要であろう。
以上を考慮して、残りの世代数を決定する。ここで解析を続行する場合、 "Continue with analysis? (yes/no):" との問いに対して、"yes" と入力する。 次に "Additional number of generations:" に対して追加の世代数を入力する。
（以下に、どのように世代数を決定するかについて、私的な見解を紹介したい）
まず、解析を一晩行うことを考えて、何世代までできるかを計算する。
例えば一晩を 14 時間（≒ 50,000 秒）とするなら、
50,000 ÷（10,000 世代の解析に要した秒数）× 10,000 世代を合計で行えばよい。
一方でデータ行列が大きい場合などは、一晩では収束に至らない可能性も十分考えられる。 この場合、適当な世代数を追加して収束するまで同様の作業を繰り返すという手もあるが、 ASDSF が一定値を切るまで解析を続ける設定で、はじめから解析をやり直す方法がある。 この場合は、下記の MrBayes Block を用いる（mcmcB.bay http://www2.tba.t-com.ne.jp/nakada/takashi/bayes/mcmcB.bay からダウンロードできる）。
begin mrbayes;
log start filename=Gene1.log replace;
mcmcp ngen=50000000 printfreq=10000 samplefreq=100
mcmcdiagn=yes diagnfreq=500000 stoprule=yes stopval=0.01
nchains=4 savebrlens=yes filename=Gene1;
mcmc;
log stop;
end;
編集する場合の参考に、簡単な解説を加えたものも示しておく。
begin mrbayes;
log start filename=Gene1.log replace;
mcmcp ngen=50000000 printfreq=10000 samplefreq=100 [ngen= で、この場合は解析の最大世代数を指定。printfreq= でデータの表示頻度を指定]
mcmcdiagn=yes diagnfreq=500000 stoprule=yes stopval=0.01[diagfreq= で、 収束の程度の評価の頻度を指定。stoprule=yes で、ASDSF が一定値を切った時点で解析を終了するように設定。 その値を stopval= で指定]
nchains=4 savebrlens=yes filename=Gene1;
mcmc;
log stop;
end;
一晩寝て待つ（一晩以上かかる場合には、きちんと起床すること）。
解析の結果、Gene1.run1.p、Gene1.run2.p、Gene1.run1.t、Gene1.run2.t、 Gene1.log の 5 個のファイルが作成される。
結果の信頼性を確認するためサンプル数が統計的に充分であったことを確認する。 以下の命令を入力して、結果を要約させる。
sump burnin=*****
ここで、burnin は樹形数の 1/4 の値を入力する（初期設定の場合）。 上記の場合には 100 世代につき 1 樹形を得ているので、 burnin の値には（最終的な世代数）÷ 100 ÷ 4 の値を入力すればよい。
上記の命令を入力すると、まずグラフが現れる。これはサンプリングした範囲での尤度のグラフで、 連鎖が正しく収束していればプロットに増減の傾向はなく、グラフ上に点（1，2 または * の文字）が散らばるはずである。
出力の一番下には、パラメータの表が出てくる。この内重要なのは一番右列の PSRF の値で、 これは連鎖が収束していてサンプル数が充分であれば 1.0 に近づくと考えられる。 この値が 0.9 〜 1.1 程度に入っていることを確認すればよい。
尤度または PSRF の値に異常があった場合は、世代数をさらに追加することが推奨される。 それでも問題が解決しない場合には別の対策を講じる必要があるが、それについては MrBayes の説明書または Bayes 法（ベイズ法）の原理 http://www2.tba.t-com.ne.jp/nakada/takashi/bayes/idea.html の 「問題点と今後の展望」を参照のこと。
以下の命令を入力して、樹形を出力させる。
sumt burnin=*****
ここで、burnin の値は上記の sump の時と同じ値を入力すること。
この結果、Gene1.con、Gene1.parts、Gene1.trprobs の 3 つのファイルが作成される。
合意樹が Gene1.con に記述されているので、これを TreeView で開く。
Gene1.con には 2 つの系統樹が保存されており、最初のものに事後確率が書き込まれている。
事後確率の表示は TreeView の Tree→
MrModelblock、mrmodeltest2.exe、mrmod.bat および mrbayes.exe を削除する。
必要に応じてファイルを削除する。Gene1.bay、MrMod.out、 mcmc.bay、Gene1.log、 Gene1.con，Gene1.parts，Gene1.trprobs は残しておくと良いだろう（その他も必要性に応じて）。
Hall, B. G. Phylogenetic Trees Made Easy: a How-To Manual 2nd ed. (Sinauer Associates, Sunderland, 2004).
Life is fifthdimension
他、プログラムのダウンロード元、引用先についてはページ上部のリンク先と、各プログラムの解説などを参照のこと。
| トップへ | 生物の起源 | 雑記 | 新刊紹介 | 地質年代表 | 生物分類表 |
| Genera Cellulatum | 系統解析 | リンク集 | プロフィール | 掲示板 | 系統解析トップへ |