授業/H17/系統学特論/PAUP*使用法/最尤法

目次

ModeltestとPAUP*による最尤系統樹の推定 †

1. Modeltest3.7の準備 †

【以下の操作は説明を簡単にするための処理。OSにおけるフォルダ（ディレクトリ）やパスなどを自分で把握している人は、自分が使いやすいように変更してもらってかまわない】

ダウンロードしたmodeltest3.7.zipをデスクトップで展開すると、modeltest3.7というフォルダができる

modeltest3.7の中にはさらに、Modeltest3.7 folderというフォルダが入っている

Modeltest3.7 folderの中に入っている全てのファイルとフォルダを、ひとつ上の階層のフォルダ（デスクトップのModeltest3.7）に移す。

2. PAUP*：modelblockPAUPb10の実行 †

PAUPを起動し、PEDIC.NEXをExecuteする。

次にFile > Open から、Modeltest のpaupblockフォルダに入っているmodelblockPAUPb10をExecuteモードで開く（Modeltestのフォルダはデスクトップにインストールしておく）

modelblockPAUPb10と同じModeltest3.7フォルダのpaupblockフォルダにmodel.scoresができる

3. PAUP*：modelblockPAUPb10の内容説明 †

modelblockPAUPb10の中身は次のようになっている：

ログファイルは、この例ではmodelfit.logという名前で、元データ〔この場合はPEDIC.NEX）と同じフォルダに保存される。NJでガイドになる樹形を推定している。

lscores 1: メモリ中のtree 1の尤度を計算
nst: 塩基置換の種類の数（1:  全て同じ、2: TvとTrで違う、6: 全て違う）
rmat: 塩基間の変化行列
base[basefreq]: 塩基組成
rates: サイトごとの進化速度の違い
shape: ガンマ分布の設定
pinv: 変化しないサイトの割合

上のようなパラメータの組み合わせを変えることで、56種類のモデルについて、NJ法で得られた樹形の尤度を求めている

4. コマンドプロンプトの使い方 †

スタートアップメニュー>全てのプログラムを表示> アクセサリ >コマンドプロンプト でコマンドプロンプトを起動する。

5. コマンドプロンプトの説明 †

ここに現れたシンプルな画面がコマンドプロンプトという、コンピュータへの命令を文字で入力して与えるところ。

ためしに、notepad と入力してみるとメモ帳が起動するはず。つまり、

「Windowsからメモ帳のアイコンをクリックすること」　＝　「コマンドプロンプトにnotepadと入力すること」

画面に表示されている、C:\Documents and Settings\tkaji>は現在参照しているディレクトリのパス（アドレス）を示す。

たとえば、デスクトップに保存されているModeltest3.7フォルダのプログラムを動かすためには、>の後に、

デスクトップ\modeltest3.7\bin\modeltest.3.7.win

と入力する。

6. パス（アドレス）の指定方法 †

コマンドプロンプトからファイルを指定するにはパス（通り道）の指定が必要。

パスが分からないときはエクスプローラー（注：IEでは無い）を立ち上げて、アドレスという部分を見る。入力の手間を省くには、パスの文字列をコピーする。

コマンドプロンプトで、

C:\Documents and Settings\tkaji\デスクトップ\modeltest3.7\bin\modeltest3.7.win

と入力してもダメ。パスの指定に空白文字（スペース）は使えない。スペースの入った文字はクォートする（" "で挟む）必要がある。

C:\"Documents and Settings"\tkaji\デスクトップ\modeltest3.7\bin\modeltest3.7.win

長いパスを入力するのが面倒ならば、

cd　でディレクトリを移動]
modeltestをスタートアップメニューに登録する
ショートカットを作る
環境変数 path に登録する

などの方法がある。

今回はcd でディレクトリを移動して操作を行う。cd の場合はクォートは必要ない

cd C:\Documents and Settings\tkaji\デスクトップ\modeltest3.7

7. Modeltestの実行の準備 †

エクスプローラでModeltest3.7のフォルダを開き、パスをコピーする（編集>コピー）。

＞の後に cd　と入力した後、コマンドプロンプトのアイコン　をクリックし、編集>貼り付け　あるいは、>　の後に

cd C:\Documents and Settings\tkaji\デスクトップ\modeltest3.7

と入力する

8. Modeltestの実行 †

model.scoresの確認: dir paupblock でpaupblockフォルダ内のファイルを一覧できる（model.scoresができていることを確認する）。

modeltestの実行: ＞の後に

bin\modeltest3.7.win < paupblock\model.scores > mydata.modeltest

