如何利用SNP信息计算亲缘关系G矩阵

生成SNP文件信息

# create marker data for 9 SNPs and 10 homozygous individuals
snp9 <- matrix(c(
  "AA",   "AA",   "AA",   "BB",   "AA",   "AA",   "AA",   "AA",  NA,
  "AA",   "AA",   "BB",   "BB",   "AA",   "AA",   "BB",   "AA",  NA,
  "AA",   "AA",   "AB",   "BB",   "AB",   "AA",   "AA",   "BB",  NA,
  "AA",   "AA",   "BB",   "BB",   "AA",   "AA",   "AA",   "AA",  NA,
  "AA",   "AA",   "BB",   "AB",   "AA",   "BB",   "BB",   "BB",  "AB",
  "AA",   "AA",   "BB",   "BB",   "AA",   NA,     "BB",   "AA",  NA,
  "AB",   "AA",   "BB",   "BB",   "BB",   "AA",   "BB",   "BB",  NA,
  "AA",   "AA",    NA,    "BB",    NA,    "AA",   "AA",   "AA",  "AA",
  "AA",    NA,     NA,    "BB",   "BB",   "BB",   "BB",   "BB",  "AA",
  "AA",    NA,    "AA",   "BB",   "BB",   "BB",   "AA",   "AA",  NA),
  ncol=9,byrow=TRUE)

# set names for markers and individuals
colnames(snp9) <- paste("SNP",1:9,sep="")
rownames(snp9) <- paste("ID",1:10+100,sep="")
str(snp9)

 chr [1:10, 1:9] "AA" "AA" "AA" "AA" "AA" "AA" "AB" "AA" ...
 - attr(*, "dimnames")=List of 2
  ..$ : chr [1:10] "ID101" "ID102" "ID103" "ID104" ...
  ..$ : chr [1:9] "SNP1" "SNP2" "SNP3" "SNP4" ...

将SNP原始文件，转化为0,1,2的格式

利用最小等位基因频率，将major为0，杂合为1，minor为2。
将缺失值随机补全。

library(synbreed)
gp <- create.gpData(geno=snp9)
gp.coded <- codeGeno(gp,impute=TRUE,impute.type="random")
geno <- gp.coded$geno
geno[1:5,1:5]

     Summary of imputation 
    total number of missing values                : 13 
    number of random imputations                  : 13 

1，利用软件包计算亲缘关系矩阵

Gmatrix <-kin(gp.coded, ret="realized")
Gmatrix[1:5,1:5]

2，利用编程手动计算亲缘关系矩阵

M1=geno
M= M1[,1:ncol(M1)]-1
p1=round((apply(M,2,sum)+nrow(M))/(nrow(M)*2),3)
p=2*(p1-.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[1:5,1:5]

由此可以看出来，两者结果是一致的