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

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

{{category 系統解析}}

*PAUP*による最尤系統推定

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

*塩基置換モデルの選択

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

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

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

　ここで使用するPAUP*で扱うことのできる塩基置換モデルの中から最も良いモデルを選ぶための自作Perlスクリプトと使用方法を[[Site-Specific Rateなモデルも含めたモデル選択]]に記述しておきましたのでモデル選択に関してはこちらを参照願います。

*解析の実行

NEXUS形式で用意してあるデータと、モデル選択の結果得られる設定ファイル(modelXXX.nex)を

 paup> Execute data.nex
 paup> Execute modelAIC.nex

などとして読み込んだ後、

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

として以下に進む。

**全樹形を探索

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

**ブートストラップ解析によるクレード支持率の推定

 paup> Bootstrap NReps=1000 GrpFreq=Yes TreeFile=bootreps.tre Search=Heuristic / MulTrees=No
 paup> SaveTrees From=1 To=1 File=bootcon.tre SaveBootP=NodeLabels
**特定の系統樹の周辺樹形を探索

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

とすれば近隣結合系統樹を初期系統樹としてTBRアルゴリズムを用いて変形した樹形を探索します。最節約系統樹を求めてstart=に適宜指定すれば最節約系統樹を初期系統樹とすることもできる。下記のようにしてランダム生成した樹形のうち、現在有効になっているモデルで最尤な樹形を初期系統樹として探索することもできます。

 paup> GenerateTrees Random NTrees=10
 paup> HSearch Swap=TBR Start=current
 paup> SaveTrees BrLens=Yes File=tree.tre

　これは、以下のように実行するのと同義です。

 paup> GenerateTrees Random NTrees=10
 paup> SortTrees
 paup> HSearch Swap=TBR Start=1
 paup> SaveTrees BrLens=Yes File=tree.tre

　この例ではブートストラップ確率は別途求める必要があります。
*配列付加順の無作為化による樹形探索範囲拡大

　後述のratchet法より効果は低いもののPAUP*のみで実行可能な樹形探索範囲拡大法。

　通常、PAUP*では初期系統樹はNJ法ではなく、3OTU〜全OTUへ1つずつ配列を追加しつつ作成する(Stepwise法)という方法を使っています。この時、配列の追加順を無作為化することで初期系統樹を複数作成することができます。その周辺探索を繰り返すことで樹形探索範囲を拡げることができます。

 paup> HSearch Swap=TBR Start=Stepwise AddSeq=Random NReps=10
 paup> SaveTrees BrLens=Yes File=tree.tre

とすることで上記の方法が実行されます。NRepsの値を変えることで反復数を増減できます。

*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";

のダブルクォート内は解析に使うモデルに応じて適宜書き換えます。modelXXX.nexファイルからLSetコマンドの内容を切り貼りして下さい。pauprat領域にはコメントがきちんと書かれていますので、適宜オプションを変更しても良いでしょう(私は反復数をデフォルトの200回から500〜5000に増やし、場合によっては1度の探索に割り当てる時間もデフォルトの300秒からより長くするか、そもそも時間制限を無くしたりしています)。

　次に、

 % pauprat data.nex

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

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

 paup> Execute data.nex
 paup> Execute lratchet.nex

　これで解析が実行され、枝長を含まない系統樹(反復した分を含む)がlratchet.treに出力されます。このファイル内の樹形は尤度によるソートがなされているわけではありませんし枝長も記録されていませんので、一旦PAUP*を終了した後、

 % paup
 paup> Execute data.nex
 paup> modelAIC.nex
 paup> GetTrees File=lratchet.tre DupTrees=Eliminate
 paup> Set Criterion=Likelihood
 paup> SortTrees
 paup> SaveTrees BrLens=Yes File=lratchet.tre Replace=Yes

とすることで重複樹形を取り除き、尤度に基づいて並び替え、枝長を記録した樹形ファイルができあがります。さらに

 paup> Defaults LScores Longfmt=Yes
 paup> LScores all / ScoreFile=scoreandparameters.txt

とすると各樹形の尤度とパラメータの推定値を出力します。

*高速化するには

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

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

　データサイズ(特にOTU数)が小さければ全樹形探索、ある程度大きければ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

*番外編

**PAUPRatのインストール

　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
 % cd ..

　./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

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