000_InoueSlidesSlide_5NovAll2013.indd

Size: px

Start display at page:

Download "000_InoueSlidesSlide_5NovAll2013.indd"

きのこいとえ
7 years ago
Views:

1 トピック 2014 年 5 月井上潤 I. 実データを用いた研究従来の分岐年代推定ベイズ分岐年代推定 II. MCMCTREE の理論的背景と実習理論的背景実習 23

2 ベイズ年代推定プログラム MCMCTREE の特徴コマンドラインで操作. (Windows: コマンドプロンプト ; Mac: ターミナル ). 樹形は固定. 化石記録は事前に調査. 二段階の解析からなる. (1: 樹長の分散共分散行列推定, 2: ベイズ解析 ). 塩基, アミノ酸配列の両方で解析が可能 ( 塩基で開発 ). ゲノムレベルの巨大データにも対応. 配列データにある程度の変異が必要. 24

3 統計的データ解析における二つの考え方古典的あるいは頻度主義者 (classical or frequentist) 確率 : 真のあるいは仮想的な母集団から繰り返し抽出する中で, その事象が生起する期待頻度. パラメータ : 未知の定数であり分布をもたない. ベイズ主義者 (Bayesian) 確率 : 母集団からのランダムな抽出における頻度とは解釈されず, パラメータの不確実性を表す. パラメータ : 確率分布する.-> 事前分布, 事後分布 25 Yang (2006) P 138.

4 事前分布と事後分布とは ( 例 : MCMCTREE の年代制約 ) パラメータがどのような値をとるかという見込みの大きさを確率分布として表したもの. 事前分布 (Prior) 解析前に設定した一様分布ヒト配列を用いない解析で得られた分布ゼブライトヨフグ事後分布 (Posterior) 配列を用いた解析で得られた分布 (Ma) ベイズ推定では事後分布の形を推定するのが目的.

5 MCMC アルゴリズムによるベイズ法の変革事後確率の計算に高次元の積分が要求されるため, 近年までベイズ法は用いられなかった 1990 年代に入り,Markov chain Monte Carlo (MCMC) アルゴリズムの導入とコンピュータの発達により解決. これによって, 最尤法では取り扱うことの出来ない複雑でパラメータの多いモデルを扱うことが可能になった. パタメータの多い年代推定の解析に適している. 27 Yang (2006)

6 分岐年代の事後分布を推定する式 f (t, r, θ X) パラメータの事後分布 = 尤度 X: 配列データ. t: 分岐年代の集合. r: 速度の集合. 事前分布 f (X t, r, θ) f (r t, θ) f (t θ) f (θ) f (X) データの周辺確率 : 事後分布の全積分が 1 になるための定数. t, r, θ についての高次元積分から得られる. 直接計算は困難. θ : 置換モデル中のパラメータや, t や r に関する事前分布のパラメータの集合. 分母の計算は困難 => 何百万回の MCMC サンプリングを行い, 同時事後分布 ( 左辺 ) から標本を生成 ( 次のページ ). 28 Yang (2006) P235

7 MCMC アルゴリズムによる計算の高速化 f (t, r, θ X) パラメータの事後確率 = 尤度 f (X t, r, θ) f (r t, θ) f (t θ) f (θ) f (X) データの周辺確率事前分布 29 MCMC アルゴリズム R = 採択比 : 新しい状態を採択するか否か決める min 1, 新しい状態 f (X t*, r*) f (r* t*, θ*) f (t* θ*) f (θ*) f (X t, r) f (r t, θ) f (t θ) f (θ) 現在の状態分母の計算は行わず, 繰り返し計算ごとに得られた分子の値だけを比較して事後確率の増減を見る. 分子だけを比較 q (θ, t, r θ*, t*, r*) q (θ*, t*, r* θ, t, r) 提案比 : 新しい状態の提案に偏りがあってもこれで修正. Yang (2006) P155

8 MCMCTREE の近似尤度計算 (Approximate Likelihood Calculation) f (t, r, θ X) パラメータの事後確率 = 尤度 => 近似尤度計算事前分布 f (X t, r, θ) f (r t, θ) f (t θ) f (θ) f (X) データの周辺確率 MCMCTREE は多変量正規分布を用いた近似尤度計算によって尤度推定をスピード化. 30 このため MCMCTREE の解析は, 1. 最尤推定値における尤度関数の勾配ベクトルと分散共分散行列の推定 2. 多変量正規分布を用いた MCMC アルゴリズムの実行という 2 段階からなる. Thorne et al. (1998) Yang (2006) PAML マニュアル

