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

遺伝疫学

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

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

ランダー・グリーン・アルゴリズムの前段階である、エルストン・スチュワート・アルゴリズム*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

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

遺伝疫学

小沢新党なみの純化路線をたどっているブクマ数が気になる！

家系の尤度

家系の「尤度（ゆうど、likelihood）」というお話です。

まずは「尤度」について説明しようと思って色々書いたり消したりしていたのですがとりあえずこっち見てもらったほうが速いのかもと思います。http://www.genstat.net/statistics.html　というか遺伝関連もそっち見てもらったほうが速いかもしれませんが・・・

さて、そういうわけなので、われわれは今、何らかの観察データがある元で、組換え価やら何やらのパラメータ推定をして最終的に遺伝的マッピングをするというランダー・グリーン・アルゴリズムに辿り着こうと努力していますので、ここでお話しているのは尤度に関してのことであります。めんどくさければ尤度という言葉は頭のなかで「確率」に置き換えても、実践上はそんなに問題がないような、いや、あるそうです。頑張りましょう。

ある人の遺伝型尤度をP(G)と書きます。Lかlで書かないといけないような気もするがあまりそうやって書いてる人はいないような。通常、何も他に情報がない場合、これはHardy-Weinberg平衡というものを仮定して得ます。例えばこんな感じで。ABO血液型について、Aの頻度をpA、BをpB、OをpOとすると

P(G=AA) = p_A^2
P(G=AB) = 2 p_A p_B
P(G=AO) = 2 p_A p_O
P(G=BB) = p_B^2
P(G=BO) = 2 p_B p_O
P(G=OO) = p_O^2

次に、表現型尤度はP(Y)です。表現型というのは、一般に文字通り表に出ているデータを意味しますので、これが基本的な観察データになります。例えばy=1なら「癌である」y=0なら「癌ではない」とします。

最後に浸透度のモデルP(Y|G)を思い出しましょう。例えば血液型についていくつか例を挙げると、

P(Y=A|G=AA) = 1.0
P(Y=A|G=AO) = 1.0
P(Y=A|G=OO) = 0.0

これは、血液型が優性遺伝形式をとっていることを示しています。浸透度を使えば、様々な遺伝形式を表せることがわかります。

ここで個人の観察の尤度P(Y)は、基本的な条件付き確率の式変形からP(Y)=P(Y|G)P(G)というふうに浸透度と遺伝型尤度の積に分解することができ、

P(Y) = \sum_{g \in G} P(Y|g)P(g)

と書くことができます。Gというのは、この人に関して可能なデータの組み合わせのことで、前回記事から大体イメージつかめると思います。gはその組み合わせの中から一つ取り出したものです。ここまでわかっているP(G)とP(Y|G)の記述の仕方を用いて、P(Y)が表せることがわかったと思います。

次にこの尤度を、個人のものとしてではなく家系、つまり人の集まりにおいて考えます。

通常このような遺伝解析で用いるような病気という表現型においては、互いに独立であるとすることができます*1。統計学的に独立であるということは、同時確率をそれぞれの確率の積で表すことができるということです。なので、ある家系の尤度は

L = \prod_i P(Y_i) =\sum_{g_1 \in G_1} \sum_{g_2 \in G_2} \cdots \sum_{g_m \in G_m} \prod_i P(Y_i|g_i)P(g_i)

として表せます。全ての可能な遺伝型についての全ての組み合わせを挙げています。

これが家系の全確率としての尤度です。

親データが存在する場合

次に、その個体に親データがある場合には、メンデルの法則にもとづき、子の遺伝子型は親の遺伝子型によって影響されます。単純に言えば、父親が血液型AA、母親がBBだったとしたら、OOの子供は生まれません*2。P(G=OO|G_p=AA,G_m=BB)=0ということです。前段でのP(G=OO)とは異なった値をとっており、つまり、尤度が親のデータによって条件付けされています。一般に

P(G_o | G_p, G_m)

と書きます。oというのはoffspringのoです。pが父親、mが母親。例えば、AOの父親とBOの母親から生まれる子供は

