参考文献：http://articles.extension.org/pages/68019/genomic-relationships-and-gblup

生成模拟的基因型原始数据

snp <- c("AA Ct GG Ag AA Ct Ga AA tt CC GG AA")
ss <- strsplit(snp,split = " ")
snpM <- matrix(ss[[1]],3,4,byrow = T)
row.names(snpM) <- c("ID1","ID2","ID3")
colnames(snpM) <- c("snp1","snp2","snp3","snp4")
# Genotype:AA Aa aa
snpM

将基因型原始数据转化为01的数据类型

library(synbreed)
gp <- create.gpData(geno = snpM)
gpb <- codeGeno(gp,label.heter = "alleleCoding", maf = 0.01, nmiss = 0.1,
                impute = TRUE, impute.type = "random", verbose = TRUE)
agenotype <- gpb$geno
# Genotype 0 1 2
MAF <- agenotype
MAF

# Genotype -1 0 2
M <- MAF-1
M

M*M’: 对角线为每个个体含有纯合的位点，比如ID1有2个纯合位点。

M%*%t(M)

M’*M：位点纯合的个体数，比如第一个位点snp1有3个个体是纯合的

t(M)%*%M

P：P矩阵是位点频率与0.5之差的2倍

p1=round((apply(M,2,sum)+nrow(M))/(nrow(M)*2),3)

set the loc MAF is …

p1 <- c(0.383,0.244,0.167,0.067)
names(p1) <- paste("snp",1:4,sep = "")
p1
p=2*(p1-0.5)
P=matrix(p,byrow=T,nrow=nrow(M),ncol=ncol(M))
P

Z矩阵，是基因型矩阵减去P矩阵

Z=as.matrix(M-P)
Z

GOF1：G1矩阵，是根据等位基因频率计算的，默认的方法

b=1-p1
c=p1*b
d=2*(sum(c))
ZZt=Z%*%t(Z)
G1=(ZZt/d)
G1 = round(G1,3)
G1

GD2：G2矩阵，根据倒数进行一些调整

D=1/(ncol(M)*(2*p1*(1-p1))) 
G2 = Z%*%(D*t(Z))
G2 =round(G2,3)
G2

G05, let the loc MAF is 0.5：G3矩阵，定义基因频率为0.5

p1=array(0.5,ncol(M))
p=p1
P=matrix(p,byrow=T,nrow=nrow(M),ncol=ncol(M))
Z=as.matrix(M-P)
b=1-p1
c=p1*b
d=2*(sum(c))
ZZt=Z%*%t(Z)
G3=(ZZt/d)
G3

GMF:G4矩阵，使用基因型的平均值作为调整

p1=round((apply(M,2,mean)),3) 
p=2*(p1-0.5) 
P = matrix(p,byrow=T,nrow=nrow(M),ncol=ncol(M)) 
Z = as.matrix(M-P)  
b=1-p1 
c=p1*b 
d=2*(sum(c)) 
ZZt = Z %*% t(Z) 
G = (ZZt/d) 
G