【原】PCA分析 | 不同品種的基因型數(shù)據(jù)繪制2D和3D的PCA圖

育種數(shù)據(jù)分析 2021-11-18

展開全文

PCA是降維的一種方法。

很多軟件可以分析PCA，這里介紹一下使用plink軟件和R語言，進(jìn)行PCA分析，并且使用ggplot2繪制2D和3D的PCA圖。

繪制后的圖如下：

2-D PCA圖：

圖片解釋，將每個(gè)品種用不同的顏色表示，同時(shí)繪制置信區(qū)間圓圈，X坐標(biāo)是PC1，解釋24.9%的變異，Y坐標(biāo)是PC2，解釋10.61%的變異?？梢钥吹剑齻€(gè)品種在PCA圖里面分的比較開，C品種的有兩個(gè)A和B的點(diǎn)，應(yīng)該是異常數(shù)據(jù)。

3-D PCA圖：

圖片解釋，將每個(gè)品種用不同的顏色表示，X坐標(biāo)是PC1，解釋24.9%的變異，Y坐標(biāo)是PC2，解釋10.61%的變異，Z坐標(biāo)是PC3，解釋1.02%的變異。可以看到，三個(gè)品種在PCA圖里面分的比較開，C品種的有兩個(gè)A和B的點(diǎn)，應(yīng)該是異常數(shù)據(jù)。

基因型數(shù)據(jù)：

共有3個(gè)品種A，B，C，共有412個(gè)個(gè)體。其中：

A品種有200個(gè)體
B品種有100個(gè)體
C品種有112個(gè)體

$ wc -l re2.ped412 re2.ped

SNP個(gè)數(shù)為：41013

$ wc -l re2.map41013 re2.map

計(jì)算思路：

1，對數(shù)據(jù)進(jìn)行清洗，將其轉(zhuǎn)化為0,1,2的形式

2，計(jì)算G矩陣

3，計(jì)算PCA的特征向量和特征值

4，根據(jù)特征值計(jì)算解釋百分比

5，根據(jù)特征向量和品種標(biāo)簽，進(jìn)行PCA的繪制

繪制代碼如下：

首先，使用plink命令，將基因型數(shù)據(jù)轉(zhuǎn)化為012的raw格式：

plink --file re2 --recodeA

結(jié)果生成plink.raw文件。

然后使用R語言，計(jì)算PCA，并繪制PCA圖。

讀取數(shù)據(jù)m012 = fread("plink.raw")

# 保留FID，IID和基因型數(shù)據(jù)g012 = m012[,-c(3:6)]dim(g012)fid = g012$FIDiid = g012$IIDlibrary(sommer)

# 整理格式，計(jì)算G矩陣setDF(g012)rownames(g012) = g012$IIDg012$IID = NULLg012$FID = NULLGmat = A.mat(g012-1)

# 計(jì)算特征值和特征向量re = eigen(Gmat)

# 計(jì)算解釋百分比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解釋百分比，命名為相應(yīng)的軸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坐標(biāo) geom_vline(xintercept = 0) + # 添加y坐標(biāo) 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"))

另外，我決定每個(gè)月開一次贊賞！如果您感覺我的文章對您有所幫助，還請贊賞支持一下，一元也是情，留個(gè)頭像行不行。您的鼓勵(lì)是我最大的動力，十分感謝。