P(G_o=AO|G_p=AO, G_m=BO) = 0.25
P(G_o=BO|G_p=AO, G_m=BO) = 0.25
P(G_o=AB|G_p=AO, G_m=BO) = 0.25
P(G_o=OO|G_p=AO, G_m=BO) = 0.25

となります。これはメンデルの分離の法則そのままです。

これを導入すると、先ほどの全家系の尤度は次のように書き直せます。ここで、観察範囲のデータ内に父親か母親がいるような個体をnon-founder、父親も母親もいないものをfounderと呼ぶことにします。

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)


以上の考え方はそのまま、遺伝型ベクトルG= {g1, ... gn}について適用できます。するとここに組換え価を入れることができるようになり、遺伝的マッピングに応用出来るようになります。

終わりに

間違っていないか不安ですが、誰か読んでくれてるんだろうか・・・一応お一方は読んでくださっているようだ！

*1:「病気の表現型が、家系内で独立でない」というのは、父親が病気になったら母親や子供も病気になるという干渉を意味していて、つまり結核のような感染症では明らかに独立ではありません。が、感染症以外のほとんどの病気に関しては独立と考えることができます

*2:突然変異が起きた場合を除く

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

遺伝疫学

あと3歩くらいでたどり着きそうかも・・・

遺伝型の相（フェイズ）

前々回、遺伝型について説明しました。これはヒトのもつ二つの染色体上のアレルを組み合わせたもので、遺伝的データのうち病気などといった表現型へと直接的な関連を示すもので、浸透度モデルに組み込まれていることを前回お話しています。常染色体上にある場合、ペアの染色体に優劣はなく、これは順序を入れ替えても特に意味の変わらない組み合わせです。例えば、ABO遺伝子座上のABO*AアレルとABO*Bアレルからなる、AB型という血液型とBA型という血液型・・・という違いは特にありません。

ところが、家系データにおいては意味が発生します。父親由来のアレルと母親由来のアレルに分けることができるからです。これをする意味は、これをすると家系データにおいてアレルの流れ（gene flow）がわかるということです。由来親をたどっていけばいいわけですから。さらに遺伝的マッピングにおいては、後述するように組換え価の推定に使えることが本質です（組換え価さえわかればマッピング自体は終わったようなものです）。このように父親・母親由来のアレルにそれぞれわけて扱っている場合、遺伝型（ジェノタイプ）の相（フェイズ）を見ていると言います。

ここでは次世代シークエンサーの標準ファイルフォーマットであるvcfを踏襲して、フェイズ付き遺伝型を|で区切って表現します。例えば父親由来を左に書き、母親由来を右に書いて、間を|で仕切ります。この書き方はvcfができるずっと前から、家系解析をしていた人たちの一部は使っていました。ご存じない方のために注釈しますと、フェイズなしの場合vcfは/で区切ります。例えば血液型ABの遺伝子型は、フェイズがわかっていなければA/B、わかればA|Bとします。フェイズがあるときは、順序に意味があります。

さて血液型についての次のような家系を考えましょう。

父親がAA、母親がAO、子がAOの遺伝型であるとき子のフェイズはわかるでしょうか。簡単にわかりますね。（突然変異が起こったのでなければ）母親からしかOをもらうことはできないので、子のフェイズはA|Oです。


しかしこのような家系、父親がAB、母親がAB、子もABの場合、子のフェイズはわかりません。A|BまたはB|Aであって決定できません。

また、いずれの場合も親のフェイズはわかりません。

二遺伝子座のフェイズ


同じ染色体上にある二つの遺伝子座上のアレルは、それぞれ二つの染色体ペアの片側に割り振ることができるので、二遺伝子座のアレルの組み合わせ、つまりハプロタイプがどちらの染色体上にあるか、つまりフェイズがどうであるかを考えることができます。

同じ染色体上にある二つの遺伝子座A、Bを考え、それぞれアレルがAとa、Bとbであるとします。


このような家系、父親がAaとBb、母親がAAとBBで、子がAaとBbであるようなものを考えましょう。母親のフェイズは祖父・祖母由来は不明であるものの、AB|ABであることが明らかです。父親は、AB|ab or ab|AB、またはAb|aB or aB|Abです。子を考えると、AB|ab, ab|AB, Ab|aB, aB|Abという可能性が挙げられます。突然変異が起こっていないなら、母親からの由来ハプロタイプがABしかありえないことを考えると、可能性はab|ABしか残りません。したがって子のフェイズが決定出来ました。

