R語言實現PCOA分析

  • 2019 年 12 月 19 日
  • 筆記

大家對主成分分析(principal components analysis, PCA) 都很熟悉，但是今天我們來介紹下主坐標分析(principal coordinate analysis, PCoA)。那麼這兩個差了個o字母具體有什麼區別？首先PCA是常用的降維演算法；利用線性變換，將數據變換到一個新的坐標系統中;然後再利用降維的思想，使得任何數據投影的第一大方差在第一個坐標(稱為第一主成分)上，第二大方差在第二個坐標(第二主成分)上。這種降維的思想首先減少數據集的維數，同時還保持數據集的對方差貢獻最大的特徵，最終使數據直觀呈現在二維坐標系。PCoA主要是探索數據相似度或者相異度可視化方法。可呈現研究數據相似性或差異性的可視化坐標，是一種非約束性的數據降維分析方法，可用來研究樣本群落組成的相似性或相異性。其實通俗的講，PCA主要是基於原始數據矩陣的降維；PCoA主要是基於樣本的原始數據計算出來的距離矩陣的降維。如果樣本數目比較多，而物種數目比較少，那肯定首選PCA;如果樣本數目比較少，而物種數目比較多，那肯定首選PCoA。

接下來我們看下在R中如何去實現，首先安裝ape包和vegan包，聯合使用才能達到最終的目的。包的安裝我們就不贅述了，其在CRAN平台，直接install.packages()。

首先是數據的導入，我們利用vegan自帶的數據dune。具體的數據集的構成大家可以直接在包的資訊中去看。接下來我們首先基於dune數據構造距離矩陣，需要用到的函數vegdist。

vegdist(x, method="bray",binary=FALSE, diag=FALSE, upper=FALSE,

na.rm = FALSE, ...)

其中主要的參數：

Method 支援了目前大部分的距離計算函數。"manhattan", "euclidean", "canberra","clark", "bray", "kulczynski","jaccard", "gower", "altGower","morisita", "horn", "mountford", "raup","binomial", "chao", "cao" or"mahalanobis"。

Diag 是否顯示對角線的值。

Upper 是否顯示對角線以上的值

library(vegan)  data(dune)  data(dune.env)  dune.dist <- vegdist(dune)#構造距離矩陣。

接下來就是利用ape中的pcoa函數獲取PCOA分析結果。當然也可以應用我們R自帶的函數cmdscale。

pcoa(D, correction="none",rn=NULL)

其中主要參數：

D 不用多說就是距離矩陣

Correction 主要指需不需進行校正。

Rn 主要是指的每個樣本的名稱標籤，可以進行自定義

接下來我們看下實例：

res=pcoa(dune.dist,correction="none")    head(res$value)

其中主要的值是特徵值的一些相關的轉換的值。

head(res$vectors)

其中主要是和PCA中主成分類似的柱坐標的值，進行了排序展示，一般選擇前兩個繪製二維可視散點圖。

biplot(res)#可視化PCOA 的結果

為了進一步完善我們的可靠性，我們還可以利用vegan中的ANOSIM相似性分析是一種非參數檢驗，用來檢驗組間（兩組或多組）差異是否顯著大於組內差異，從而判斷分組是否有意義。

anosim(x, grouping, permutations = 999, distance = "bray", strata = NULL, parallel = getOption("mc.cores"))
其中的參數我們不再贅述，直接看下實例：
dune.ano <- with(dune.env, anosim(dune.dist, Management))  summary(dune.ano)
plot(dune.ano)
至此，我們的PCOA的分析過程可以實現，那麼如何優化我們輸出的可視化影像，我們需要用到ggplot2這個包可以對我們的值進行更加友好的可視化。直接看實例：
  library(ggplot2)  P= dune.ano$signif  R=round( dune.ano$statistic,3)  data=res$vectors[,1:2]  data=data.frame(data)  colnames(data)=c('x','y')  eig=as.numeric(res$value[,1])    type=dune.env$Management    data=cbind(data,type)    p =ggplot(data, aes(x=x, y=y, color=type))+    geom_point(alpha=.7, size=2) +    labs(x=paste("PCoA 1 (",format(100* eig[1] / sum(eig), digits=4), "%)", sep=""),    y=paste("PCoA 2 (", format(100*eig[2] / sum(eig), digits=4), "%)", sep=""),    title="PCoA")    p +theme_classic()+stat_ellipse(type ="t", linetype = 2)+    annotate("text",x=0.004,y=-0.15,parse=TRUE,size=4,label=paste('R:',R),family="serif",fontface="italic",colour="darkred")+    annotate("text",x=0,y=-0.17,parse=TRUE,size=4,label=paste('p:',P),family="serif",fontface="italic",colour="darkred")

		  				

		  				

													

							

																
Tags:  

																

								

							

						

						
					

					
									

													
													
											


			



	

		
VirMach 便宜 VPS