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 (アミノ酸組成を推定しない)」と設定することで，どのようなデータを用いても同じ値が得られます．

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