codemlのCodonFreqの選択について

選択圧(淘汰圧)を解析するためにPAMLのcodemlがよく使われる。

いくつか日本語でも解説が書かれており、それに従って使っていることが多かったが、CodonFrequencyについての記述は日英共に解説が少ない。どうやって選択すべきか調べていたが、灯台もと暗しで知人がきちんと統計的なモデル選択の方法を論文に記載していたので書き記しておく。

Zhang, Z. and Nikaido, M. (2020) Inactivation of ancV1R as a predictive signature for the loss of vomeronasal system in mammals. Genome Biology and Evolution16(6), pp. 766-778

CodonFreqのパラメータ数

まず codeml では4つのパラメータを選択することができる。コドン置換モデルにおける平衡コドン頻度をどれにするかを選ぶ。

https://snoweye.github.io/phyclust/document/pamlDOC.pdf

  • 0. 頻度が等しいと仮定する場合(標準遺伝暗号の場合、1/61)
  • 1. 平均塩基頻度から計算
  • 2. コドンの3つの位置での平均塩基頻度から計算
  • 3. 自由パラメータ

がある。それぞれパラメータ数は0, 3, 9, 60である。

機械学習や統計をかじったことがある人は掴めるだろうが、パラメータ数は多ければ多いほど Overfitting が起きやすくなる。だから、パラメータ数が小さく高い尤度を出すモデルを設定してやることが良いモデル選択となる。

AIC (赤池情報量規準) を用いたモデル選択

それぞれ、0-3の4つのモデルを設定したcodeml.ctlファイルを入れたディレクトリを作成し、model0で実行する。

      seqfile = ../day1.fa
     treefile = ../primate_lysozyme_spe.tre
      outfile = ./model0_CF0_result.txt

        noisy = 3   * 0,1,2,3,9: how much rubbish on the screen
      verbose = 0   * 1: detailed output, 0: concise output
      runmode = 0   * 0: user tree;  1: semi-automatic;  2: automatic
                    * 3: StepwiseAddition; (4,5):PerturbationNNI

      seqtype = 1   * 1:codons; 2:AAs; 3:codons-->AAs
    CodonFreq = 0   * 0:1/61 each, 1:F1X4, 2:F3X4, 3:codon table
        clock = 0   * 0: no clock, unrooted tree, 1: clock, rooted tree
        model = 0    
                    * models for codons:
                        * 0:one, 1:b, 2:2 or more dN/dS ratios for branches

上の例はCodonFreqが0の場合。

CodonFreqを0-4に変えてそれぞれ実行し、Outputfileの中のlnLの値を取り出す。

最後に AIC (赤池情報量規準) を求めてやる。AICとは、統計モデルの良し悪しを評価するための指標の1つである。これが先程の話に出てくる、「パラメータ数を抑えながらも尤度が高いモデル」を選択する指標となる。求め方は以下の通り。

AIC = -2lnL + 2k

lnLに先程のOutputfileのlnLを、kに前項のパラメータ数を代入して求める。この値が4つの中で最も小さなものが今回の配列における良いモデルであることが推定できる。

参考文献

Zhang, Z. and Nikaido, M. (2020) Inactivation of ancV1R as a predictive signature for the loss of vomeronasal system in mammals. Genome Biology and Evolution16(6), pp. 766-778

今回の記事の内容はZhangさんにご指導いただきました。ありがとうございます。

コメントを残す

メールアドレスが公開されることはありません。 が付いている欄は必須項目です