参考文献：

Vanraden P M. Efficient methods to compute genomic predictions.[J]. Journal of Dairy Science, 2008, 91(11):4414-4423.

Forni S, Aguilar I, Misztal I. Different genomic relationship matrices for single-step analysis using phenotypic, pedigree and genomic information.[J]. Genetics, selection, evolution : GSE, 2011, 43(1):1.

http://articles.extension.org/pages/68019/genomic-relationships-and-gblup

常用的G矩阵公式

$$

• Z 表示标记矩阵的发生率（incidence）

• 表示测序个体考虑近交和选择后的MAF

计算方法

1，我们假定现在有3个二倍体的个体，每个个体有4个SNP位点信息，基因型如下

Genotype <- data.frame(row.names = paste("Ind",1:3,sep=""),
                       snp1=c("AA","AA","TT"),
                       snp2=c("CG","CG","CC"),
                       snp3=c("GG","GC","GG"),
                       snp4=c("AT","AA","AA"))

Genotype

2，现在讲基因型进行变为0,1,2的形式，0表示Major frequent（大写字母）的纯合子，2表示Minor frequent（小写字母）的纯合子，1表示杂合子。转化的文件命名为MAF矩阵

MAF <-  data.frame(row.names = paste("Ind",1:3,sep=""),
                   snp1=c(0,0,2),
                   snp2=c(1,1,0),
                   snp3=c(0,1,0),
                   snp4=c(1,0,0))

MAF

3，M矩阵，将MAF进一步转化为-1，0,1的形式

M <- MAF -1 M

MM’矩阵有一些特性，M’是M的转置矩阵

M <- as.matrix(M)

M%*%t(M)

MM’的对角线，表示个体纯合子的个数，比如个体1有两个纯合子，个体2有两个纯合子，个体3有4个纯合子

M’M的特性

t(M)%*%M

M’M的对角线表示每个位点纯合子的个体数，比如第一个位点，有3个个体是纯合的，第二个位点，有1个位点是纯合的。

4，P矩阵，P矩阵是Minor 配子的频率所计算的矩阵，P=2*(p-0.5)，用于M矩阵的中心化

比如，第一个位点，p为0.3333，第二个位点p为0.3333，第三个位点p为0.1667，第四个位点p为0.1667，则P矩阵为

p1 <- 0.3333;

p2 <- 0.3333;

p3 <- 0.1667p4 <- 0.1667;

pp <- data.frame(snp1 = rep(p1,3),snp2 = rep(p2,3),snp3=rep(p3,3),snpp4=rep(p4,3),row.names = paste("Ind",1:3,sep=""))

P <- (pp-0.5)*2P <- as.matrix(P)

P

5，Z矩阵，为M矩阵的中心化

Z <- M - P

Z

6，G矩阵的计算

tZZ <- Z%*%t(Z)

tZZ

fenmu <- 2*sum(p1*(1-p1) + p2*(1-p2) + p3*(1-p3) + p4*(1-p4))

fenmu

1.44448888

G <- tZZ/fenmu

G

7，G矩阵使用sommer软件包计算结果

sommer::A.mat(M)

结论：两者结果一致，通过手动计算和软件包计算一致。

8，自己编写函数，计算G矩阵(数据格式-1,0,1）

G_matrix <- function(data){
    M=data 
    # 计算每个SNP minor allet的频率p值
    p1=round((apply(M,2,sum)+nrow(M))/(nrow(M)*2),3) 
    # 计算P矩阵
    p=2*(p1-0.5)
    P=matrix(p,byrow=T,nrow=nrow(M),ncol=ncol(M)) 
    # 就是那Z矩阵
    Z=as.matrix(M-P)    # 计算G矩阵公式的分母
    c=p1*(1-p1)
    d=2*(sum(c))    # 计算G矩阵公式的分子
    ZZt=Z%*%t(Z)    # 计算G矩阵
    G=(ZZt/d)    return(G)
}

G_matrix(M)

9 全部R代码

Genotype <- data.frame(row.names = paste("Ind",1:3,sep=""),snp1=c("AA","AA","TT"),snp2=c("CG","CG","CC"),snp3=c("GG","GC","GG"),snp4=c("AT","AA","AA"))

Genotype

MAF <-  data.frame(row.names = paste("Ind",1:3,sep=""),snp1=c(0,0,2),snp2=c(1,1,0),snp3=c(0,1,0),snp4=c(1,0,0))

MAF

M <- MAF -1 ; M

M <- as.matrix(M)

M%*%t(M)

t(M)%*%M

p1 <- 0.3333;
p2 <- 0.3333;
p3 <- 0.1667
p4 <- 0.1667;
pp <- data.frame(snp1 = rep(p1,3),snp2 = rep(p2,3),snp3=rep(p3,3),snpp4=rep(p4,3),row.names = paste("Ind",1:3,sep=""))
P <- (pp-0.5)*2
P <- as.matrix(P);P

Z <- M - P; Z

tZZ <- Z%*%t(Z); tZZ

fenmu <- 2*sum(p1*(1-p1) + p2*(1-p2) + p3*(1-p3) + p4*(1-p4)); fenmu

G <- tZZ/fenmu; G

sommer::A.mat(M)

G_matrix <- function(data){
    M=data 
    # 计算每个SNP minor allet的频率p值
    p1=round((apply(M,2,sum)+nrow(M))/(nrow(M)*2),3) 
    # 计算P矩阵
    p=2*(p1-0.5)
    P=matrix(p,byrow=T,nrow=nrow(M),ncol=ncol(M)) 
    # 就是那Z矩阵
    Z=as.matrix(M-P)
    # 计算G矩阵公式的分母
    c=p1*(1-p1)
    d=2*(sum(c))
    # 计算G矩阵公式的分子
    ZZt=Z%*%t(Z)
    # 计算G矩阵
    G=(ZZt/d)
    return(G)
}

G_matrix(M)

​