ランダー・グリーン・アルゴリズム（５）

こういうのはですね、途中で間をあけると飽きて終わりになっちゃうのでですね、書ききってしまおうかと。

エルストン・スチュワート・アルゴリズム

ランダー・グリーン・アルゴリズムの前段階である、エルストン・スチュワート・アルゴリズム *1について説明します。

前回、全家系の尤度を得ました。

$L = \prod P(Y) =\sum_{g_1 \in G_1} \sum_{g_2 \in G_2} \cdots \sum_{g_m \in G_m} \prod P(Y_i|g_i) \prod_{founder} P(g_i) \prod_{non founder} P(g_o|g_p, g_m)$

しかしこれは、各個人すべての可能な遺伝型について全ての組み合わせを見ており、大変なだけではなく無駄があります。これを、ループのない（近親婚のない）家系においては情報量を失わずに計算量を大幅に減らすのがエルストン・スチュワート・アルゴリズムです。peeling algorithmとも呼ばれます。

この図の家系で考えてみます。

この図について、前回ご紹介した尤度を用いると次のような計算になります。

$L=\sum_{g_1 \in G_1} \sum_{g_2 \in G_2} \sum_{g_3 \in G_3} \sum_{g_4 \in G_4} \sum_{g_5 \in G_5} \sum_{g_6 \in G_6} \sum_{g_7 \in G_7} \sum_{g_8 \in G_8} \sum_{g_9 \in G_9} \sum_{g_{10} \in G_{10}} P(Y_1|g_1)\times$
$P(g_1)P(Y_2|g_2)P(g_2)P(Y_3|g_3)P(g_3)P(Y_4|g_4)P(g_4|g_1,g_2)P(Y_5|g_5)P(g_5|g_1,g_2)\times$
$P(Y_6|g_6)P(g_6)P(Y_7|g_7)|P(g_7|g_3,g_4)P(Y_8|g_8)P(g_8|g_5,g_6)P(Y_9|g_9)P(g_9|g_5,g_6)\times$
$P(Y_{10}|g_{10})P(g_{10}|g_5,g_6)$

総当りしてるわけですね。大学入試の確率の問題を力技で解く要領です。式の改行は私の環境に合わせてやりました。

ここで、和の数字の順番から見て、「上から下に向かって尤度を書き下している」ことにご注意下さい。図でも矢印でこれを表現しました。本質的にはこの点だけ把握していれば十分。

もうちょいうまくやれないものか、ということをやるのがエルストン・スチュワート・アルゴリズムです。

まず核家族を考えます。次の図のように、家系を核家族に分解します。

核家族では、親と子の二世代しかいません。

ここで、核家族間をつなぐ個体linking individualに着目します。なぜなら条件分岐が、これらつなぐ個体に関して絞れる可能性が高いからです。このサンプルについて遺伝型を固定し、それを親とする核家族での尤度を検討します。そしてこれを、家系図の下から数え上げます。

ここでは核家族NF2について、つなぐ個体4の、特定の遺伝型G4の条件下での尤度を見ましょう。すでに見たように、子の遺伝型は親の遺伝型のもとで周辺分布をとるので

$P(Y_3, Y_7|G_4) = \sum_{g_3 \in G_3} \sum_{g_7 \in G_7} P(Y_3|g_3) P(Y_7|g_7) P(g_3) P(g_7|g_3, G_4)$

となります。同様に核家族NF3では

$P(Y_6, Y_8, Y_9, Y_{10}|G_5) = \sum_{g_6 \in G_6} \sum_{g_8 \in G_8} \sum_{g_9 \in G_9} \sum_{g_{10} in G_{10}} P(Y_6|g_6) \times$
$P(Y_8|g_8) P(Y_9|g_9) P(Y_{10}|g_{10}) P(g_6) P(g_8|G_5, g_6) P(g_9|G_5, g_6) P(g_{10}|G_5, g_6)$

