R語言做主坐標(biāo)舉例分析

使用R語言vegan包里的varespec數(shù)據(jù)集

首先加載vegan這個包

library(vegan)

然后通過data()函數(shù)來獲得這個數(shù)據(jù)

data("varespec")

這個數(shù)據(jù)集是24行44列，每行是一個樣地，每列是樣地里的植物種類。數(shù)值具體代表什么我還沒有看明白，現(xiàn)在的理解就是度量這種植物在這個樣地豐富度的一個標(biāo)準(zhǔn)吧。

• 首先是將數(shù)據(jù)集賦值給新的變量，并以site1-24對新的數(shù)據(jù)集的行進行命名

df<-varespec
rownames(df)<-paste0("site",1:24)

• 計算距離

bray_dist<-vegdist(df,method = "bray")

• 使用ape這個包中的pcoa()函數(shù)做PCoA分析

library(ape)
df.pcoa<-pcoa(bray_dist,correction = "cailliez")

df.pcoa$vectors能夠獲得用于畫圖的數(shù)據(jù)

 

df.pcoa$values可以獲得坐標(biāo)軸上顯示的百分比

 

df.plot<-data.frame(df.pcoa$vectors)
head(df.plot)
library(ggplot2)
x_label<-round(df.pcoa$values$Rel_corr_eig[1]*100,2)
y_label<-round(df.pcoa$values$Rel_corr_eig[2]*100,2)
x_label
y_label
ggplot(data=df.plot,aes(x=Axis.1,y=Axis.2))+
  geom_point()+
  theme_bw()+
  theme(panel.grid = element_blank())+
  geom_vline(xintercept = 0,lty="dashed")+
  geom_hline(yintercept = 0,lty="dashed")+
  labs(x=paste0("PCoA1 ",x_label,"%"),
       y=paste0("PCoA2 ",y_label,"%"))

 

通過上圖我們可以看到這些樣地大體上可以分為兩組，如果自己手頭有樣地的分組數(shù)據(jù)就可以看看這個結(jié)果是不是和自己的分組數(shù)據(jù)一致。

下面人為的給他分個組，然后添加一個表示分組的橢圓

df.plot$group<-ifelse(df.plot$Axis.1<0,"AAA","BBB")
ggplot(data=df.plot,aes(x=Axis.1,y=Axis.2,
                        color=group,shape=group))+
  geom_point(size=5)+
  theme_bw()+
  theme(panel.grid = element_blank())+
  geom_vline(xintercept = 0,lty="dashed")+
  geom_hline(yintercept = 0,lty="dashed")+
  labs(x=paste0("PCoA1 ",x_label,"%"),
       y=paste0("PCoA2 ",y_label,"%"))+
  stat_ellipse(data=df.plot,
               geom = "polygon",
               aes(fill=group),
               alpha=0.3)+
  scale_fill_manual(values = c("#e31a1c","#1f78b4"))