GCTA计算PCA时，要两步走：

• 第一步：构建kinship矩阵（有两种方法）
• 第一种：Yang的方法--make-grm 0
• 第二种：Van的方法--make-grm-alg 1
• 第二步：PCA计算

1. GCTA构建kinship的方法

使用--make-grm-alg进行设置，0位Yang的方法，1位Van的方法。「Yang的方法：」

GRM = sum{[(xij- 2pi)*(xik- 2pi)] / [2pi(1-pi)]}/N

「Van的方法：」

GRM = sum[(xij- 2pi)(xik- 2pi)] / sum[2pi(1-pi)]

2. GCTA计算PCA

2.1 Yang的方法

gcta64 --bfile ../test --make-grm --make-grm-alg 0 --out kinship_yang
gcta64 --grm kinship_yang --pca 20 --out pca_re

结果生成：

pca_re.eigenval  pca_re.eigenvec  pca_re.log

2.2 Van的方法

这里，将0变为1.

gcta64 --bfile ../test --make-grm --make-grm-alg 1 --out kinship_yang
gcta64 --grm kinship_yang --pca 20 --out pca_re

结果生成：

pca_re.eigenval  pca_re.eigenvec  pca_re.log

3. PCA可视化

这里，先对数据进行处理，计算每个主成分解释百分比，以及前几个PCA的累计百分比。

pcaal = fread("pca_re.eigenval")
head(pcaal)
pcaal$por = pcaal$V1/sum(pcaal$V1)
pcaal$cumula = cumsum(pcaal$por)
pcaal$index = as.factor(1:dim(pcaal)[1])
head(pcaal)

3.1 碎石图（折线图）

这里，选择前10个主成分。

# 选择最佳的PCA个数：碎石折线图
pcaal[1:10,] %>%
  ggplot(aes(x=index,y=por, group=1))+
  geom_point(size=4)+
  geom_line()+
  labs(title="Scree plot: PCA")

3.2 碎石图（条形图）

这里，选择前10个主成分。

# 选择最佳的PCA个数：碎石条形图
pcaal[1:10,] %>%
  ggplot(aes(x=index,y=por))+
  geom_col()+
  labs(title="Scree plot: PCA on scaled data")

3.3 PCA可视化

「代码：」

ggplot(pcaec,aes(x = V3,y = V4)) + geom_point() + 
  xlab(paste0("PC1 ",round(pcaal$por[1],4),"%")) + 
  ylab(paste0("PC2 ",round(pcaal$por[2],4),"%"))

4. PCA分析拓展

1，PCA分析，可以根据分组，绘制置信区间

2，PCA分析中，可以将PCA的百分比和累计百分比绘制到一张图上面。

❝

欢迎关注我的公众号：育种数据分析之放飞自我。主要分享R语言，Python，育种数据分析，生物统计，数量遗传学，混合线性模型，GWAS和GS相关的知识。

❞