R-RDA

AOA=read.table("16AOA.csv",header=T,row.names=1,sep=",")
AOB=read.table("16AOB.csv",header=T,row.names=1,sep=",")
envi=read.table("16env.csv",header=T,sep=",",row.names=1)
name <- read.table("name.csv", header=T,sep=",",colClasses=c("character","character"))#用作產生圖例，對樣本進行分類
Aspe<-read.table("Aspe.csv",sep=",",header=T,row.names=1)
Bspe<-read.table("Bspe.csv",sep=",",header=T,row.names=1)
library(vegan)
decorana(Aspe)
decorana(Bspe)#DCA函數，用於決定選擇進行RDA或者CCA
#如果DCA排序前4個軸中最大值超過4，選擇單峰模型排序更合適。如果是小於3，則選擇線性模型更好(Lepx & Smilauer 2003)。如果介於3-4之間，單峰模型和線性模型都可行
#AOA RDA plot
A.rda=rda(Aspe,envi,center=T,scale=T)
A.rda.data <- data.frame(A.rda$CCA$u[,c(1,2)], group$V2)  #提取樣本得分
colnames(A.rda.data) <- c("RDA1","RDA2","group")
A.spe <- A.rda$CCA$v[,c(1,2)]  #物種得分
A.spe<- as.data.frame(A.spe)
env <- A.rda$CCA$biplot[,c(1,2)] #環境因子得分
env <- as.data.frame(env)
A.rda1 =round(A.rda$CCA$eig[1]/sum(A.rda$CCA$eig)*100,2) #第一軸標簽
A.rda2 =round(A.rda$CCA$eig[2]/sum(A.rda$CCA$eig)*100,2) #第二軸標簽
#繪圖對象的創建
#顏色自定義，用於對shape的線條和實心進行上色
cols=c("CK"="red","N"="grey","W"="blue","-W"="green","NW"="black","N-W"="darkgreen")
Treat=name$treat
Layer=name$layer
Aplot<- ggplot(data=A.rda.data,aes(RDA1,RDA2))
#幾何對象
#fill填充顏色，21和24分別為空心圓與空心三角形，需要用對應的顏色填充
Aplot<- Aplot+ geom_point(aes(colour=Treat,shape=Layer,fill=Treat),size=3) + 
    scale_shape_manual(values=c(21,24))+scale_colour_manual(values=cols,breaks=c("CK","N","W","-W","NW","N-W"))+theme(legend.title=element_blank())+scale_fill_manual(values=cols,breaks=c("CK","N","W","-W","NW","N-W"))+
    labs(title="AOA RDA Plot",x=paste("RDA1 ",A.rda1," %"),y=paste("RDA2 ",A.rda2," %")) +
    theme(text=element_text(family="serif"))
#去除背景以及網格線
Aplot=Aplot+theme_bw() + theme(panel.grid=element_blank()) 
Aplot=Aplot + geom_hline(yintercept=0) + geom_vline(xintercept=0)+ geom_segment(data = env,aes(x=0,y=0,xend = env[,1], yend = env[,2]), colour="black", size=1, arrow=arrow(angle=35, length=unit(0.3, "cm"))) +
    geom_text(data=env, aes(x=env[,1], y=env[,2], label=rownames(env)),
              size=3.5, colour="black",hjust = (1 - 2 * sign(env[ ,1])) / 2,
              angle = (180/pi) * atan(env[ ,2]/env[ ,1]))
Aplot
#RDA更詳細的分析
A.sum=summary(A.rda)
A.sum
#檢驗環境因子相關顯著性（Monte Carlo permutation test）
A.perm=permutest(A.rda,permu=999) # permu=999是表示置換循環的次數
A.perm
A.ef=envfit(A.rda,envi,permu=999)#每個環境因子顯著性檢驗
A.ef
#CCA1 和 CCA2 兩列所對應的值是環境因子箭頭與排序軸夾角的余弦值，表示環境因子與排序軸的相關性。r2表示環境因子對物種分布的決定系數，r2越小，表示該環境因子對物種分布影響越小。Pr表示相關性的顯著性檢驗
#第二張圖
AOB RDA plot
B.rda=rda(Bspe,envi,center=T,scale=T)
B.rda.data <- data.frame(B.rda$CCA$u[,c(1,2)], group$V2)  #提取樣本得分
colnames(B.rda.data) <- c("RDA1","RDA2","group")
B.spe <- B.rda$CCA$v[,c(1,2)]  #物種得分
B.spe<- as.data.frame(B.spe)
B.env <- B.rda$CCA$biplot[,c(1,2)] #環境因子得分
B.env <- as.data.frame(B.env)
B.rda1 =round(B.rda$CCA$eig[1]/sum(B.rda$CCA$eig)*100,2) #第一軸標簽
B.rda2 =round(B.rda$CCA$eig[2]/sum(B.rda$CCA$eig)*100,2) #第二軸標簽
Bplot<- ggplot(data=B.rda.data,aes(RDA1,RDA2))
#幾何對象
Bplot<- Bplot+ geom_point(aes(colour=Treat,shape=Layer,fill=Treat),size=3) + 
    scale_shape_manual(values=c(21,24))+scale_colour_manual(values=cols,breaks=c("CK","N","W","-W","NW","N-W"))+theme(legend.title=element_blank())+scale_fill_manual(values=cols,breaks=c("CK","N","W","-W","NW","N-W"))+
    labs(title="AOB RDA Plot",x=paste("RDA1 ",B.rda1," %"),y=paste("RDA2 ",B.rda2," %")) +
    theme(text=element_text(family="serif"))
