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
注：测试数据小，结果有点失真。
