PAUP*による最尤系統推定の変更点

• 追加された行はこのように表示されます。
• 削除された行はこのように表示されます。

{{category 系統解析}}

*概要

　最尤系統樹の探索には様々なソフトが利用されるが、ここでは最もよく使用されているPAUP* 4.0 beta版での探索方法を解説する。他のソフトウェアでは利用できる塩基置換モデルに制限があったりするため、最適な塩基置換モデルの選択を行ってから、その塩基置換モデルに基づいた最尤推定を行うという至極当然のことができない。このため、論文にする場合など、厳密な解析が要求される場合には以下の方法を利用することを強く推奨する。なお、以下の説明ではWindows/Mac/UNIX版で共通に利用できるコマンドラインインターフェイスを想定しており、Mac版でのみ利用できるメニューからの操作については触れない。また、アミノ酸配列の解析も触れない。

*塩基置換モデルの選択
　最尤推定のキモはモデルの選択にこそあると言っても過言ではない。適したモデルを利用していない最尤推定では最尤推定の意味が無いと言ってもいい。

　そもそも、最尤推定が何故良いかと言うと、尤度比やAIC(赤池情報量基準)やBIC(ベイズ情報量基準)を用いて「モデルの評価ができる」ということにある。これによって考えられ得る限りの中では最も良いモデルを選択することができる。(注: 正確にはモデル選択と最尤推定は別のものであるが、現在では切っても切り離せない密接な関係にあると言っていいと思う)

　ただし、これもあくまで「考えられ得る限りでは」の話であり、「とても考えられないモデル」がより良いモデルである可能性は残る。つまり、「最も良いモデルを用いて最尤推定を行う」ためには、そもそも「ありそうなモデル」を全てリストアップした上で、「最も良いモデルを選択」しなければならない。

　しかしここで、塩基配列やアミノ酸配列の置換モデルは過去の研究によりほぼ「考えられ得る限り」のモデルは出尽くしていると思われる。このため、系統樹の最尤推定に限って言えば既存のモデルから最も良いモデルを選択すればよい。

　モデル選択用の自作Perlスクリプトと使用方法を[[Site-Specific Rateなモデルも含めたモデル選択]]に記述したのでこちらも併せてご覧下さい。

　系統樹の最尤推定に最適なモデルを選択するには、以下のソフトを利用する。

-PAUP* (以下ではVer.4.0beta10を利用)
-Modeltest (以下ではVer.3.5を利用)

　なお、筆者の環境ではmodeltest 3.5の配布ファイルがソフトによってはうまく展開できないことがあった。もし展開できない場合には別のソフトで展開してみるといいだろう。

　まず、解析対象の塩基配列をPAUP*が利用するNEXUS形式で用意し(試すだけならmodeltestに付属のデータでもよい)、PAUP*を適当な作業ディレクトリをカレントとして起動し、

 paup> Execute sample.nex

などとしてデータを読み込む。

　次に、modeltestに付属しているmodelblockPAUPb10というファイル(これにはNEXUS形式で実行コマンドが記述してある)を同様に読み込む。

 paup> Execute modelblockPAUPb10

　すると、ファイル内のコマンドを次々に実行して全56モデルでの対数尤度が計算され、上記ファイルを編集していない限りはカレントディレクトリに「model.scores」というファイルが出力される。

　出力された「model.scores」を以下のようにしてmodeltestに読み込む。

 % modeltest < model.scores

　結果は標準出力(通常は画面)に出力されるので、

 % modeltest < model.scores > output.txt

などとしてファイルに保存することもできる。modeltestにはいくつかオプションがあるのでドキュメントやオプション無しでの起動時のメッセージを参照して適宜設定する。

　modeltestの出力は、尤度比検定によって最適モデルを検討した結果とAICによって最適モデルを検討した結果に分けられる。また、PAUP*においてそれぞれの最適モデルで最尤系統樹を探索するためのコマンドも出力されるので、適宜そのコマンドをデータのNEXUSファイルにペーストしてexecuteすればデータと解析の設定が同時に読み込まれる。別のNEXUSファイルとして保存してexecuteしたり、直接コマンドラインから入力してもよい(その場合コマンド末尾のセミコロンは不要)。ちなみに、コマンドの領域は以下の部分。

 BEGIN PAUP;
 Lset hogehoge;
 END;

