PCA分析
GCTA
gcta64 --bfile bestqc --make-grm --make-grm-alg 0 --out kinship0 ;
gcta64 --grm kinship0 --pca 3 --out gcta
Plink
plink --file bestqc --pca 5 header tabs
具体可见我之前总结的：CSDN
6.1 利用plink计算pca-----直接出结果
plink --file plink --pca 3 header --out plink_pca
#如
plink --file plink --pca 3 header --out plink_pca3 --chr-set 30
#默认是提取前20个pca
# header 结果文件加表头，默认没有表头
out.file 特征值.eigenval .eigenvec
#解释的方差 （PCA1/（PCA1 + PCA2 +PCA3））
结果文件具体内容
6.2 利用GCTA计算pca-----分2步；先计算亲缘关系矩阵---再计算PCA
1.转bed格式
# 因为我没有fam文件，用plink转换二进制文件
# 省略 plink --file plink --make-bed --out plink --chr-set 30
2.计算亲缘关系
gcta_1.94.0beta/gcta64 --bfile plink --make-grm --make-grm-alg 0 --out plink_GCTA_Yang
# GCTA主要是分析人的数据，如果分析其他动植物需要对染色体体进行重新编号，不然会报错：
Error: Line 56287 of [plink.bim] contains illegal chr number, please check
主要做人的 要求为数字染色体:
#用plink提取染色体 plink --file plink --make-bed --out plink --chr-set 30 --chr 1-13
3.pca
gcta_1.94.0beta/gcta64 --grm plink_GCTA_Yang --pca 3 --out plink_GCTA_Yang_pca3
PCA可视化做主成分分析加"置信区间
head(iris) ord <-prcomp(iris[,1:4]) #R语言中，PCA的函数 princomp ， prcomp ord <-princomp(iris[,1:4]) summary(ord) summ <- summary(ord) #PCA解释度 xlab <- paste0("PC1(",round(summ$importance[2,1]*100,2),"%)") ylab <- paste0("PC2(",round(summ$importance[2,2]*100,2),"%)") #我们看到PC1的方差解释率达92.46%，PC2的方差解释率为5.31%， dt<-ord$x df<-data.frame(dt,iris$Species) head(df) library(ggplot2) ggplot(data = df,aes(x=PC1,y=PC2,color=iris.Species))+ # 绘制置信椭圆： stat_ellipse(aes(fill=iris.Species), type = "norm", geom = 'polygon',alpha=0.2,color=NA,linetype = 1)+ #geom = 'polygon' 填充 geom_point()+labs(x=xlab,y=ylab,color="")+ guides(fill=F) ggplot(data = df,aes(x=PC1,y=PC2,color=iris.Species))+ # 绘制置信椭圆： stat_ellipse(aes(fill=iris.Species), type = "norm", geom = 'polygon',alpha=0.2,color=NA,linetype = 2, show.legend = FALSE)+ #geom = 'polygon' 填充 #默认的“t”假设多元t分布，而“norm”假设多元正态分布 geom_point()+labs(x=xlab,y=ylab,color="")+ guides(fill=F) p1
#参考 R语言教程|如何做主成分分析加"置信区间"？ - 知乎
#每个群组加上置信椭圆
stat_ellipse( mapping = NULL, data = NULL, geom = "path", position = "identity", ..., type = "t", level = 0.95, segments = 51, na.rm = FALSE, show.legend = NA, inherit.aes = TRUE ) #stat_ellipse 函数默认对每个类别的数据计算自己的置信区间。 #stat_ellipse 函数默认会使用 ggplot 中的 aes 设置， #如果希望 stat_ellipse 使用自己的 aes 设置，需要将参数 inherit.aes 设置为 FALSE。
level:置信区间水平，一般0.95/ 0.99
根据0.9的置信区间水平添加的点圈，可以根据自己的需求调节置信水平，将尽量多的点囊括进来
PCA 置信区间
我们需要根据点的分布计算他们的置信区间。
type = "t" 假设数据服从T分布
type = "norm" 假设数据服从正态分布（常用）