このとき、父親由来のハプロタイプabについて、二通りの考えができます。

  1. 父親のハプロタイプ・フェイズはAB|ab or ab|ABであったため、連鎖したハプロタイプabが伝達された。
  2. 父親のハプロタイプ・フェイズはAb|aB or aB|Abであったのだが、組換えが起こってabが伝達された。

このどちらであるかはこの情報だけでは決定はできません。しかし、AB遺伝子座間の組換え価というパラメータの元に尤度を立てることができます。これをもとに最尤法で、伝達されたハプロタイプがどれであるかとか、組換え価がどれくらいであるかを推定していきます。次回以降見て行きたいと思います。これだけではあまりに情報が少ないので推定は不定になるでしょうが、家系を大きくしたり、増やしたり、周辺の座位数を増やしたりして情報量を大きくしていくことで、推定の正確性を増すことができます。

横道にそれる

ちなみに「親のフェイズは決定できない」と書いたのは、その個体の親データがなければ今回のエントリに書いたような考察はできないということですが、親データのない個体の集団データがある場合、集団におけるハプロタイプ頻度というものに着目してEMアルゴリズム*1やGibbsサンプラー*2を用いてフェイズを取ることはでき、後者はゴールドスタンダードとしてHapMap計画におけるphased haplotype推定に使われました。

*1:Excoffier L and Slatkin M. Mol Biol Evol 1995

*2:Stephens M, Smith NJ, Donnely P. AJHG 2001

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

遺伝疫学

ランダー・グリーン・アルゴリズムの紹介、第二回です。いつになったら辿り着くのやら！

メンデルの法則

まず第一に、これはメンデルの法則にもとづくモデルです。メンデルの優性、分離、独立の法則を思い浮かべましょう。

優性の法則 law of dominance

優性の法則は、ある遺伝子座上にA、aのアレルがあり、Aが優性ならば、AAとAaはどちらもAの表現型を表し、aaだけがaの表現型を表し、これらの遺伝型が1:2:1の比率で生まれるので、優性形質Aが3:1の比率で生まれる。というふうなものですが、これはさすがにプリミティブすぎて現在の遺伝学には適用不可能です。

ここは、メンデルの優性の法則を、浸透度penetranceの定義で置き換えます。これは、ある遺伝型のもとで、表現型が一定の確率で発生するとするものです。遺伝型G、表現型Yについて、P(Y|G)を浸透度と呼びます。P(Y|G=AA)=1、P(Y|G=Aa)=1、P(Y|G=aa)=0とすれば先ほどの優性遺伝を表します。このように浸透度は、メンデルの優性の法則を包含する概念です。

日本では三つの法則として習ったと思うのですが、今手元の欧米の教科書は"Mendel's Two Laws"と言ってこの優性の法則は挙げられていません。

分離の法則 law of segregation（第一法則）

分離の法則は、親が同じ遺伝子座上に祖父由来A、祖母由来aの二つのタイプのアレルをそれぞれの常染色体に持っている時、子にAが伝達される確率とaが伝達される確率は等しい(p=0.5)とするものです。これはまあいいですね。遺伝統計学理論のほとんどは、この確率0.5を付与のものとしています。現在でも行われるTransmission Disequilibrium Testという遺伝的関連解析法は、ほとんどこれだけを用いて解析を構築しています。

独立の法則 law of independent segregation（第二法則）

独立の法則は、別々の遺伝子座A、Bがあり、そのそれぞれにアレルA、aと、B、bがあるとき、その組み合わせ、AB、Ab、aB、abが子に伝達される確率は等しい（p=0.25）するものです。これは分離の法則のもとで、すべての遺伝子座が互いに独立だと考えるなら正しいと思われます。


