{{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=100 TreeFile=bootreps.tre Search=Heuristic paup> SaveTrees From=1 To=1 File=bootcon.tre SaveBootP=NodeLabels 　この例ではStepwise法で生成した初期系統樹の周囲をTBRによる枝交換で探索する。OTU数が多いときにはStepwise法による初期系統樹生成に時間がかかるため、以下のように初期系統樹を近隣結合樹にし、さらに枝交換をNNIやSPRにして1レプリケート当たりの樹形探索範囲を絞る。 paup> Bootstrap NReps=100 TreeFile=bootreps.tre Search=Heuristic / Start=NJ Swap=NNI 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 　以上と同様に、以下のように実行することで最節約系統樹を初期系統樹として樹形探索することもできます。 paup> Set Criterion=Parsimony paup> HSearch AddSeq=Random NReps=1000 paup> Set Criterion=Likelihood paup> HSearch Swap=TBR Start=current 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の値を変えることで反復数を増減できます。 　ただ、この方法ではOTUの枝を付加する場所を最尤法で決めているため、OTU数が多い場合には初期系統樹の生成自体に時間がかかってしまい、肝心の樹形探索に割く時間が相対的に減ってしまいます。そのような場合には下記のratchet法を使う方が良いでしょう。Stepwise法による初期系統樹生成の際のみ最節約法を使えるオプションの追加が望まれます。 *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 以上。