R語言繪圖：PCA分析和散點(diǎn)圖(共7頁)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上PCA分析和散點(diǎn)圖gaom 今天主要跟大家演示一下簡單的PCA分析，并且以散點(diǎn)圖的形式將結(jié)果展示出來。 首先在進(jìn)行PCA分析之前，先跟大家稍微討論下什么是PCA分析。PCA分析又叫主成分分析，其實(shí)從字面上來理解我們可以發(fā)現(xiàn)它其實(shí)是和樣品分組相關(guān)的。舉個簡單的例子，我們觀察了某種植物的株高、葉片大小、果實(shí)大小等等多種性狀，并記錄每種性狀對應(yīng)的數(shù)值。這時候我們想看看根據(jù)這些性狀信息看看我們觀察的樣本是否明顯的分組現(xiàn)象。每一種性狀相當(dāng)于一個維度。利用PCA分析可以將結(jié)果投影到一個低維的向量空間（具體計算就不詳述了）。類似的比如我們多個樣本的表達(dá)譜數(shù)據(jù)，每個基因在各個樣品的表

2、達(dá)情況就可以算作一個維度。如果大家對PCA算法感興趣的話，可以自行百度，在這里就不進(jìn)行太多的描述了。畢竟今天主要是教大家怎么利用R進(jìn)行PCA分析和結(jié)果展示。 還是第一步，我們先準(zhǔn)備好我們用來分析的數(shù)據(jù)。setwd("C:/Users/gaom/Desktop")#打開文件所在路徑，并將文件所在目錄作為工作目錄data<-read.table(file = "test_data.txt",header = T,sep = "t")#讀取數(shù)據(jù)，并將首行作為列名dim(data)# 1 2999 13head(data)# ID_RE

3、F T01 T02 T03 T04 T05 T06# 1 1007_s_at 10. 11. 10. 11. 12. 12.# 2 1053_at 9. 8. 9. 9. 9. 9.# 3 117_at 8. 8. 8. 8. 8. 8.# 4 121_at 12. 12. 12. 12. 12. 12.# 5 1255_g_at 7. 7. 7. 7. 7. 7.# 6 1294_at 8. 9. 9. 8. 8. 9.# T07 T08 T09 T10 T11 T12# 1 10. 10. 10. 11. 10. 11.# 2 9. 9. 8. 8. 9. 8.# 3 8. 8. 8. 8

4、. 8. 8.# 4 13. 13. 13. 13. 13. 13.# 5 7. 7. 7. 7. 7. 7.# 6 9. 9. 9. 9. 9. 9. 上述數(shù)據(jù)為從GEO數(shù)據(jù)庫隨意找的基因表達(dá)。其中第一列為基因探針號，后續(xù)幾列則為T01到T12的12個樣品對應(yīng)的表達(dá)量數(shù)據(jù)，每三個樣品為一組。因?yàn)閿?shù)據(jù)是拼湊的，所以這里不關(guān)注探針具體信息了。 準(zhǔn)備好數(shù)據(jù)之后我們就開始進(jìn)行PCA計算了。其實(shí)代碼非常簡單。pca<- prcomp(t(data,-1), scale=T)head(pca$x)# PC1 PC2 PC3 PC4 PC5 PC6# T01 -43. -44. 8. 3. -7.

5、14.# T02 42. -19. -25. 21. -5. 14.# T03 -2. -21. -11. 17. 15. -34.# T04 8. -4. 22. -11. 8. -6.# T05 18. -5. 18. -10. -0. -5.# T06 27. -5. 26. -11. -4. 5.# PC7 PC8 PC9 PC10 PC11 PC12# T01 -1. 9. -1. -1. 3. 2.e-13# T02 5. 8. 9. 1. 1. 2.e-13# T03 -5. -9. -11. 2. -4. 2.e-13# T04 27. 7. 9. -3. -21. 2.e-1

6、3# T05 7. -8. -2. -17. 29. 2.e-13# T06 -30. -0. -5. 12. -8. 2.e-13summary(pca)# Importance of components:# PC1 PC2 PC3 PC4 PC5 PC6# Standard deviation 21.9980 21.7992 18.5932 16.67518 16.1346 15.16897# Proportion of Variance 0.1614 0.1585 0.1153 0.09272 0.0868 0.07672# Cumulative Proportion 0.1614 0

