進化系統樹推定

生命情報学（３，４）
進化系統樹推定
阿久津達也
京都大学化学研究所
バイオインフォマティクスセンター
進化系統樹
進化系統樹



種間（もしくは遺伝子間）の進化の関係を表す木
以前は形態的特徴をもとに構成
現在は配列情報をもとに構成
有根系統樹と無根系統樹



有根系統樹: 根（共通の祖先に対応）がある系統樹
無根系統樹: 根のない系統樹
いずれも葉にのみラベル（種に対応）がつく
有根系統樹
無根系統樹
系統樹を扱う際は、頂点でなく節点という用語を使う
進化系統樹の個数
系統樹の同型性



同型な系統樹：本質
的に同じつながり方
（グラフとして同型）
をした系統樹
(a)と(b)：同型
(a)と(c)：非同型
(d)と(e)：同型
(d)と(f)：非同型
系統樹の個数

葉の個数と節点の個数の関係（葉の個数＝ n）
有根系統樹： 2n-1 個の節点、 2n-2 本の枝
無根系統樹： 2n-2 個の節点、 2n-3 本の枝



無根系統樹と有根系統樹の関係


根は無根系統樹の任意の枝に設定可能
⇒ 有根系統樹の個数＝(2n-3)×無根系統樹の個数
枝が 2n-3 本の無根系統樹に1本の枝を追加
⇒ 2n-3 種類の異なる無根系統樹
#葉
#枝
（無根系統樹）
#無根系統樹
#有根系統樹
3
3
1
3
4
5
3
3×5
5
7
3×5
3×5×7
6
9
3×5×7
3×5×7×9
…
…
…
…
n
2n-3
(2n-5)!!
(2n-3)!!
系統樹
の個数
（例）


値は無根系統樹の任意の枝に設定可能
⇒ 有根系統樹の個数＝(2n-3)×無根系統樹の個数
枝が 2n-3 本の無根系統樹に1本の枝を追加
⇒ 2n-3 種類の異なる無根系統樹
UPGMA法
UPGMA法 (Unweighted pair group method using arithmetic averages)

距離行列法の一つ


距離行列法：葉のペアの距離データから系統樹を構成
アルゴリズム


各配列のみからなるクラスタを作成
クラスタが２個になるまで以下を繰り返す



クラスタ間距離が最小のクラスタどうしを併合
新しくできたクラスタと他のクラスタ間の距離を計算
クラスタ間距離
1
Dkl 
dij

| Ck || Cl | iCk , jCl

距離の性質
Dkl 
Dil | Ci |  D jl | C j |
| Ci |  | C j |
UPGMA法の例
C6  C4  C5  {4,5}  C7  C1  C2  {1,2} 
C8  C7  C3  {1,2,3}  C9  C6  C8  {1,2,3,4,5}
UPGMA法の正当性
分子時計＝進化速度一定性の仮定
 枝長＝分子時計により刻まれた時間
 分子時計仮説が成立 
「任意の葉までの枝長の和が等しい」が任意節点について成立
⇒UPGMA法は系統樹を正しく再構成
以下の例で、(a) は仮説を満たさないが、(b)は満たす

近隣結合法
近隣結合法（１）
UPGMA法では以下の例(a)で、2と3が最初に選ばれたため、
正しい系統樹を再構成できなかった
 最初に、3と4、もしくは、1と2が選ばれれば良い
 つまり、兄弟関係（近隣関係）にある葉が選ばれれば良い
 近隣結合法では（距離が加法性を満たすという仮定のもとで）
常に近隣関係にある葉を選ぶ

近隣結合法（2）

近隣結合法



距離が加法性を満たせば系統樹を正しく再構成
ただし、計算されるのは無根系統樹
加法性：節点間の距離＝節点間を結ぶパスの枝長の和
近隣結合法（3）


加法性：節点間の距離＝節点間を結ぶパスの枝長の和
i と j の親を k とすると、葉 m に対し以下が成立
d km  12 (d im  d jm  d ij )
近隣結合法：アルゴリズム


M を葉集合とし、L=M とする
|L|>2の間、以下を繰り返す





L の中から Dij が最小となるペア (i,j) を選ぶ
新しい節点 k を作り，L 中のすべての m について
d km  12 (d im  d jm  d ij ) とする
k を M に追加し，k と i, j を枝で結び，枝長を
d ik  12 (d ij  ri  rj ), d jk  d ij  d ik とする
i, j を L から削除し，k を L に追加する
i, j 間に枝をはり，その長さを dij とする
ただし、 Dij  d ij  (ri  rj ), ri  |L1|2  dik
kL
アウトグループ


近隣結合法：無根系統樹しか計算されない
無根系統樹における根の決定法