これが成立しない場合があります。つまり遺伝子座間は必ずしも統計的独立ではなかったのです。この、独立の法則の例外は、ベイトソンにより1905年、モルガンにより1911年に報告されました。これは、二つの遺伝子座が非常に近くにある場合に起こります。その時、同じ染色体上の二つの遺伝子座は、前回の減数分裂の図を見るとわかると思いますが一緒になって伝達されます。つまり、ABとabしか伝達されず（それぞれp=0.5）、AbやaBとして伝達される可能性がほとんどなくなります（いずれもp=0）。右図の右矢印のような感じです。このような状況を完全連鎖complete linkと言います。

ところが、ここに前回紹介した乗換え現象が起こるため、同じ染色体上だけどやや遠くにある遺伝子座間ではその間で乗換えが起こり、AbやaBとして伝達されることも起きるのです。これが上図の下側にしめされています。しかし、毎回必ず乗換えが起こるばかりではないので、この場合の伝達確率は、0~0.5の値をとることなるでしょう。すなわち、結論としてはAB、Ab、aB、abの伝達確率が0.25からずれることになります。統計屋さんなら感覚的にわかると思いますが、ここが解析のキモです。0.25からずれている場合、完全ではないものの「連鎖」と呼びます。モルガンはこの連鎖の法則を「遺伝の第三法則」とも呼びました。

遺伝統計学的には、ABが乗換えcross overの影響でAbに変わった場合を組換えrecombinationと呼ぶようにしているようです。すると、A遺伝子座とB遺伝子座の間で二回乗り換えが起こると、並びはABのままなので組換えは起こっていないことになります。乗換えは、生物学的現象なので、遺伝子座間で2回以上起きることができますが、組換えは遺伝子座間で起きたか起きないかのバイナリの現象を表します。この定義は後の組換え価の定義に影響してきます。

「近くなら組換えは起こりづらい」「遠くなら組換えは起こりやすい」というのが遺伝的マッピングのファンダメンタルな考えです。これをもとに、全ての理論が組み立てられています。

後のために用語をひとつ導入しますと、同じ染色体上のアレルの並びをハプロタイプと言います。アレルと同様、ハプロタイプも、人は二個ずつ持ちます。分離の法則がアレルについての法則であるのに対し、独立の法則はハプロタイプについての法則であると考えることもできます。

組み換え価recombination rate

独立の法則において現れた「連鎖」と「組換え」の現象についてもう少し話を進めます。

親が祖父由来AB、祖母由来abを持っていたのに、連鎖が不完全であったため、AbやaBと組み変わって伝達されるようになる可能性を組み換え価として表します。ある個体のある染色体上の遺伝子座iにおいて、それが祖父由来なら0、祖母由来なら1と記述し、Siで表すとすると、二つの遺伝子座間の組み換え価は

\theta_{ij} = P(S_i \ne S_j)

と定義できます。前述したように、これは二遺伝子座間が近いなら小さく、遠いなら大きいと考えられます。

このθは最大値が0.5です。なぜなら生物学的に、非常に近傍にあって乗り換えが起こらない(θ=0)から、非常に遠くにあったりべつべつの染色体上にあってアレルの組み合わせがほぼランダムになる(θ=0.5)までの値しか取り得ないからです。従ってこれは確率の要件を満たさず、recombination rateと呼ばれ、「組換え価」の訳語ならよいのですが、「組換え確率」と呼んでしまうとやや問題があります(これは1以上を取りうるのに罹患率とか発生率と訳されてしまったincidence rateの訳語問題と同じです)。

ここで、もう一つの遺伝子座kを考え、
\theta_{ij} \lt \theta{ik}

であったとしましょう。組換え価の大きさが遺伝子座間の距離を反映するとするなら、これは染色体上の遺伝子座の位置の並びがi, j, kであることを示唆するかもしれません。もちろん互いに別方向の可能性もありますが、ここは近傍の場合は組換え価に加法性があるので、

\theta_{ij} + \theta_{jk} \approx \theta_{ik}

であるなら、少なくともi, j, kもしくはk, j, iの並びである蓋然性が高いです。これが最も初期に提唱された遺伝的マッピングの方法です。このようにして、最初にゲノム上の遺伝子の並びについての考察を行ったのはスターツバントで、1913年のことだそうです*1 。

マップ関数

すぐその後、ホールデンがこの組換え価を遺伝子座間の距離の関数と捉え、遺伝的距離を定義しました*2 。ホールデンのマップ関数と呼ばれます。