#去除背景以及網格線
Bplot=Bplot+theme_bw() + theme(panel.grid=element_blank()) 
Bplot=Bplot + geom_hline(yintercept=0) + geom_vline(xintercept=0)+ geom_segment(data = env,aes(x=0,y=0,xend = env[,1], yend = env[,2]), colour="black", size=1, arrow=arrow(angle=35, length=unit(0.3, "cm"))) +
    geom_text(data=env, aes(x=env[,1], y=env[,2], label=rownames(env)),
              size=3.5, colour="black",hjust = (1 - 2 * sign(env[ ,1])) / 2,
              angle = (180/pi) * atan(env[ ,2]/env[ ,1]))
Bplot			  
B.sum=summary(B.rda)
B.sum
#檢驗環境因子相關顯著性（Monte Carlo permutation test）
B.perm=permutest(B.rda,permu=999) # permu=999是表示置換循環的次數
B.perm
B.ef=envfit(B.rda,envi,permu=999)#每個環境因子顯著性檢驗
B.ef

x

45

 

1

AOA=read.table("16AOA.csv",header=T,row.names=1,sep=",")

2

AOB=read.table("16AOB.csv",header=T,row.names=1,sep=",")

3

envi=read.table("16env.csv",header=T,sep=",",row.names=1)

4

name <- read.table("name.csv", header=T,sep=",",colClasses=c("character","character"))#用作產生圖例，對樣本進行分類

5

Aspe<-read.table("Aspe.csv",sep=",",header=T,row.names=1)

6

Bspe<-read.table("Bspe.csv",sep=",",header=T,row.names=1)

7

library(vegan)

8

decorana(Aspe)

9

decorana(Bspe)#DCA函數，用於決定選擇進行RDA或者CCA

10

#如果DCA排序前4個軸中最大值超過4，選擇單峰模型排序更合適。如果是小於3，則選擇線性模型更好(Lepx & Smilauer 2003)。如果介於3-4之間，單峰模型和線性模型都可行

11

#AOA RDA plot

12

A.rda=rda(Aspe,envi,center=T,scale=T)

13

A.rda.data <- data.frame(A.rda$CCA$u[,c(1,2)], group$V2)  #提取樣本得分

14

colnames(A.rda.data) <- c("RDA1","RDA2","group")

15

A.spe <- A.rda$CCA$v[,c(1,2)]  #物種得分

16

A.spe<- as.data.frame(A.spe)

17

env <- A.rda$CCA$biplot[,c(1,2)] #環境因子得分

18

env <- as.data.frame(env)

19

A.rda1 =round(A.rda$CCA$eig[1]/sum(A.rda$CCA$eig)*100,2) #第一軸標簽

20

A.rda2 =round(A.rda$CCA$eig[2]/sum(A.rda$CCA$eig)*100,2) #第二軸標簽

21

#繪圖對象的創建

22

#顏色自定義，用於對shape的線條和實心進行上色

23

cols=c("CK"="red","N"="grey","W"="blue","-W"="green","NW"="black","N-W"="darkgreen")

24

Treat=name$treat

25

Layer=name$layer

26

Aplot<- ggplot(data=A.rda.data,aes(RDA1,RDA2))

27

#幾何對象

28

#fill填充顏色，21和24分別為空心圓與空心三角形，需要用對應的顏色填充

29

Aplot<- Aplot+ geom_point(aes(colour=Treat,shape=Layer,fill=Treat),size=3) + 

30

    scale_shape_manual(values=c(21,24))+scale_colour_manual(values=cols,breaks=c("CK","N","W","-W","NW","N-W"))+theme(legend.title=element_blank())+scale_fill_manual(values=cols,breaks=c("CK","N","W","-W","NW","N-W"))+