最も長い葉の間の中点を根とする
対象とする種と離れていることが明らかな種（アウトグループ）を追加して
系統樹を作成し、アウトグループに接続する接点を根とする
最節約法
最節約法：問題設定


最小コストの置換で各配列を生成する系統樹を計算
２種類の探索が必要


系統樹の形状（トポロジー)
⇒ 難しい
内部節点への配列の割り当て ⇒ 動的計画法で計算可能
最節約法：アルゴリズム




Sk(a) : 節点 k に塩基 a を割り当てた時のコスト
S(a,b) : a を b に置換するためのコスト
動的計画法により葉から根へという順で Sk(a) を計算
配列割り当ては各位置ごとに独立に計算可能
k が内部節点の場合（i,j は k の子節点）
S k (a)  min ( Si (b)  S (a, b))  min ( S j (b)  S (a, b))
b
k が葉の場合
 0 if a  sk [1]
S k (a)  
  otherwise
b
最節約
法：
実行例
S k (a)  min ( Si (b)  S (a, b))  min ( S j (b)  S (a, b))
b
b
進化の確率モデル
進化の確率モデル

P(b|a,t) : 長さ t の枝をたどることにより塩基 a が塩基 b
に置換される確率

P(y|x,t) : 長さ t の枝をたどることにより配列 x が配列 y
に置換される確率
P( y | x, t )   P( yu | xu , t )
u
置換行列（置換確率行列）
 P ( A | A, t ) P ( C | A , t ) P ( G | A, t ) P ( T | A, t ) 


 P(A | C, t ) P(C | C, t ) P(G | C, t ) P(T | C, t ) 
S (t )  
P ( A | G , t ) P (C | G , t ) P (G | G , t ) P ( T | G , t ) 


 P ( A | T, t ) P (C | T , t ) P (G | T, t ) P ( T | T, t ) 


乗法性
S (t ) S ( s )  S (t  s )
Jukes-Cantor行列（１）



置換速度行列を R とし、 S (t )  I  Rt を仮定
S(t) の乗法性より
S (t  t )  S (t )  S (t ) S (t )  S (t )
 S (t )( I  Rt )  S (t )  S (t ) Rt
よって

dS (t )
 S (t ) R
dt
A,C,G,T が対等であるとして

対称性より

よって
dr (t )
 3 r (t )  3 s (t )
dt
  3

 
R


 

 r (t )

 s (t )
S (t )  
s (t )

 s (t )



 

 3

 

 3
 



 3 
s(t ) s (t ) s (t ) 

r (t ) s (t ) s (t ) 
s(t ) r (t ) s (t ) 

s(t ) s (t ) r (t ) 
ds(t )
  r (t )   s (t )
dt
Jukes-Cantor行列（２）
  3

 
R


 

を


 

 3

 

 3
 



 3 
dS (t )
 S (t ) R
dt
 r (t )

 s (t )
S (t )  
s (t )

 s (t )

に代入して、
dr (t )
 3 r (t )  3 s (t )
dt
これを解いて
s(t ) s (t ) s (t ) 

r (t ) s (t ) s (t ) 
s(t ) r (t ) s (t ) 

s(t ) s (t ) r (t ) 
r (t )  14 (1  3e 4 t )
ds(t )
  r (t )   s (t )
dt
s(t )  14 (1  e 4 t )
この S(t) が Jukes-Cantor行列
r (t )  s(t )  14 if t   は十分時間が経つと A,C,G,T の割
合いが等しくなることを意味する
Jukes-Cantor行列の応用例

TATAT(=x)が t 時間後にTTAAA(=y)に置
換する確率は
P(T | A, t )  P(A | T, t )  14 (1  e 4 t )
より
P(T | T, t )  P(A | A, t )  14 (1  3e  4 t )
1
P( y | x, t )  5 (1  3e  4 t ) 2 (1  e  4 t )3
4
最尤法による系統樹推定

系統樹の尤度（確率）
P( s1 ,, sn | T , t ) 
 P( s
2n2
2 n 1
sn1 ,, s2 n1

配列 s1,…,sn に対し、尤度最
大の系統樹のトポロジー T
と、枝長 ti を計算


ただし、p(i) は i の親節点
この計算は難しく、様々な手
法が提案されている
)  P( si | s p (i ) , ti )
i 1
尤度計算の例（１）

1塩基あたりの尤度
P( xu1 , xu2 | T , t1 , t2 )   qa P( xu1 | a, t1 ) P( xu2 | a, t2 )
a

配列全体の尤度
N
P( x1 , x 2 | T , t1 , t 2 )   P( xu1 , xu2 | T , t1 , t 2 )
i 1
尤度計算の例（2）

CGのみからなる配列