この二つをまとめます。今度はそのさらに親である個体1、個体2の条件下で

$P(Y_3, \ldots, Y_{10}|G_1, G_2) = ( \sum_{g_4 \in G_4} P(Y_3, Y_7|g_4) P(Y_4|g_4) P(g_4|G_1,G_2) ) \times$
$( \sum_{g_5 \in G_5} P(Y_6, Y_8, Y_9, Y_{10}|g_5) P(Y_5|g_5) P(g_5|G_1,G_2) )$

最後にこれを折り畳みます。

$P(Y_1,\ldots,Y_10) = \sum_{g_1 \in G_1} \sum_{g_2 \in G_2} P(Y_1|g_1) P(Y_2|g_2) P(Y_3, \ldots, Y_{10}|g_1, g_2) P(g_1) P(g_2)$

うん、綺麗になりました！

開発者の一人ロバート・エルストンは、発表当初、先祖から子孫に受け継がれていくという遺伝因子の性質から考えると、家系を下から畳み込むという方法で大丈夫なのかという懸念が表明されたが、すぐに情報量のロスがないことが確かめられたと述懐していました。

応用

ここに述べたように、エルストン・スチュワート・アルゴリズムは、遺伝学上の新しい発見というよりは、計算量を大きく軽減し、計算機負荷を軽減するものでした。これが開発された1970年代のコンピュータの性能を考えれば、今以上に重要な性質です。

このアルゴリズムを実装した、fortranで書かれたプログラムLIPED、LINKAGEが使われ、数々のメンデル遺伝性疾患の原因遺伝子座位が同定されました。ハンチントン病の原因遺伝子座位を同定し、後のHuntingtinの発見につながった歴史的な論文*2における連鎖解析は、LIPEDを用いて行われました。また、若年性乳癌のBRCA1の遺伝的座位を同定した論文*3 においては、LIPEDとLINKAGEが併用されています。

さらにLINKAGEは遺伝学データをコンピュータプログラムが扱うための標準フォーマットを提供しました。GWASの時代は完全に標準でしたし（それはplinkのおかげという面もあるが）、今、次世代シークエンサーの遺伝的多型データはvcf形式が標準であるものの、それでもなお、linkage形式に変換するプログラムがすぐに用意されるくらい、まだまだ使われています。

これら二つのソフトウェア開発に主導的な役割を果たしたロックフェラー大学名誉教授ユルグ・オットーは、今中国にいるはずです。やるなぁ中国。

計算機負荷上の限界

エルストン・スチュワート・アルゴリズムの計算量は、人数に対して線形、座位数について指数的です。

これが、分子生物学方面でのゲノム学の発展にあたって大きな困難となりました。ゲノム学が発展して、遺伝型を得られる人数も増えましたが、それ以上に爆発的に増えたのは遺伝型を取得可能な座位数だったからです。エルストン・スチュワート・アルゴリズムが発表された1971年というのは、年代から考えて、サンガーやギルバートのDNA配列決定法よりも前の時代です。これが使われ始めた頃、ヒトゲノム上の観測可能な座位数は数個とかだったらしいです！前述のハンチントン病論文で使われていたマーカー数が、全染色体数すらカバーできていない11個だったということが知られています。それが、ものの数年から十年ちょっとで数百、数千、そして万を越えていきました。

そこに登場したのがランダー・グリーン・アルゴリズムだったわけです。

ちなみに2013年現在の観測可能なヒト遺伝的多型座位数は、一塩基多型に限っても一人あたり1200万を超えています。今や個人個人の全ゲノムシークエンスを観測している時代なので、これで打ち止めでこれ以上大きく増えることはなさそうです。こっから先はエピゲノムなどを解析に統合していく方向性は現れるでしょうが、エピゲノムマーカーが数千万もあるということはなさそうに思います。

*1:Elston RC and Stewart J. Hum Hered 1971.

*2:Guesella et al. Nature 1983.

*3:Hall et al. Science 1990