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