n1個は同じで n2個は変化（この例では、n1=8, n2=3）
CCGGCCGCGCG
CGGGCCGGCC G
P(C, C | T , t1 , t 2 )  qC rt1 rt2  qG st1 st2  qA st1 st2  qT st1 st2
 14 (rt1 rt2  3st1 st2 )  161 (1  3e -4 (t1 t2 ) )
P(C, G | T , t1 , t2 )  (rt1 st2  st1 rt2  2st1 st2 ) 
1
4
1
16
(1  e
- 4 ( t1  t 2 )
)
よって、系統樹全体の尤度は、
1
- 4 ( t1  t 2 ) n1
- 4 ( t1  t 2 ) n2
1
2
P( x , x | T , t1 , t 2 )  n1  n2 (1  3e
) (1  e
)
16
尤度は、t1+t2 が一定なら変わらないことに注意
⇒ 根の位置の任意性
ブートストラップ法

系統樹の信頼性




推定された系統樹の信頼性は不明
⇒信頼性評価の必要性
必ずしも系統樹全体でなくても、大きな分岐などの特定の特
徴の信頼性を評価できれば良い
⇒ブートストラップ法による評価
ある確率モデルのもとで、
ブートストラップ法による頻度≈ 特徴の事後確率
系統樹に対するブートストラップ法
1.
2.
3.
もとのアラインメントから各カラムを重複を許してランダムに抽出
ステップ 1 で作ったアラインメントから系統樹を推定
ステップ 1,2 を数多く（例えば1000回）繰り返し、各特徴の頻度分布を
計算し、頻度により信頼性を評価
最大合致部分系統樹
最大合致部分系統樹 (maximum agreement subtree)


系統樹推定には様々な方法 ⇒ 複数の系統樹を統合
一部の種だけを考えることにより、複数の系統樹から共通の
系統樹を得る（⇒種の数の最大化）
入力: 同じラベル集合{1,2,…,n}を持つ有根系統樹 T1,T2,…,TN
出力: Ti|B がすべて同型となる要素数最大の B
Ti|B: B に含まれない葉とそれに接続する枝を削除し、さらに子が一つになっ
た内部節点を繰り返し縮約して得る
⇒ B={1,3,5,6,8}
最小共通祖先 (LCA: Least Common Ancestor)


lcai(a,b): a,b をラベルにもつ葉の共通の祖先で最も
根から遠いもの
半順序: 節点 u が節点 v の子孫 ⇔ u  v (u  v)
LCAの例
□: lcai(2,3)
○: lcai(4,7)
△： lcai(5,7)
動的計画法アルゴリズム
a≠b として L(a,b), R(a,b) を次のように定義
L(a, b)  { (a, x) | Ti (lca i (a, x)  lca i (a, b)) }
R(a, b)  { ( x, b) | Ti (lca i ( x, b)  lca i (a, b)) }
すべての a≠b に対して、W(a,b) を動的計画法で計算
W ( a, a )  1
W (a, b)  max {W (a, c)}  max {W (c, b)}
( a ,c )L ( a ,b )
( c ,b )R ( a ,b )
ただし、 L(a,b), R(a,b) のいずれかが空の時は、W(a,b)=-∞
すると、 | B | max {W (a, b)}
a b
B はトレースバックで
計算可能
動的計画法の例
下図の例に対する計算例（一部分）
L(1,3)  {(1,1)}
R (1,3)  {( 3,3)}
L(5,8)  {( 5,5), (5,6)} R (5,8)  {(8,8)}
W (1,3)  W (1,1)  W (3,3)  2
W (5,8)  W (5,6)  W (8,8)  3
W (1,8)  W (1,3)  W (5,8)  5
B={1,3,5,6,8}
時間計算量の解析
定理：最大合致部分系統樹問題は O(n3N) 時間で計算可能
証明: 正当性の証明は省略。
O(n2) 個の (a,b) と N 個の系統樹について、LCAを計算するの
にかかる時間の合計は O(n2N) 。
L(a,b)（およびR(a,b)）１個あたりの要素数は O(n)であり、あらかじ
め lcai(a,b)間の半順序関係の表を作っておけば、全体で O(n3N)
時間。
W(a,b) は全体で O(n3) 時間で計算可能。
子頂点の数が制約されない一般の場合は NP困難
まとめ

系統樹の個数


距離行列法




UPGMA法：最小距離のクラスタを結合
近隣結合法：近隣の節点を結合
最節約法: 置換コストが最小の系統樹を計算
進化の確率モデル:




無根系統樹で (2n-5)!! 有根系統樹で (2n-3)!!
Jukes-Cantor行列を用いて配列間の置換確率を計算
最尤法：確率が最大となる系統樹を推定
ブートストラップ法：系統樹の特徴の信頼性評価
最大合致部分系統樹：複数の系統樹の比較

Download Report