選択圧(淘汰圧)を解析するために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 Evolution, 16(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 Evolution, 16(6), pp. 766-778
今回の記事の内容はZhangさんにご指導いただきました。ありがとうございます。
