Microbial ecological data analysis - redundancy analysis

Microbial ecological data analysis —— Redundancy analysis

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"))# Used to generate legend , Classify the samples 
library(vegan)
decorana(sa)#DCA function , Used to decide whether to RDA perhaps CCA
# If DCA Before ordering 4 The maximum of axes exceeds 4, It is more appropriate to select unimodal model for sorting . If it's less than 3, The linear model is better (Lepx & Smilauer 2003). If it's between 3-4 Between , Both unimodal model and linear model are feasible 
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")  Show species and environmental factors . If quadrat and environmental factors are displayed , It can be expressed as display=c("si","bp"), If the species 、 Quadrat and environmental factors are displayed at the same time , You can set 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# Classified documents 
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) # First axis label , Show the degree of explanation 
rda2 =round(sam$CCA$eig[2]/sum(sam$CCA$eig)*100,2) # Second axis label , Show the degree of explanation 
line_x = c(0,envi[1,2],0,envi[2,2],0,envi[3,2])  # Lines and ling_g The number is the same ,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
# Start redrawing RDA chart 
# Fill in the sample data , Respectively by RDA1,RDA2 by X,Y Axis , Different samples are distinguished by color 
ggplot(data=samples,aes(RDA1,RDA2)) + geom_point(aes(color=sample),size=2)# The data used to generate the legend is sample, No samples！, Otherwise, it will prompt an error : Aesthetics must be either length 1 or the same as the data (6): colour, x, y
#geom_point(aes(color=name),size=2
# Populate microbial species data , Different species are distinguished graphically ,seq Increase the number of shapes 
 + geom_point(data=species,aes(shape=spece),size=2) + scale_shape_manual(values=seq(0,15))+
# Fill in environmental factor data , Show directly 
geom_text(data=envi,aes(label=en),color="blue") + 
# add to 0 Scale vertical and horizontal lines 
geom_hline(yintercept=0) + geom_vline(xintercept=0)+ 
# Add a straight line with an arrow whose origin points to the environment factor ,size The thickness of the straight line can be adjusted 
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"))) +
# Add horizontal and vertical axis labels 
labs(title="RDA Plot",x=paste("RDA1 ",rda1," %"),y=paste("RDA2 ",rda2," %")) +
# Title font formatting 
#theme(text=element_text(family="serif"))+
# Remove background colors and excess gridlines 
theme_bw() + theme(panel.grid=element_blank()) 
# Be accomplished , Save as vector graph, etc 
ggsave("NRRDA2.PDF")      
#RDA More detailed analysis ,
summary(sam)
# Test the significance of environmental factors （Monte Carlo permutation test）
permutest(sam,permu=999) # permu=999 Is the number of permutation cycles 
ef=envfit(sam,env,permu=999)# Significance test of each environmental factor 
ef

I see that many students ask questions about the imported data format , I wrote a detailed ranking analysis article again , It was posted on the official account “ When scientific research is in progress ” On , Article title ：R mapping -RDA Sorting analysis , There are screenshots of data format , You can scan the QR code and pay attention to view .