|
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个个体。其中：
$ wc -l re2.ped412 re2.pedSNP个数为：41013 $ wc -l re2.map41013 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$FIDiid = g012$IIDlibrary(sommer)# 整理格式，计算G矩阵setDF(g012)rownames(g012) = g012$IIDg012$IID = NULLg012$FID = NULLGmat = 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"))另外，我决定每个月开一次赞赏！如果您感觉我的文章对您有所帮助，还请赞赏支持一下，一元也是情，留个头像行不行 |
|
|
。您的鼓励是我最大的动力，十分感谢
。