7、.3198 0.4351 0.52780 0.6146 0.69133# PC7 PC8 PC9 PC10 PC11# Standard deviation 14.48695 14.01978 13.4814 13.09112 12.8896# Proportion of Variance 0.06998 0.06554 0.0606 0.05714 0.0554# Cumulative Proportion 0.76131 0.82685 0.8875 0.94460 1.0000# PC12# Standard deviation 2.859e-13# Proportion of Vari

8、ance 0.000e+00# Cumulative Proportion 1.000e+00 上述數(shù)據(jù)中，pca$x就是后面我們畫pca圖要用的數(shù)據(jù)。而在summary(pca)中我們看到的Proportion of Variance就是各個主成分的方差占所有方差的比值，即對應(yīng)的貢獻(xiàn)率。而Cumulative Proportion則對應(yīng)的百分比累積值。從上述結(jié)果看這組數(shù)據(jù)pca結(jié)果并不是很好，所以應(yīng)該肯定會有一些分組的結(jié)果不太好。不過我們今天主要是展示結(jié)果，就不在意這些細(xì)節(jié)了。 做完上述的計算，下面就進(jìn)入我們的結(jié)果展示階段。 首先用基本畫圖函數(shù)展示。plot(pca$x,1:2 )group

9、 <- factor(c(rep("A1", 3), rep("A2", 3),rep("B1", 3),rep("B2", 3)#這里我們添加分組信息colour_group<-rainbow(length(unique(group)#利用rainbow函數(shù)選擇顏色colour<-colour_groupas.numeric(factor(group)#創(chuàng)建顏色向量colour# 1 "#FF0000FF" "#FF0000FF" "#FF000

10、0FF" "#80FF00FF" "#80FF00FF"# 6 "#80FF00FF" "#00FFFFFF" "#00FFFFFF" "#00FFFFFF" "#8000FFFF"# 11 "#8000FFFF" "#8000FFFF"plot(pca$x,1:2,col = colour ,pch = c(21,22,23,24)group)#在plot函數(shù)中我們把分組信息和顏色方案添加進(jìn)去legend

11、("topleft", legend = levels(group),col = colour_group, pch = c(21,22,23,24)#添加legendtitle("test") 這是我們用基本函數(shù)對pca分析結(jié)果的展示。除此外我們也可以利用ggplot2包進(jìn)行相同的圖片繪制。示例如下：library(ggplot2)group2<-data.frame(group)pca_reuslt<-as.data.frame(pca$x)pca_reuslt<-cbind(pca_reuslt,group2)p<-ggpl

12、ot(pca_reuslt)+geom_point(aes(x=pca_reuslt,1,y=pca_reuslt,2,color=pca_reuslt$group,shape = pca_reuslt$group),size=5)p<-p+theme(legend.title =element_blank()+labs(x="PCA1",y="PCA2")p 好了，上面那些基本的結(jié)果展示我們已經(jīng)結(jié)束了。下面我們開始把這個圖的檔次再提高一點(diǎn)。比如，我們畫了二維的，現(xiàn)在我們畫個三維的PCA結(jié)果吧。library(scatterplot3d)par(m

13、ar=c(5.1, 4.1, 4.1, 8.1), xpd=TRUE)scatterplot3d(pca_reuslt,1:3,pch=20 ,color=colour,angle=45,main = "test_3D",cex.symbols = 2,mar=c(5.1, 4.1, 4.1, 8.1)legend("right",legend = group,col = colour,pch = 20,bg="white",xpd=TRUE,inset= -0.5)#設(shè)置位置為right后，可以用inset來移到legend位置。 除此之外，我們可以考慮把相同的組進(jìn)行一個圈定，方便我們更好的觀察結(jié)果。library(ggfortify)#使用這個包時可能要注意R的版本，我剛開始用較老的版本就用不了這個包。data2<-cbind(t(data,-1),group2)a