標準化ジェノタイプ

• 0,1,2のジェノタイプデータは平均はアレル頻度に依存するし、分散もアレル頻度に依存する
• 今、遺伝率を計算するときには、全座位のもたらす「分散」の和が欲しいので、すべての座位を分散的に同じように扱えるとよい
• そうすると、線形回帰してやったとき、遺伝要因の分散は、線形回帰係数の二乗の和になるから
• Rでやってみる

# 平均0,分散1に変換する関数
my.standardize <- function(x){
	m <- mean(x)
	v <- 1/length(x) * sum((x-m)^2)
	ret <- (x-m)/sqrt(v)
	ret
}

# No. loci
n.marker <- 100
# 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))
}
# Risk of each marker's allele
r <- rnorm(n.marker)
# 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

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

g.st <- g
for(i in 1:n.marker){
	g.st[,i] <- my.standardize(g[,i])
}
lm.out <- lm(pheno~g)
pheno. <- lm.out$fitted
lm.out2 <- lm(pheno~g.st)
pheno.2 <- lm.out2$fitted
var(pheno.)
var(pheno.2)
sum(lm.out2[[1]][-1]^2)
# narrow-sense heritability
h1 <- var(R)/var(pheno)
h1
sum(lm.out2[[1]][-1]^2)/var(pheno)

> var(pheno.)
[1] 37.50983
> var(pheno.2)
[1] 37.50983
> sum(lm.out2[[1]][-1]^2)
[1] 37.53385
> # narrow-sense heritability
> h1 <- var(R)/var(pheno)
> h1
[1] 0.6962973
> sum(lm.out2[[1]][-1]^2)/var(pheno)
[1] 0.6990415