GCTA计算PCA时,要两步走:
- 第一步:构建kinship矩阵(有两种方法)
- 第一种:Yang的方法
--make-grm 0
- 第二种:Van的方法
--make-grm-alg 1
- 第一种:Yang的方法
- 第二步: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相关的知识。
GCTA 系列教程
GCTA学习2 | 软件下载安装–windows和Linux
GCTA学习3 | GCTA的两篇NG:fast-LMM和fast-GLMM