# Makes and plots GNMDS for plant and fungal communities, # as well as procrustes comparison between the two. #Set working directory getwd() setwd("WD") list.files() # Load required libraries. library(MASS) library(vegan) # Load the vegetation table, with weighted values for each species. veg<-read.table("vegetationdata", header=TRUE) # Transpose the matrix: veg_t<-t(veg) # Make a Bary-Curtis distance matrix: dist.m.veg<-vegdist(veg_t, method="bray") # define an empty object called mds.v: mds.v<-NULL # making 100 "mds"s from initial starting configurations, allocating them into the mds object. Makes the three first axis, 1000 iterations. for(i in 1:100) {mds.v[[i]]<-isoMDS(dist.m.veg, initMDS(dist.m.veg, k=2), k=2, maxit=1000, tol=1e-7)} # distance for each ordination: mds.stress.v<-unlist(lapply(mds.v,function(v){v[[2]]})) # Order the stress values from lowest to highest order(mds.stress.v) # Choose the best stress values. mds.stress.v[23] #Scaling of axis to half change units and varimax rotation. mds.best.veg<-postMDS(mds.v[[23]], dist.m.veg) gnmds1.veg<-mds.best.veg$points[,1] gnmds2.veg<-mds.best.veg$points[,2] # Load table with factors enf<-read.table("variables", header=TRUE) names(enf) # Make mean values for each plot (we only have 10 entries for vegetation data. mean.P<-c(mean(enf$P[1:5]), mean(enf$P[6:10]), mean(enf$P[11:15]), mean(enf$P[16:20]), mean(enf$P[21:25]), mean(enf$P[26:30]), mean(enf$P[31:35]), mean(enf$P[36:40]), mean(enf$P[41:45]), mean(enf$P[46:50])) mean.N<-c(mean(enf$N[1:5]), mean(enf$N[6:10]), mean(enf$N[11:15]), mean(enf$N[16:20]), mean(enf$N[21:25]), mean(enf$N[26:30]), mean(enf$N[31:35]), mean(enf$N[36:40]), mean(enf$N[41:45]), mean(enf$N[46:50])) mean.C<-c(mean(enf$C[1:5], na.rm=TRUE), mean(enf$C[6:10]), mean(enf$C[11:15]), mean(enf$C[16:20]), mean(enf$C[21:25]), mean(enf$C[26:30]), mean(enf$C[31:35]), mean(enf$C[36:40]), mean(enf$C[41:45]), mean(enf$C[46:50])) mean.Wr<-c(mean(enf$Wr[1:5], na.rm=TRUE), mean(enf$Wr[6:10]), mean(enf$Wr[11:15]), mean(enf$Wr[16:20]), mean(enf$Wr[21:25]), mean(enf$Wr[26:30]), mean(enf$Wr[31:35]), mean(enf$Wr[36:40]), mean(enf$Wr[41:45]), mean(enf$Wr[46:50], na.rm=TRUE)) mean.RVl<-c(mean(enf$RVl[1:5], na.rm=TRUE), mean(enf$RVl[6:10]), mean(enf$RVl[11:15]), mean(enf$RVl[16:20]), mean(enf$RVl[21:25]), mean(enf$RVl[26:30]), mean(enf$RVl[31:35]), mean(enf$RVl[36:40]), mean(enf$RVl[41:45]), mean(enf$RVl[46:50])) mean.RHl<-c(mean(enf$RHl[1:5], na.rm=TRUE), mean(enf$RHl[6:10]), mean(enf$RHl[11:15]), mean(enf$RHl[16:20]), mean(enf$RHl[21:25]), mean(enf$RHl[26:30]), mean(enf$RHl[31:35]), mean(enf$RHl[36:40]), mean(enf$RHl[41:45]), mean(enf$RHl[46:50])) mean.Rt<-c(mean(enf$Rt[1:5]), mean(enf$Rt[6:10]), mean(enf$Rt[11:15]), mean(enf$Rt[16:20]), mean(enf$Rt[21:25]), mean(enf$Rt[26:30]), mean(enf$Rt[31:35]), mean(enf$Rt[36:40]), mean(enf$Rt[41:45]), mean(enf$Rt[46:50])) mean.V<-c(rep("R",5 ), rep("S", 5)) mean.V=as.factor(mean.V) # Make a matrix of the mean values. enf.10<-cbind(mean.C, mean.N, mean.P, mean.RHl, mean.RVl, mean.Wr) enf.10<-as.data.frame(enf.10) enf.10<-cbind(enf.10, mean.V) enf.10.<-apply(enf.10, 2, rev) enf.10.<-as.data.frame(enf.10.) # Make sure the vectors ar numeric. enf.10.$mean.C<-as.numeric(enf.10.$mean.C) enf.10.$mean.N<-as.numeric(enf.10.$mean.N) enf.10.$mean.P<-as.numeric(enf.10.$mean.P) enf.10.$mean.RHl<-as.numeric(enf.10.$mean.RHl) enf.10.$mean.RVl<-as.numeric(enf.10.$mean.RVl) enf.10.$mean.Wr<-as.numeric(enf.10.$mean.Wr) # Envfit the matrices. gnmds.xy.veg <- envfit(cbind(mds.best.veg$points.1, mds.best.veg$points.2), enf.10.[1:3], perm=999) # Plot the vegetation GNMDS. # Define the size of the plot window. quartz(width=6.771654, height=3.38582677) # Define the parameters of the plot. par(mai=c(0.6,0.75,0.1,0.15), lwd=0.5, mar=c(3, 3, 1, 1), mgp=c(2,1,0), cex=0.7, mfrow=c(1,2)) # Plot an empty plot. plot(mds.best.veg$points.1, mds.best.veg$points.2, xlab="GNMDS1", ylab="GNMDS2", type="n", cex.axis=0.5, cex.lab=0.7) # Plot the GNMDS points for each sapmpled plot. points(mds.best.veg$points.1[1],mds.best.veg$points.2[1], pch=21, bg="#74C476", col="#74C476", cex=0.7) points(mds.best.veg$points.1[2],mds.best.veg$points.2[2], pch=22, bg="#74C476", col="#74C476", cex=0.7) points(mds.best.veg$points.1[3],mds.best.veg$points.2[3], pch=23, bg="#74C476", col="#74C476", cex=0.7) points(mds.best.veg$points.1[4],mds.best.veg$points.2[4], pch=24, bg="#74C476", col="#74C476", cex=0.7) points(mds.best.veg$points.1[5],mds.best.veg$points.2[5], pch=25, bg="#74C476", col="#74C476", cex=0.7) points(mds.best.veg$points.1[6],mds.best.veg$points.2[6], pch=21, col="#005A32", cex=0.7, lwd=1.2) points(mds.best.veg$points.1[7],mds.best.veg$points.2[7], pch=22, col="#005A32", cex=0.7, lwd=1.2) points(mds.best.veg$points.1[8],mds.best.veg$points.2[8], pch=23, col="#005A32", cex=0.7, lwd=1.2) points(mds.best.veg$points.1[9],mds.best.veg$points.2[9], pch=24, col="#005A32", cex=0.7, lwd=1.2) points(mds.best.veg$points.1[10],mds.best.veg$points.2[10], pch=25, col="#005A32", cex=0.7, lwd=1.2) # Add legend legend("topleft", "a)", bty="n", cex=0.7) # Add environmental parameters plot(gnmds.xy.veg, arrow.mul=0.3, add=TRUE, col="black", cex=0.6) ####### GNMDS FOR FUNGAL OTUS ############### #Load the FINAL fungal otu table OTU<-read.table("OTUtable", header=TRUE) #Make a presence/absence matrix. Replaces values greater than zero with one. PA<-replace(OTU, OTU>0, 1) dim(PA) names(PA) #Transpose your matrix: PA_t<-t(PA) dim(PA_t) # Make a distance matrix, using bray-curtis dist.m<-vegdist(PA_t, method="bray") # define an empty object called mds: mds.otu<-NULL # making 100 "mds"s from initial starting configurations, allocating them into the mds object. Makes the three first axis, 1000 iterations. for(i in 1:100) {mds.otu[[i]]<-isoMDS(dist.m, initMDS(dist.m, k=2), k=2, maxit=1000, tol=1e-7)} # distance for each ordination: mds.stress<-unlist(lapply(mds.otu,function(v){v[[2]]})) # Find the mds with lowest stress value. order(mds.stress) mds.stress[54] mds.best.f<-postMDS(mds.otu[[54]], dist.m) #Make GMNDS axis 1 and 2 mds.best.f gnmds1.<-mds.best.f$points.1 gnmds2.<-mds.best.f$points.2 # envfit the data gnmds.xy <- envfit(cbind(mds.best.f$points.1, mds.best.f$points.2), enf[2:9], perm=999, na.rm=TRUE) #Plot your data. This will add the plot to the exsisting plot. #par(mai=c(0.6,0.75,0.1,0.15), lwd=0.5, mar=c(3, 3, 1, 1), mgp=c(2,1,0), cex=0.7) # plot the points plot(gnmds1., gnmds2., xlab="GNMDS1", ylab="GNMDS2", type="n", cex.axis=0.5, cex.lab=0.7, lwd=1.2) points(mds.best.f$points.1[1:5],mds.best.f$points.2[1:5], pch=21, col="#005A32", cex=0.7, lwd=1.2) points(mds.best.f$points.1[6:10],mds.best.f$points.2[6:10], pch=22, col="#005A32", cex=0.7, lwd=1.2) points(mds.best.f$points.1[11:15],mds.best.f$points.2[11:15], pch=23, col="#005A32", cex=0.7, lwd=1.2) points(mds.best.f$points.1[16:20],mds.best.f$points.2[16:20], pch=24, col="#005A32", cex=0.7, lwd=1.2) points(mds.best.f$points.1[21:25],mds.best.f$points.2[21:25], pch=25, col="#005A32", cex=0.7, lwd=1.2) points(mds.best.f$points.1[26:30],mds.best.f$points.2[26:30], pch=21, bg="#74C476", col="#74C476", cex=0.7) points(mds.best.f$points.1[31:35],mds.best.f$points.2[31:35], pch=22, bg="#74C476", col="#74C476", cex=0.7) points(mds.best.f$points.1[36:40],mds.best.f$points.2[36:40], pch=23, bg="#74C476", col="#74C476", cex=0.7) points(mds.best.f$points.1[41:45],mds.best.f$points.2[41:45], pch=24, bg="#74C476", col="#74C476", cex=0.7) points(mds.best.f$points.1[46:50],mds.best.f$points.2[46:50], pch=25, bg="#74C476", col="#74C476", cex=0.7) # Add legends. legend("topleft", "b)", cex=0.7, bty="n") legend("topright", c("Ridge1", "Ridge2", "Ridge3", "Ridge4", "Ridge5", "Snowbed1","Snowbed2", "Snowbed3", "Snowbed4","Snowbed5"), col=c("#74C476","#74C476","#74C476","#74C476","#74C476", "#005A32", "#005A32", "#005A32", "#005A32", "#005A32"), pt.bg=c("#74C476","#74C476","#74C476","#74C476","#74C476", "white", "white", "white", "white", "white"), pch=c(21,22,23,24,25,21,22,23,24,25), cex=0.5, bty="n") # Add environmental data. plot(gnmds.xy, arrow.mul=0.6, add=TRUE, col="black", cex=0.6) mds.best.veg.r<-mds.best.veg[1:5, ] mds.best.veg.s<-mds.best.veg[6:10,] mds.best.veg.s mds.veg<-rbind(mds.best.veg.s, mds.best.veg.r) mds.veg gnmds1.m<-c(mean(gnmds1.[1:5]), mean(gnmds1.[6:10]), mean(gnmds1.[11:15]),mean(gnmds1.[16:20]), mean(gnmds1.[21:25]), mean(gnmds1.[26:30]), mean(gnmds1.[31:35]), mean(gnmds1.[36:40]), mean(gnmds1.[41:45]), mean(gnmds1.[46:50])) gnmds1.m gnmds2.m<-c(mean(gnmds2.[1:5]), mean(gnmds2.[6:10]), mean(gnmds2.[11:15]),mean(gnmds2.[16:20]), mean(gnmds2.[21:25]), mean(gnmds2.[26:30]), mean(gnmds2.[31:35]), mean(gnmds2.[36:40]), mean(gnmds2.[41:45]), mean(gnmds2.[46:50])) gnmds2.m scores.mean<-cbind(gnmds1.m, gnmds2.m) #Procrustes comparison #procrustes comparisons.md ?procrustes scores.mean hei<-procrustes(mds.veg[,1:2], scores.mean, permutations=999) protest(mds.veg[,1:2], scores.mean, permutations=999) quartz(width=3.1496063, height=3.1496063) par(mai=c(0.6,0.75,0.1,0.15), lwd=0.5, mar=c(3, 3, 1, 1), mgp=c(2,1,0), cex=0.7) plot(procrustes(mds.veg, scores.mean, permutations=999), cex.axis=0.5, cex.lab=1, main="") points(mds.veg, col="darkgray", pch=16) points(hei$Yrot, col="gray12", pch=15) legend("topleft", pch=c(16, 15), c("Fungal communities GNMDS", "Vegetation GNMDS"), col=c("darkgray", "gray12"), bty="n") protest(mds.veg, scores.mean, permutations=999) #To visualize the comparisons. plot(procrustes(mds.best, mds.secbest, permutations=999)) cor.test(mds.veg[,2],scores.mean[,1],method="k") cor.test(mds.veg[,2],scores.mean[,2],method="k")