邓飞
GCTA学习5 | GCTA计算PCA及可视化
2022-1-17 19:02
阅读:1528

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相关的知识。


转载本文请联系原作者获取授权,同时请注明本文来自邓飞科学网博客。

链接地址:https://wap.sciencenet.cn/blog-2577109-1321475.html?mobile=1

收藏

分享到:

当前推荐数:0
推荐到博客首页
网友评论0 条评论
确定删除指定的回复吗?
确定删除本博文吗?