R語言繪製精美PCoA圖
阿新 • • 發佈:2019-02-14
什麼是PCoA?principal coordinate analysis
微生物群落結構受多種因素影響,例如光照、溫度、人群性別、年齡等。
以最大限度地保留原始樣本的距離關係,使相似的樣本在圖形中的距離更為接近,相異的樣本距離更遠。
輸入OTU表格之後,執行上面程式碼,就可以出來圖形(當然結果的資料輸入是經過一定修改,自己根據需求定義樣本數目,點的形狀和顏色就可以)。
微生物群落結構受多種因素影響,例如光照、溫度、人群性別、年齡等。
要了解目的分組是否與某種因素存在聯絡,我們常常會用到PCA、PCoA等排序方法。
PCoA能夠將樣本之間的相似性距離(虛擬距離),經過投影后,在低維度空間進行歐幾里德距離展示,以最大限度地保留原始樣本的距離關係,使相似的樣本在圖形中的距離更為接近,相異的樣本距離更遠。
因此相比於PCA,PCoA以樣本距離為整體考慮,更符合生態學資料特徵,應用也更為廣泛。
雖然一般的16S或者巨集基因組等分析流程當中都會包含PCoA分析,但如果自己想要更改分組的形狀,或者挑選特定的OTU進行分析,那麼自己進行操作會高效很多。
PCoA的作圖主要分為三個步驟:
選擇特定的相似性距離並計算距離矩陣。距離的選擇可以有Bray-curits、Unifrac等,不同的距離有不同的作用和意義(具體可以參考 微生物β多樣性常用計算方法比較)。相似性距離可以利用R的GUniFrac和vegan等包計算,也可以利用QIIME計算。
進行PCoA分析,也就是利用表徵分析選擇最能表示樣本距離的座標軸。這個可以利用R的ape包的pcoa()命令完成。
PCoA圖形展示。圖形可以用ordiplot()命令展示,但如果需要比較美觀的圖形,建議用ggplot來畫。
下面我們以R為基礎,展示如何根據Unweighted Unifrac距離來畫PCoA圖:
otu_table
tree檔案包括以下值(value)
edge.length
tip.label
Nnode
edge
###匯入需要的R包 library(GUniFrac) #用於計算Unifrac距離 library(ape) # 用於pcoa分析 library(ggplot2) #用於畫圖 ##讀檔案 Otu_tab <- read.table("otu_table",row.names=1,header=T,sep="\t",check.names=FALSE) Tree <- read.tree("Muscle.align.tree.nhx") #輸入OTU表格 Otu_tab <- as.data.frame(t(Otu_tab)) Otu_tab_rff <- Rarefy(Otu_tab)$otu.tab.rff unifracs <- GUniFrac(Otu_tab_rff,Tree,alpha=c(0, 0.5, 1)) du <- unifracs$unifracs[, , "d_UW"] # 計算Unweighted UniFrac距離 Group <- c('A','B','C') #按照目的輸入樣本 shape <- c("A" =16,"B" =17,"C" =16) #定義點形狀 color <- c("A" ='#CCFF33',"B" ='#CCFF33',"C" ='#CCFF33') #定義點顏色 PCOA <- pcoa(du, correction="none", rn=NULL) #利用PCOA()指令做pcoa分析 result <-PCOA$values[,"Relative_eig"] pro1 = as.numeric(sprintf("%.3f",result[1]))*100 pro2 = as.numeric(sprintf("%.3f",result[2]))*100 x = PCOA$vectors sample_names = rownames(x) pc = as.data.frame(PCOA$vectors) pc$names = sample_names legend_title = "" group = Group pc$group = group xlab=paste("PCOA1(",pro1,"%)",sep="") ylab=paste("PCOA2(",pro2,"%)",sep="") pca=ggplot(pc,aes(Axis.1,Axis.2)) + #用ggplot作圖 geom_point(size=3,aes(color=group,shape=group)) + # geom_text(aes(label=names),size=4,vjust=-1) + labs(x=xlab,y=ylab,title="PCOA",color=legend_title,shape=legend_title) + geom_hline(yintercept=0,linetype=4,color="grey") + geom_vline(xintercept=0,linetype=4,color="grey") + scale_shape_manual(values=shape) + scale_color_manual(values=color) + theme_bw()
輸入OTU表格之後,執行上面程式碼,就可以出來圖形(當然結果的資料輸入是經過一定修改,自己根據需求定義樣本數目,點的形狀和顏色就可以)。