\theta = M(d)

ここで、乗換えは完全にランダムにおけると仮定しましょう。ゲノム上どの場所でも同様に起こるというものです。また、その確率は距離に比例して大きくなると考えます*3。乗換えの起こる回数は整数値で、ほとんど起こらないのでポアソン分布で表せると考えられます。このポアソン分布のパラメータ（平均値かつ分散）をxで書くなら、乗換えがc回起こる確率は

P(C=c) = \frac{x^c e^{-x}}{c!}

です。組換えという観点からすると、前述したように、乗換え回数が1回でも3回でも5回でも現象としては同じ物を観察します。従って、組換え価は

\theta = P(c=1|x) + P(c=3|x) + P(c=5|x) + \cdots
= e^{-x} \left( x + x^3/3! + x^5/5! + \cdots \right)
= e^{-x} \sinh(x) = \left(1-e^{-2x}\right)/2

sinhのテイラー展開を利用してます。この式中のxを遺伝的距離と呼び、これが組換え価を用いて次のように求めることができることがわかります。

x = -\frac{1}{2}\log\left(1-\theta\right)

遺伝的距離の単位としては、100回の減数分裂で1回の組換えが起こるような距離を1cM(モルガン)と呼びました。この単位は現在でも組換え価の表現において用いられており、現在では1cMはだいたい1Mbp、つまり物理的な100万塩基対に対応すると推定されています。ヒトゲノム長が30億塩基なので、全ゲノムの遺伝的長さは30M、すなわち1回の減数分裂で全ゲノムに30の乗り換えが起こっている計算です。

先述のスターツバント論文から教科書に引用されている数値を示すと、線虫の遺伝子y、w、miにおいて、
\theta_{y,w} = 1.3%
\theta_{w,mi} = 32.6%
\theta_{y,mi} = 33.8%
だそうなので、遺伝的距離は、y, w間が1.3cM、w, mi間が52.2cM、y, mi間が55.8cMと計算出来ます。それで、y, w, miの順にだと言う彼の考察はあてはまります。ちなみに組換え価の加法性は近傍でしか成立しづらいと述べており、ここでもなんだか足し算が合わないように思われますが、これはかなり遠くにある遺伝子座間は二回以上の乗換えを経ている可能性があるためです。

さきほど言ったように組換え価は実際には距離の単純な関数ではなく、例えば女性の組換え価は男性より高いので、同じ遺伝子座間でも女性の方が遺伝的距離が長いという結果を得ます。改良されたマップ関数も多数存在し、コサンビのマップ関数が最も有名です。次のようなものです。

x = \frac{1}{4}\log\frac{1+2\theta}{1-2\theta}

これで先程の遺伝的距離を計算し直すと、1.3cM、38.9cM、41.1cMとなります。


ここで組換え価は、動物においては実験的に、二項分布パラメータとして最尤法で推定することができます。ヒトの場合はそのように単純にはわかりません。

*1:Sturtevent AH. J Exp Zool. 1913 さすがに読んでない

*2:Haldane JBS. J Genet. 1919

*3:後述する遺伝的距離と物理的距離にはややズレがありますが、このズレの原因はこのような仮定の現実とのズレにあります

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

遺伝疫学

お気に入り経由で以下の記事を読んでとてもおもしろかったのですが、

http://ama.an-pan-man.com/archives/814

原文も読んで確認もしたんですけど、エリック・ランダーが自己紹介するにあたって「ランダー・グリーンアルゴリズム*1」について一言も触れてない。

これはとても美しいアルゴリズムで、隠れマルコフモデルの遺伝家系データへの適用で、最初にこれを思いついたのは天才だと私は思う。今でもこのアルゴリズムは応用されており、たとえば次世代シーケンサーのデータにおいても行われるハプロタイプの相決定や、imputationアルゴリズムと言ってヒトゲノム上数十万のマーカーデータから、相関構造を利用して一千万以上のマーカーデータを推定するというのに使われます。

わたしは数学出身ではありません。EMアルゴリズムのDempster論文も、なんとか読み終えることは出来ましたというレベル。私などがこれを正しく紹介できるかわかりませんが、他のお気に入り経由で読んだ（やや古い記事ではありますが）