31

    labs(title="AOA RDA Plot",x=paste("RDA1 ",A.rda1," %"),y=paste("RDA2 ",A.rda2," %")) +

32

    theme(text=element_text(family="serif"))

33

#去除背景以及網格線

34

Aplot=Aplot+theme_bw() + theme(panel.grid=element_blank()) 

35

Aplot=Aplot + geom_hline(yintercept=0) + geom_vline(xintercept=0)+ geom_segment(data = env,aes(x=0,y=0,xend = env[,1], yend = env[,2]), colour="black", size=1, arrow=arrow(angle=35, length=unit(0.3, "cm"))) +

36

    geom_text(data=env, aes(x=env[,1], y=env[,2], label=rownames(env)),

37

              size=3.5, colour="black",hjust = (1 - 2 * sign(env[ ,1])) / 2,

38

              angle = (180/pi) * atan(env[ ,2]/env[ ,1]))

39

Aplot

40

#RDA更詳細的分析

41

A.sum=summary(A.rda)

42

A.sum

43

#檢驗環境因子相關顯著性（Monte Carlo permutation test）

44

A.perm=permutest(A.rda,permu=999) # permu=999是表示置換循環的次數

45

A.perm

46

A.ef=envfit(A.rda,envi,permu=999)#每個環境因子顯著性檢驗

47

A.ef

48

#CCA1 和 CCA2 兩列所對應的值是環境因子箭頭與排序軸夾角的余弦值，表示環境因子與排序軸的相關性。r2表示環境因子對物種分布的決定系數，r2越小，表示該環境因子對物種分布影響越小。Pr表示相關性的顯著性檢驗

49

#第二張圖

50

AOB RDA plot

51

B.rda=rda(Bspe,envi,center=T,scale=T)

52

B.rda.data <- data.frame(B.rda$CCA$u[,c(1,2)], group$V2)  #提取樣本得分

53

colnames(B.rda.data) <- c("RDA1","RDA2","group")

54

B.spe <- B.rda$CCA$v[,c(1,2)]  #物種得分

55

B.spe<- as.data.frame(B.spe)

56

B.env <- B.rda$CCA$biplot[,c(1,2)] #環境因子得分

57

B.env <- as.data.frame(B.env)

58

B.rda1 =round(B.rda$CCA$eig[1]/sum(B.rda$CCA$eig)*100,2) #第一軸標簽

59

B.rda2 =round(B.rda$CCA$eig[2]/sum(B.rda$CCA$eig)*100,2) #第二軸標簽

60

Bplot<- ggplot(data=B.rda.data,aes(RDA1,RDA2))

61

#幾何對象

62

Bplot<- Bplot+ geom_point(aes(colour=Treat,shape=Layer,fill=Treat),size=3) + 

63

    scale_shape_manual(values=c(21,24))+scale_colour_manual(values=cols,breaks=c("CK","N","W","-W","NW","N-W"))+theme(legend.title=element_blank())+scale_fill_manual(values=cols,breaks=c("CK","N","W","-W","NW","N-W"))+

64

    labs(title="AOB RDA Plot",x=paste("RDA1 ",B.rda1," %"),y=paste("RDA2 ",B.rda2," %")) +

65

    theme(text=element_text(family="serif"))

66

#去除背景以及網格線

67

Bplot=Bplot+theme_bw() + theme(panel.grid=element_blank()) 

68

Bplot=Bplot + geom_hline(yintercept=0) + geom_vline(xintercept=0)+ geom_segment(data = env,aes(x=0,y=0,xend = env[,1], yend = env[,2]), colour="black", size=1, arrow=arrow(angle=35, length=unit(0.3, "cm"))) +

69

    geom_text(data=env, aes(x=env[,1], y=env[,2], label=rownames(env)),

70

              size=3.5, colour="black",hjust = (1 - 2 * sign(env[ ,1])) / 2,

71

              angle = (180/pi) * atan(env[ ,2]/env[ ,1]))

72

Bplot     

73

B.sum=summary(B.rda)

74

B.sum

75

#檢驗環境因子相關顯著性（Monte Carlo permutation test）

76

B.perm=permutest(B.rda,permu=999) # permu=999是表示置換循環的次數

77

B.perm

78

B.ef=envfit(B.rda,envi,permu=999)#每個環境因子顯著性檢驗

79

B.ef

參考文獻：

圖例： http://blog.sina.com.cn/s/blog_dc25c5920102xj9e.html

https://www.cnblogs.com/yangk/p/4761994.html

https://blog.csdn.net/tanzuozhev/article/details/51108040