plink计算的PCA为什么和GCTA计算的不一样？

今天度过了求知的一天，求知的快乐就是这么朴实无华且枯燥。

今天同事问了我一个问题，为什么plink计算的pca和GCTA计算得不一样？然后就引出的今天的查看说明文档，也证明了世界上就怕认真二字。

我们的发现：

1，GCTA的说明文档中，有bug，公式没有写全：

第一个公式还要除以N。给出的2010 NG上有写，但是软件的说明文档中不完整。

2，GCTA计算PCA时，中间要构建G矩阵，G矩阵构建的方法有两种：

• yang # 作者的方法，默认的方法
• VanRaden #GS中GBLUP构建的G矩阵方法

两种方法计算PCA的代码：

system("gcta --bfile aaa --autosome-num 29 --make-grm --make-grm-alg 1 --out kinship")
system("gcta --grm kinship --pca 3  --out gcta")

system("gcta --bfile aaa --autosome-num 29 --make-grm --make-grm-alg 0 --out kinship0")
system("gcta --grm kinship0 --pca 3  --out gcta0")

3，plink计算PCA时，用的是yang的方法

所以，如果如果plink的PCA和GCTA的VanRaden方法相遇时，结果就不一致了。

主要是因为计算PCA时G矩阵的方法不一致。

4，怎么证明呢？

• 手动证明（自己编写代码验证）
• 使用R包的函数证明

有一个包叫AGHmatrix包，里面有个Gmatrix，它构建矩阵时可以选择构建的方法：

结果证明了两者确实不一样。

5，自己构建G矩阵，手动计算PCA

# 计算特征值和特征向量
re = eigen(Gmat)

# 计算解释百分比
por = re$values/sum(re$values)

# 整理格式
pca_re1 = re$vectors[,1:3]
pca_re2 = data.frame(pca_re1,Ind = iid)

6，plink的--kinship构建的G矩阵

plink的--kinship构建的G矩阵不是VanRaden的矩阵，而是yang的矩阵，所以很少用于GBLUP的分析

7，pca用什么方法？

推荐用Yang的方法构建G矩阵，得到的PCA结果。也就是plink的--pca的结果，同样也是gcta默认的计算PCA的参数，--make-grm-alg 0。

8，为什么要用GCTA计算PCA？

因为GCTA给出了每个PCA的特征值，可以用于计算PCA解释的百分比。plink默认没有给出所有的（应该也可以指定PCA的个数，然后手动计算，待验证）。

也可以用plink的--kinship计算yang的G矩阵，然后手动计算PCA，这样就可以计算百分比了，计算代码：

# 计算特征值和特征向量
re = eigen(Gmat)

# 计算解释百分比
por = re$values/sum(re$values)

特别致谢

同事的一个问题，让我理清了很多东西，也让我知道了“世界上就怕认真二字”，也更让我坚信“说明文档是最香的”，多查资料，相互佐证，自己验证，然后豁然开朗。