satomacoto: Pythonで隠れマルコフモデルのFilteringの例

を読みまして、例示というのもそれなりの意味はあるものと思いました。

そこで、自分の理解の確認がてら、ランダー・グリーンアルゴリズムをご紹介してみようと思います。わかりやすさのため、元論文そのままの表現ではなく、その後の拡張であるKruglyak論文*2や教科書*3 *4などでの表現も必要に応じて流用しています。

問題設定

ランダー・グリーン・アルゴリズムが解こうとする問題は、

「ある家系があって、その家系内で疾患が多発している。また、その中の十分な数の方々からDNAの提供を受け、DNAデータがある。この病気の原因遺伝子座位を知りたい」

です。こんな感じの家系があるわけです。

2013年現在、ここ数年の進歩を元にすると、百万円くらい出せばもう一人分の十分な精度のヒトゲノム全長を得られてしまうので、この程度の大きさの家系なら全例シークエンスして、単純に比較してあとは機能予測すりゃいいじゃんてなもんですが（そーゆー論文がいっぱい出ています）、なにしろ25年ほど前に考えだされたアルゴリズムです。この時点では、ヒトゲノム上でタイピング可能な多型はたかだか数百〜数千程度であったとご想像下さい。問いかけをより正確に書けば、これらをマーカーとして、どれが原因遺伝子座に近いか？ということになります。

歴史的意味しかないと思ってしまうかもしれないけれど、こんな25年前のアルゴリズムが今の次世代シーケンサー解析で応用されている事自体が、この進歩の速い学問分野においては、かなり先駆的な業績だったと思います。

生物学的な前提

遺伝統計学を正しく行うには、生物学の知識もしっかりしていなければなりません。

まず、DNAがヒトなど生物の基本的設計図になっていて、それをもとに転写産物やタンパク質がつくられ生物がつくられます。セントラルドグマと言いますが、その例外の存在もすでにわかっています。

DNAは、22本の常染色体と1本の性染色体からなっており、そのそれぞれが一対のペアからなっていて、ただ男性の性染色体だけはXとYという別々のものからなります*5。

常染色体上のある一点を取ると、そこにはペアの染色体に由来する二つの情報があります。情報、というのは、ヌクレオチド（アデニンA、シトシンC、グアニンG、チミンT）でもいいし、挿入欠失（片側はACTGTCなのにもう片側はAC--TCとなっている、だとか、集団間での同様の違いとか）でもいいし、リピート配列の違いでも言いし、HLAやCYPのように遺伝子そのものでもいいし、とても長いコピー数多型でもいいです。一つの情報のことをアレル（対立遺伝子）と言います。すると、一つの場所には一対二つのアレルがあることになります。このような場所のことを遺伝的座位（ローカス）と呼びます。二つのアレルの組み合わせを遺伝型（ジェノタイプ）と呼びます。

一番わかり易い例は血液型です。血液型はA, B, Oに対応したアレルが3種類あり*6、可能なジェノタイプは、

AA
AB
AO
BB
BO
OO

です。それで、最終的に人間の特徴として現れるもの、表現型（フェノタイプ）は、上からA, AB, A, B, B, Oとなります。

このペアの染色体のうち片方だけが、親から子に伝達されます。片親から子という状況を考えていますが、ここで、親の染色体ペアのそれぞれが「祖父由来」「祖母由来」であることに注意しましょう。生物学的に重要なのは、この伝達がされるとき、真核生物では、ですけど、減数分裂meiosisといって、受け渡される染色体がどちらか選ばれますが、この時「祖父由来」の染色体がまるごと子に伝達されるわけではありません。乗り換えcross overという現象が起こるため、祖父由来と祖母由来がモザイク状になった染色体が、子に受け継がれることになります。

*1:Lander ES and Green P. PNAS 1987

*2:http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1915045/pdf/ajhg00019-0253.pdf:Kryglyak L et al. AJHG 1996.

*3:Balding DJ, Bishop M, Cannings C. Handbook of Statistical Genetics. John Wiley & Sons 2007.

*4:Thomas DC. Statistical Methods in Genetic Epidemiology. Oxford University Press 2004.

