DNAシーケンスをアミノ酸シーケンスに変換してアライメントしたwork-conc.ptnファイルから進化系統樹を推定する作業を行う．系統樹推定のプログラムとしてはMolphyを用いる．~shimo/class/enshu200209ディレクトリにある実例をこのページでもクリックして見れるようにしてあるので，それを参考にすること．

１．まずwork-conc.ptnを覗いてみよう．このままでは見づらいので，

protst -a work-conc.ptn > work-conc-ptn.txt [enter]

と形式を変換してから見ると良い．work-conc-ptn.txtにある"CONSENSUS"というのは入力した配列の各列で最も頻度の多かったアミノ酸を並べたシーケンスである（多数決で決めたということ）．このCONSENSUSと同じアミノ酸はドットで表示され，それと異なる場合だけアミノ酸のアルファベットが書かれている．このように共通祖先から進化するにしたがってDNAが少しずつ変化する．比較的最近枝分かれしたばかりの生物（例えばヒトとチンパンジー）ではDNAの違いはごくわずかであるが，ずっと昔に枝分かれした生物（例えばヒトとイヌ）ではDNAの違いも大きい．逆にこの違いの情報を利用して，生物がどのような順序でいつころ分岐したかを推定することが可能である．

２．最初に近隣結合法（neighbor-joining, NJ法）を用いて推定を行う．コマンドは

cmd311 [enter]

である．このシェルスクリプトの内容は

#!/usr/local/bin/bash

# molphy nj analysis

outgroup=24

file=work-conc

protml -mfD $file.ptn > $file.dis

njdist -o ${outgroup} -w -t$file-nj -T$file-nj $file.dis > $file.nj

mv njdist.eps $file-nj.eps

repnamegb $file-nj.eps > $file-nj-name.eps

bbfix $file-nj-name.eps

である．ここでoutgroup=24というのは後で述べる方法によって自分で変更する必要がある．NJ法ではまずDNAシーケンス間の「距離」をもとめる．これを行っているのが "protml -mfD work-conc.ptn > work-conc.dis"である．ここでwork-conc.disはシーケンス間の距離行列である．距離行列からNJ法で系統樹を推定しているのが"njdost -o 24 -w -twork-conc-nj -Twork-conc-nj work-conc.dis > work-conc.nj"である．結果のグラフィックスはEPSファイルとしてnjdist.epsに書かれるが，namegb.txtにかかれた情報をもとに生物の名前を書き込んだものがwork-conc-nj-name.epsである．

　さて，自分でNJ法を行う場合はまずoutgroup=1などと適当な値を設定してcmd311を実行してみる．出力されたwork-conc.njをみると系統樹がテキスト形式で入っている．このなかで他の生物と最も遠い親戚にあたる生物が直接ついている枝の番号を読み取る．ここではNC_000891 (学名Ornithorhynchus anatinus，通称名platypus，カモノハシ）は，他の哺乳類とは最も遠い親戚なのでこれがついている枝の番号２４を読み取る．そしてこの番号をcmd311のoutgroup=24に書き込んで再びcmd311を実行する．

　outgroupは系統樹の根を指定している．DNAシーケンスを比較して得られる系統樹は本来根の位置が分からない木（unrooted tree)である．そこであらかじめ最も遠縁と分かっている生物（すなわちアウトグループ）を含めて系統樹を推定し，このアウトグループのついている枝が根の位置であると考える．

　NJ法による推定結果はwork-conc-nj.treにもテキスト形式で書かれている．このカッコ形式の意味は先ほどのグラフィックスによる結果と見比べれば理解できるだろう．実数は枝の長さである．この枝の長さの情報を取り除いて系統樹の形だけをカッコ式で表したのがwork-conc-nj.tplである．

３．NJ法で求めた系統樹は必ずしも信頼性が高いとはいえない．これに対して最尤法(Maximum-likelihood method, ML法）はよりよい方法と言われている．ML法では候補となる系統樹毎にそのもっともらしさの指標である尤度を進化確率モデルに基づいて計算する．そして最も尤度が高い系統樹を最尤系統樹として推定結果とする．実際には尤度の対数をとった対数尤度を計算して比較する．ML法の問題点は尤度の計算に時間がかかることである．そこでNJ法で求めたwork-conc-nj.tplから出発して，尤度を大きくすうるように系統樹を少しずつ変更して行き，もうそれ以上尤度を大きく出来なくなったら最尤推定の近似とみなすという方法が考えられる．これを実行するコマンドは

cmd331 [enter]

である．結果はwork-conc.mlr, work-conc-ml.tpl, work-conc-ml.tre, work-conc-ml.eps, work-conc-ml-name.epsなどに書かれているので内容を確認する．

　得られた系統樹work-conc-ml-name.epsを見ると，枝に数値がついている．これはローカルブートストラップ確率（LBP)と呼ばれて，その枝の信頼性を測る目安である．これが小さい枝はあまり信頼できない．ここでは根のほうからLBP=74, 86の枝はあまり信用できない．次の97は一見すると十分高いのだが，それもあわせてこれら３個の枝の信頼性を次に調べる．

４．work-conc-ml.tplを編集してwork-conc-constraint.tplを作成する．これは先ほどの３個の枝の長さをゼロに縮め，不確定な部分のカッコを { } に変更したものである．たとえばもとの系統樹が (((((A,B),C),D),E)F)）という構造の場合にA,B,C,D,Eの順序が不確定の場合には({A,B,C,D,E},F)とする．こうするとA,B,C,D,Eから作られる部分系統樹のあらゆる組み合わせに展開される．たとえば，(((((A,B),C),D),E)F)）や((((A,B),(C,D),E)F)）などの組み合わせ１０５通りに展開される．この展開した系統樹一覧を書き出すためのコマンドが

cmd351 [enter]

である．結果はwork-conc-user.tplに入る．これは近似的に尤度を計算していてその大きい順に展開されている．さらに正確な尤度を系統樹一覧のすべてに対して計算して尤度を最大にする系統樹を得るには

cmd371 [enter]

を実行する．この計算には比較的時間がかかるので，バックグラウンドジョブとして

cmd371 & [enter]

として実行したほうが良いかもしれない．（あやまって何度も同じコマンドを実行しないこと！　バックグラウンドジョブのことが良く理解できていないならばやらないこと．）　実行結果はwork-conc-user.mlrとwork-conc-user.treに書かれている．尤度を最大にするのは７番目の系統樹でグラフィックはwork-conc-user-name.epsである．

５．後で用いるために，以下のファイルを作成しておく．まずwork-conc-constraint.tplを作成するときに出てきたA,B,C,D,E,Fなどの部分系統樹をwork-conc-clade.txtというファイルに書いておく．そしてwork-conc-user.tplを編集して最初の行をwork-conc-user.tplのようにしておく．つまり一行目は系統樹の個数をあらわす整数だけ残し，あとは＃をつけてコメントにする．