

GCTA学习5 | GCTA计算PCA及可视化

育种数据分析之放飞自我

2022-01-12

导读：GCTA绘制PCA。

GCTA学习5 | GCTA计算PCA及可视化 _#2022.01.12

微信公众号真的要凉了吗，昨天的阅读量才290，是时间亮出我的镇楼图了（上图就是）。

GCTA这款软件，写了几篇了，后面将介绍单性状遗传力评估，以及多性状遗传力和遗传相关评估，感觉它与传统的评估软件，比如ASReml，DMU比较像，但是使用范围上更偏向医学。它的显著特征是速度快，里面还有很多GWAS方面不同模型的参数，真是一款强大的软件啊。“取法于上，仅得为中，取法于中，故为其下。”我学习好的软件，希望掌握个中不溜，就很不错了。

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),"%"))