Branch モデル
旧世界ザルの共通祖先の枝を対象に、 正の自然選択 (d_\text{N} / d_\text{S} > 0) を検出する。
まず対立仮説として、注目する枝で d_\text{N} / d_\text{S} が異なると仮定して codeml を動かす。
系統樹 (lysozymeSmall.fa.treefile) をエディタで編集して、 注目する枝 (i.e. d_\text{N} / d_\text{S} > 0 を想定する枝) に #1 を振る:
(Hsa_Human:0.0087447996,Hla_gibbon:0.0130303081,(((Cgu/Can_colobus:0.0156654142,Pne_langur:0.0181835564)#1:0.0280807319,Mmu_rhesus:0.0070703774):0.0147707727,(Ssc_squirrelM:0.0140657314,Cja_marmoset:0.0076397044):0.0477319527):0.0256856586);
Warning in fortify(data, ...): Arguments in `...` must be used.
✖ Problematic arguments:
• as.Date = as.Date
• yscale_mapping = yscale_mapping
• hang = hang
ℹ Did you misspell an argument name?
#1 を振った枝とその他の枝で d_\text{N} / d_\text{S} が異なる、 というシナリオの .ctl を書く:
outfile = branch_alt * 何でもいいけど branch モデルの対立仮説と分かるように
model = 2 * 記号の有無で異なる ω を推定
NSsites = 0 * サイト間では ω は一定
fix_omega = 0 * ω の値を配列から推定
omega = 1 * 推定は ω=1 からスタート
* 他のパラメータは変更なし
この .ctl を指定して codeml を実行:
指定した出力ファイル branch_alt を見てみる。 まずは最後の方に書かれている各枝の \omega の値:
dS tree:
(Hsa_Human: 0.011030, Hla_gibbon: 0.016794, (((Cgu/Can_colobus: 0.018929, Pne_langur: 0.022497): 0.009367, Mmu_rhesus: 0.008406): 0.018929, (Ssc_squirrelM: 0.017701, Cja_marmoset: 0.010258): 0.052331): 0.030198);
dN tree:
(Hsa_Human: 0.007565, Hla_gibbon: 0.011517, (((Cgu/Can_colobus: 0.012981, Pne_langur: 0.015429): 0.032838, Mmu_rhesus: 0.005765): 0.012982, (Ssc_squirrelM: 0.012140, Cja_marmoset: 0.007035): 0.035889): 0.020710);
w ratios as labels for TreeView:
(Hsa_Human #0.68581 , Hla_gibbon #0.68581 , (((Cgu/Can_colobus #0.68581 , Pne_langur #0.68581 ) #3.50573 , Mmu_rhesus #0.68581 ) #0.68581 , (Ssc_squirrelM #0.68581 , Cja_marmoset #0.68581 ) #0.68581 ) #0.68581 );
Time used: 0:01
期待通り、#1 を振った枝とそれ以外の枝で異なる \omega が計算されている。 (旧世界ザルの共通祖先で \omega = 3.50573, そのほかの枝で \omega = 0.68581)
もう一か所は、パラメータ数と対数尤度が記述された行:
lnL(ntime: 11 np: 14): -904.636553 +0.000000
パラメータ数が14、対数尤度が -904.636553 であることを示している。
次に、帰無仮説として指定した枝の \omega を1で固定したモデルの 当てはまりを計算する。
.ctl の以下のパラメータを変更する:
outfile = branch_null * 何でもいいけど branch モデルの帰無仮説と分かるように
fix_omega = 1 * ω の値を固定
omega = 1 * ω=1
* 他のパラメータは変更なし
この .ctl を指定して codeml を実行:
同様に出力ファイル branch_null を見てみる。 各枝の \omega の値は:
w ratios as labels for TreeView:
(Hsa_Human #0.685577 , Hla_gibbon #0.685577 , (((Cgu/Can_colobus #0.685577 , Pne_langur #0.685577 ) #1 , Mmu_rhesus #0.685577 ) #0.685577 , (Ssc_squirrelM #0.685577 , Cja_marmoset #0.685577 ) #0.685577 ) #0.685577 );
パラメータ数と対数尤度は:
$ grep "lnL" branch_null
lnL(ntime: 11 np: 13): -905.484183 +0.000000
まとめると次のようになった:
| 対立仮説 |
3.50573 |
14 |
-904.636553 |
| 帰無仮説 |
1 |
13 |
-905.484183 |
最後に尤度比検定を行って、 「注目する枝の \omega が他の枝より高いようだけど、 これは選択の緩和 (\omega=1) じゃなくて 正の自然選択 (\omega>1) だ。」 ということを統計的に主張できるかどうか確かめる。
尤度比検定の方法はいくつかあるけど、R とか Python 使うのがいいんじゃないだろうか。 (エクセルとかでもできるらしい。)
Python の尤度比検定の関数を使って検定する:
from scipy.stats import chi2
alt_lnL = -904.636553 # 対立仮説の対数尤度
null_lnL = -905.484183 # 帰無の対数尤度
lr_stat = 2 * (alt_lnL - null_lnL)
alt_np = 14 # 対立仮説のパラメータ数
null_np = 13 # 帰無仮説のパラメータ数
df = alt_np - null_np
p_val = chi2.sf(lr_stat, df)
print(p_val)
$ python3 lrt.py
0.19290903422911437
結果、尤度比検定の p-value は 0.1929 > 0.05 となった。 これを解釈すると、帰無仮説を棄却することができない = 正の自然選択とは主張できないとなる。
Branch-Site モデル
Site モデルや Branch モデルは枝全体やサイト全体で \omega を平均するため、 時に検出力が弱くなる。
つまり、ある枝で特定のサイトに正の自然選択が本当に働いていたとしても、 他のサイトや枝の \omega が小さければそれに引っ張られて \omega > 1 を検出できない。
そこで、Branch-Site モデルは特定の枝の特定のサイトに働いた自然選択を検出する。
使う系統樹は Branch モデルと同じ。
まず対立仮説として、#1 を振った枝で \omega > 1 のサイトがある、 というシナリオの .ctl を書く:
outfile = bs_alt * 何でもいいけど branch-site モデルの対立仮説と分かるように
model = 2 * 記号の有無で異なる ω を推定
NSsites = 2 * ω > 1 のサイトを仮定
fix_omega = 0 * ω の値を配列から推定
omega = 1 * 推定は ω=1 からスタート
* 他のパラメータは変更なし
この .ctl を指定して codeml を実行:
出力ファイル bs_alt を見てみる。 以下のような記述があるはず:
MLEs of dN/dS (w) for site classes (K=4)
site class 0 1 2a 2b
proportion 0.29433 0.35964 0.15574 0.19029
background w 0.00000 1.00000 0.00000 1.00000
foreground w 0.00000 1.00000 5.78754 5.78754
各サイトクラスは次のように解釈する:
- site class 0
-
#1 を振った枝もそれ以外も \omega < 1 であるサイト
-
今回は全サイトのうち約29%がこれにあたり、\omega = 0 の強い純化選択をうけている。
- site class 1
-
#1 を振った枝もそれ以外も \omega = 1 であるサイト
-
今回は全サイトのうち約36%がこれにあたる。
- site class 2a
-
#1 を振った枝で \omega > 1、 それ以外で \omega < 1 であるサイト
-
今回は全サイトのうち約16%がこれにあたり、
#1 の枝では \omega = 5.78754 の強い正の自然選択が働いている。
- site class 2b
-
#1 を振った枝で \omega > 1、 それ以外で \omega = 1 であるサイト
-
今回は全サイトのうち約19%がこれにあたる。
次は Branch モデルと同じく、パラメータ数と対数尤度が記述された行:
$ grep "lnL" bs_alt
lnL(ntime: 11 np: 16): -901.562791 +0.000000
さらに、Branch-Site モデルではどのサイトが正の自然選択を受けているかを示す記述がある:
Bayes Empirical Bayes (BEB) analysis (Yang, Wong & Nielsen 2005. Mol. Biol. Evol. 22:1107-1118)
Positive sites for foreground lineages Prob(w>1):
14 R 0.795
21 R 0.798
23 I 0.799
37 G 0.585
41 R 0.713
50 R 0.707
62 R 0.583
87 D 0.796
126 Q 0.699
BEB 法 (Bayes empirical Bayes法) により求められた、 そのサイトが \omega > 1 で正の選択下にある事後確率を示す。 この事後確率が 0.95 や 0.99 を超えていた場合に正の自然選択が働いたサイトとする論文をよく見る。
今回はそういうサイトはなさそう。
次に、帰無仮説として指定した枝のサイトの \omega を1で固定したモデルの当てはまりを計算する。
.ctl のパラメータを変更する:
outfile = bs_null * 何でもいいけど branch-site モデルの帰無仮説と分かるように
fix_omega = 1 * ω の値を固定
omega = 1 * ω=1
* 他のパラメータは変更なし
この .ctl を指定して codeml を実行:
同様に出力ファイル bs_null を見てみる。 サイトクラスを見ると #1 の枝の \omega が1になっている:
MLEs of dN/dS (w) for site classes (K=4)
site class 0 1 2a 2b
proportion 0.26177 0.32442 0.18479 0.22902
background w 0.00000 1.00000 0.00000 1.00000
foreground w 0.00000 1.00000 1.00000 1.00000
パラメータ数と対数尤度は:
$ grep "lnL" bs_null
lnL(ntime: 11 np: 15): -902.301501 +0.000000
最後に、同じく尤度比検定を行って、対立仮説が採択されるか確かめる:
from scipy.stats import chi2
alt_lnL = -901.562791 # 対立仮説の対数尤度
null_lnL = -902.301501 # 帰無の対数尤度
lr_stat = 2 * (alt_lnL - null_lnL)
alt_np = 16 # 対立仮説のパラメータ数
null_np = 15 # 帰無仮説のパラメータ数
df = alt_np - null_np
p_val = chi2.sf(lr_stat, df)
print(p_val)
$ python3 lrt.py
0.22417862391381319
結果、尤度比検定の p-value は 0.2242 > 0.05 となった。 Branch-Site モデルでも、 帰無仮説を棄却することができない = 正の自然選択とは主張できない となった。
