首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >PCA分析 | 不同品种的基因型数据绘制2D和3D的PCA图

PCA分析 | 不同品种的基因型数据绘制2D和3D的PCA图

作者头像
邓飞
发布2021-03-30 10:50:26
2K0
发布2021-03-30 10:50:26
举报

PCA是降维的一种方法。

很多软件可以分析PCA,这里介绍一下使用plink软件和R语言,进行PCA分析,并且使用ggplot2绘制2D和3D的PCA图。

绘制后的图如下:

2-D PCA图:

图片解释,将每个品种用不同的颜色表示,同时绘制置信区间圆圈,X坐标是PC1,解释24.9%的变异,Y坐标是PC2,解释10.61%的变异。可以看到,三个品种在PCA图里面分的比较开,C品种的有两个A和B的点,应该是异常数据。

3-D PCA图:

图片解释,将每个品种用不同的颜色表示,X坐标是PC1,解释24.9%的变异,Y坐标是PC2,解释10.61%的变异,Z坐标是PC3,解释1.02%的变异。可以看到,三个品种在PCA图里面分的比较开,C品种的有两个A和B的点,应该是异常数据。

基因型数据:

共有3个品种A,B,C,共有412个个体。其中:

  • A品种有200个体
  • B品种有100个体
  • C品种有112个体
$ wc -l re2.ped
412 re2.ped

SNP个数为:41013

$ wc -l re2.map
41013 re2.map

计算思路:

1,对数据进行清洗,将其转化为0,1,2的形式

2,计算G矩阵

3,计算PCA的特征向量和特征值

4,根据特征值计算解释百分比

5,根据特征向量和品种标签,进行PCA的绘制

绘制代码如下:

首先,使用plink命令,将基因型数据转化为012的raw格式:

plink --file re2 --recodeA

结果生成plink.raw文件。

然后使用R语言,计算PCA,并绘制PCA图。

读取数据
m012 = fread("plink.raw")

# 保留FID,IID和基因型数据
g012 = m012[,-c(3:6)]
dim(g012)
fid = g012$FID
iid = g012$IID
library(sommer)

# 整理格式,计算G矩阵
setDF(g012)
rownames(g012) = g012$IID
g012$IID = NULL
g012$FID = NULL
Gmat = A.mat(g012-1)

# 计算特征值和特征向量
re = eigen(Gmat)

# 计算解释百分比
por = re$values/sum(re$values)

# 整理格式
pca_re1 = re$vectors[,1:3]
pca_re2 = data.frame(pca_re1,Ind = iid)
pca_re2$Gen = fid

# 把PCA1,PC2,PC2的j解释百分比,命名为相应的轴
xlab = paste0("PC1(",round(por[1]*100,2),"%)")
ylab = paste0("PC2(",round(por[2]*100,2),"%)")
zlab = paste0("PC3(",round(por[3]*100,2),"%)")

# 绘制2-D PCA图
ggplot(pca_re2, aes(x=X1, y=X2,color=Gen))  + geom_point(size=2) +
  # stat_ellipse(level = 0.95, size = 1) + 
  stat_ellipse(aes(fill=Gen),
               type ="norm", geom ="polygon",alpha=0.2,color=NA)+
  geom_hline(yintercept = 0)  + # 添加x坐标
  geom_vline(xintercept = 0) + # 添加y坐标
  labs(x = xlab,y = ylab,color="")+
  guides(fill=F)+
  theme_bw() # 主题

# 绘制3-D PCA图
library(scatterplot3d)
pca_re2 = pca_re2 %>% mutate(colour = case_when(
  Gen == "A" ~ "red",
  Gen == "B" ~ "green",
  Gen == "C" ~ "blue",
))

scatterplot3d(pca_re2[,1:3],color=pca_re2$colour,
              pch = 16,angle=30,
              box=T,type="p",
              xlab = xlab, ylab = ylab, zlab = zlab,main = "3D PCA Plot",
              lty.hide=2,lty.grid = 2)

legend("topright",c("A","B","C"),fill=c('red','green',"blue"))
本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2021-03-18,如有侵权请联系 cloudcommunity@tencent.com 删除

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

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

本文参与 腾讯云自媒体分享计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档