*解析の実行

上記のデータと設定を読み込んだ後、

 paup> Set Criterion=Likelihood [系統推定法を最尤法に設定]

として以下に進む。

**全樹形を探索

 paup> BandB MulTrees=No [全樹形から最尤系統樹を探索]
 paup> SaveTrees BrLens=Yes File=tree.tre [枝長も含めてファイルに保存]

で樹形探索が実行される。樹形の探索が終わったら

 paup> Bootstrap NReps=1000 GrpFreq=Yes TreeFile=boot.tre Search=Heuristic / Swap=NNI MulTrees=No

でブートストラップ解析を実行できる。

**特定の系統樹の周辺樹形を探索

 paup> HSearch Swap=TBR Start=NJ MulTrees=No
 paup> SaveTrees BrLens=Yes File=tree.tre

とすれば近隣結合系統樹を初期系統樹としてTBRアルゴリズムを用いて変形した樹形を探索します。最節約系統樹を求めてstart=に適宜指定すれば最節約系統樹を初期系統樹とすることもできる。下記のようにしてランダムに選んだ樹形のいくつかを初期系統樹として探索することもできます。

 paup> RandTrees NReps=10
 paup> HSearch Swap=TBR Start=current MulTrees=No
 paup> SaveTrees BrLens=Yes File=tree.tre

　この例ではブートストラップ確率は別途求める必要がある。

　また、上記の代わりに

 paup> Puzzle NPuzzles=1000000 GrpFreq=Yes UseApproxL=No TreeFile=boot.tre MinMemReq=no
 paup> SaveTrees BrLens=Yes File=tree.tre

とすることでquartet puzzling法を用いて近隣結合系統樹の周辺の樹形探索とブートストラップ確率もどき(quartet puzzling support value)の計算を一度に行うこともできる。NPuzzles=の値を変更することでステップ数を変えられます。

ブートストラップ確率を枝にラベリングした上で枝長も含めて系統樹をファイルに出力することはPAUP*ではできない?ようです。もしできるのでしたらやり方を教えて下さい。

*ratchet法による樹形探索範囲の拡大

　OTU(系統樹の枝の末端数)が増えてくると上記のように初期系統樹の周辺だけを探すなどしなければ樹形探索に非常に時間がかかります。しかし、樹形空間上で初期系統樹に近い樹形だけを探索するだけだと、初期系統樹(多くの場合近隣結合樹)からあまり近くないけど実は最尤な樹形を見逃してしまう可能性があります。そこで、ランダムに重み付けをしたデータから作成した近隣結合樹の周辺を探しまくるratchet法という手法があります。ここではそれをPAUP*を使用して行う方法について述べます。

　最尤法による系統推定でratchet法を適用するには以下のソフトを用いる。

-PAUPRat

　また、当然PAUP*も必要です。

