生命情報学 （３，４）
進化系統樹推定
阿久津 達也
京都大学 化学研究所
進化系統樹
種間（もしくは遺伝子間）の進化の関係を表す木
以前は形態的特徴をもとに構成
有根系統樹と無根系統樹
有根系統樹: 根（共通の祖先に対応）がある系統樹 無根系統樹: 根のない系統樹 いずれも葉にのみラベル（種に対応）がつく 有根系統樹 無根系統樹 系統樹を扱う際は、頂点でなく節点という用語を使う系統樹の同型性
同型な系統樹：本質 的に同じつながり方 （グラフとして同型） をした系統樹 (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法
(Unweighted pair group method using arithmetic averages)
l k j C C i ij l k kld
C
C
D
,|
||
|
1
| | | | | | | | j i j jl i il kl C C C D C D D UPGMA法の例
}
5
,
4
,
3
,
2
,
1
{
}
3
,
2
,
1
{
}
2
,
1
{
}
5
,
4
{
8 6 9 3 7 8 2 1 7 5 4 6
C
C
C
C
C
C
C
C
C
C
C
C
UPGMA法の正当性
分子時計＝進化速度一定性の仮定 枝長＝分子時計により刻まれた時間 分子時計仮説が成立 「任意の葉までの枝長の和が等しい」が任意節点について成立 ⇒UPGMA法は系統樹を正しく再構成 以下の例で、(a) は仮説を満たさないが、(b)は満たす近隣結合法（１）
UPGMA法では以下の例(a)で、2と3が最初に選ばれたため、 正しい系統樹を再構成できなかった 最初に、3と4、もしくは、1と2が選ばれれば良い つまり、兄弟関係（近隣関係）にある葉が選ばれれば良い 近隣結合法では（距離が加法性を満たすという仮定のもとで） 常に近隣関係にある葉を選ぶ近隣結合法（
2）
近隣結合法
距離が加法性を満たせば系統樹を正しく再構成
ただし、計算されるのは無根系統樹
近隣結合法（
3）
加法性： 節点間の距離 ＝ 節点間を結ぶパスの枝長の和 i と j の親を k とすると、葉 m に対し以下が成立)
(
2 1 ij jm im kmd
d
d
d
近隣結合法： アルゴリズム
M を葉集合
とし、
L=M
とする
|L|>2の間
、以下を繰り返す
L の中から
D
ijが最小となるペア (i,j) を選ぶ
新しい節点 k
を作り，L 中のすべての m について
とする
k を M に追加し，k と i, j を枝で結び
，枝長を
とする
i, j を L から削除し，k を L に追加
する
i, j 間に枝をはり
，その長さを d
ijとする
)
(
2 1 ij jm im kmd
d
d
d
ik ij jk j i ij ikd
r
r
d
d
d
d
21(
),
L k ik L i j i ij ijd
r
r
r
d
D
(
),
| |1 2ただし、
アウトグループ
近隣結合法： 無根系統樹しか計算されない 無根系統樹における根の決定法 最も長い葉の間の中点を根とする 対象とする種と離れていることが明らかな種（アウトグループ）を追加して 系統樹を作成し、アウトグループに接続する節点を根とする最節約法： 問題設定
最小コストの置換
で各配列を生成する系統樹を計算
２種類の探索が必要
系統樹の形状（トポロジー) ⇒ 難しい
最節約法： アルゴリズム
S
k(a)
:
節点 k に塩基 a を割り当てた時のコスト S(a,b)
:
a を b に置換するためのコスト 動的計画法により葉から根へという順で Sk(a) を計算 配列割り当ては各位置ごとに独立に計算可能))
,
(
)
(
(
min
))
,
(
)
(
(
min
)
(
a
S
b
S
a
b
S
b
S
a
b
S
j b i b k
otherwise
]
1
[
if
0
)
(
k ks
a
a
S
k
が内部節点の場合（i,j は k の子節点）k
が葉の場合最節約
法：
実行例
))
,
(
)
(
(
min
))
,
(
)
(
(
min
)
(
a
S
b
S
a
b
S
b
S
a
b
S
j b i b k
進化の確率モデル
P(b|a,t) :
長さ t の枝をたどることにより塩基 a が塩基 b に置換される確率 P(y|x,t) :
長さ t の枝をたどることにより配列 x が配列 y に置換される確率
u u ux
t
y
P
t
x
y
P
(
|
,
)
(
|
,
)
置換行列 （置換確率行列） ) , T | T ( ) , T | G ( ) , T | C ( ) , T | A ( ) , G | T ( ) , G | G ( ) , G | C ( ) , G | A ( ) , C | T ( ) , C | G ( ) , C | C ( ) , C | A ( ) , A | T ( ) , A | G ( ) , A | C ( ) , A | A ( ) ( t P t P t P t P t P t P t P t P t P t P t P t P t P t P t P t P t S 乗法性)
(
)
(
)
(
t
S
s
S
t
s
S
Jukes-Cantor行列（１）
置換速度行列を R とし、 を仮定 S(t) の乗法性より よって A,C,G,T が対等であるとして 対称性より よってt
R
I
t
S
(
)
3 3 3 3 R t R t S t S t R I t S t S t S t S t S t t S ) ( ) ( ) )( ( ) ( ) ( ) ( ) ( ) ( R t S dt t dS ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( t r t s t s t s t s t r t s t s t s t s t r t s t s t s t s t r t S ) ( ) ( ) ( ) ( 3 ) ( 3 ) ( t r t s dt t ds t s t r dt t dr Jukes-Cantor行列（２）
を に代入して、 これを解いて この S(t) が Jukes-Cantor行列 は十分時間が経つと A,C,G,T の割 合いが等しくなることを意味する 3 3 3 3 R R t S dt t dS ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( t r t s t s t s t s t r t s t s t s t s t r t s t s t s t s t r t S ) ( ) ( ) ( ) ( 3 ) ( 3 ) ( t r t s dt t ds t s t r dt t dr ) 1 ( ) ( ) 3 1 ( ) (t 14 e 4 t s t 14 e 4 t r
s
t
t
t
r
(
)
(
)
14if
TATAT(=x)が t 時間後にTTAAA(=y)に置
換する確率は
より
Jukes-Cantor行列の応用例
)
3
1
(
)
A,
|
A
(
)
T,
|
T
(
)
1
(
)
T,
|
A
(
)
A,
|
T
(
4 4 1 4 4 1 t te
t
P
t
P
e
t
P
t
P
3 4 2 4 5(
1
3
)
(
1
)
4
1
)
,
|
(
y
x
t
e
te
tP
系統樹の尤度（確率）
配列 s
1,…,s
nに対し、尤度最
大の
系統樹のトポロジー
T
と、枝長
t
iを計算
ただし、p(i) は i の親節点 この計算は難しく、様々な手
法が提案されている
最尤法による系統樹推定
2 2 1 ) ( , , 1 2 1,
,
|
,
)
(
)
(
|
,
)
(
1 2 1 n i i i p i s s n nT
t
P
s
P
s
s
t
s
s
P
n n
1塩基あたりの尤度
配列全体の尤度
尤度計算の例（１）
a u u a u ux
T
t
t
q
P
x
a
t
P
x
a
t
x
P
(
1,
2|
,
1,
2)
(
1|
,
1)
(
2|
,
2)
N i u ux
T
t
t
x
P
t
t
T
x
x
P
1 2 1 2 1 2 1 2 1)
,
,
|
,
(
)
,
,
|
,
(
CGのみからなる配列
n
1個は同じで n
2個は変化 （この例では、n
1=8, n
2=3）
尤度計算の例（
2）
G
CGGGCCGGCC
G
CCGGCCGCGC
)
e
1
(
)
2
(
)
,
,
|
G
C,
(
)
e
3
1
(
)
3
(
)
,
,
|
C
C,
(
) ( 4 -16 1 4 1 2 1 ) ( 4 -16 1 4 1 T A G C 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 2 1 t t t t t t t t t t t t t t t t t t t t t ts
s
r
s
s
r
t
t
T
P
s
s
r
r
s
s
q
s
s
q
s
s
q
r
r
q
t
t
T
P
2 2 1 1 2 1 2 1(
1
3
e
)
(
1
e
)
16
1
)
,
,
|
,
(
x
1x
2T
t
1t
2 n n -4 (t t ) n -4 (t t ) nP
よって、系統樹全体の尤度は、
尤度は、t1+t2 が一定なら変わらないことに注意 ⇒ 根の位置の任意性
系統樹の信頼性
推定された系統樹の信頼性は不明 ⇒信頼性評価の必要性 必ずしも系統樹全体でなくても、大きな分岐などの特定の特 徴の信頼性を評価できれば良い ⇒ブートストラップ法による評価 ある確率モデルのもとで、 ブートストラップ法による頻度≈ 特徴の事後確率 系統樹に対するブートストラップ法
1. もとのアラインメントから各カラムを重複を許してランダムに抽出 2. ステップ 1 で作ったアラインメントから系統樹を推定 3. ステップ 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)
lca
i(a,b)
: a,b をラベルにもつ葉の
共通の祖先
で最も
根から遠いもの
半順序
: 節点 u が節点 v の子孫 ⇔
u
v
(
u
v
)
LCAの例
動的計画法アルゴリズム
}
))
,
(
lca
)
,
(
(lca
|
)
,
(
{
)
,
(
}
))
,
(
lca
)
,
(
(lca
|
)
,
(
{
)
,
(
b
a
b
x
T
b
x
b
a
R
b
a
x
a
T
x
a
b
a
L
i i i i i i
)}
,
(
{
max
)}
,
(
{
max
)
,
(
1
)
,
(
) , ( ) , ( ) , ( ) , (W
a
c
W
c
b
b
a
W
a
a
W
b a R b c b a L c a
a≠b として L(a,b), R(a,b) を次のように定義
すべての a≠b に対して、W(a,b) を動的計画法で計算
ただし、 L(a,b), R(a,b) のいずれかが空の時は、W(a,b)=-∞
)}
,
(
{
max
|
|
B
W
a
b
b a
すると、
B はトレースバックで 計算可能動的計画法の例
B={1,3,5,6,8}5
)
8
,
5
(
)
3
,
1
(
)
8
,
1
(
3
)
8
,
8
(
)
6
,
5
(
)
8
,
5
(
2
)
3
,
3
(
)
1
,
1
(
)
3
,
1
(
)}
8
,
8
{(
)
8
,
5
(
)}
6
,
5
(
),
5
,
5
{(
)
8
,
5
(
)}
3
,
3
{(
)
3
,
1
(
)}
1
,
1
{(
)
3
,
1
(
W
W
W
W
W
W
W
W
W
R
L
R
L
下図の例に対する計算例（一部分）
時間計算量の解析
定理： 最大合致部分系統樹問題は 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困難
