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

邓飞

发布于 2022-02-09 09:43:12

1.5K0

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

4. PCA分析拓展

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

分组PCA

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

相关系列阅读：

第一篇：GCTA学习1 | 抛砖引玉--初步介绍

第二篇：GCTA学习2 | 软件下载安装--windows和Linux

第三篇：GCTA学习3 | GCTA的两篇NG：fast-LMM和fast-GLMM

第四篇：GCTA学习4 | GCTA说明文档--功能分类及常见问题

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2022-01-12，如有侵权请联系 cloudcommunity@tencent.com 删除

数据分析

本文分享自育种数据分析之放飞自我微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

数据分析

登录后参与评论

0 条评论

热度