はじめに
私の所属する生物測定学研究室では、輪読の一環としてFalconerの「Introduction to Quantitative Genetics (量的遺伝学入門)」を読んでいます。
この量的遺伝学とは、簡単にいえばメンデル遺伝学を多くの遺伝子が関与できるように拡張した理論であり、花弁の色や豆のしわといった質的な形質しか扱えなかったメンデル遺伝学と違い、伸長や体重といった量的に変化する形質の遺伝法則を扱うことが可能です。現在栽培されている作物における多くの形質は、ほとんどが量的に変化する形質であり、この量的遺伝学が現在の品種改良を大きく支えてきたことは簡単に想像が付くかと思います。
さて、この量的遺伝学では、まずはじめに平均効果と呼ばれる量が導入されます。これはAA,Aa,aaという型を持つ遺伝子について、それぞれどれくらいの効果を持っているかを表す指標であり、ホモ接合基準モデルと呼ばれる別の定義に対して新しく導入される量となっています。ホモ接合基準モデルの方が直観的であるため、習った当初は、なぜわざわざ平均効果を導入しなければならないのかがわからず、非常に苦しめられました。また、この平均効果の定義にはHWE(ハーディーワインバーグ平衡)を仮定するものと仮定しないものがありますが、多くの本でHWEが仮定されているのかされていないのかがはっきりと書かれておらず、混乱を招いています。
今回、この平均効果をHWEを仮定しない形で線形回帰モデルを使いつつ考え直してみたところ、多くのことがはっきりとしましたのでこちらにまとめようと思います。平均効果(正確にはそれに付随する置換効果)は現在の研究とはあまり関係ないと考える人もいるかと思いますが、近年盛り上がりを見せているGS(ゲノミック予測)などを支えるRidge回帰やLasso回帰、BayesBは、基本的に平均効果の推定を行っています。このことからも平均効果を正しく理解する重要性が分かるかと思います。
筆者は平均効果より置換効果の方が本質だと思っていますので、以下では置換効果⇒平均効果という流れで説明を行います。
目次
ホモ接合基準モデル
まずはじめに平均効果モデルと対をなすホモ接合基準モデルについて復習をします。
遺伝子型が$aa,Aa,AA$を取り得る遺伝子について、その遺伝子型値を$-a,+d,+a$と対応させるモデルをホモ接合基準モデルと言います。なぜこのように呼ばれているかというと、$-a,+d,+a$の原点がホモ接合である$AA$と$aa$の中間に設定してあるからです(図1)。
図1. ホモ接合基準モデルの模式図
$+d$がない場合を考えると、$-a,0,+a$となりますが、これは$a$がひとつ増えると効果が$+a$増えるという線形性を表しており、このような効果は相加効果と呼ばれています(図2)。一方で、$d$は線形性からのずれを表しており、優性効果と呼ばれています(図3)。一般に、$d=0$の時は優性効果なし、$|d|<|a|$の時は不完全優性、$|d|=|a|$の時は完全優性、$|d|>|a|$の時は超優性と呼ばれます。このようなモデルは、遺伝子型値を完璧に表せますし、私たちの直観と非常に一致するわけですが、線形性だけで説明しようとすると少し問題が生じます。そしてそれを克服する過程で平均効果モデルが生まれました。以下ではこれを詳しく見ていきます。
図2. 相加効果の線形性
図3. 優性効果の非線形性
ホモ接合基準モデルの課題
先ほど説明したように、遺伝子型値は、相加効果と優性効果に分けられます。どちらも非常に重要な効果なわけですが、線形・非線形という観点から見ると、後者は$a$の増加に対して非線形であり、非線形性は多くの場合モデルを複雑にするため好まれません。従って、線形性のみを使って遺伝子型値を説明してあげたい(近似してあげたい)、というモチベーションが湧いてきます。もしこれが可能になれば、遺伝子型から遺伝子型値を線形性のみで予測することができるようになるということです。
そこで、視点を改めて、線形性のみ遺伝子型値をうまく表すモデルを考えます。先ほどの相加効果も線形性で遺伝子型値を表現するモデルの一つですが、以下の図(図4)を見るとそのモデル化が、あまりよくないことが分かるかと思います。なぜなら、$a$の数が$0$の時のみ予測値とのずれが生じ、他ではずれが生じていないからです(アンバランス)。
図4. 相加効果による遺伝子型値のモデル化
このようなアンバランスさ予測精度の$a$の個数依存性を生じさせるため好まれません。むしろ以下のように(図5)、どれも$a$の個数に関わらず同程度のずれを含むように考えるべきです。
図5. 不正確さのバランスの良いモデル
こうすることで、どの遺伝子型値も同程度の正確さ(あるいは不正確さ)で予測することができるようになります。このようなモデルをどのように求めればよいかが次の課題となるわけですが、これは単純に$y=\alpha x+\beta$といった線形回帰モデルを当てはめ、最小二乗法などで解けば求まります。せっかくなので解いてみましょう。今回のモデルは次のように表すことができます。
\begin{pmatrix} +a \\ +d \\ -a \end{pmatrix} = \begin{pmatrix} 1 & 2 \\ 1 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} \beta \\ \alpha \end{pmatrix} + \begin{pmatrix} e_{1} \\ e_{2} \\ e_{3} \end{pmatrix}
ここで、$X,\mathbf{y}$を次のように定義します。
\begin{align}
X&=\begin{pmatrix} 1 & 2 \\ 1 & 1 \\ 1 & 0 \end{pmatrix} \\
\mathbf{y}&=\begin{pmatrix} +a \\ +d \\ -a \end{pmatrix}
\end{align}
この時、最小二乗解は次のように求まります。
\begin{pmatrix} \beta \\ \alpha \end{pmatrix} = (X^TX)^{-}X^T\mathbf{y}
それぞれ求めてあげると次が求まります。
\begin{pmatrix} \hat{\beta} \\ \hat{\alpha} \end{pmatrix} = \begin{pmatrix} -5/2a+1/3d \\ a \end{pmatrix}
図6は$a=5,d=3$として実際に当てはめた例になります。不正確さが均等に分配されているのが分かるかと思います。
図6. 最小二乗法による線形近似
さて、これで線形近似が終わったかというとまだ終わってはいません。一般に量的遺伝学では、サンプル数($aa,Aa,AA$の数)は$1:1:1$ではないからです。さらに厄介なことに、量的遺伝学ではこのサンプル数は数ではなく頻度によって与えられています。この時どのようにすれば線形近似ができるでしょうか。以下では、サンプル数が頻度によって定められているということを活用し、線形近似をする方法を考えていきます。
線形回帰によるモデル化と置換効果との関連
サンプル数が$N_p,N_h,N_q$と与えられている場合は先ほどのように最小二乗法で解を求めればよいわけですが、量的遺伝学では、$aa,Aa,AA$の数は$P,H,Q\space (s.t.P+H+Q=1)$といったように頻度(確率)で与えられるため考え方がやや特殊です。もちろん全体の数をNと仮定することによりそれぞれのサンプル数を$PN,HN,QN$と置くことができます。それに基づいて最小二乗法を行うことにより、求めることができるのですが(そしてその解は$N$に依らない)、優性偏差の期待値が$0$になることを自然な形で示すことができないのでちょっと都合がよくありません。
ではどうすればよいかですが、実は線形回帰の最小二乗法については$\alpha = \frac{Cov(\mathbf{y},\mathbf{x})}{Var(\mathbf{x})}$,$\beta=E[\mathbf{y}]-\frac{Cov(\mathbf{y},\mathbf{x})}{Var(\mathbf{x})}E[\mathbf{x}]$となることが知られています。これは確率論的で、確率さえわかっていればサンプル数関係なしに計算ができるため、今回の問題設定と非常に相性が良いと言えます。そこで、今回はこの公式を用いて解いていくことにします。
まず、$\mathbf{x}$についてですが、これは$a$の数を表しており$0,1,2$を取る確率変数です。ここではその確率を$P,H,Q$とします。この確率変数の平均は$E[\mathbf{x}]=H+2Q$、分散は$V[\mathbf{x}]=H+4Q-(H+2Q)^2$となります。また、$\mathbf{y}$は$-a,+d,+a$を取る確率変数で、その確率は$\mathbf{x}$と同じ$P,H,Q$を取ります。平均は$E[\mathbf{y}]=a(Q-P)+dH$、分散は$V[\mathbf{y}]=a^2(P+Q)+d^2H-(a(Q-P)+dH)^2$となります。また、今更ではあるのですが、確率変数$\mathbf{y}$については期待値が$0$になるように中心化してから計算をするのが一般的です。そこで$\mathbf{y'}=\mathbf{y}-E[\mathbf{y}]$として新しく確率変数を定義し計算を行います。この$\mathbf{y'}$は期待値$E[\mathbf{y}]=0$、分散$V[\mathbf{y'}]=a^2(P+Q)+d^2H-(a(Q-P)+dH)^2$を取ります。
また、この計算方法では$Cov(\mathbf{y'},\mathbf{x})$を計算する都合で、$\mathbf{y'},\mathbf{x}$の同時確率$p(\mathbf{y'},\mathbf{x})$が必要です。ここで$\mathbf{y'},\mathbf{x}$は一切独立性を持ちませんので次のように仮定するのが自然です。
p(\mathbf{y'},\mathbf{x})=
\begin{cases}
P \space :\space\mathbf{x}=0, \mathbf{y'}=-a-E[\mathbf{y}]\\
H \space :\space\mathbf{x}=1, \mathbf{y'}=d-E[\mathbf{y}] \\
Q \space :\space\mathbf{x}=2, \mathbf{y'}=a-E[\mathbf{y}] \\
0 \space :\space Other \\
\end{cases}
それではこれらの情報を用いながら、具体的に傾き$\alpha$を求めていこうと思います。まず共分散$Cov(\mathbf{y'},\mathbf{x})$は次のように計算されます。
\begin{align}
Cov(\mathbf{y'},\mathbf{x})&=Cov(\mathbf{y},\mathbf{x})\\&=E[\mathbf{y}\mathbf{x}]-E[\mathbf{y}]E[\mathbf{x}]
\end{align}
一式目から二式目の変形では、確率変数に定数を加えても、共分散は変わらないことを利用しました。ここで$E[\mathbf{y}]E[\mathbf{x}]=(a(Q-P)+dH)(H+2Q)$であり、$E[\mathbf{y}\mathbf{x}]$は次のように求まります。
\begin{align}
E[\mathbf{y}\mathbf{x}]&=0*-a*P+1*d*H+2*a*Q\\&=dH+2aQ
\end{align}
よって、$Cov(\mathbf{y'},\mathbf{x})=dH+2aQ-(a(Q-P)+dH)(H+2Q)$が求まります。$V[\mathbf{x}]$の分散は$V[\mathbf{x}]=H+4Q-(H+2Q)^2$ですので、傾き$\alpha$が次のように求まります。
\alpha=\frac{dH+2aQ-(a(Q-P)+dH)(H+2Q)}{H+4Q-(H+2Q)^2}
次に切片$\beta$は$\beta=E[\mathbf{y'}]-\alpha E[\mathbf{x}]$ですので、次のように求まります。
\begin{align}
\beta&=0-\frac{dH+2aQ-(a(Q-P)+dH)(H+2Q)}{H+4Q-(H+2Q)^2}(H+2Q)\\&=
-\frac{dH+2aQ-(a(Q-P)+dH)(H+2Q)}{H+4Q-(H+2Q)^2}(H+2Q)
\end{align}
これで直線回帰の式$y=\alpha x+ \beta$が求まりました。$\alpha$は$aa\rightarrow Aa$、$Aa\rightarrow AA$という風に、アリルをひとつ置換したときの集団平均の上昇度を表しているため、置換効果と呼ばれます。このように線形回帰でモデル化すると、置換効果がなぜ置換効果と呼ばれているかがすっきりわかるかと思います。量的遺伝学の枠組みでは、この置換効果と同等に重要な量として平均効果$\alpha_1,\alpha_2$と呼ばれる量が定義されます。この平均効果も置換効果から直ちに定義ができるのですが、ここで少し立ち止まって、推定された直線の挙動と$P,H,Q$との関連を探ってみます。
次のグラフ図7は横軸が$\mathbf{x}$、縦軸が$\mathbf{y'}$を表しています。$a=5,d=3$の条件は変わりません。赤い軸は推定された回帰直線を表しており、遺伝子型頻度は$P=0.5,H=0.0,Q=0.5$となっています。
図7. $P=0.5,H=0.0,Q=0.5$での回帰
$H=0.0$はヘテロ接合が無いことと同値ですので、直線はホモ接合の二点を通ります。
次のグラフ図8は$P=0.8,H=0.0,Q=0.2$の条件のもとでの回帰になります。先ほどと同様ホモ接合の二点を通る一方、$\mathbf{y'}$軸の中心点($y=0$の場所)が大きく下がっているのが分かるかと思います。
図8. $P=0.8,H=0.0,Q=0.2$での回帰
次のグラフ図9は少しヘテロ接合の割合を増やし$P=0.8,H=0.1,Q=0.1$としたものです。$\mathbf{y'}$軸の中心点は以前低いままですが、ヘテロ接合の影響を受け、傾きがやや大きくなっているのが分かります。
図9. $P=0.8,H=0.1,Q=0.1$での回帰
このように、$P,H,Q$は傾きの大きさだけでなく、$\mathbf{y'}$軸の中心点の位置を大きく変えるということを頭に入れておいてください。
平均効果
次に平均効果という量を考えます。まず推定された直線において、$x=0$の値を$2\alpha_1$、$x=2$の値を$2\alpha_2$とします。線形性が成り立ちますのでその間の$x=1$の値は$(2\alpha_1+2\alpha_2)/2=\alpha_1+\alpha_2$となることがわかります。また傾きは$y$の増加量/$x$の増加量ですので、$\alpha=(\alpha_2-\alpha_1)$の関係が成り立ちます。概略図を次のグラフ図10に示してみました。
図10. 平均効果の概略図
ここで$x=0$は$aa$、$x=1$は$Aa$、$x=2$は$AA$に対応しますので、$\alpha_1$は$a$の効果、$\alpha_2$は$A$の効果と考えることができそうです。この$\alpha_1,\alpha_2$を$a,A$に対する平均効果と呼びます。
通常、回帰は傾きと切片のパラメータで表されますが、もちろん、直線が通るべき2つの点でも表現可能です。この平均効果によるパラメトリゼーションは、傾きと切片は直線が通るべき2点で置き換えられること、また直線回帰$f(x)$においては$f((x+y)/2)=(f(x)+f(y))/2$が成り立つことをうまく利用した表現方法と言えるでしょう。なお、$2\alpha_1=-\frac{dH+2aQ-(a(Q-P)+dH)(H+2Q)}{H+4Q-(H+2Q)^2}(H+2Q)$、$2\alpha_2=\frac{dH+2aQ-(a(Q-P)+dH)(H+2Q)}{H+4Q-(H+2Q)^2}(2-H+2Q)$となることが直ちに求まります。
優性偏差
さて、最初の話に立ち返ります。今回は、線形の関数でなるべくうまく$\mathbf{y}$を表すことを目標にしてきました。しかし$\mathbf{y}$は優性効果$d$がある限り線形ではなく、したがって直線回帰で完璧に$\mathbf{y}$を表すことはできません。このモデル式と実際のデータとのずれは一般に誤差と呼ばれますが、この優性効果$d$に由来するずれのことを量的遺伝学では優性偏差と呼んでいます。今回は線形回帰$y=\alpha x+\beta$のパラメータを最小二乗法で求めましたが、最小二乗法で求まったパラメータによる回帰は誤差の期待値を$0$にするという著しい性質を持っています。従ってこの優性偏差も期待値は0となります。優性偏差の期待値が0になることはかなり一般的な事実として教科書に書かれていますが、これはHWEなどの仮定から来るものではなく、線形回帰からのずれとして優性偏差を定義したことに由来することが、ここからもわかるかと思います。これが線形回帰として置換効果を定義したかった理由の一つです。
ハーディーワインバーグ平衡(HWE)との関連
最後にHWEと置換効果の関連を見ていきます。HWEは、「任意交配のもとで遺伝子や遺伝子型の頻度が変化しない」ことを主張した理論です(厳密にはそれに加え多くの条件が存在する)。この結論より$a,A$の遺伝子型頻度$p,q$と遺伝子頻度$P,H,Q$の間に$P=p^2,H=2pq,Q=q^2$という関係性が生じます。今回は、線形の関数で$\mathbf{y}$をうまく表すことを目標にしたわけですが、そこで予測されたパラメータは$P,H,Q$によって完全に決定されることがわかりました。逆にいうと、HWEのもと、$P,H,Q$が変わらない限り、予測関数$y=\alpha x+\beta$は常に正しく機能するということができます。これがHWEのありがたい点です。
また、先ほどもとまった置換効果$\alpha$ですが、$P=p^2,H=2pq,Q=q^2$を代入すると$\alpha=a+d(q-p)$が求まります。これはよく教科書に書かれている置換効果に一致します。従って上記の置換効果導出方法はHWEに依らない一般的な導出方法であることがわかります。また、HWEのもと導出される平均効果$\alpha_1=q\alpha,\alpha_2=p\alpha$も直ちに求まります。やや計算が厄介ですが、お時間のある方はぜひ試してみてください。