plink + R构建G矩阵

1、测试数据

root@PC1:/home/test# ls
outcome.map  outcome.ped
root@PC1:/home/test# cat outcome.map
1       snp1    0       55910
1       snp2    0       85204
1       snp3    0       122948
1       snp4    0       203750
1       snp5    0       312707
1       snp6    0       356863
1       snp7    0       400518
1       snp8    0       487423
1       snp9    0       578716
1       snp10   0       639876
root@PC1:/home/test# cat outcome.ped
DOR 1 0 0 0 -9 C C C C T T G G A G A A G G G C A G T T
DOR 2 0 0 0 -9 G G G G T G G G G G A A A G C C A G T T
DOR 3 0 0 0 -9 G G G G T G G G G G A A A G G C G G T A
DOR 4 0 0 0 -9 G G C C G G G G G G C C G G G G G G A A
DOR 5 0 0 0 -9 G G C C G G G A G G C C A G G C G G A A
DOR 6 0 0 0 -9 G G C C G G A A G G A A A A C C G G A A
DOR 7 0 0 0 -9 G G C C G G A G A A A A G G C C A A A A
DOR 9 0 0 0 -9 G G C C G G A G A A A A G G C C A A A A

2、使用plink软件转换为0、1、2次等位基因个数的形式

root@PC1:/home/test# ls
outcome.map  outcome.ped
root@PC1:/home/test# plink --file outcome --recode A --out temp > /dev/null; rm *log *nosex
root@PC1:/home/test# ls
outcome.map  outcome.ped  temp.raw
root@PC1:/home/test# cat temp.raw
FID IID PAT MAT SEX PHENOTYPE snp1_C snp2_G snp3_T snp4_A snp5_A snp6_C snp7_A snp8_G snp9_A snp10_T
DOR 1 0 0 0 -9 2 0 2 0 1 0 0 1 1 2
DOR 2 0 0 0 -9 0 2 1 0 0 0 1 0 1 2
DOR 3 0 0 0 -9 0 2 1 0 0 0 1 1 0 1
DOR 4 0 0 0 -9 0 0 0 0 0 2 0 2 0 0
DOR 5 0 0 0 -9 0 0 0 1 0 2 1 1 0 0
DOR 6 0 0 0 -9 0 0 0 2 0 0 2 0 0 0
DOR 7 0 0 0 -9 0 0 0 1 2 0 0 0 2 0
DOR 9 0 0 0 -9 0 0 0 1 2 0 0 0 2 0

3、使用R计算G矩阵

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

G <- read.table("temp.raw", header = T)
id <- G$IID
G <- G[,-(1:6)]
rownames(G) <- id

maf <- apply(G, 2, sum)/(2 * nrow(G))
sum2pq <- sum(2 * maf * (1 - maf))      ## 分母部分

temp1 <- 2 * (maf - 0.5)
temp2 <- matrix(temp1, byrow = T, nrow = nrow(G), ncol = ncol(G))
Z <- as.matrix(G - 1 - temp2)           ## 分子部分

Gmatrix <- Z %*% t(Z) / sum2pq          ## 计算G矩阵
Gmatrix

注：测试数据小，结果有点失真。