と入力する。処理が正常に終われば（メッセージも無しに）プロンプト( > )が表示される（図）。同じフォルダにできたmydata.modeltestというファイルに結果が書き込まれている。

9. Modeltestの結果 †

これまでの手順どおりに処理を行うと、デスクトップのmodeltest3.7というフォルダに、mydata.modeltestという名前のテキストファイルができている

mydata.modeltestをPAUP*やメモ帳(notepad)などで開く

10. 尤度比検定 (hLRTs)の結果 †

尤度比検定は、最尤推定に用いたモデル２つを比較して、２つのモデルに明らかに差があるかどうかを統計的に検定する方法（この説明はかなり直感的です）。例えば、JCは全ての塩基間の置換率が同じであるとするモデルだが、K2Pはトランジションとトランスバージョンでは置換率が異なるとしている。K2Pの方がパラメータ数が多いので、現実データに対する説明能力はJCより高くなるが、その２つから得られる尤度の違いを比較してみて、そんなに違いがなければわざわざ複雑なモデルを使わなくてもよいだろうという考え方。

次のAICについても、2004年度の集中講義で来られた三中さんが、千葉大での講義ノートをホームページに載せています。Modeltestのデータの見方については、そちらを参照して下さい（http://cse.niaes.affrc.go.jp/minaka/cladist/NOTES/modeltest.html）

ひとつのモデルが選択されて（この例ではF81）

BEGIN PAUP;
Lset  Base=(0.3853 0.1549 0.1458)  Nst=1  Rates=equal  Pinvar=0;
END;

というPAUP blockが作られた。

11. AIC †

AIC(Akaike Information Criterion)は、

AIC(M1) = -2 x MLL(M1) + 2 x k1
但しMLL(M1)はモデルM1のもとでの最大対数尤度、k1はモデルM1に含まれるパラメータ数）

AICは最尤推定につかうモデルを選ぶ基準になるもので、AICの値が小さいほど、真のモデルに近いと考える。やっていることは、パラメータの数が多いモデルは、いくら尤度が高くてもだめですよということ。現実のデータにモデルを当てはめるとき、パラメータを多くすれば、全てのデータをうまく説明するようなモデルを作ることはできる。でも、それは現実には即していないので、なるべくパラメータ数を少なくおさえつつも尤度が高くなるモデルを選びましょうというもの。

Modeltestで最小のAICを与えるモデルが選択されて（この例ではK81uf）

BEGIN PAUP;
Lset  Base=(0.3021 0.1755 0.1489)  Nst=6  Rmat=(0.9583 1.7489 0.0929 0.5681 1.7489)  Rates=equal  Pinvar=0.4237;
END;

というPAUP blockが作られた。

12. PAUP*でMLによる系統解析 †

系統学特論：PAUPブロックから、ML, non-coding ModeltestのPAUP block（下に示したもの）を貼り付ける

begin paup;
set autoclose=yes warnreset=no increase=auto;
set criterion=parsimony; [最尤法で探索する系統樹の最初の樹形を節約法で推定。NJを使っても構わない]
hsearch; [推定には発見探索法を用いる]
set criterion=likelihood; [系統樹推定の方法を最尤法に変更]
[**ここにModeltestの結果のLsetから始まる行を挿入**
Lsetは最尤推定を行うときにどんなモデルを使うかが書かれている]
hsearch start=1; [さっき節約法で推定した最初の樹形から発見探索法で最尤系統樹を探す]
savetrees brlens=yes maxDecimals=4 file=ML_out.tre replace=yes; [得られた系統樹を保存]
end;

上のコメント部分にModeltestの結果得られたlsetの行（ここでは　Lset Base=(0.3853 0.1549 0.1458) Nst=1 Rates=equal Pinvar=0;）を貼り付ける

貼り付け後

begin paup;
set autoclose=yes warnreset=no increase=auto;
set criterion=parsimony; 
hsearch;
set criterion=likelihood;
Lset  Base=(0.3021 0.1755 0.1489)  Nst=6  Rmat=(0.9583 1.7489 0.0929 0.5681 1.7489)  Rates=equal  Pinvar=0.4237;
hsearch start=1;
savetrees brlens=yes maxDecimals=4 file=ML_out.tre replace=yes;
end;

Executeする。解析結果はML_out.treに保存される。

13. 解析結果：PAUP*でdescribeして結果確認 †

コマンドラインで describe と入力すると、解析の詳細と樹形が表示される

14. Tree View で結果を表示 †

ML_out.treをTreeViewで開いてPhylogram表示