9 化石制約の選び方条件 : 化石の系統的位置がわかっている. 地層 1 地層 2 地層 3 地層 4 現在分岐 A の上限 : 化石なし A 分岐 B の上限 : 姉妹群最古の化石 Clade c Clade b 年代推定に使える化石その他の化石 B C 分岐 C の上限 : ステム最古の化石 Clade a Modified from Benton and Donoghue (2007) 分岐 A の下限分岐 B の下限分岐 C の下限 31 - 下限制約の設定は難しくない. - 上限制約の設定には, より詳細な取り決めが必要.

10 33.7 分子年代推定用の化石記録 Sponge Anemone Hydra Nematode Daphnia Louse Aphid Rhodnius Benton et al. (2009) Beetle Mosquito Benton, M. J., P. C. J. Donoghue, and R. J. Asher Calibrating and constraining molecular clocks. Pages in Timetree of Life (S. B. Hedges, and S. Kumar, eds.). Oxford University Press. Fruit Fly Wasp Honey bee Limpet Sea Hare Biomphalaria Capitella Leech Sea Urchin Amphioxus Tunicates Vertebrates 脊椎動物根幹の制約は化石記録が整備されつつある. 他は依然 Fossil Records II (Benton 1993) に頼らざるを得ない. Lizard Zebra Fish Zebrafinch Chicken Crocodile Toad Emu Medaka Stickleback Tetraodon Opossum Platypus Takifugo Elephant Shark Lamprey Kangaroo Tenrec 61.5 Elephant 239 Armadillo 48.4 Invertebrates Cow Sheep Dolphin Pig 62.5 Horse Cat Dog Bat 61.5 Shrew 61.5 Hedgehog Rabbit Pika Human Squirrel Guinea Pig 10.4 Rat Mouse Tree Shrew Bushbaby Macaque Orangutan Chimp Neanderthal Lemur Marmoset Gibbon 32

11 トピック I. 実データを用いた研究従来の分岐年代推定ベイズ分岐年代推定 II. MCMCTREE の理論的背景と実習理論的背景実習 33

12 関連ファイルのダウンロード mcmctree 簡単な例題で Google 検索 34 デスクトップに例題とスライドをダウンロード

13 解析プログラムの準備エディタ Win:? Mac: TextWrangler etc. 系統樹の図示 FigTree 分岐年代推定 MCMCTREE (PAML で検索 : paml4.x.tgz) Mac: terminal に入る. cd into paml4.x/src 以下をタイプ. make -f Makefile プログラムをコピー. baseml mcmctree Win: 次のスライド (2 枚 ) を参照. => make がない場合は, Xcode をインストール. => mcmctreeexamplevert6/examples ディレクトリにコピー. 35

14 ディレクトリ設定 -1 (Windows).tar 解凍プログラムのインストール 1. 解凍プログラム Lhaplus ( フリー ) を Google 検索して, ダウンロード & インストール. 2. ダウンロードした mcmctreeexamplevert6.tar を,Lhaplus エイリアス ( デスクトップ ) にドラッグ & ドロップ. 得られたファイルにまだ.tar が着いていたら, この.tar ファイルを再度 Lhaplus にドラッグ & ドロップ. 拡張子の表示以下のチェックを外す整理 > フォルダと検索のオプション > 表示 > 登録されている拡張子は表示しないファイルのアクセス権を変更 mcmctreeexamplevert6/examples ディレクトリを右クリックし以下のチェックを外す. プロパティ > 読み取り専用 36

15 ディレクトリ設定 -2 (Windows) 試しに変更エディタで開く試しにエディタで開く.baseml.ctl を右クリックして以下を選択. SAKURA ( など ) で開くコマンドプロント : ディレクトリの移動 mcmctreeexamplevert6/examples ディレクトリに入る. コマンドプロンプト ( プログラム > アクセサリ > コマンドプロンプト ) に以下のコマンドを入力. cd Desktop cd mcmctreeexamplevert6/examples dir 37 baseml を走らせる 1. pamlx.x/bin/ にある以下のプログラムを mcmctreeexamplevert6/examples にコピー baseml.exe, mcmctree.exe 2. コマンドプロンプトに以下のコマンドを入力. baseml

