ABO遺伝子型の次世代分布をマルコフ推移でみる
ABO遺伝子型の次世代分布をマルコフ推移でみる
ABO遺伝子型(AA, AO, BB, BO, AB, OO)の次世代分布をマルコフ推移させて, マルコフ定常するのか発散するのかをシミュレーションしてみる.
きっかけとしては,日本だとBB型の人は3%しかいなく,単純に世代を繰り返していくと消失するのではと思ったので確認してみる.
オチとしては,ハーディ・ワインベルクの法則というものがあるので定常はする.
他に似たような計算をしている方(特に2番目の論文は同じ計算をしている.):
ABO血液型と遺伝子型
A, B, AB, O型と血液型を分ける一番有名な方法をABO血液型という. さらに遺伝子型に分けると,
- A型(AA, AO)
- B型(BB, BO)
- O型(OO)
- AB型(AB)
の6状態に分けることができる.
例として,親1(AA), 親2(AO)とすると, 子の遺伝子型は,親1親2のもつ遺伝子型の片方をランダムに遺伝するので$AA, AO$のどちらかになる. それ以外の遺伝子型は基本ならない(突然変異除く).
子の遺伝子型パターンを表で表すと,以下になる
親1 \ 親2 | A | O |
---|---|---|
A | AA | AO |
A | AA | AO |
AA型になる割合は$\frac{2}{4}=\frac{1}{2}$, AO型になる割合も同様に$\frac{2}{4}=\frac{1}{2}$ になる.
遺伝子型のマルコフ推移を考える
上記を踏まえて遺伝子型のマルコフ推移を考える.
遺伝子型の状態ベクトル
遺伝子型の状態ベクトルを以下のように定義する.
- AAである:$(1, 0, 0, 0,0,0)$
- AOである:$(0, 1, 0, 0,0,0)$
- BOである:$(0, 0, 1, 0,0,0)$
- BBである:$(0, 0, 0, 1,0,0)$
- ABである:$(0, 0, 0, 0,1, 0)$
- OOである:$(0, 0, 0, 0,0, 1)$
なお,すべての要素を足して1になるよう状態確率ベクトルとして拡張すると, $(\frac{1}{2}, \frac{1}{2}, 0, 0, 0, 0)$ はAAになる確率$\frac{1}{2}$, AOになる確率$\frac{1}{2}$と表現することができる.
状態遷移図を考える
まず,状態が多いので,各遺伝子型を起点にどの遺伝子型に遷移するのかを確認する. 上記にもあった遷移表を書くことで確認できる.
AAを起点にした場合
行を起点にした遺伝子型(AA),列をペアになる遺伝子型として遷移表を書くと以下になる.
\ | A | A | . | A | O | . | B | O | . | B | B | . | A | B | . | O | O |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
A | AA | AA | . | AA | AO | . | AB | AO | . | AB | AB | . | AA | AB | . | AO | AO |
A | AA | AA | . | AA | AO | . | AB | AO | . | AB | AB | . | AA | AB | . | AO | AO |
AAを起点にして,子がAA型になる割合を考える. AAという記号をAA型の割合とすれば,indexに世代を与えるとt+1世代目のAA型の割合は以下の式になる.
$$ AA_{t+1} = 1 \cdot AA_{t} + \frac{1}{2} AO_{t}+ \frac{1}{2} AB_{t} $$
状態遷移図で表すと,
他の遷移状態を式にすると,
これらの式を状態ベクトルと推移行列を使って表す.
遷移元がAAのときの推移行列は,6x6行列で以下になる. 各行がペア相手に該当, 各列が推移先t+1に該当している. ペア相手がAAで子がAAになる推移確率は(1行, 1列)に該当している. 各行の総和は1になる.
$$ T_{ \cdot | AA} = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ \frac{1}{2} & \frac{1}{2} & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{2} & 0 & 0 & \frac{1}{2} & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 \\ \frac{1}{2} & 0 & 0 & 0 & \frac{1}{2} & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 \end{pmatrix} $$
相手がAOの場合の子の状態ベクトルは以下の式で求まる.
AAを起点にマルコフ推移の極限をみる. (上記の状態遷移図の推移を繰り返したときをみる) チャップマン・コルモゴロフにより, n回推移したときの推移行列は,推移行列の累乗により求まる. この推移行列は,AAから全ペアをかけ合わせて1回遷移してn-1回目までAAをかけ合わせてどの状態に推移するかを表す行列になる.各行は最初にかけあわせたペアを表し,(2,1)要素はAAに最初にAOをかけあわせてn回AAをかけあわせて最終的にAAになる確率になる.
$$ T_{ n | AA } = T_{ \cdot | AA }^{ n } $$
pythonで20回推移した推移行列を計算する.
import scipy as sp from scipy import linalg T_AA = sp.matrix([ [1, 1/2, 0, 0, 1/2, 0], [0, 1/2, 1/2, 0, 0, 1], [0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 0, 0], [0, 0, 1/2, 1, 1/2, 0], [0, 0, 0, 0, 0, 0], ]).T linalg.fractional_matrix_power(T_AA, 20) > array([[1.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00], [9.99999046e-01, 9.53674316e-07, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00], [9.99998093e-01, 9.53674316e-07, 0.00000000e+00, 0.00000000e+00, 9.53674316e-07, 0.00000000e+00], [9.99998093e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 1.90734863e-06, 0.00000000e+00], [9.99999046e-01, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 9.53674316e-07, 0.00000000e+00], [9.99998093e-01, 1.90734863e-06, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00]])
1列目が全て1になる.つまり,AAを起点にしてAAをかけ合わせ続けるとAAばかりになっていく.(当たり前)
AOを起点にした場合
\ | A | A | . | A | O | . | B | O | . | B | B | . | A | B | . | O | O |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
A | AA | AA | . | AA | AO | . | AB | AO | . | AB | AB | . | AA | AB | . | AO | AO |
O | AO | AO | . | AO | OO | . | BO | OO | . | BO | BO | . | AO | BO | . | OO | OO |
BOを起点にした場合
BOを起点に
\ | A | A | . | A | O | . | B | O | . | B | B | . | A | B | . | O | O |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
B | AB | AB | . | AB | BO | . | BB | BO | . | BB | BB | . | AB | BB | . | BO | BO |
O | AO | AO | . | AO | OO | . | BO | OO | . | BO | BO | . | AO | BO | . | OO | OO |
BBを起点にした場合
\ | A | A | . | A | O | . | B | O | . | B | B | . | A | B | . | O | O |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
B | AB | AB | . | AB | BO | . | BB | BO | . | BB | BB | . | AB | BB | . | BO | BO |
B | AB | AB | . | AB | BO | . | BB | BO | . | BB | BB | . | AB | BB | . | BO | BO |
ABを起点にした場合
\ | A | A | . | A | O | . | B | O | . | B | B | . | A | B | . | O | O |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
A | AA | AA | . | AA | AO | . | AB | AO | . | AB | AB | . | AA | AB | . | AO | AO |
B | AB | AB | . | AB | BO | . | BB | BO | . | BB | BB | . | AB | BB | . | BO | BO |
OOを起点にした場合
\ | A | A | . | A | O | . | B | O | . | B | B | . | A | B | . | O | O |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
O | AO | AO | . | AO | OO | . | BO | OO | . | BO | BO | . | AO | BO | . | OO | OO |
O | AO | AO | . | AO | OO | . | BO | OO | . | BO | BO | . | AO | BO | . | OO | OO |
推移行列を一般化
各行を現在の状態,各列を次の状態になるように推移行列を作る. (1,1)要素はAAから次の世代がAAになる確率になる. これは,AAを起点にした場合の $AA_{t+1} = 1 \cdot AA_{t} + \frac{1}{2} AO_{t} + \frac{1}{2}AB_{t} $ の式に該当する.
推移行列は以下になる.各行と各列のはじめはラベルを表している. この推移行列の各要素の値は,固定値でなく状態ベクトルの値によって変わる.
$(AA, AO, BO, BB, AB, OO) = (1, 0 , 0, 0, 0, 0)$を代入すればAAを起点とした場合の推移行列と同じになる.
pythonで実装する場合,引数を状態ベクトルにし,推移行列を返す関数にすればいい.
def T_Mat(state_vec): AA, AO, BO, BB, AB, OO = state_vec return sp.matrix([[AA+1/2*AO+1/2*AB, 1/2*AO+1/2*BO+OO, 0, 0, 1/2*BO+BB+1/2*AB, 0], [1/2*AA+1/4*AO+1/4*AB, 1/2*AA+1/2*AO+1/4*BO+1/4*AB+1/2*OO, 1/4*BO+1/2*BB+1/4*AB, 0, 1/4*BO+1/2*BB+1/4*AB, 1/4*AO+1/4*BO+1/2*OO], [0, 1/2*AA+1/4*AO+1/4*AB, 1/4*AO+1/2*BO+1/2*BB+1/4*AB+1/2*OO, 1/4*BO+1/2*BB+1/4*AB, 1/2*AA+1/4*AO+1/4*AB, 1/4*AO+1/4*BO+1/2*OO], [0, 0, 1/2*AO+1/2*BO+OO, 1/2*BO+BB+1/2*AB, AA+1/2*AO+1/2*AB,0], [1/2*AA+1/4*AO+1/4*AB, 1/4*AO+1/4*BO+1/2*OO, 1/4*AO+1/4*BO+1/2*OO, 1/4*BO+1/2*BB+1/4*AB, 1/2*AA+1/4*AO+1/4*BO+1/2*BB+1/2*AB, 0], [0, AA+1/2*AO+1/2*AB, 1/2*BO+BB+1/2*AB, 0, 0, 1/2*AO+1/2*BO+OO]]).astype("float") T_Mat(sp.array([1, 0, 0, 0, 0, 0])) > matrix([[1. , 0. , 0. , 0. , 0. , 0. ], [0.5, 0.5, 0. , 0. , 0. , 0. ], [0. , 0.5, 0. , 0. , 0.5, 0. ], [0. , 0. , 0. , 0. , 1. , 0. ], [0.5, 0. , 0. , 0. , 0.5, 0. ], [0. , 1. , 0. , 0. , 0. , 0. ]])
次世代の状態確率ベクトルとその極限を求める
次の世代の遺伝子型の状態確率ベクトル$\vec{x}_{i+1}$を求める. 推移行列により以下の式で求まる.
$$ \vec{ x_{ i+1 } } = \vec{ x_{ i } } T_{ i } \\ \vec{ x_{ i+1 } }^{T} = T_{ i }^{T} \vec{ x_{i} }^{T} $$
遺伝子型の初期分布を現在の日本の分布にする. ググってみるとこの辺の数字が出る.
AA型 8%,AO型 31%,BO型 19%,BB型 3%,AB型 10%,OO型 29%(ソースになるような出典なし)
$$ \vec{x_0} = (0.08, 0.31, 0.19, 0.03, 0.1, 0.29) $$
これを初期状態ベクトルとして遷移させて次世代がどのように変化していくかをみる. pythonで書くと以下のループを回していけばいい.
vec_x = sp.array([0.08, 0.31, 0.19, 0.03, 0.1, 0.29]) for i in sp.arange(1, 10+1): vec_x = vec_x.dot(T_Mat(vec_x)).A1 print(i, vec_x) > 1 [0.081225 0.3078 0.189 0.030625 0.09975 0.2916 ] 2 [0.081225 0.3078 0.189 0.030625 0.09975 0.2916 ] 3 [0.081225 0.3078 0.189 0.030625 0.09975 0.2916 ] 4 [0.081225 0.3078 0.189 0.030625 0.09975 0.2916 ] 5 [0.081225 0.3078 0.189 0.030625 0.09975 0.2916 ] 6 [0.081225 0.3078 0.189 0.030625 0.09975 0.2916 ] 7 [0.081225 0.3078 0.189 0.030625 0.09975 0.2916 ] 8 [0.081225 0.3078 0.189 0.030625 0.09975 0.2916 ] 9 [0.081225 0.3078 0.189 0.030625 0.09975 0.2916 ] 10 [0.081225 0.3078 0.189 0.030625 0.09975 0.2916 ]
結果として,1世代推移させるだけほぼ変わらず定常になる.
ハーディー・ワインベルクの法則(平衡)
血液型のような親から2者択一のように決まる遺伝子の分布は, ハーディー・ワインベルクの法則が成り立つ. これは次のような条件で今の世代と次世代の分布が変わらないことを示している.
法則が成立するための条件 ハーディー・ワインベルクの法則が成立するためには、当該個体群が以下の条件を全て満たしている必要がある。
分布が変わらないことは以下の式で示される. wikipediaの内容を丁寧に展開してみる.
$AA:Aa:aa=p2:2pp':p'^2$となる。遺伝子型がAAとなるのは、遺伝子プールから集めた2個の対立遺伝子が両方ともAだった場合で、そうなる確率はp×p=p2 この次世代集団のA遺伝子の遺伝子頻度をq、a遺伝子の遺伝子頻度をq’とする。q算出の分母となる遺伝子プール内の遺伝子総数は、各個体がAあるいはa遺伝子のいずれかを合計2個ずつ持つので、 $(p2+2pp'+p'^2)\times 2$ となる。分子を与える遺伝子プール内のA遺伝子の総数は、遺伝子型AAの個体が2個ずつ、Aaの個体が1個ずつのA遺伝子を持つので、$p2×2+2pp'$となる。 従って、$q=(p2×2+2pp')/{(p2+2pp'+p'^2)×2}=p$ となる(q'についても同様にq’=p'となる)。
i+1世代の遺伝子総数(各A, aを足した総数)は,i世代の遺伝子総数をN_iとしたとき $2N_{i}(p2+2pp'+p'^2)$
Aの総数は,$p^{2}$はAを2個持つので,$(p^{ 2 } \times 2 + 2pp')N_{ i }$
次世代がAを持つ割合(確率)は,
前世代のAをもつ割合と同じになる.