【R语言】绘制PCA之二:PCA加置信度

最新推荐文章于 2024-05-16 20:03:02 发布

梦仔生信进阶

最新推荐文章于 2024-05-16 20:03:02 发布

阅读量1.5k

点赞数

分类专栏：初探-----R命令文章标签：机器学习人工智能

本文链接：https://blog.csdn.net/weixin_55649639/article/details/131764546

版权

初探-----R命令专栏收录该内容

14 篇文章 0 订阅

订阅专栏

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。