*5:ただしX、Y間にはpseudoautosomal regionというペアと考えられる領域もあります

*6:ほんとはもっと多くて、かなり大雑把に言ってもAとOは二つずつにわかれて合計5個ですが

環境省ゲノムのやつ

福島県住民のゲノム検査断念の理由を

ゲノムの解読過程で機械的な誤りが生じるため、専門家からは「親子間に違いがあったとしても放射線の影響なのか、他の要因なのかも区別できない」と妥当性を疑問視する声が続出。さらに、実子ではなかったことが解析で明らかになった場合の倫理的な課題も指摘されていた。

http://mainichi.jp/select/news/20130312k0000m040022000c.html

と書かれていますが、ちょっとよくわからないです。

まず第一に後者の倫理面の問題ですが、これは確かに重大な課題ですが、これ自体が計画全体を中止させるものとまでは言えないと思います。とはいえこの問題点の指摘は正当です。

最初のやつ『ゲノムの解読過程で機械的な誤りが生じるため、専門家からは「親子間に違いがあったとしても放射線の影響なのか、他の要因なのかも区別できない」』というやつなんですけど。これよくわからない。

まず第一に『ゲノムの解読過程で機械的な誤りが生じるため』・・・他の要因かどうか区別できない、って言いますけど、これいい始めたら今、毎週トップジャーナルに報告が続いている希少疾患のシークエンス研究報告は一体どうなるんじゃ・・・全部エラーだったとか！！？

実際にはかなりシークエンスの精度は上がってきており、その上どうしても確認したい場合はサンガーという確認法があるので（精度が上がりエラーがそれほど多くない状況下では）あまり問題ないはずです。これは現在の国際的コンセンサスだと思いますけど。


次に『親子間に違いがあったとしても放射線の影響なのか、他の要因なのかも区別できない』というやつですが、これは複雑です。

  1. そもそも、血液採取して細胞成分を取った時に見られる親と子の遺伝的変異の違いは、生殖細胞ゲノムの違いです。放射線は関係ない。「放射線か、他の要因なのか」というレベルの問題ではありません。放射線を含め、出生後に生じた変異の影響は入っていないもの（とされるもの）を見ているのです。原発事故以前に生まれていた子を見るのなら、原理的に言って何もわかりません。
  2. 原発事故中・事故後に受精し生まれた子については、親の生殖細胞に入った変異が子の生殖細胞ゲノムにおいて見える可能性があります。つまり血液採取によるゲノムに、です。親と子の間で、自然の状態でも突然変異が入り、この変異率については現在世界で精力的に推定が詰められていて、つまりコントロールデータがあります。この状況下では親子間の違いに放射線が影響したかどうかは、大きな違いがあったならわかりますし、統計的に判断できないくらいの小さい違いであるかどうかもわかります。（ただしこの方法は当初研究の計画には入っていないはずです）

どっちもあてはまらないよ〜。ほんとに専門家がそんなこと言ったんかいな。

実のところ、今回の環境省の計画の大問題は、福島医大にゲノムの専門家はいないのに、環境省と福島医大だけで計画を決定し、その計画が前提知識を決定的に欠いていたということでした。ただ、事故後出生児を見ることをメインにするとか、体細胞変異率を見る方向にするとかで、やりようはあります。また、計画中止は難しかろうと思いました。だって、代替案もなくこの計画を中止するというのはちょっと恥ずかしい。

でも、とくに言い訳もなく先送りということだそうです。

なんでこんなことが許されちゃうかというと、記者の実力不足ではないでしょうかね。もうちょっと、何が問題点なのか、どうしてその上でこの計画が立案されたのかを論理的に詰めて欲しい。

などと考えておりました。誰にも専門分野はあります。環境省は素晴らしいなあと思うこと多いです。ですが専門外に手を出すときは、専門家にちゃんと意見聞いたほうがいいと思います。厚労省が金だした東北大ですら（ゲノムの専門家がいないので？）手こずってると聞きますのに。

※2013-03-12　一部過剰な表現を訂正しました。

「ラウンダバウト（ロン・ポワン）」について

http://www.yomiuri.co.jp/national/news/20130223-OYT1T00652.htmについて

