はじめに
私が博士課程で研究をしている量的遺伝学と呼ばれる分野は、ゲノム情報が簡単に読めるようになる前となった後でモデルが大きく変化しました。特に混合モデルと呼ばれるモデルについては、ゲノム情報が読まれる前はABLUPと呼ばれる系譜情報を元にしたモデルが主流でしたが、ゲノム情報が読まれるようになってからはGBLUPと呼ばれるゲノム情報を元にしたモデルが取って変わりました。どちらも基本的にモデルの形は一緒で、それぞれ次のように表されます。
ABLUP
\mathbf{y}=\mathbf{X}\mathbf{b}+\mathbf{Z}\mathbf{u}+\mathbf{e}\space where\space \mathbf{u}\sim MVN(\mathbf{0},\mathbf{A}\sigma_a^2)
GBLUP
\mathbf{y}=\mathbf{X}\mathbf{b}+\mathbf{Z}\mathbf{u}+\mathbf{e}\space where\space \mathbf{u}\sim MVN(\mathbf{0},\mathbf{G}\sigma_g^2)
要するに背景に仮定してる分散共分散構造が$A$なのか$G$なのかという違いしかありません。ここで出てくる$A$という行列は、系譜情報に基づいて、ゲノムがランダムに落ちてきたときに仮定される分散共分散構造で、詳しくはこちらの記事にまとめてあります。一方の$G$行列はゲノム情報を元に計算される行列で、系譜に基づいてゲノムがランダムに落ちてきたある試行に基づいて計算される行列です。つまり$A$行列の実現値であり、なぜそれが分散共分散構造を持つのか、についてはなかなかイメージがつきづらいかもしれません。
よくある説明が、マーカー効果自体が確率的に発生しており、その効果を引きずっているからというものです。つまり、たとえどんな効果を持つ生物種であろうと期待される分散・共分散が$G$行列に反映されているという主張です。言葉ではなかなかイメージがつきずらいと思いますので今回はシミュレーションを組んでみました。本当に簡単なものですが参考になれば幸いです。
シミュレーション実装内容
こちらが今回書いたシミュレーション内容です。個体数を多くし過ぎると計算が重くなるので、今回は2個体だけを選んできています。やっていることは単純で、毎回マーカー効果をある分布に従って生成し、遺伝子型値を生成し、1000回の反復を元に、分散共分散行列を推定するというものです。
library(gaston)
genome <- read.vcf("xxx.vcf") # すみません、ゲノム情報は非公開です
# 二個体だけ選び名前を付ける
genomeSel <- genome[1:2, ]
gMat <- as.matrix(genomeSel)
rownames(gMat) <- c("Indiv1", "Indiv2")
# マーカー効果を発生させそれを元に遺伝子型値を決めることを1000回反復する
gEst <- sapply(1:1000, function(i) {
effect <- rnorm(ncol(gMat), 0, 1)
g <- gMat %*% effect
return(g)
})
v11 <- character(0)
v22 <- character(0)
v12 <- character(0)
# 累積で分散を取っていく(sapplyだと重いのでfor文にしました)
for (i in 2:1000) {
print(i)
V <- var(t(gEst[, 1:i]))
v11 <- c(v11, V[1, 1])
v22 <- c(v22, V[2, 2])
v12 <- c(v12, V[1, 2])
}
# こちらは理論値。
GRM <- (gMat %*% t(gMat))
# プロット
plot(
v11,
pch = 19,
cex = 0.5,
col = 2,
ylim = c(min(GRM) - 10000, max(GRM) + 10000),
ylab = "Variance"
)
points(v22, pch = 19, cex = 0.5, col = 3)
points(v12, pch = 19, cex = 0.5, col = 4)
abline(h = c(GRM[1, 1], GRM[2, 2], GRM[1, 2]), col = 1, lwd = 1.5)
legend(
x = "topright",
legend = c("GRM11", "GRM22", "GRM12"),
col = 1:3,
pch = 19
)
ということで実際にこのコードを回してみた結果がこちらです(legendは適当に移動させましたm(__)m)。黒の線でかかれているところが理論値で、どの線もそこに向かって収束していることが分かります。ということで、$G$行列の分散共分散構造は、どんなマーカー効果を持っているかに依らず、マーカー遺伝子型によって決まってくる分散共分散構造である、ということが確かに確かめられました。めでたしめでたし。
