R語言作圖之PCA作圖和散點圖

PCA分析和散點圖今天主要跟大家演示一下簡單的PCA分析，并且以散點圖的形式將結(jié)果展示出來。首先在進(jìn)行PCA分析之前，先跟大家稍微討論下什么是PCA分析。PCA分析又叫主成分分析，其實從字面上來理解我們可以發(fā)現(xiàn)它其實是和樣品分組相關(guān)的。舉個簡單的例子，我們觀察了某種植物的株高、葉片大小、果實大小等等多種性狀，并記錄每種性狀對應(yīng)的數(shù)值。這時候我們想看看根據(jù)這些性狀信息看看我們觀察的樣本是否明顯的分組現(xiàn)象。每一種性狀相當(dāng)于一個維度。利用PCA分析可以將結(jié)果投影到一個低維的向量空間(具體計算就不詳述了)。類似的比如我們多個樣本的表達(dá)譜數(shù)據(jù)，每個基因在各個樣品的表達(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_REFT01T02T03T04T05T06##11007_s_at10.19858611.80567610.86795311.76366012.07223212.108312##21053_at9.5940748.7131089.2470969.4332659.0923299.005518##3117_at8.5817638.6036808.8044258.6617008.6349798.606976##4121_at12.02231512.65532912.62733412.79139012.96176112.885307##51255_g_at7.2285697.2146007.2371317.2934177.2767997.268233##61294_at8.8284879.3802779.2979898.8589858.9957729.126825##T07T08T09T10T11T12##110.64686810.85274410.67589811.13766310.79673711.102408##29.0876819.0272088.9652838.9583099.2750108.940965##38.6258388.5772448.6467518.6258438.6251648.522129##413.40204413.24012613.08888313.23409913.38290313.472223

#5 7.197440 7.262662 7.289796 7.232249 7.202364 7.306229#6 9.002385 9.003561 9.006278 9.006721 9.018183 9.164313上述數(shù)據(jù)為從GEO數(shù)據(jù)庫隨意找的基因表達(dá)。其中第一列為基因探針號，后續(xù)幾列則為T01到T12的12個樣品對應(yīng)的表達(dá)量數(shù)據(jù)，每三個樣品為一組。因為數(shù)據(jù)是拼湊的，所以這里不關(guān)注探針具體信息了。準(zhǔn)備好數(shù)據(jù)之后我們就開始進(jìn)行PCA計算了。其實代碼非常簡單。pca<-prcomp(t(data[,-1]),scale=T)head(pca$x)##PC1PC2PC3PC4PC5PC6##T01-43.457435-44.9500318.3055713.210563-7.428048114.818150##T0242.067255-19.142248-25.57404121.120294-5.793099014.702922##T03-2.123455-21.512488-11.19247417.58300615.2149034-34.730308##T048.166077-4.77481422.837578-11.3641288.4021038-6.921738##T0518.214073-5.83680718.522768-10.941626-0.6183613-5.548845##T0627.219529-5.51932826.649872-11.054961-4.14804135.097715##PC7PC8PC9PC10PC11PC12##T01-1.9663429.2181269-1.520882-1.0608353.0484982.731227e-13##T025.8321978.97930189.3861871.6687611.7054742.674666e-13##T03-5.168168-9.7483411-11.5703202.618203-4.2214562.738955e-13##T0427.7829867.58290079.726761-3.391763-21.9004852.730871e-13##T057.039535-8.9173716-2.239005-17.51443329.7009062.736544e-13##T06-30.026232-0.8253129-5.20703712.349414-8.9006762.681674e-13summary(pca)#Importanceofcomponents:# PC1PC2PC3PC4PC5PC6#Standarddeviation21.998021.799218.593216.6751816.134615.16897#ProportionofVariance0.16140.15850.11530.092720.08680.07672

#CumulativeProportion0.16140.31980.43510.527800.61460.69133# PC7 PC8 PC9 PC10 PC11#Standarddeviation 14.48695 14.01978 13.4814 13.09112 12.8896#Proportion ofVariance 0.06998 0.06554 0.0606 0.05714 0.0554#Cumulative Proportion 0.76131 0.82685 0.8875 0.94460 1.0000# PC12#Standarddeviation2.859e-13#ProportionofVariance0.000e+00#CumulativeProportion1.000e+00上述數(shù)據(jù)中，pca$x就是后面我們畫pca圖要用的數(shù)據(jù)。而在summary(pca)中我們看到的ProportionofVaria^各個主成分的方差占所有方差的比值，即對應(yīng)的貢獻(xiàn)率。而CumulativeProportio則對應(yīng)的百分比累積值。從上述結(jié)果看這組數(shù)據(jù)pca結(jié)果并不是很好，所以應(yīng)該肯定會有一些分組的結(jié)果不太好。不過我們今天主要是展示結(jié)果，就不在意這些細(xì)節(jié)了。做完上述的計算，下面就進(jìn)入我們的結(jié)果展示階段。首先用基本畫圖函數(shù)展示。plot(pca$x[,1:2])oOJozI-40-2020oOJozI-40-202040PC1group<-factor(c(rep("A1",3),rep("A2",3),rep("B1",3),rep("B2",3)))這里我們添加分組信息colour_group<-rainbow(length(unique(group)))#利用rainbow函數(shù)選擇顏色colour<-colour_group[as.numeric(factor(group))]#倉ll建顏色向量colour# [1] "#FF0000FF" ”#FF0000FF" ”#FF0000FF" "#80FF00FF" "#80FF00FF”# [6] "#80FF00FF" "#00FFFFFF" "#00FFFFFF" "#00FFFFFF" "#8000FFFF"# [1l] "#8000FFFF" "#8000FFFF"plot(pca$x[,1:2],col=colour,pch=c(21,22,23,24)[group])#在plot函數(shù)中我們把分組信息和顏色方案添加進(jìn)去legend("topleft",legend=levels(group),col=colour_group,pch=c(21,22,23,24))#添加legendtitle("test")test口A2△B2-40PC1這是我們用基本函數(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<-ggplot(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")20-oO-2e<o20-oO-2e<oCL知A2B1%?▲■+-40-.I I I-25 0 25PCA1好了，上面那些基本的結(jié)果展示我們已經(jīng)結(jié)束了。下面我們開始把這個圖的檔次再提高一點。比如，我們畫了二維的，現(xiàn)在我們畫個三維的PCA結(jié)果吧。library(scatterplot3d)par(mar=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位置。test3DA1A1test3DA1A1A1A2A2A2B1B1B1B2B2B2PC1除此之外，我們可以考慮把相同的組進(jìn)行一個圈定，方便我們更好的觀察結(jié)果。library(ggfortify)#使用這個包時可能要注意的版本，我剛開始用較老的版本就用不了這個包。dat