はじめに

量的遺伝学を勉強すると、ゲノム置換効果という効果が現れるが、これの解釈はなかなか難しい。結論を述べるとこの置換効果とは、線形の相加効果と非線形の優性効果を無理やり線形で近似したもので、近似で表しきれない部分は優性偏差として説明される。一方、この置換は無理やり近似したことにより、後代に遺伝子しない優性効果を含んでしまったり、集団の遺伝子頻度に依存するような問題を抱えている。

ゲノム解読技術がまだ未発達で、ゲノム情報にアクセスすることは困難であった時代において(量的遺伝学が作られた時代を指す)、推定可能な遺伝的効果は置換効果だけであり、近似であったとしてもこの効果が広く用いられていたのには納得がいく一方、ゲノム情報へのアクセスが簡単になった現代において、この置換効果がまだ使われているのは少しばかり不思議である。特に選抜により集団構造が劇的に変化する育種において、置換効果が良い指標であるとは考えにくい。

置換効果や優性偏差とは異なる表現方法として、aとdで表す方法が知られている。これは名前が特についていないためここでは相加効果a、優性効果dということにする。このaとdは遺伝子頻度によらず一意に決まるものであるためこちらを使うほうが合理的である。

そこで、ここでは育種において置換効果を使うとよくないことを示すため、以下のような設定を考え、予測精度を比較した。設定： 1. ヘテロ性の高い集団においてモデルを構築する。 2. 選抜などにより遺伝子頻度が大きく変わり、かつ個体が固定された集団において、予測を行う。

必要なパッケージと関数定義(無視してよい)

# 一般化逆行列の推定に必要
library( MASS )

# 遺伝子頻度に基づいてマーカー遺伝子型を生成する関数
func <- function( p, nInd ){
  p2 <- p^2
  q2 <- ( 1-p )^2
  h <- 2 * (1-p) * p
  nP <- round( nInd*p2 )
  nQ <- round( nInd*q2 )
  nH <- nInd - nP - nQ
  vec <- sample( c( rep(-1,nP), rep(0,nH), rep(1,nQ) ) )
  return( vec )
}

# マーカー遺伝子型から優性偏差の推定用の行列を作成する関数
Devfunc <- function( i, pVec, Mat ){
  p <- p[i]
  q <- (1-p)
  vec <- Mat[,i]
  vec[vec==0] <- -2*q^2
  vec[vec==1] <- 2*p*q
  vec[vec==2] <- -2*p^2
  return(vec)
}

# マーカー遺伝子型の個体を固定する関数
fix <- function(i,Mat){
  vec <- Mat[,i]
  vec[vec==1] <- sample(c(0,2), sum(vec==1), replace = TRUE )
  return(vec)
}

データの生成

まずは表現型データを疑似的に生成する。遺伝率は仮に0.8としておく。データ生成ではフィッシャーによる無限対立遺伝子仮説を仮定し、正規分布に従って効果を生成した。

h <- 0.8 # 遺伝率
r <- ( 1-h )/h　# 誤差分散の割合

nInd <- 100 # 個体数
nMar <- 200 # マーカー数

# 遺伝子頻度の決定
p <- sample(seq(0,0.5,0.01),nMar,replace = T)

# マーカー遺伝子型の生成 
M1 <- lapply( p, func, nInd = nInd ) 
M2 <- matrix( unlist( M1 ), ncol = nMar, byrow = FALSE )

# 相加効果aと優性効果d用のデザイン行列
Madd <- M2 # 相加効果a用
Mdom <- ( Madd + 1 ) %% 2 # 優性効果d用

# 置換効果と優性偏差用のデザイン行列
M <- M2 + 1 # 置換効果用
Mdev1 <- lapply( 1:nMar, Devfunc, pVec = p, Mat = M )
Mdev2 <- matrix( unlist( Mdev1 ), ncol = nMar, byrow = FALSE )
Mdev <- Mdev2 # 優性偏差

# 無限対立遺伝子仮説を仮定し効果を生成する。
add <- rnorm( nMar, 0, 1 ) # 相加効果a
dom <- rnorm( nMar, 0, 1 ) # 優性効果d
Inter <- matrix( rep( 1, nInd ) ) # 切片

# データ生成
y0 <- Inter  + Madd %*% add + Mdom %*% dom # 遺伝効果+切片
Vg <- var(y0) # 遺伝分散の計算
Ve <- r*Vg # 誤差分散の計算
y <- y0 + rnorm( nInd, 0, sqrt(Ve) ) # 誤差の付与