16 MCMCTREE 解析の全体像 1. 勾配ベクトルと分散共分散行列の推定 Vertebrate6.phy 塩基配列データ MCMCTREE mcmctree.ctl Vertebrate6.tre 種の系統樹 BASEML out.bv 樹長の勾配と分散共分散行列 2. MCMC アルゴリズムの実行 MCMCTREE Vertebrate6.phy 塩基配列データ ( 実際は使われない ) mcmctree.ctl Vertebrate6.tre 種の系統樹 in.bv (out.bv) 勾配ベクトル分散共分散行列 38 out_usedata2 分岐ごとの進化速度と分岐年代

17 MCMCTREE を用いた分岐年代推定 ( 近似尤度計算 ) 0. 準備 => 進化速度の事前確率を推定 1. 勾配ベクトルと分散共分散行列の推定 (usedata=3) 2. MCMC アルゴリズムの実行 (usedata=2) 39

18 0: mcmctree.ctl ファイル進化速度 ( 遺伝子 [partition] 全体 ) の事前確率 seed = -1 seqfile = Vertebrate6.phy treefile = Vertebrate6.tre outfile = out_usedata3 ndata = 1 usedata = 3 * 0: no data; 1:seq like; 2:normal a clock = 2 * 1: global clock; 2: independent ra RootAge = '<10' * constraint on root age, used if no model = 7 * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85 alpha = 0.5 * alpha for gamma rates at sites ncatg = 5 * No. categories in discrete gamma cleandata = 0 * remove sites with ambiguity data (1: BDparas = kappa_gamma = 6 2 alpha_gamma = 1 1 * birth, death, sampling * gamma prior for kappa * gamma prior for alpha rgene_gamma = * gamma prior for rate for genes sigma2_gamma = * gamma prior for sigma^2 (for finetune = 0: * times, rat * finetune = * times, rates print = 1 burnin = sampfreq = 50 nsample = *** Note: Make your window wider (100 columns) when running 40

19 0: 進化速度 ( 全体 ) の事前分布 (rgene_gamma) を計算 Vertebrate6.phy 塩基配列データ Vertebrate6PE.tre 1 較正点付き系統樹 baseml.ctl BASEML パラメータ推定 mlb 進化速度 BASEML を用いて, 全体の進化速度を概算. 41

20 0: baseml.ctl ファイル配列ファイルアウトファイル tree ファイル進化モデル ( ここでは GTR) seqfile = Vertebrate6.phy * sequence data file name outfile = mlb * main result file treefile = Vertebrate6PE.tre * tree structure file n noisy = 3 * 0,1,2,3: how much rubbish on the scr verbose = 1 * 1: detailed output, 0: concise outpu runmode = 0 * 0: user tree; 1: semi-automatic; 2 * 3: StepwiseAddition; (4,5):Perturbat model = 7 Mgene = 0 fix_kappa = 0 kappa = 2 * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85 * 0:rates, 1:separate; 2:diff pi, 3:di * initial or given kappa fix_alpha = 0 alpha = 0.5 * initial or given alpha, 0:infinity Malpha = 0 * 1: different alpha's for genes, 0: o ncatg = 5 * # of categories in the dg, AdG, or n 42 分子時計の一定性を仮定したモデルを選択 fix_rho = 1 rho = 0. * initial or given rho, 0:no correla npark = 0 * rate-class models. 1:rK, 2:rK&fK, 3: * clock = 1 * 0: no clock, unrooted tree, 1: clock nhomo = 1 * 0 & 1: homogeneous, 2: kappa's, 3: N getse = 1 * 0: don't want them, 1: want S.E.s of RateAncestor = 0 cleandata = 0 * (1/0): rates (alpha>0) or ancestral remove sites with ambiguity data (1:y

21 0: シーケンスファイル Vertebrate6.phy 配列データ Carp ---ATGGCA---AGCCTACGAAAAACACACCCTCTCATTAAAATCGCTAACGACGCACTAGTTGA Fugu ---ATGGCC---AGCCTACGCAAAACCCACCCCCTACTAAAAATCGTAAACGACATAGTAATTGA Frog ---ATGGCACCCAACATCCGTAAATCTCATCCATTAATTAAAATTATTAATAATTCTTTCATTGA Human ATGACCCCAATACGCAAAACTAACCCCCTAATAAAATTAATTAACCACTCATTCATCGA Shark ---ATGGCCACAAACATCCGAAAAACTCACCCTCTACTAAAAATTGTAAATCATGCCCTAATTGA Stingray ATGACTACTACAAACACCCGCAAAACCCATCCCCTATTCAAAATTATCAATAACTCCGTAATTGA Phylip 形式 43

22 0: tree ファイル Vertebrate6PE.tre 1 較正点付き系統樹 6 1 (((Carp,Fugu),(Frog,Human)),(Shark,Stingray))'@4.5'; Newick format FigTree で開く Carp Fugu Shark : おおよその分岐年代.Benton et al. (2009) など化石論文を参照. 注意 :4.5 は 450 Mya を意味する.MCMCTREE は 1-10 の値を用いた計算に適している. 44 全体の進化速度をラフに推定. =>MCMC アルゴリズムの事前確率に用いる.

23 0: baseml を走らせるコマンドプロンプト ( あるいはターミナルから ) baseml と入力してリターンを押す. [JunInoue:12tr]$ baseml BASEML in paml version 4.4d, March 2011 ns = 6 ls = 1146 Reading sequences, sequential format.. Reading seq # 1: Carp h-m-p C /14 30 h-m-p C Out... lnl = Calculating SE's Time used: 0:00 [JunInoue:12tr]$ 45 アウトファイル (mlb という名前になっている ). Substitution rate is per time unit

24 1: MCMCTREE を用いた分岐年代推定 ( 近似尤度計算 ) 0. 準備 1. 勾配ベクトルと分散共分散行列の推定 (usedata=3) 2. MCMC アルゴリズムの実行 (usedata=2) 46

25 1. 勾配ベクトルと分散共分散行列の推定 Vertebrate6.phy 塩基配列データ Vertebrate6.tre 制約付き系統樹 MCMCTREE mcmctree.ctl BASEML out_usedata2 out.bv 枝長の勾配と分散共分散行列枝長の最尤推定値 (maximum likelihood estimates; MLEs) を計算し, 枝長の勾配と分散共分散行列を求める. 47

26 1: mcmctree.ctl ファイル - usedata3 配列ファイル tree ファイルアウトファイル勾配, 分散共分散推定モードを指定進化モデル ( ここでは GTR) seed = -1 seqfile = Vertebrate6.phy treefile = Vertebrate6.tre outfile = out_usedata3 ndata = 1 usedata = 3 * 0: no data; 1:seq like; 2:normal a clock = 2 * 1: global clock; 2: independent ra RootAge = '<10' * constraint on root age, used if no model = 7 * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85 alpha = 0.5 * alpha for gamma rates at sites ncatg = 5 * No. categories in discrete gamma cleandata = 0 * remove sites with ambiguity data (1: BDparas = kappa_gamma = 6 2 alpha_gamma = 1 1 * birth, death, sampling * gamma prior for kappa * gamma prior for alpha rgene_gamma = * gamma prior for rate for genes sigma2_gamma = * gamma prior for sigma^2 (for finetune = 0: * times, rat * finetune = * times, rates print = 1 burnin = sampfreq = 50 nsample = *** Note: Make your window wider (100 columns) when running 48

27 1: tree ファイル Vertebrate6.tre 制約付き系統樹 6 1 Newick 形式 (((Carp,Fugu)'L(1.5,0.1,1.0,1e-300)',(Frog,Human)'B(3.3, 3.5)')'B(4.16, 4.22)',(Shark,Stingray)'L(1.9,0.1,1.0,1e-300)')'B(4.22, 4.63)'; B(4.22, 4.63) FigTree で開く Carp L(1.5,0.1,1.0,1e-300) Fugu B(4.16, 4.22) Frog B(3.3, 3.5) Human Shark L(1.9,0.1,1.0,1e-300) Stingray : 化石制約 ( 他の論文から調べる ) Carp Fugu Frog Human Shark Stingray B(4.22, 4.63) ペア制約 ( 上限 4.63, 下限 4.22) L(1.5,0.1,1.0,1e-300) 49 下限制約 (1.5, p=0.1, c=1.0, hard bound)

28 1: MCMCTREE (usedata=3) を走らせるコマンドプロンプトから mcmctree と入力してリターンを押す. [JunInoue:12tr]$ mcmctree MCMCTREE in paml version 4.4c, August 2010 Reading options from mcmctree.ctl.. Reading master tree. (((Carp, Fugu), (Frog, Human)), (Shark, Stingray)); h-m-p C /15 26 h-m-p C /15 27 h-m-p C Out... lnl = Calculating SE's Time used: 0:00 [JunInoue:12tr]$ 50 枝長, 勾配ベクトル, 分散共分散行列の書かれた out.bv ファイルが得られる ( 後述 ).

29 MCMCTREE を用いた分岐年代推定 ( 近似尤度計算 ) 0. 準備 1. 勾配ベクトルと分散共分散行列の推定 (usedata=3) 2. MCMC アルゴリズムの実行 (usedata=2) 51

30 2. MCMC アルゴリズム Vertebrate6.phy 塩基配列データ ( 実際は使われない ) Vertebrate6.tre 制約付き系統樹 MCMCTREE mcmctree.ctl out_usedata2 mcmc.out FigTree.tre in.bv (out.bv) 勾配ベクトルと分散共分散行列ベイズ法によって枝ごとの進化速度と年代を推定 52

31 2: in.bv ファイル (out.bv ファイルの名前を変更 ) 種数枝長付き最尤樹枝長ベクトル対数尤度の勾配ベクトル 1 階微分 in.bv ファイル (out.bv:usedata = 3 解析より得られる ) 対数尤度の 2 階微分ヘッセ行列 ( 枝長の分散共分散行列 ) 枝長の MLE を平均とし, 枝長の分散共分散行列を用いた多変量正規分布が得られる. 尤度を繰り返し計算ごとに近似注意 : 配列データは使われない. 53 Thorne, 1998, p1650 左下 Yang 2006 p22, MCT manual.

32 2: mcmctree.ctl ファイル - usedata2 54 MCMC アルゴリズムモードを指定進化速度 ( 遺伝子 [partition] 全体 ) の事前分布枝間の進化速度変動を調節する事前分布 finetune (1 は自動調節 ) 繰り返し計算で用いられる, 各パラメータの step length MCMC 繰り返し計算 2 回以上別々に計算して同じぐらいの推定値が出れば良い. seed = -1 seqfile = Vertebrate6.phy treefile = Vertebrate6.tre outfile = out_usedata2 ndata = 1 usedata = 2 * 0: no data; 1:seq like; 2:normal a clock = 2 * 1: global clock; 2: independent ra RootAge = '<10' * constraint on root age, used if no model = 7 * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85 alpha = 0.5 * alpha for gamma rates at sites ncatg = 5 * No. categories in discrete gamma cleandata = 0 * remove sites with ambiguity data (1: BDparas = kappa_gamma = 6 2 alpha_gamma = 1 1 * birth, death, sampling * gamma prior for kappa * gamma prior for alpha rgene_gamma = * gamma prior for rate for genes sigma2_gamma = * gamma prior for sigma^2 (for finetune = 1: * times, rat * finetune = * times, rates print = 1 burnin = sampfreq = 50 nsample = *** Note: Make your window wider (100 columns) when running

33 rgene_gamma = 遺伝子 (partition) 全体の進化速度 μ( ミュー ) の事前分布を設定するパタメータ. μ はガンマ分布すると仮定. 平均 (m), 標準偏差 (s), 形状パラメータ (α), 尺度パラメータ (β) の関係 : α = (m/s) 2 β = m/s Substitution rate is per time unit (0: 準備で得た mlb ファイル ) から m = 0.08, s=0.08 とすると, α = (0.08/0.08) 2 = 1 β = 0.08/ =12.5 MCMCTREE マニュアル

34 sigma2_gamma = 枝間における進化速度の変動を特定 σ 2 はガンマ分布すると仮定. 大きければ変動が激しい. パラメータ設定の根拠はあまり明確になっていない. =>MULTIDIVTIME のマニュアル. rttm*brownmean = 1 rttm: 根幹の分岐年代. brownmean: σ 2 の平均 (m). rttm = 4.5 とすると, m = 1/4.5 m = s= 1/4.5 とすると, α = (m/s) 2 = 1 β = m/s 2 =4.5 t2 t1 56

35 2: MCMCTREE (usedata=2) を走らせる [JunInoue:12tr]$ mcmctree MCMCTREE in paml version 4.4c, August 2010 Reading options from mcmctree.ctl.. Reading master tree. (((Carp, Fugu), (Frog, Human)), (Shark, Stingray));... r_n (0.0392, ) lnl ( , ) time prior: Birth-Death-Sampling rate prior: Log-Normal Time used: 0:14 [JunInoue:12tr]$ out_usedata2 mcmc.out FigTree.tre が得られる. 57

36 2: finetune を調節する (finetune = 0, 手動調節の場合 ) 採択率 (acceptance proportion ) はの間になるように finetune を調節 [JunInoue:12tr]$ mcmctree MCMCTREE in paml version 4.4c, August Starting MCMC (np = 17)... finetune steps (time rate mixing para RatePara?): paras: 5 times, 1 mu, 1 sigma2 (& rates, kappa, alpha) -10% % % % % % % % ^C [JunInoue:12tr]$ 大きいので, これを 0.01 にしてみる. finetune = * times, rates, mixing, para 58

37 2: mcmc.out ファイル (MCMC サンプリングが記載 ) t mcmc.out ファイル σ2 対数尤度繰り返し全体の進化速度 μ 枝ごとの進化速度 r 分岐年代 t, 置換速度 r, パラメータ θ のランダムな集合から開始. 個々の繰り返しで,t, r, θの変更を提案. 繰り返し計算の k 回ごとに標本抽出 : t, r,θを mcmc.out に保存 R = 計算終了時に平均と信頼区間を要約. 59 min ここでは 50 回ごとに標本抽出. 1, mcmctree.ctl ファイル f (X t*, r*) f (r* t*, θ*) f (t* θ*) f (θ*) f (X t, r) f (r t, θ) f (t θ) f (θ) q (θ, t, r θ*, t*, r*) q (θ*, t*, r* θ, t, r) Yang (2006)

38 2: アウトファイル (out_usedata2 ファイル ) Species tree for TreeView. Branch lengths = posterior mean (((1_Carp, 2_Fugu) 9, (3_Frog, 4_Human) 10 ) 8, (5_Shark, (((Carp: , Fugu: ): , (Frog: , (Shark: , Stingray: ): ); (((Carp: , Fugu: ) ' ': , ' ': ) ' ': , (Shark: 3. ' ': ) ' '; rategram locus 1: (((Carp: , Fugu: ): , (Frog: , (Shark: , Stingray: ): ); Posterior median mean & 95% credibility interval 推定年代と 95% 信頼区間 t_n (4.2264, ) (age of jeffnode t_n (4.1598, ) (age of jeffnode t_n (1.9057, ) (age of jeffnode t_n (3.3021, ) (age of jeffnode t_n (2.2630, ) (age of jeffnode mu (0.0631, ) sigma (0.0269, ) r_n (0.0416, ) r_n (0.0513, ) 60

39 2: FigTree.tre ファイル FigTree で作図. 61

40 参考文献岸野洋久,Jeffrey L. Thorne 分子進化速度のベイズ型階層モデル. 統計数理.50: Thorne, J. L., H. Kishino, and I. S. Painter Estimating the rate of evolution of the rate of molecular evolution. Molecular Biology and Evolution 15: Yang, Z. Yang, Z 分子系統学への統計的アプローチ計算分子進化学. 藤博幸, 加藤和貴, 大安裕美, 訳. 東京, 共立出版. 62

2011 年 6 月 9 日

2011 年 6 月 9 日簡易マニュアル : MCMCTREE (PAML) の近似尤度計算を用いた分岐年代推定 Jun Inoue, Mario dos Reis, and Z iheng Yang ( 井上訳 ) この簡易マニュアルでは,Inoue et al. (2010) のデータを例として年代推定プログラム MCMCTREE の解析手順を説明します.MCMCTREE で解析を行うには,Windows