##### Inquiries about this R code should be directed to ## Bucky Gates, tagates@ncsu.edu library("vegan") FaunaData<-read.csv("~/Desktop/COI Pie Charts/Goldberg database FamilyHigher Taxonomic Data.csv", header=T, row.names=1) TaxaColOnly<-(ncol(FaunaData)-1) TaxaOnlyRows<-(nrow(FaunaData)-2) ############################################## #Pie Charts################################### ############################################## pienames<-colnames(FaunaData) piecolors<-c("antiquewhite1", "aliceblue", "aquamarine2", "bisque2", "cadetblue2", "azure3", "lightpink1", "lightgoldenrod1", "mistyrose1", "slategray1") for (i in 1:(ncol(FaunaData)-1)){ datarows<-nrow(FaunaData)-2 pieFileName<-paste("COI pie charts", pienames[i], sep=" ") postscript(file= pieFileName) pieslices<-FaunaData[c(1:(nrow(FaunaData)-2)),i] pietitle<-paste("Pie Chart for", pienames[i], sep=" ") lbls<- paste(row.names(FaunaData[1:datarows]), round((100*pieslices/sum(pieslices)), 1), "%", sep=" ") pie(pieslices, labels=lbls, main=pienames[i], cex=0.8, col=piecolors) legend("topright", row.names(FaunaData[1:datarows,]), cex=0.8, fill=piecolors) dev.off() } # Pie chart for combined taxa RowTotal<-(1:TaxaOnlyRows) names(RowTotal)<-rownames(FaunaData[1: TaxaOnlyRows,]) for (w in 1: TaxaOnlyRows) { RowTotal[w]<-sum(FaunaData[w, 1: TaxaColOnly]) } postscript(file= "Total Pie Chart") lbls<- paste(row.names(FaunaData[1:length(RowTotal)]), round((100*RowTotal/sum(RowTotal)), 1), "%", sep=" ") pie(RowTotal, labels=lbls, main="Total Proportions of All Specimens", cex=0.8, col=piecolors) legend("topright", row.names(FaunaData[1:length(RowTotal),]), cex=0.8, fill=piecolors) dev.off() ################################################################ ##Summary Histograms############################################ ################################################################ SW<-(1: TaxaColOnly) names(SW)<-colnames(FaunaData[1: TaxaColOnly]) FossilDens<-(1: TaxaColOnly) names(FossilDens)<-colnames(FaunaData[1: TaxaColOnly]) SiteProp<-(1: TaxaColOnly) names(SiteProp)<-colnames(FaunaData[1: TaxaColOnly]) barVec<-(1: (TaxaColOnly*3)) for (l in 1: TaxaColOnly) { FossilDens[l]<-round(sum(FaunaData[1: TaxaOnlyRows,l]) /FaunaData[nrow(FaunaData), l], digits=3) barVec[2*l+(l-2)]<-FossilDens[l] SiteProp[l]<-round((sum(FaunaData[1: TaxaOnlyRows, l])/sum(FaunaData[1: TaxaOnlyRows, 1: TaxaColOnly])), digits=3) barVec[2*l+(l-1)]<-SiteProp[l] SW[l]<- round(diversity(FaunaData[1: TaxaOnlyRows,l], index = "shannon", MARGIN = 1, base = exp(1)), digits=3) barVec[3*l]<-SW[l] } barMatrix<-rbind(FossilDens, SiteProp, SW) barMatrix<-as.matrix(barMatrix) write.csv(barMatrix, "COI Complete Site Biodiversity Data Table") barMatrix[1,1]<-(barMatrix[1,1]/11) ## This is because the value is so large that it makes the other values extremely small on the graph...the true value is listed above the bar and in the accompanying data table barMatrix1<-cbind(barMatrix[,1:4]) barMatrix3<-cbind(barMatrix[,9:13]) barMatrix2<-cbind(barMatrix[,5:8]) barVec1<-barVec[1:12] barVec2<-barVec[13:24] barVec3<-barVec[25:39] postscript(file="COI barchart1") CompGraph<-barplot(barMatrix1, beside=T, names.arg=colnames(barMatrix1), axes=F, axisnames=T, cex.axis=1) text(CompGraph, barVec1+0.1, paste(barVec1), cex = 1, col = "red") dev.off() postscript(file="COI barchart2") CompGraph<-barplot(barMatrix2, beside=T, names.arg=colnames(barMatrix2), axes=F, axisnames=T, cex.axis=1) text(CompGraph, barVec2+0.1, paste(barVec2), cex = 1, col = "red") dev.off() postscript(file="COI barchart3") CompGraph<-barplot(barMatrix3, beside=T, names.arg=colnames(barMatrix3), axes=F, axisnames=T, cex.axis=1) text(CompGraph, barVec3+0.1, paste(barVec3), cex = 1, col = "red") dev.off() ##X-Y of Fossil abundance vs biodiversity (Shannon-Wiener) AbundVec<-(1: TaxaColOnly) names(AbundVec)<-colnames(FaunaData[1: TaxaColOnly]) for (q in 1: TaxaColOnly) { AbundVec[q]<-sum(FaunaData[1: TaxaOnlyRows, q]) } AbundSW<-lm(SW~AbundVec) AbundSW ################################################################ ##Rarefaction curves############################################ ################################################################ rareMatrix<-FaunaData[-c(nrow(FaunaData)-1, nrow(FaunaData)),-ncol(FaunaData)] rareMatrix<-t(rareMatrix) postscript(file="COI Rarefaction curves color") RCurve<-rarecurve(rareMatrix, step=1, col=(1:nrow(rareMatrix))) dev.off() ################################################################ ## Shannon-Weiner Subsampling ############################ ################################################################ library("vegan") library("fossil") FaunaData<-read.csv("~/Desktop/COI Pie Charts/Goldberg database FamilyHigher Taxonomic Data.csv", header=T, row.names=1) FaunaData<-FaunaData[,-ncol(FaunaData)] FaunaData<-FaunaData[-c((nrow(FaunaData)-1):(nrow(FaunaData))), ] #This deletes the rows about Total specimens and total sediment processed TotalRuns<-100 ## This number designates the total number of replicates run for subsampling Results.Matrix<-matrix(nrow=TotalRuns, ncol=((ncol(FaunaData)-1)*2)) #This matrix will house the actual subsampled biodiversity for the COI site as well as the difference between this number and that of the comparable site. Matrix.names<-colnames(FaunaData) Matrix.names<-Matrix.names[-1] Results.Matrix.names<-vector(length=(length(Matrix.names)*2)) for (i in 1:length(Matrix.names)){ SW.names<-paste(Matrix.names[i], "SW", sep=" ") M.H.names<-paste(Matrix.names[i], "Mor.Horn", sep=" ") Results.Matrix.names[i]<-SW.names Results.Matrix.names[i+(length(Results.Matrix.names)/2)]<-M.H.names } colnames(Results.Matrix)<-Results.Matrix.names COI.Sample.Vectors<-matrix(nrow=nrow(FaunaData), ncol=TotalRuns) #This will recordeach vector of subsampled COI species COI.Sample.List<-list() #This will record each of the Vectors for each site. for (c in 1:(ncol(FaunaData)-1)) { ## Loop goes across the columns of the Data chart starting with the second column print(paste(c, "/", ncol(FaunaData)-1, sep="")) SWa<- diversity(FaunaData[,(c+1)], index = "shannon", MARGIN = 1, base = exp(1)) # This is the SW of the OMNH site to compare COI for (r in 1:TotalRuns) { ## Loop goes down the rows of the Data chart sampData=vector(length=nrow(FaunaData)) #This creates the vector we will accumulate sample data names(sampData)=row.names(FaunaData) BinVec<-1:nrow(FaunaData) #This creates a vector of a single number for each animal taxon that we will use to taxa sample from. Even though all taxa are not present in every site, we use a filled vector here because sampling weights will dictate the proportion of taxa that are samples, meaning that those with a 0 weight will not be sampled. WeightVec<-FaunaData[,1]/sum(FaunaData[,1]) #This is placing weights on the random sampling of equal to the current proportion of fauna discovered from the COI site sample.vec<-rep(0, nrow(FaunaData)) #This is the initial sampling vector with all values set to 0. for (g in 1:sum(FaunaData[,(c+1)])){ #In short, this loop is randomly drawing a number from vector BinVec that corresponds to an animal group in our analysis. Once chosen that rgoup gets 1 added. Once completed we will have a randomly drawn biodiversity sample. q<-sample(BinVec, size=1, prob=WeightVec) sample.vec[q]<-sample.vec[q]+1 } COI.Sample.Vectors[,r]<-sample.vec #saving the smapling vector in data matrix SWcoi<- diversity(sample.vec, index = "shannon", MARGIN = 1, base = exp(1)) #This is the subsampled COI SW index Results.Matrix[r,c]<-SWcoi-SWa Results.Matrix[r,(c+(ncol(Results.Matrix)/2))]<-morisita.horn(sample.vec, FaunaData[,c]) } COI.Sample.List[[c]]<-COI.Sample.Vectors } save(COI.Sample.List, file="COI Subsampling Data R List.RData") save(Results.Matrix, file="COI Subsampling Results Matrix.RData") write.csv(Results.Matrix, file="COI Subsampling Results Matrix.csv") ## Plotting the SW Subsamples postscript(file="COI diversity subsample plot.ps") plot.new() plot.window(xlim=c(0, 13), ylim=c(-1.75, 1.25)) title(xlab="COI Compared to OMNH Localities", ylab="H' difference", cex=0.5) title(main="Variation in Subsampled Shannon-Wiener Index between COI and OMNH localities") axis(side=1, at=1:12, labels=Matrix.names, tick=T, cex.axis=0.5) axis(2) for (z in 1:(ncol(Results.Matrix)/2)){ ## The following vectors are placing the SW values and Simpson values into corret spots in the plot xPosition.SW<-rep(z, nrow(Results.Matrix)) xPosition.SW<-xPosition.SW-0.2 xPosition.M.H.<-rep(z, nrow(Results.Matrix)) xPosition.M.H.<-xPosition.M.H.+0.2 ### Circle.size<-3 points(xPosition.SW, Results.Matrix[,z], cex= Circle.size, col="salmon") points(xPosition.M.H., Results.Matrix[,(z+(ncol(Results.Matrix)/2))], cex= Circle.size, col="royalblue") } #below is text for the M.H. values value<-round(mean(Results.Matrix[,(z+(ncol(Results.Matrix)/2))]), digits=2) text(xPosition.M.H., max(Results.Matrix[,(z+(ncol(Results.Matrix)/2))])+.15, labels=value, cex=1) #below is point and text for the original COI SW value minus the selected OMNH site coi.SW.orig<-diversity(FaunaData[,1], index="shannon", base=exp(1)) other.SW<-diversity(FaunaData[,z+1], index="shannon", base=exp(1)) b<-coi.SW.orig-other.SW points(xPosition.SW[1], b, col="black", cex= Circle.size) #text(xPosition.SW[1], b+.1, labels=round(b, digits=2), cex=0.5) #below is the text of the mean value of SW for each site. value2<-round(mean(Results.Matrix[,z]), digits=2) text(xPosition.SW[1], max(Results.Matrix[,z])+0.15, labels=value2, cex=1) dev.off()