# 遺伝分散に対する優性効果の割合の計算
Ratio <- var(Mdom %*% dom) / var(Madd %*% add+Mdom %*% dom)

y0とyを合わせてプロットすると次のような図が表示される。遺伝率0.8だとかなりきれいにy=xの軸に乗る。

plot( y0, y, xlab = "Genotypic_Value + Intercept", ylab = "Phenotype_Value", pch = 19 )
abline(0,1)

また、相加効果、優性効果、誤差分散の割合を可視化すると次のようになる。

barplot( t( as.matrix( c(Vg*Ratio, Vg*(1-Ratio), Ve ) ) ), 
         names.arg = c("Dominance", "Additive", "Residual") )

(addとdomの間に共分散がないと仮定していることに注意)

パラメータ推定

それではこの生成されたデータをもとに相加効果a、優性効果d、置換効果、優性偏差の推定を行う。正確に言うと推定されたパラメータ自体は遺伝効果ではなくマーカー遺伝子型と掛け合わされて遺伝効果となるが、言葉がややこしいので上記のように書くこととする。

# 相加効果aと優性効果d
X1 <- cbind( Inter, Madd, Mdom )
b1 <- ginv( t(X1) %*% X1 ) %*% t(X1) %*% y

# 置換効果と優性偏差
X2 <- cbind( Inter, M, Mdev )
b2 <- ginv( t(X2) %*% X2 ) %*% t(X2) %*% y

ここでb1とb2がそれぞれ推定されたパラメータであるが、二つの効果が混ざっているので分離する。

aPos <- 1:nMar + 1 
dPos <- 1:nMar + nMar + 1

# 相加効果aと優性効果dの分離
add1 <- X1[,aPos] %*% b1[aPos] # 相加効果
dom1 <- X1[,dPos] %*% b1[dPos] # 優性効果

# 置換効果と優性偏差の分離
add2 <- X2[,aPos] %*% b2[aPos] # 置換効果
dom2 <- X2[,dPos] %*% b2[dPos] # 優性偏差

相加効果aと優性効果d、置換効果と優性偏差はそれぞれの効果は共分散を持たないのが好ましい。どれくらい共分散をもって推定されているかをここで調べておくことにする。

barplot( c(cov(add1, dom1),cov(add2, dom2)), names.arg = c("相加効果aと優性効果d", "置換効果と優性偏差" ) )

ここでは後者の推定のほうが共分散が大きくなっているようである。

集団の更新と固定

今回の目的は集団の遺伝子頻度が大きく変動し、そのうえで固定された場合に予測精度がどうなるかを調べることである。ここはかなり適当であるが、遺伝子頻度を新たに定義し、ヘテロの遺伝子(1で符号化されている)を0か1にランダムに振り分けることでそれを実現する。

# 遺伝子頻度の更新
p2 <- sample(seq(0,0.5,0.01),nMar,replace = T)

# 遺伝子頻度に基づいてマーカー遺伝子型を生成
M1.2 <- lapply( p2, func, nInd = nInd )
M2.2 <- matrix( unlist( M1.2 ), ncol = nMar, byrow = FALSE )

# ヘテロになっている部分を無理やり固定
Mfix1 <- lapply( 1:nMar, fix, Mat = (M2.2 + 1) )
Mfix2 <- matrix( unlist( Mfix1 ), ncol = nMar, byrow = FALSE )

# 相加効果aと優性効果d用のデザイン行列
Mfadd <- Mfix2 - 1
Mfdom <- Mfix2 %% 2

# 優性偏差用のデザイン行列
Mfdev1 <- lapply( 1:nMar, Devfunc, pVec = p2, Mat = Mfix2 )
Mfdev2 <- matrix( unlist( Mfdev1 ), ncol = nMar, byrow = FALSE )

新しい集団での推定

今作り出した集団で新しく効果を推定する。置換効果はdを含んでしまうが、新しい集団では固定化によりdの効果はなくなっているため、置換効果での推定値は真値とずれることが予想される。

# 相加効果aと優性効果d
X3 <- cbind( Inter, Mfadd, Mfdom )
y3 <- X3%*%b1

# 置換効果と優性偏差
X4 <- cbind( Inter, Mfix2, Mfdev2 )
y4 <- X4%*%b2

# 真値
ytrue <- Inter  + Mfadd %*% add + Mfdom %*% dom

真値との相関を調べたのが次の図である。

barplot( c(cor(ytrue, y3),cor(ytrue, y4)), names.arg = c("相加効果aと優性効果d", "置換効果と優性偏差" ) )

反復を設けて統計的に調査