　PAUPRatは元々最節約法による系統推定でratchet法を適用するためのソフトですが、[R. A. Vos氏のページで説明されている|http://www.sfu.ca/~rvosa/likelihoodratchet/]ように、最尤法にも転用できます。

　まず、解析対象のデータをNEXUS形式で用意し、[Vos氏が用意されているスクリプト|http://www.sfu.ca/~rvosa/likelihoodratchet/input.nex]の

 begin pauprat;
 ...
 end;

という部分を用意したNEXUS形式データの末尾にエディタなどで書き加えます。この時、

 dimensions nchar=;

とある=の後ろに、

 begin data;
 ...
 end;

で囲まれたデータ領域内にあるnchar=の後ろの数字(塩基配列のconsensus length?)を書き加えておきます。また、

 startcmd "lset nst=1  base=equal  rates=equal  pinv=0";

のダブルクォート内は解析に使うモデルに応じて適宜書き換えます。上記のようにmodeltestの出力を使う場合は

 Lset ...;

という部分の末尾のセミコロンを除いて記述しておけばいいでしょう。pauprat領域にはコメントがきちんと書かれていますので、適宜オプションを変更しても良いでしょう(私は反復数をデフォルトの200回から500〜5000に増やし、場合によっては1度の探索に割り当てる時間もデフォルトの300秒からより長くしたりしています)。

　次に、

 % pauprat hogehoge.nex

などとしてlratchet.nexを作成します。このファイルに解析に使う実行コマンドが出力されます。

　lratchet.nexができたら、PAUP*を起動してデータファイルと実行スクリプトファイルをexecuteします。

 paup> Execute hogehoge.nex
 paup> Execute lratchet.nex

　これで解析が実行され、枝長を含まない系統樹(反復した分を含む)がlratchet.treに出力されます。枝長を含めて出力するなら

 paup> SaveTrees BrLens=Yes File=hogehoge.tre

としましょう。

*高速化するには

　画面出力やディスクアクセスを抑制したりするオプションを利用することで多少は高速化することができます。詳しくはコマンドリファレンスのSETコマンドの項などを参照して下さい。正式版のリリース時に公開されるというソースコードが手に入れば、コンパイル時にCPU最適化を行うことでさらに高速化も可能でしょう。

*で、結局どれが良いの?

　データサイズが小さければ全樹形探索、ある程度大きければratchet、あまりにも大きくてブートストラップ確率の計算に死ぬほど時間がかかってしまう場合はquartet puzzling、でしょうか。近隣結合系統樹の周辺だけっていうのはかなりヤバいと思います。quartet puzzlingもそれに該当するので使う場合はnpuzzlesの値をかなり大きくした方が良いでしょう。また、ratchetの結果と比較してみるのもいいでしょう(ratchetの反復数も大きく取ること)。

*参考URL

-http://cse.niaes.affrc.go.jp/minaka/files/notes/modeltest.html
-http://cse.niaes.affrc.go.jp/minaka/files/notes/phylostat1.html
-http://paup.csit.fsu.edu/Quick_start_v1.pdf
-http://paup.csit.fsu.edu/Cmd_ref_v2.pdf

*番外編

**Modeltestのインストール

　Windows / MacOS Xなら配布ファイルに実行可能バイナリが収録されているのでそれを利用すればよい。展開したらPATHの通ったディレクトリに置いておく。その他の場合も基本的には配布ファイルを展開してmakeするだけ。以下FreeBSDの場合。

 % unzip modeltest3.5.zip
 % cd Modeltest3.5/source
 % make
 % su
 # cp modeltest3.5 /usr/local/bin/modeltest
 # exit

　Modeltest3.5/source/Makefileはコンパイラなどに応じて適宜書き換える。gccならLinuxでも書き換える必要は無いと思う。Linuxの場合は/usr/binへインストールした方が良いかもしれない。

**PAUPRatのインストール

　Modeltest同様Windows / MacOS X用バイナリがあるのでWindows / MacOS X上で実行するならそちらを利用する。Mac環境向けにはGUIも備えた専用のバイナリがある。その他の環境ではコンパイルする必要がある。PAUPRatのコンパイルには[NEXUS Class Library|http://hydrodictyon.eeb.uconn.edu/ncl/]が必要なのでまずはこれをコンパイルする(インストールは必要無い、というか、インストールしてもヘッダがインストールされないんですね。タコconfigureのせいで。インストールする場合は別途ncl.hを適当なディレクトリにコピーする)。

 % tar xvzf ncl-2.0.tar.gz
 % cd ncl-2.0
 % ./configure
 % make

　./configureでは--helpを付けるとオプションが表示されるので適当に指定する。大抵の環境では何も指定しなくてもいいと思いますが。場合によってはMakefileを編集する必要もあるかもしれません。

　ということでいよいよPAUPRatのコンパイルとインストールです。

 % tar xvzf pauprat.tgz
 % cd pauprat

　ここで、Makefileを以下のように編集(NCLをインストールしてncl.hも適当なディレクトリにコピーしてある場合は不要)。

 LIBDIR		= ../ncl-2.0/src
 LIBS		= ../ncl-2.0/src/libncl.a

　これはpauprat.tgzとncl-2.0.tar.gzの展開ディレクトリが同じと仮定してますので注意。んで、編集が終わったら以下のようにする。

 % make
 % su
 # cp ./pauprat /usr/local/bin

以上。
*PAUP*による最尤系統推定
　内容が古くなったので削除します。新しく[分子系統学演習 - データセットの作成から仮説検定まで|http://www.fifthdimension.jp/documents/molphytextbook/]を書きましたのでこちらをご覧下さい。ただしPAUP*については扱っていません。