### Diversity calculations rm(list=ls()) rull <- read.csv("dataset1894Rull1987_red.csv", sep=",") #read in dataset with taxa columns in numeric format #### (removed 4 rows from original file: analysisUnitName, Thickness, Sample name, radiocarbon dates, ) str(rull) # 1. step: I need to subsample the dataset in a way that I can process with ordinations (remove all non-count parts of data # and remove all taxa that are not within the pollen sum (for this I need to keep the "group" column and taxa column)) rull2<-rull[-2,-c(3:5)] #set all na values to zero values in the counts rull2[is.na(rull2)]<-0 str(rull2) #to calculate percentages, I need the whole terrestrial pollen dataset rull_terr<-with(rull2,rull2[rull2$group == "TRSH" | rull2$group == "UPHE" ,]) rull_terr_nogroup<-rull_terr[,-2] depth<-rull[1,c(1,6:ncol(rull))] rull_terr_nogroup_depth<-rbind(rull[1,c(1,6:ncol(rull))],rull_terr_nogroup) trull<-t(rull_terr_nogroup_depth) trull<-as.data.frame(trull) str(trull) write.table(trull,"trull_terr.csv",col.names = FALSE,sep=",") rull_div<-read.csv("trull_terr.csv") str(rull_div) #Remove rows with pollen sum below 100 rull.pollen.sum<-rowSums(rull_div[,-c(1:2)]) rull_div_sum<-cbind(rull_div,rull.pollen.sum) rull_div_def<-with(rull_div_sum,rull_div_sum[rull_div_sum$rull.pollen.sum > 100,]) rull_div_def_def<-rull_div_def[,-ncol(rull_div_def)] #diversity calculations: library (vegan) library(ggplot2) library(xlsx) library(gridExtra) pollen.data<-rull_div_def_def pollen.data[is.na(pollen.data)]<-0 # get pollen levels and taxa on separate sheets pollen.counts<-pollen.data [,-1] pollen.level<-pollen.data [,c(2)] # calculate pollen sums pollen.sum<-rowSums(pollen.counts) #### I need to remove pollen sums of less than 100 # create two 'lists' of same dimension as pollen.counts pollen.percentages<-pollen.counts pollen.proportions<-pollen.counts # loop in each row and in each column to get proportions and % for(j in 1:nrow(pollen.counts)) { for(k in 1:ncol(pollen.counts)) { pollen.proportions[j,k]=pollen.counts[j,k]/pollen.sum[j] pollen.percentages[j,k]=(pollen.counts[j,k]/pollen.sum[j])*100 } } # check if sum of % and proportions is 100 and 1... print(rowSums(pollen.percentages)) print(rowSums(pollen.proportions)) ###### # Rarefaction for pollen richness ###### # vector's instance pollen.rarefaction<-c(0) # minimum pollen sum min.pollensum<-100 # rarefaction pollen.counts<-sapply(pollen.counts,as.integer) pollen.rare.end<-round((rarefy(pollen.counts,min.pollensum, se=TRUE)),2) # transpose data, and put it in a dataframe, this is easier to handle pollen.rare.t<-t(pollen.rare.end) pollen.rare<-as.data.frame(pollen.rare.t) # get S=rarefaction (column=1), se=standard error (column=2) rare.index<-pollen.rare[,1] rare.error<-pollen.rare[,2] # get min/max values based on the Standard Error of the Mean rare.error.plus<-rare.index+rare.error rare.error.minus<-rare.index-rare.error write.xlsx(data.frame(pollen.level,rare.index), "biodiversity_rull.xlsx", sheetName="richness")