round about、フランスではrond point（ロン・ポワン）。自分の経験について少しメモしてみます。ラウンダバウト議論はほとんどイギリスを参考にしてるっぽい（そもそも名前からしてそうか）ので、フランスの話はちょっと見当違いかもしれないけど。。。

田舎では

  • 欧州の田舎では、一般道を100km/hオーバーで走る車がいるので、適切な間隔で配置されたロン・ポワンは速度抑制に役立っているっぽい。
    • 信号だと加減速が過剰になってしまいそう。
    • では、そんな速度で走らない日本での有用性は？とは思う。
  • 90km/h（法定速度）から急減速し、10Rくらい？のロン・ポワンを同乗者に迷惑かけずに走り切るにはマニュアル車でないと厳しい気がする。日本と違って欧州はマニュアル車だらけである。
  • 信号で止まらず走り続けられるというのはストレスは少ないかも。
    • 危険性に関して言うと、特に大きめのロン・ポワンがあるような道の場合、歩行者はそもそもいない。
    • 村などに入って歩行者がいるようになると、やはり信号が出現したり、また最近では速度表示警告付きのオービスがどんどん設置されてる。

都会（パリ）では

  • 田舎はほとんどがロン・ポワンで信号が数えるほどだが、逆にパリにはそんなに多くはなく、やはり田舎型の交通状況に適していることを反映しているのでは。
  • 信号と組み合わされたロン・ポワンも稀ではない。この場合「どの信号を見ればいいの？」は、初期にはわりとむずかしい。
  • これは各国違うだろうが、フランスではロン・ポワンの交通上暗黙のルールがあり、それは「絶対的な左優先（日本で言う右）」。ま、そうじゃないと出れなくなるから、5〜6車線ロン・ポワンの凱旋門とか。
    • 普段の道路では「右優先」なので、ロン・ポワン内でルールを入れ替えていることになる。
    • 車は右（日本で言う左）側からロン・ポワンに進入することになるので、つまり進入側が待つことになるということで、読売紙でも「過去の日本のロータリーとの違い」として紹介されたルールと同じものだと思います。
    • こちらのページにもあるけれど、法律とは関係なく、そこの住民が完全に同意する暗黙のルールがない限り、ロン・ポワン的な構造は維持が難しいのではないかと思う。
    • というわけで、事故ったら右の車が必ず立場が悪くなる
    • ただし、凱旋門の巨大ロン・ポワンだけは例外扱いで、責任関係立証不可能とみなされ、どんな場合でも50:50の過失割合になるそうである（経験者談・・・僕じゃないよ！）。
    • 自信がないので外周をぐるぐる周る・・・はルール違反らしく、怒られる。侵入してから出口までの距離が遠いなら、内側に入る。マナー問題に過ぎないけど。
  • わりと大きなスペースに設置され、かつ交通量がそれほど多くない場合、ロン・ポワンの外周に沿って路面駐車場が設置してあることもある。駐車後出発する際に行き先を自由に決められるというのは便利
  • 歩行者は、ロン・ポワンは大回りして向こうに行かなければならず、単純にめんどうである。
    • 歩道橋は見たことない。
    • 交通量が少ない時は、歩行者は余裕でロン・ポワンを横断する。フランスは、絶対的な「歩行者優先」。たとえどんなに車側がルールを守っていて、歩行者がルールを破っていても、それでも歩行者優先。これもロン・ポワンを使用するとしたら重要なこととは思う
  • ちなみにフランスを含め欧州では、日本と比べ車は赤信号を厳密に遵守する。それどころか黄色でも大体の車は止まる。かなりの急ブレーキをかけてでも止まる。あの適当で怠惰なフランス人なのに！
    • どうもその理由は、少なくともフランスでは、速度超過や駐車違反は単にお金払えばいいが、信号無視についてはかなり厳しくて出頭しなければならなくなるから、らしい。
    • ロン・ポワンが好まれる理由の一つにそれもあるかも。日本も信号無視もっと厳しくして、ロン・ポワンどんどん入れる、というのはアリかも

感想

慣れてみると、特に田舎道でのロン・ポワンは大変気持ちいいです。都会には、なくていいです。