multidistirbute の使い方: mtREV + F + Γ モデルを
用いたアミノ酸配列解析
2008 年 6 月 30 日 改訂
井上 潤

〜〜〜このページはまだ作成途中です〜〜〜



multidistribute を初めて使う方は,まずはこちらを参考に塩基配列の解析をやってみて下さい.

アミノ酸を用いた最尤法は,mtREV などで知られる既存の rate matrix を元に解析を行うのが一般的です.multidistribute のアミノ酸解析でも,既存の rate matrix を用いて年代推定を行います.このため,塩基配列データのように rate matrix に用いるパラメータを推定する必要がありません.ただ,F モデル (データからアミノ酸組成を計算して mtREV モデルに組み込む) や Γ モデル (サイトごとの進化速度差異を補正する) を用いた解析を行うことは可能です.ここでは mtREV matrixを用いて F 補正およびガンマ補正を施した解析を紹介します.multidistirubte の HP にもありますが,熊澤先生が作成されたデータセットを用いて解説します.英文の付属マニュアルがあるので,そちらも参考にして下さい.



フローチャート


codeml.c の書きかえ
F + Γ モデルを用いない場合は,codeml の解析を行う必要はありません.

F + Γ モデルを用いる場合は,codeml のソースコードを少しだけ書きかえてアミノ酸組成とΓ補正のパラメータを得る必要があります.フローチャートにあるように,codeml のoutfile に従って 4 種類のファイルを作成します.1 パーティションの解析に 20 ファイル程度必要で,これを例えば 12 パーティションの解析であれば 12 × 20 ファイル作成することになります.煩雑になるので,パーティションごとにフォルダを作成し,分けて保存した方が良いです.





codeml.c (src というフォルダーに入っている) を書き換えます.「for (itree)」で検索すると,paml3.14beta であれば752 行目あたりに (paml3.15 であれば 772 行目あたりに),

} /* for(itree) */

が表れます.この直前に,熊澤先生が書かれたマニュアルに指定されている 7 行のコード (下にもあります) をコピー&ペーストして保存します.同じバージョンであってもダウンロードによって「 } /* for(itree) */」直前の行が異なっており,上の図とまったく同じではないようなので気を付けてください.

そして,paml のマニュアルに従って通常通りのコンパイルを行います.コンパイル後は「codeml_F」のようにファインダや mv コマンドを使ってプログラムの名前を変更しても,問題なく作動します.

fprintf(fout, "\nRoot, U, V, freq\n");
FOR(i,com.ncode) Root[i]=exp(Root[i]*.01);
matout(fout, Root, 1, com.ncode);
matout(fout, U, com.ncode, com.ncode);
matout(fout, V, com.ncode, com.ncode);
matout(fout, com.pi, 1, com.ncode);
exit(0);



メモ

  • codeml の解析は baseml に比べると非常に早く終わります.40 OTU, 100 残基でも,おそらく 10 分はかからないと思います.これ以上時間がかかった場合は,infile 各種を見直した方が良いです.一つ一つ infile をうまく動いた解析のファイルと交換すると良い確認になります.

  • paml3.15.OSX_G5.tar.gz は G4 ではうまくコンパイルできないようなので,old version のページにある v3.14 の PAML をコンパイルして使いました.

  • estbranches は Win で解析をしましたが,こちらも長くても 10 程度で 1 パーティションの解析が終了しました.どうも Win での解析は,改行コードは Mac は読みませんが,Unix は問題なく読み込むようです.

  • 小さいデータセット (e.g., <50 残基ぐらいでしょうか) では,codeml の推定値が 0 になることがあります.この場合,estbranches (アミノ酸バージョン) は動き出しますが,解析が終わらなくなります.0 の代わりに非常に小さい値を入れるなど (mtREV の原著論文 [Adachi and Hasegawa,1996] で,モデルを作成する際に matrix の要素が 0 になるのを避けるために同様の操作をしています.) の工夫がいるかも知れません.




codeml の output file を estbranches の infile に書きかえる

上の図に従って,codeml の output file から必要なデータを抜き取り,estbranches のインファイルを作成します.Excel などを使って工夫して作業にあたる必要があります.




ガンマ補正について


gamma 補正を行う場合は,mtREV などの matrix にカテゴリーごとの rate の差を組み込まなければなりません.このために matrix を一度,固有値と固有ベクトルに分解して,固有値 (20×1 の行列) の値をそれぞれ rate で累乗し,固有値・固有ベクトルを別々に最尤法の計算 (ここでは estbranches) にもっていくことになります.


**** 以下のように書きましたが,そうでない場合も出てきまし
た.まだよく理解していません.
****

つまり matrix から最初に得られる固有値と固有ベクトルは,用いるデータには無関係なものです.実際に T3 で配布されている mtmam の固有値,固有ベクトルは,codeml で 「aaRatefile = mtmam.dat, 「model=2 (アミノ酸組成を推定しない)」と設定することで,どのようなデータを用いても同じ値が得られます.

****************