2016-03-16

ゲノムワイド情報を使って未検出座位の代用をする

前の記事では既知のリスク座位に加えて未知のリスク座位を想定した
未知のリスク座位があって、それに関する情報が皆無であれば、なすすべもない
しかしながら、個人のジェノタイプの遠近関係はゲノムワイド情報があれば定量できる
個人間の遺伝的遠近が、フェノタイプ値の遠近と関連があれば、それによって遺伝率が計算できる、というのは、血縁関係の遠近を使って遺伝率を算出した、20世紀の方法でもわかる
それを考えてみる
既知リスク座位と、未知多数座位とを発生し、人工的に構造化を入れる
そのリスクを計算し、既知座位のみでの線形回帰、既知・未知すべての線形回帰をしてみる
もちろん、構造化情報を使っていないから、遺伝率は低くなるままである

## 一部の座位はわかっている
## 多数の座位は弱い寄与だが不明
# No. loci
n.marker <- 10
# population size
N <- 10^5
# Genotype matrix
g <- matrix(0,N,n.marker)
# allele freq of markers
f <- 0.1 + runif(n.marker) * 0.4
# HWE and random genotype generation
for(i in 1:n.marker){
	g[,i] <- sample(0:2,N,replace=TRUE,prob=c(f[i]^2,2*f[i]*(1-f[i]),(1-f[i])^2))
}

# Unknown loci
n.marker.unknown <- 100
# Unknown Genotype matrix
g.unknown <- matrix(0,N,n.marker.unknown)
# allele freq of markers
f <- 0.1 + runif(n.marker) * 0.4
# HWE and random genotype generation
for(i in 1:n.marker){
	g.unknown[,i] <- sample(0:2,N,replace=TRUE,prob=c(f[i]^2,2*f[i]*(1-f[i]),(1-f[i])^2))
}
# Population structured
# 何度かサンプルの一部を取り出し、
# そのサンプルについて、すべての既知リスク座位と一部の未知座位との
# ジェノタイプをそろえることで、構造化のまねごとをする
g.st <- g
g.st.unknown <- g.unknown
n.st <- 10
for(i in 1:n.st){
	s <- sample(1:N,N/5)
	t <- sample(1:n.marker.unknown,n.marker.unknown/20)
	for(j in 1:length(t)){
		g.st.unknown[s,t[j]] <- g.unknown[s[1],t[1]]
	}
	# main loci should be also structured
	for(j in 1:n.marker){
		g.st[s,j] <- g[s[1],j]
	}
}

# Risk of each marker's allele
r <- rnorm(n.marker)

# Risk of each unknown marker's allele
# only a part of markers are truely risky
n.unknown.true <- 10
r.unknown <- rep(0,n.marker.unknown)
r.unknown[1:n.unknown.true] <- rnorm(n.unknown.true)

# Total genetic risk when addtive model
R <- rep(0,N)
for(i in 1:n.marker){
	R[which(g[,i]==1)] <- R[which(g[,i]==1)] + r[i]
	R[which(g[,i]==2)] <- R[which(g[,i]==2)] + r[i]*2

}

# Total genetic risk when addtive model
R.unknown <- rep(0,N)
for(i in 1:n.marker){
	R.unknown[which(g.unknown[,i]==1)] <- R.unknown[which(g.unknown[,i]==1)] + r.unknown[i]
	R.unknown[which(g.unknown[,i]==2)] <- R.unknown[which(g.unknown[,i]==2)] + r.unknown[i]*2
}

# Total risk
R.total <- R + R.unknown

# narrow-sense heritability
h <- 0.7
# environmental variance 
eV <- (1-h)/h * var(R.total)
# sum of genetic and environmental risks
R. <- R.total + rnorm(N,0,sqrt(eV))
# phenotype is R. itself
pheno <- R.

# narrow-sense heritability of whole loci
lm.out.whole <- lm(pheno~cbind(g,g.unknown))
pheno.whole <- lm.out.whole$fitted

hwhole <- var(pheno.whole)/var(pheno)
hwhole
# narrow-sense heritability of known loci only
lm.out.partial <- lm(pheno~g)
pheno.partial <- lm.out.partial$fitted

hpartial <- var(pheno.partial)/var(pheno)
hpartial

結果

> hwhole
[1] 0.7017883
> # narrow-sense heritability of known loci only
> lm.out.partial <- lm(pheno~g)
> pheno.partial <- lm.out.partial$fitted
> 
> hpartial <- var(pheno.partial)/var(pheno)
> hpartial
[1] 0.159487

さて集団構造化・サンプル間の遺伝的遠近の情報を使うにはどうしたらよいだろうか
ここで登場するのが、線形混合モデルである
なるモデルを立てる
- $X\beta$ は既知座位での線形和の項(固定効果項)
- はともに乱雑項で
  - $E$ はすべての個人に関して独立な正規乱数項
  - $G$ は個人の間でバラツキの出方が大きめだったり小さめだったりの凸凹を組み込んだ乱数項(人数ｘ人数の分散共分散行列)を持つ、原点を中心とした人数次元の正規乱数
- ここで $G$ の分散共分散行列を指定するにあたり、ゲノムワイドSNPによる個人間関係分散共分散行列を指定すればよい。特に、分散・共分散の値を標準化しておけば、そのスケーリング係数推定をするだけで狭義遺伝率の計算が容易になる
- 分散共分散行列は、対称行列でpositive definiteであるので、この行列をこれスキー分解したもののように式表示してあることなども気にしておくと文書を眺めるときに楽になる
実際GWASデータからこの分散・共分散行列を作ったり、後述するように連鎖不平衡を考慮して補正したりする方法については、テクニカルなことになるのでこのペイパーのOnline methodを参照するのがよい