微生物生态数据分析——冗余分析

微生物生态数据分析——冗余分析

sa=read.table("NRRDA.csv",header=T,row.names=1,sep=",")
env=read.table("NRenv.csv",header=T,sep=",",row.names=1)
name <- read.table("name.csv", header=F,sep=",",colClasses=c("character","character"))#用作产生图例，对样本进行分类
library(vegan)
decorana(sa)#DCA函数，用于决定选择进行RDA或者CCA
#如果DCA排序前4个轴中最大值超过4，选择单峰模型排序更合适。如果是小于3，则选择线性模型更好(Lepx & Smilauer 2003)。如果介于3-4之间，单峰模型和线性模型都可行
sam=rda(sa,env,center=T,scale=T)
plot(sam,scaling=3)#plot(gts.rda,display=c("sp","bp"),scaling=3)
# display=c("sp","bp") 表示显示物种与环境因子。如果显示样方与环境因子，可以表示为display=c("si","bp")，如果物种、样方和环境因子同时显示，可以设定display=c("sp","si","bp")。
library(ggplot2)
new=sam$CCA
samples<-data.frame(sample=row.names(new$u),RDA1=new$u[,1],RDA2=new$u[,2]) 
name=name$V2#分类文件
species<-data.frame(spece=row.names(new$v),RDA1=new$v[,1],RDA2=new$v[,2])  
envi<-data.frame(en=row.names(new$biplot),RDA1=new$biplot[,1],RDA2=new$biplot[,2]) 
rda1 =round(sam$CCA$eig[1]/sum(sam$CCA$eig)*100,2) #第一轴标签,显示解释度
rda2 =round(sam$CCA$eig[2]/sum(sam$CCA$eig)*100,2) #第二轴标签,显示解释度
line_x = c(0,envi[1,2],0,envi[2,2],0,envi[3,2])  #行数与ling_g数量一致，envi[line_g,2]
line_x  
line_y = c(0,envi[1,3],0,envi[2,3],0,envi[3,3])  
line_y  
line_g = c("grazing","grazing","Tot bio","Tot bio","ND","ND")  
line_g  
line_data = data.frame(x=line_x,y=line_y,group=line_g)  
line_data
#开始重绘RDA图
#填充样本数据，分别以RDA1,RDA2为X,Y轴，不同样本以颜色区分
ggplot(data=samples,aes(RDA1,RDA2)) + geom_point(aes(color=sample),size=2)#生成图例的数据使用的是sample，不是samples！，不然会提示错误: Aesthetics must be either length 1 or the same as the data (6): colour, x, y
#geom_point(aes(color=name),size=2
#填充微生物物种数据，不同物种以图形区分,seq增加形状数量
 + geom_point(data=species,aes(shape=spece),size=2) + scale_shape_manual(values=seq(0,15))+
#填充环境因子数据，直接展示
geom_text(data=envi,aes(label=en),color="blue") + 
#添加0刻度纵横线
geom_hline(yintercept=0) + geom_vline(xintercept=0)+ 
#添加原点指向环境因子的带箭头的直线，size可以调节直线粗细
geom_line(data=line_data,aes(x=x,y=y,group=group),color="green") + 
geom_segment(data=line_data,aes(x=x,y=y,xend = line_data[,1], yend = line_data[,2],group=group),color="black",size=1.5,arrow=arrow(angle=35, length=unit(0.3, "cm"))) +
#添加横纵坐标轴标签
labs(title="RDA Plot",x=paste("RDA1 ",rda1," %"),y=paste("RDA2 ",rda2," %")) +
#标题字体格式设置
#theme(text=element_text(family="serif"))+
#去除背景颜色及多余网格线
theme_bw() + theme(panel.grid=element_blank()) 
#大功告成，保存为矢量图等等
ggsave("NRRDA2.PDF")      
#RDA更详细的分析，
summary(sam)
#检验环境因子相关显著性（Monte Carlo permutation test）
permutest(sam,permu=999) # permu=999是表示置换循环的次数
ef=envfit(sam,env,permu=999)#每个环境因子显著性检验
ef

看到有很多同学问导入的数据格式问题，我重新写了一份详细的排序分析文章，发在了公众号“科学研究进行时”上，文章名称：R绘图-RDA排序分析，里面有数据格式的截图，大家可以扫描二维码关注查看。