最後にこれまでの施行を1000回行い、どちらで遺伝効果を調べたほうが良いかを調査する。関数の中身はこれまでとほとんど同じである。

nRep <- 10000
store0 <- rep(NA,nRep)
store1 <- rep(NA,nRep)
store2 <- rep(NA,nRep)
store3 <- rep(NA,nRep)
store4 <- rep(NA,nRep)

for( i in 1:nRep ){
  
  nInd <- 100
  nMar <- 200
  
  aPos <- 1:nMar + 1
  dPos <- 1:nMar + nMar + 1
  
  p <- sample(seq(0,0.5,0.01),nMar,replace = T)
  
  M1 <- lapply( p, func, nInd = nInd )
  M2 <- matrix( unlist( M1 ), ncol = nMar, byrow = FALSE )
  
  Madd <- M2
  Mdom <- ( Madd + 1 ) %% 2
  
  M <- M2 + 1
  Mdev1 <- lapply( 1:nMar, Devfunc, pVec = p, Mat = M )
  Mdev2 <- matrix( unlist( Mdev1 ), ncol = nMar, byrow = FALSE )
  Mdev <- Mdev2
  
  add <- rnorm( nMar, 0, 1 )
  dom <- rnorm( nMar, 0, 1 )
  Inter <- matrix( rep( 1, nInd ) )
  
  y0 <- Inter  + Madd %*% add + Mdom %*% dom
  Vg <- var(y0)
  Ve <- r*Vg
  y <- y0 + rnorm( nInd, 0, sqrt(Ve) )
  
  Ratio <- var(Mdom %*% dom) / var(Madd %*% add+Mdom %*% dom)
  store0[i] <- Ratio
  
  X1 <- cbind( Inter, Madd, Mdom )
  b1 <- ginv( t(X1) %*% X1 ) %*% t(X1) %*% y
  
  X2 <- cbind( Inter, M, Mdev )
  b2 <- ginv( t(X2) %*% X2 ) %*% t(X2) %*% y

  add1 <- X1[,aPos] %*% b1[aPos]
  dom1 <- X1[,dPos] %*% b1[dPos]
  
  add2 <- X2[,aPos] %*% b2[aPos]
  dom2 <- X2[,dPos] %*% b2[dPos]
  
  store1[i] <- cov( add1, dom1 )
  store2[i] <- cov( add2, dom2 )
  
  p2 <- sample(seq(0,0.5,0.01),nMar,replace = T)
  
  M1.2 <- lapply( p2, func, nInd = nInd )
  M2.2 <- matrix( unlist( M1.2 ), ncol = nMar, byrow = FALSE )
  
  Mfix1 <- lapply( 1:nMar, fix, Mat = (M2.2 + 1) )
  Mfix2 <- matrix( unlist( Mfix1 ), ncol = nMar, byrow = FALSE )
  
  Mfadd <- Mfix2 - 1
  Mfdom <- Mfix2 %% 2
  
  Mfdev1 <- lapply( 1:nMar, Devfunc, pVec = p2, Mat = Mfix2 )
  Mfdev2 <- matrix( unlist( Mfdev1 ), ncol = nMar, byrow = FALSE )
  
  X3 <- cbind( Inter, Mfadd, Mfdom )
  y3 <- X3%*%b1
  
  X4 <- cbind( Inter, Mfix2, Mfdev2 )
  y4 <- X4%*%b2
  
  ytrue <- Inter  + Mfadd %*% add + Mfdom %*% dom
  
  store3[i] <- cor( ytrue, y3 )
  store4[i] <- cor( ytrue, y4 )
}

結果の可視化

上記の計算結果を可視化する。色は遺伝効果のうち優性効果が占める割合(Ratio)をに基づいてついている。

plot(store1, store2, xlab = "相加効果aと優性効果d", ylab = "置換効果と優性偏差", main = "Covariance", col = rainbow(nMar)[order(store0)], pch = 19 )
abline(0,1)

plot(store3, store4, xlab = "相加効果aと優性効果d", ylab = "置換効果と優性偏差", col = rainbow(nMar)[order(store0)], pch = 19 )
abline(0,1)

思いのほか大きな差は現れなかったが、相加効果aと優性効果dを用いたほうが置換効果と優性偏差を用いた場合よりもわずかに推定精度が上昇した。おそらくマーカー効果の推定精度がそもそもあまり高くないのが理由であろう。また、面白いことに、置換効果と優性偏差を用いるよりも相加効果aと優性効果dを用いたほう共分散が小さくなることが判明した。なぜその違いが生じているかは自明ではなく、今後考えていく必要がある。

マーカーの符号化に関して

Soh Kimura

2025-09-29