HWE とLD

R 数理モデル 統計 遺伝医学 医学

アレルA とアレルa がある。
アレルA の存在確率がpのとき、アレルa の存在確率は1-pである。
1世代交配して得られる遺伝子型は(AA, Aa, aa)=(p^2,2p(1-p),(1-p)^2)である。
これは、アレルA とアレルA を持つ者同士は強制的に子孫残してね、というような神の見えざる手が働かない、つまり、完全にランダムなときにこの比率になる。
この場合にHWE という。
実際には、完全にランダム、とはならず、とある係数F(近交係数と呼ばれる)により
\begin{pmatrix}AA\\Aa\\ aa \end{pmatrix}=\begin{pmatrix}p^2+p(1-p)F\\2p(1-p)(1-F)\\(1-p)^2+p(1-p)F\end{pmatrix}
となる。
AA とaa からFが計算できるので、ゴリ押しすると

N1 <- c(36, 48, 16) # p = 0.6, F = 0
N2 <- c(42, 36, 22) # p = 0.6, F = 0.25
HWE_F <- function(vec){
  p <- ((vec[1]-vec[3])/sum(vec)+1)/2
  F <- 1-vec[2]/sum(vec) /2/p/(1-p)
  return(list(p=p, F=F, 
              HWE=c(p^2, 2*p*(1-p), (1-p)^2),
              HWD=c(p^2+p*(1-p)*F, 2*p*(1-p)*(1-F), (1-p)^2+p*(1-p)*F)))
}

HWE_F(N1) # F = 0

$p
[1] 0.6

$F
[1] 0

$HWE
[1] 0.36 0.48 0.16

$HWD
[1] 0.36 0.48 0.16

HWE_F(N2) # F > 0

$p
[1] 0.6

$F
[1] 0.25

$HWE
[1] 0.36 0.48 0.16

$HWD
[1] 0.42 0.36 0.22

　
アレルA/a, B/b の組を考える。確率はそれぞれ、p_a, 1-p_a, p_b, 1-p_bである。
これがお互いになんの関係もなく、交配したとしたら、次世代の遺伝子型の確率は次のようになる。
\begin{matrix} & A:p_a & a:1-p_a & \\ B:p_b & p_ap_b & (1-p_a)p_b & p_b \\ b:1-p_b & p_a(1-p_b) & (1-p_a)(1-p_b)&1-p_b \\ & p_a & 1-p_a & 1 \end{matrix}
いま、A/a, B/b が独立ではなく、連鎖不平衡という現象によって偏って次世代に交配すると
\begin{pmatrix}AB\\Ab\\aB\\ab\end{pmatrix}=\begin{pmatrix}p_ap_b+\delta\\p_a(1-p_b)-\delta\\(1-p_a)p_b-\delta\\(1-p_a)(1-p_b)+\delta\end{pmatrix}
ただし、\delta=r\sqrt{p_a(1-p_a)p_b(1-p_b)}
AB とab から\deltaを消去すると
\begin{pmatrix}p_a+p_b\\p_a-p_b\end{pmatrix}=\begin{pmatrix}1+AB-ab\\Ab-aB\end{pmatrix}
\begin{pmatrix}1&1\\1&-1\end{pmatrix}\begin{pmatrix}p_a\\p_b\end{pmatrix}=\begin{pmatrix}1+AB-ab\\Ab-aB\end{pmatrix}
となる二次方程式をsolve で解けばよい。

M1 <- matrix(c(42, 18, 28, 12), 2) # p1 = 0.6, p2 = 0.7, r = 0
M2 <- matrix(c(48, 12, 22, 18), 2) # r ~ 0.25
LD <- function(mat){
  a <- matrix(c(1, 1, 1, -1), 2)
  b <- c(1+(mat[1,1]-mat[2,2])/sum(mat), (mat[2,1]-mat[1,2])/sum(mat))
  p12 <- solve(a, b)
  delta <- mat[1,1]/sum(mat) - prod(p12)
  LE <- matrix(c(prod(p12), p12[1]*(1-p12[2]), (1-p12[1])*p12[2], prod(1-p12)), 2)
  return(list(p = p12, r = delta/sqrt(prod(c(p12, 1-p12))),
              LE = LE, LD = LE + matrix(c(1, -1, -1, 1), 2)*delta))
}

LD(M1) # No disequilibrium

$p
[1] 0.6 0.7

$r
[1] -2.472663e-16

$LE
     [,1] [,2]
[1,] 0.42 0.28
[2,] 0.18 0.12

$LD
     [,1] [,2]
[1,] 0.42 0.28
[2,] 0.18 0.12

LD(M2) # with disequilibrium

$p
[1] 0.6 0.7

$r
[1] 0.2672612

$LE
     [,1] [,2]
[1,] 0.42 0.28
[2,] 0.18 0.12

$LD
     [,1] [,2]
[1,] 0.48 0.22
[2,] 0.12 0.18