血液型から遺伝子型を推定する

このページでは、血液型の基本と表現型の情報を使った遺伝型の確率的な推定方法について考える。

血液型についての基本

血液型にA, B, AB, O 型があることは誰でも知っている(ABO式の血液型）。しばしば性格などとの関連性が疑似科学の例として取り上げられたりもする。 ある調べによると、日本の場合は人口のうち A型が約 38%、B型が 22%、O型が 31%、AB型が 9% を占めている。

血液型は、両親から受け継いだ A, B, Oの３つの遺伝子型の組み合わせによって決まることが分かっている（下表）。 ただし、血液の抗体検査によって分かるのは、血液型（A,B,AB,O）であって、遺伝子型の組み合わせの情報は完全には分からない。

遺伝子型	表現型（血液型）
AA	A
AO	A
BB	B
BO	B
AB	AB
OO	O

$X, Y, x, y \in \{A, B, O\}$ として、 父親の遺伝子型が $XY$、母親の遺伝子型が $xy$とすると、遺伝の法則が教えるところによれば、 子供には $Xx$, $Xy$, $Yx$, $Yy$ のひとつが等しい確率で受け継がれる。

ある集団（例えば日本）で、A, B, Oの遺伝子型の存在確率を、それぞれ、$p_a, p_b, p_o$ と表すことにする。 ここで、その定義から $p_a + p_b + p_o = 1$ である。 これらの遺伝子型が一様に分布しているとすると、例えば、平均的にAB型の子が生まれる確率は、Aを持つ父親とBを持つ母親の場合と、Bを持つ父親とAを持つ母親の２通りを考えて、$ 2 p_a p_b $となる。 すべての場合を表にまとめると

血液型	遺伝子型から生じる確率	日本での人口比率
A	${p_a}^2 + 2p_a p_o$	0.38
B	${p_b}^2 + 2p_b p_o$	0.22
AB	$2 p_a p_b$	0.09
O	${p_o}^2$	0.31

実際に得られているは血液型の割合であって、遺伝子型のそれではないので、データから$p_a, p_b, p_o$を推定するうまい方法は無いものだろうか。

ここで、 $$ {p_a}^2 + 2p_a p_o = 0.38 $$ などのように置いて、$p_a, p_b, p_o$の連立方程式を解こうとしても駄目である。 そもそも、未知数よりも条件式が多いので解は見つからない。

EMアルゴリズム

隠れた変数の導入

遺伝子型のパターン毎の人口を、AAの人が$n_{aa}$、 以下同様に $n_{ao}$、$n_{bb}$、 $n_{bo}$、$n_{ab}$、 $n_{oo}$ と表すことにしよう。 そして、全人口を $N = n_{aa} + n_{ao} + n_{bb} + n_{bo} + n_{ab} + n_{oo}$ とおく。 実際に知ることができるのはA,B,AB,Oの血液型ごとの人口であるので、遺伝子型ごとの人口は全く未知の（隠れた）変数と言える （ただし、この中で、$n_{aa}, n_{oo}$ は、血液型から知ることができるので、本当に未知なのは $n_{aa}, n_{ao}, n_{bb}, n_{bo}$）。

仮に、遺伝子型の比率 $p_a, p_b, p_o$ が分かっていたとして、そのとき パターン毎の人口が　$n_{aa}, n_{ao}, n_{bb}, n_{bo}, n_{ab}, n_{oo}$ となる確率を考えてみる。 Aの遺伝子は全体として$\eta_a = 2 n_{aa} + n_{ao} + n_{ab}$、Bは $\eta_b = 2 n_{bb} + n_{bo} + n_{ab}$、 Oは $\eta_o = 2 n_{oo} + n_{ao} + n_{bo}$ であること等々を考慮すると、 $$ P(n_{aa}, n_{ao}, n_{bb}, n_{bo}, n_{ab}, n_{oo} \,|\, p_a, p_b, p_o ) = {p_a}^{\eta_a} \times {p_b}^{\eta_b} \times {p_o}^{\eta_o} $$ となる。

隠れた変数が既知と仮定し尤度を最大化する

逆に、もしも $n_{aa}, n_{ao}, n_{bb}, n_{bo}, n_{ab}, n_{oo}$ が分かっていたときに いちばん尤もらしい遺伝子型の分布は、上記の $P$ を最大化するような $p_a, p_b, p_o$ である、と考えることができる。 確率の代わりにその対数（対数尤度）を考えて、 $$ \begin{eqnarray} L(p_a, p_b, p_o) & = & \log \, P(n_{aa}, n_{ao}, n_{bb}, n_{bo}, n_{ab}, n_{oo} \,|\, p_a, p_b, p_o ) \\ & = & \eta_a \log p_a + \eta_b \log p_b + \eta_o \log p_o \end{eqnarray} $$ を最大化してみよう。

$G(p_a, p_b, p_o) = p_a + p_b + p_o - 1 = 0$という拘束条件のもとで、$L$を最大化したいのであるから、ラグランジェの未定乗数法により $$ \begin{eqnarray} \frac{\partial}{\partial p_a} \left( L + \lambda G \right) = 0 \\ \frac{\partial}{\partial p_b} \left( L + \lambda G \right) = 0 \\ \frac{\partial}{\partial p_o} \left( L + \lambda G \right) = 0 \end{eqnarray} $$ となるような$p_a^*, p_b^*, p_o^*$を見つければよい。 実際に計算してみると、$\lambda=-N$ であり、 $$ p_a^* = \eta_a /N, \; p_b^* = \eta_b /N, \; p_o^* = \eta_o /N $$ が得られる。

隠れた変数を予想する

上記のように、隠れた変数（遺伝子型の分布）が既知の場合には、知りたい量（$p_a, p_b, p_o$）を最尤推定できるが、 そもそも遺伝子型の分布についての計測データは無いのだから、ここで話は終わりのように思える。 そこを何とかするために、既知の量である血液型の分布を使って、遺伝子型の分布の予想を立ててみる。

ある $p_a, p_b, p_o$ のときに、A, B, Oの血液型が生じる確率は表のとおりであった。 例えば、血液型がA型という条件の下で遺伝子型AAの生じる確率は、 $$ P(aa|A) = \frac{P(aa,A)}{P(A)} $$ となり、その数は、 $$ n_{aa} = n_A P(aa|A) = \frac{n_A {p_a}^2}{{p_a}^2 + 2p_a p_o} $$ と見積もれる。

同様に、 $$ \begin{eqnarray} n_{ao} = n_A P(ao|A) & = & \frac{2 n_A p_a p_o}{{p_a}^2 + 2p_a p_o} \\ n_{bb} = n_B P(bb|B) & = & \frac{n_B {p_b}^2}{{p_b}^2 + 2p_b p_o} \\ n_{bo} = n_B P(bo|B) & = & \frac{2 n_B p_b p_o}{{p_b}^2 + 2p_b p_o} \\ n_{ab} & = & n_{AB} \\ n_{oo} & = & n_{O} \end{eqnarray} $$ となる。このように、もし、$p_a, p_b, p_o$ がわかっていれば、既知の量 $n_A, n_B, n_O, n_{AB}$ を使って隠れた変数の期待値は求めることができる。

堂々巡りの解消：EMステップ

以上を要約すると、「もし遺伝子型の分布がわかっていれば、知りたいパラメータが最尤推定できる」、 こと、そして「もし知りたいパラメータが分かっていたら（既知の量も使いながら）遺伝子型の分布の期待値を計算できる」 という堂々巡りの構造になっている。卵が先か鶏が先か、といった感じである。

そこで、以下のようなアルゴリズムが考案されている：

1. 知りたいパラメータ $p_a, p_b, p_o$ に仮の値を設定する
2. 既知の量($n_A, n_B, n_O, n_{AB}$)を使って、隠れ変数の期待値を求める
3. 隠れ変数を使って、尤度を最大化するパラメータ $p_a^*, p_b^*, p_o^*$ を求める
4. $p_a, p_b, p_o$ と $p_a^*, p_b^*, p_o^*$ の差が大きければ、 $p_a^*, p_b^*, p_o^*$ を $p_a, p_b, p_o$ と置き直して2に戻る。 そうでなければ、繰り返しを終了する。

こうした反復を繰り返すと、パラメータの最尤推定値に収束することがわかっている（EMアルゴリズム）。 上記の2は予測ステップ（Eステップ）、3は最大化ステップ（Mステップ）と呼ばれている。

　練習：遺伝子型の比率の推定

このページで説明したEMアルゴリズムに従って、A,B,O遺伝子型の割合 $p_a, p_b, p_o$ を計算するPythonコードを作成し、 日本人の $p_a, p_b, p_o$ を推定してみなさい。