setwd("~/SWITCHdrive/3_Projects/EUALGAE/Yaiza Tejido/Dryad") library(dichromat) library(ggplot2) mycolors <- c(colorschemes$Categorical.12[8],colorschemes$Categorical.12[2]) ## Add an alpha value to a colour add.alpha <- function(col, alpha=1){ if(missing(col)) stop("Please provide a vector of colours.") apply(sapply(col, col2rgb)/255, 2, function(x) rgb(x[1], x[2], x[3], alpha=alpha)) } # Multiple plot function # # ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects) # - cols: Number of columns in layout # - layout: A matrix specifying the layout. If present, 'cols' is ignored. # # If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE), # then plot 1 will go in the upper left, 2 will go in the upper right, and # 3 will go all the way across the bottom. # multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) { library(grid) # Make a list from the ... arguments and plotlist plots <- c(list(...), plotlist) numPlots = length(plots) # If layout is NULL, then use 'cols' to determine layout if (is.null(layout)) { # Make the panel # ncol: Number of columns of plots # nrow: Number of rows needed, calculated from # of cols layout <- matrix(seq(1, cols * ceiling(numPlots/cols)), ncol = cols, nrow = ceiling(numPlots/cols)) } if (numPlots==1) { print(plots[[1]]) } else { # Set up the page grid.newpage() pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout)))) # Make each plot, in the correct location for (i in 1:numPlots) { # Get the i,j matrix positions of the regions that contain this subplot matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE)) print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row, layout.pos.col = matchidx$col)) } } } df <- read.csv("Tejio_Nunez_2019_Ecol_Eng.csv", header = TRUE, sep = ",", quote="\"", dec=".", fill = TRUE, comment.char="", na.strings = "NA", as.is = F) #Select only columns that are needed and give the short names names(df) names(df)[c(2,4,5,6,7,8,9,11,12,13,14,15,16,17,18)] <- c("day", "algae", "treatment", "replicate", "sterile", "abs", "fluo", "cult.dens", "protozoa", "dw", "bacteria", "no3.hl", "po4.hl", "no3.ic", "po4.ic") df <- df[,c(2,4,5,6,7,8,9,11,13,12,14,15,16,17,18)] names(df) meta.inf <- data.frame( col.nos = 6:15, col.names = names(df)[6:15], dat.class = c(rep("growth",4),rep("biology",2),rep("nutrient",4)), axis.labels = c("Absorbance (750 nm)", "Fluorescence", "Culture density (cells/ml)", "Dry weight (g/l)", "Bacteria (cfu/ml)", "Protozoa (per ml)", "NO3-N (mg/l) (HL)", "PO4-P (mg/l) (HL)", "NO3-N (mg/l) (IC)", "PO4-P (mg/l) (IC)")) #Define factors and reorder levels df$treatment <- factor(df$treatment,levels(df$treatment)[c(5,3,4,2,1,6)]) df$algae <- factor(df$algae,levels(df$algae)[c(2,1)]) df$replicate <- factor(df$replicate) df$sterile <- factor(df$sterile) df <- df[order(df$algae, df$day, df$treatment, df$replicate),] #Calculate daily averages and std errors for growth, biology and nutrients daily.means <- aggregate(df$abs, by = list(df$day, df$treatment, df$algae), FUN = mean)[c(1:3)] names(daily.means) <- c("day","treatment","algae") for(i in 1:nrow(meta.inf)){ variable <- meta.inf$col.names[i] j <- meta.inf$col.nos[i] daily.means <- cbind(daily.means,aggregate(df[,j], by = list(df$day, df$treatment, df$algae), FUN = mean)[4]) daily.means <- cbind(daily.means, aggregate(df[,j], by = list(df$day, df$treatment, df$algae), FUN = sd)[4]) names(daily.means)[c(2*i-1,2*i)+3] <- c(paste(variable,"mean",sep = "."),paste(variable,"sd",sep = ".")) } daily.means$algae.long <- ifelse(daily.means$algae == "Cv", "Chlorella vulgaris", "Acutodesmus obliquus") daily.means$treatment.long <- ifelse(daily.means$treatment == "MM", "Mineral medium", ifelse(daily.means$treatment == "FT 0.2", "Fish tank (0.2 micron)", ifelse(daily.means$treatment == "FT", "Fish tank", ifelse(daily.means$treatment == "DF", "after drum filter", ifelse(daily.means$treatment == "UV", "after UV lamp", "Fish tank (4-7 micron)"))))) write.csv(daily.means, file = "daily_means.csv") #Calculate removal efficiencies days <- c(13,18) #Build new dataframe. Take all relevant factors. nutrients <- df[which(df$day == 0),c(2:15)] nutrients[,(1:4)+length(nutrients)] <- NA #For every algae, take data from last day and place them next to nutrient levels of day 0 for(i in 1:nlevels(df$algae)){ temp1 <- df[which(df$algae == levels(df$algae)[i]),] temp1 <- temp1[which(temp1$day == days[i]),] nutrients[(i-1)*24+1:24,length(nutrients)-(3:0)] <- temp1[,length(temp1)-(3:0)] nutrients[(i-1)*24+1:24,5:10] <- temp1[,6:11] } nutrients[,(1:4)+length(nutrients)] <- 100 - (nutrients[,length(nutrients)-(3:0)] / nutrients[,length(nutrients)-(7:4)] * 100) names(nutrients)[length(nutrients)-(11:0)] <- c("no3.hl.before","po4.hl.before","no3.ic.before","po4.ic.before","no3.hl.after","po4.hl.after","no3.ic.after","po4.ic.after","no3.hl.re","po4.hl.re","no3.ic.re","po4.ic.re") #Calculate averages and errors nutrients.avg <- aggregate(nutrients$no3.hl.re, by = list(treatment = nutrients$treatment, algae = nutrients$algae), FUN = mean)[,c(1:2)] for(i in length(nutrients)-3:0){ nutrients.avg <- cbind(nutrients.avg, aggregate(nutrients[,i], by = list(treatment = nutrients$treatment, algae = nutrients$algae), FUN = mean)[,3]) nutrients.avg <- cbind(nutrients.avg, aggregate(nutrients[,i], by = list(treatment = nutrients$treatment, algae = nutrients$algae), FUN = sd)[,3]/2) names(nutrients.avg)[length(nutrients.avg)-1:0] <- c(paste(names(nutrients)[i],"avg",sep = "."),paste(names(nutrients)[i],"sem",sep = ".")) } write.csv(nutrients.avg, file = "nutrient_removal_means.csv") #Start plotting figures continue <- T while(continue){ choice <- readline("\nWhich figure should be produced? \ta: Fig 1, 2 (growth) \tb: Fig 3 (Protozoa) \tc: not used \td: Fig 4 (nutrients) \te: Fig 5 (Correlations) \tf: not used \tg: presentation \ts: Stats \tt: Growth rates \tx: Exit\n") #Growth curves for all growth data, one page per algae if(choice == "a"){ yliml <- c(0.05,50,200000,0,0.05,40,1000000,0) ylimu <- c(2.5,50000,2000000000,5,2.5,50000,150000000,8) xlimu <- c(13,18) logax <- c("y","y","y","") for(i in 1:nlevels(daily.means$algae)){ filename <- paste("figure",i,".pdf",sep = "") pdf(file=paste("figures/",filename,sep = ""), width = 20*0.394, height = 16*0.394, colormodel = "rgb", pointsize = 11) par(mfcol = c(nlevels(df$treatment)-2,length(which(meta.inf$dat.class == "growth"))-1), mar = c(0,4.5,0,0.5), oma = c(3.5,2,0.5,0)) temp1 <- daily.means[which(daily.means$algae == levels(daily.means$algae)[i]),] for(j in 1:length(which(meta.inf$dat.class == "growth"))){ if(j == 2) next # skip fluorescence data col <- 2*j+2 for(k in 3:nlevels(temp1$treatment)){ temp2 <- temp1[which(temp1$treatment == levels(temp1$treatment)[k]),] temp3 <- temp2[which(!is.na(temp2[,col])),] plot(temp3$day,temp3[,col],type = "n", xlab = "", ylab = "", log = logax[j], xlim = c(0,xlimu[i]), ylim = c(yliml[(i-1)*4+j],ylimu[(i-1)*4+j]) , las = 1, yaxt = "n", xaxt = "n") if(j == 1){ axis(2,at = c(0,0.05,0.1,0.2,0.5,1,2), labels = c(0,0.05,0.1,0.2,0.5,1,2), las = 1) } if(j == 3){ axis(2,at = c(1e5,1e6,1e7,1e8,1e9), labels = c(expression(10^5),expression(10^6),expression(10^7),expression(10^8),expression(10^9)), las = 1) } if(j == 4){ axis(2,at = 0:10, labels = 0:10, las = 1) } if(k == nlevels(temp1$treatment)){ axis(1, at = 2*(0:(ceiling(xlimu[i]/2))), labels = 2*(0:(ceiling(xlimu[i]/2)))) par(xpd = NA) title(xlab = "Time (days)", line = 2) par(xpd = F) }else{ axis(1, at = 2*(0:(ceiling(xlimu[i]/2))), labels = FALSE) } title(ylab = meta.inf$axis.labels[j], line = 3) lower <- c(temp3[,col]-temp3[,col+1]) lower[lower < 0] <- 0.001 upper <- c(temp3[,col])+c(temp3[,col+1]) lines(temp3$day,temp3[,col], col = mycolors[2]) #plot control curves (MM and FT0.2) for(l in 1:2){ ctrl <- temp1[which(temp1$treatment == levels(temp1$treatment)[l]),] ctrl <- ctrl[which(!is.na(ctrl[,col])),] lctrl <- c(ctrl[,col]-ctrl[,col+1]) lctrl[lctrl < 0] <- 0.001 uctrl <- c(ctrl[,col])+c(ctrl[,col+1]) lines(ctrl$day,ctrl[,col], lty = l+1, col = mycolors[1]) points(ctrl$day,ctrl[,col], pch = 20, col = mycolors[1]) for(m in 1:nrow(ctrl)){ lines(c(ctrl$day[m],ctrl$day[m]),c(uctrl[m],lctrl[m]), col = mycolors[1]) } } for(l in 1:nrow(temp3)){ lines(c(temp3$day[l],temp3$day[l]),c(upper[l],lower[l]), col = mycolors[2]) } points(temp3$day,temp3[,col], pch = 20, col = mycolors[2]) if(j == 1){ if(temp3$treatment[1] == "FT 4-7"){ mtext(text = expression(paste("Fish tank (4-7 ",mu,"m)")), side = 2, outer = F, line = 5, cex = 0.8) }else{ mtext(text = temp3$treatment.long[1], side = 2, outer = F, line = 5, cex = 0.8) } } if(k == 3){ #mtext(text = meta.inf$axis.labels[j], side = 3, outer = F, line = 1, cex = 0.8) } } } #mtext(text = paste("Algae species:",temp1$algae.long[1]), side = 3, line = 3, adj = 0, outer = T) dev.off() } answer <- readline("Print more figures? (y/n) ") continue <- ifelse(answer == "y",T,F) } #Protozoa (growth and correlation with algae growth) if(choice == "b"){ pdf(file="figures/figure3.pdf", width = 16*0.394, height = 7*0.394, colormodel = "rgb", pointsize = 11) par(mfrow = c(2,4), mar = c(0,3.5,0,0.5), oma = c(3.5,0,0.5,0)) xlimr1 <- c(20,20) xliml2 <- c(1000,1000) xlimr2 <- c(10000000,10000000) yliml1 <- c(10,10) ylimu1 <- c(2000000,2000000) yliml2 <- c(0.01,0,1000000,0.5) ylimu2 <- c(1.3,0,200000000,2.7) day <- c(13,18) logax <- c("x","x","xy","x") for(i in 1:nlevels(daily.means$algae)){ temp1 <- daily.means[which(daily.means$algae == levels(daily.means$algae)[i]),] temp1 <- temp1[which(!is.na(temp1$protozoa.mean)),] plot(temp1$day,temp1$protozoa.mean, xlab = "", ylab = "", log = "y", type = "n", xlim = c(0,xlimr1[i]), ylim = c(yliml1[i],ylimu1[i]), xaxt = "n", yaxt = "n") title(ylab = "Protozoa (per ml)", line = 2.5) if(i == 2){ par(xpd = NA) title(xlab = "Time (days)", line = 2.5) par(xpd = F) axis(1, at = 5*(0:4), labels = 5*(0:4)) }else{ axis(1, at = 5*(0:4), labels = FALSE) } axis(2, at = c(10,100,1000,10000,100000,1000000), labels = c(10,expression(10^2),expression(10^3),expression(10^4),expression(10^5),expression(10^6)), las = 1) for(j in 3:nlevels(temp1$treatment)){ temp2 <- temp1[which(temp1$treatment == levels(temp1$treatment)[j]),] lines(temp2$day,temp2$protozoa.mean,lty = j-2, col = "black") points(temp2$day,temp2$protozoa.mean, pch = j+12, col = "black") } #mtext(text = temp1$algae.long[i], side = 2, outer = F, line = 5) temp3 <- df[which(df$algae == levels(daily.means$algae)[i]),] temp3 <- temp3[which(!is.na(temp3$protozoa)),] temp3 <- temp3[which(temp3$sterile == 0),] temp3 <- temp3[which(temp3$day == day[i]),] for(j in 1:4){ if(j == 2) next # skip fluorescence data if(j == 3){ temp3$cult.dens <- ifelse(temp3$cult.dens == 0, NA,temp3$cult.dens) } plot(temp3$protozoa,temp3[,j+5], xlab = "", ylab = "", log = logax[j], xlim = c(xliml2[i],xlimr2[i]), ylim = c(yliml2[j],ylimu2[j]), xaxt = "n", yaxt = "n", las = 1, col = "black") title(ylab = meta.inf$axis.labels[j], line = 2.5) if(i == 1){ axis(1, at = c(1000,10000,100000,1000000,10000000), labels = FALSE, las = 1) }else{ axis(1, at = c(1000,10000,100000,1000000,10000000), labels = c(expression(10^3),expression(10^4),expression(10^5),expression(10^6),expression(10^7)), las = 1) par(xpd = NA) title(xlab = "Protozoa (per ml)", line = 2.5) par(xpd = F) } if(j == 1){ axis(2, at = c(0,0.2,0.4,0.6,0.8,1,1.2), labels = c(0,0.2,0.4,0.6,0.8,1,1.2), las = 1)} if(j == 3){ axis(2, at = c(1000000,10000000,100000000), labels = c(expression(10^6),expression(10^7),expression(10^8)), las = 1) } if(j == 4){ axis(2, at = c(0.5,1,1.5,2,2.5), labels = c(0.5,1,1.5,2,2.5), las = 1) } pvalue <- round(cor.test(temp3$protozoa,temp3[,j+5], method = "spearman")$p.value, digits = 3) if(pvalue < 0.001){ mtext(" p < 0.001",side = 1, line = -1.5, adj = 0, cex = 0.5) }else{ if(pvalue > 0.01){ mtext(paste(" p =",round(pvalue,digits = 2)),side = 1, line = -1.5, adj = 0, cex = 0.5) }else{ mtext(paste(" p =",pvalue),side = 1, line = -1.5, adj = 0, cex = 0.5) } } } } dev.off() mean(daily.means$protozoa.mean[which(daily.means$algae == "Cv" & daily.means$day == 13)][3:6]) mean(daily.means$protozoa.mean[which(daily.means$algae == "Ao" & daily.means$day == 14)][3:6]) mean(daily.means$protozoa.mean[which(daily.means$algae == "Ao" & daily.means$day == 18)][3:6]) temp1 <- df[which((df$algae == "Cv" & df$day == 13) & df$sterile == 0),] cor.test(temp1$protozoa,temp1$abs, method = "spearman") cor.test(temp1$protozoa,temp1$cult.dens, method = "spearman") cor.test(temp1$protozoa,temp1$dw, method = "spearman") temp1 <- df[which((df$algae == "Ao" & df$day == 18) & df$sterile == 0),] cor.test(temp1$protozoa,temp1$abs, method = "spearman") cor.test(temp1$protozoa,temp1$cult.dens, method = "spearman") cor.test(temp1$protozoa,temp1$dw, method = "spearman") temp2 <- df[which((df$algae == "Ao" & df$day == 14) & df$sterile == 0),] temp3 <- df[which((df$algae == "Ao" & df$day == 18) & df$sterile == 0),] wilcox.test(temp1$protozoa,temp2$protozoa) wilcox.test(temp2$protozoa,temp3$protozoa) answer <- readline("Print more figures or statistics? (y/n) ") continue <- ifelse(answer == "y",T,F) } #Bacteria if(choice == "c"){ pdf(file="figures/figureBacteria.pdf", width = 11*0.394, height = 10*0.394, colormodel = "rgb", pointsize = 10) par(mfrow = c(2,2), mar = c(4.5,4.5,0.5,0.5), oma = c(0,0,0,0)) xlimr <- c(15,20) day <- c(13,18) mycolors <- c("red","red","black","black","black","black") for(i in 1:nlevels(daily.means$algae)){ temp1 <- daily.means[which(daily.means$algae == levels(daily.means$algae)[i]),] temp1 <- temp1[which(!is.na(temp1$bacteria.mean)),] plot(temp1$day,temp1$bacteria.mean, xlab = "", ylab = "Bacteria (cfu/ml)", log = "y", type = "n", xlim = c(0,xlimr[i]), ylim = c(1000,20000000), xaxt = "n", yaxt = "n") axis(2, at = c(1000,10000,100000,1000000,10000000), labels = c(expression(10^3),expression(10^4),expression(10^5),expression(10^6),expression(10^7)), las = 1) if(i == 2){ axis(1, at = 5*(0:4), labels = 5*(0:4)) title(xlab = "Time (days)", line = 2.5) }else{ axis(1, at = 5*(0:4), labels = FALSE) } for(j in 1:6){ temp2 <- temp1[which(temp1$treatment == levels(temp1$treatment)[j]),] lines(temp2$day,temp2$bacteria.mean, lty = j, col = mycolors[j]) points(temp2$day,temp2$bacteria.mean, pch = 14+j, col = mycolors[j]) } #mtext(text = temp1$algae.long[i], side = 2, outer = F, line = 5) temp1 <- df[which(df$algae == levels(df$algae)[i]),] temp1 <- temp1[which(temp1$day == day[i]),] plot(temp1$protozoa,temp1$bacteria, log = "xy", xaxt = "n", yaxt = "n", ylab = paste("Bacteria (cfu/ml)"), xlab = "", xlim = c(1000,1000000), ylim = c(10000,100000000)) axis(2, at = c(1000,10000,100000,1000000,10000000,100000000), labels = c(expression(10^3),expression(10^4),expression(10^5),expression(10^6),expression(10^7),expression(10^8)), las = 1) if(i == 2){ axis(1, at = c(1000,10000,100000,1000000), labels = c(expression(10^3),expression(10^4),expression(10^5),expression(10^6)), las = 1) title(xlab = paste("Protozoa (no/ml)"), line = 2.5) }else{ axis(1, at = c(1000,10000,100000,1000000), labels = FALSE) } pvalue <- round(cor.test(temp1$protozoa,temp1$bacteria, method = "spearman")$p.value, digits = 3) if(pvalue > 0.1){ mtext(paste(" p =",round(pvalue, digits = 2)),side = 1, line = -2, adj = 0, cex = 0.7) }else{ mtext(paste(" p =",pvalue),side = 1, line = -2, adj = 0, cex = 0.7) } } dev.off() temp1 <- df[which(df$algae == "Cv"),] temp1 <- temp1[which(temp1$day %in% c(3,11,13)),] temp2 <- temp1[which(temp1$sterile == 0),] temp3 <- temp1[which(temp1$sterile == 1),] wilcox.test(temp2$bacteria[which(temp2$day == 3)],temp3$bacteria[which(temp3$day == 3)]) wilcox.test(temp2$bacteria[which(temp2$day == 11)],temp3$bacteria[which(temp3$day == 11)]) wilcox.test(temp2$bacteria[which(temp2$day == 13)],temp3$bacteria[which(temp3$day == 13)]) cor.test(temp2$bacteria[which(temp2$day == 13)],temp2$protozoa[which(temp2$day == 13)], method = "spearman") temp1 <- df[which(df$algae == "Ao"),] temp1 <- temp1[which(temp1$day %in% c(3,11,18)),] temp2 <- temp1[which(temp1$sterile == 0),] temp3 <- temp1[which(temp1$sterile == 1),] wilcox.test(temp2$bacteria[which(temp2$day == 3)],temp3$bacteria[which(temp3$day == 3)]) wilcox.test(temp2$bacteria[which(temp2$day == 11)],temp3$bacteria[which(temp3$day == 11)]) wilcox.test(temp2$bacteria[which(temp2$day == 18)],temp3$bacteria[which(temp3$day == 18)]) cor.test(temp2$bacteria[which(temp2$day == 18)],temp2$protozoa[which(temp2$day == 18)], method = "spearman") answer <- readline("Print more figures? (y/n) ") continue <- ifelse(answer == "y",T,F) } #Nutrient removal if(choice == "d"){ nutrients.m <- rbind(t(matrix(data = nutrients.avg[,3], nrow = 6, ncol = 2)),t(matrix(data = nutrients.avg[,5], nrow = 6, ncol = 2)),t(matrix(data = nutrients.avg[,7], nrow = 6, ncol = 2)),t(matrix(data = nutrients.avg[,9], nrow = 6, ncol = 2))) pdf(file="figures/figure4.pdf", width = 20*0.394, height = 8*0.394, colormodel = "rgb", pointsize = 10) par(mfrow = c(1,2), mar = c(2,3.5,0.5,0), oma = c(0,0,0,0)) nutrients.label <- c("nitrate removal efficiency (%)","phosphate removal efficiency (%)","nitrate removal efficiency (%)","phosphate removal efficiency (%)") yliml <- c(60,90,60,90) #Plot Hach-Lange: i = 1,2; Plot IC: i = 3,4 for(i in 1:2){ barplot(nutrients.m[(i-1)*2+1:2,], ylim = c(yliml[i],100), las = 1, ylab = "", beside = TRUE, names.arg = levels(nutrients.avg$treatment), xpd = F, border = NA, col = gray.colors(2)) title(ylab = nutrients.label[i], line = 2.5) for(j in 1:2){ for(k in 1:6){ middle <- nutrients.avg[(j-1)*6+k,(i*2)+1] sem <- nutrients.avg[(j-1)*6+k,(i*2)+2] lines(c(3*k-2.5+j,3*k-2.5+j),c(middle-sem,middle+sem)) } } if(i == 1){ mtext("Hach-Lange", side = 2, line = 5) } if(i == 3){ mtext("Ion-chromatography", side = 2, line = 5) } } #legend(x = -25,y = 88.5, c("Chlorella vulgaris","Acutodesmus obliquus"), xpd = NA, ncol = 2, bty = "n", fill = gray.colors(2), border = NA) dev.off() set1 <- nutrients[which(nutrients$algae == "Cv" & nutrients$treatment != "MM"),] kruskal.test(set1$no3.hl.re ~ set1$treatment) wilcox.test(set1$no3.hl.re[which(set1$sterile == 0)], set1$no3.hl.re[which(set1$sterile == 1)]) mean(set1$no3.hl.re[which(set1$sterile == 1)]) mean(set1$no3.hl.re[which(set1$sterile == 0)]) #kruskal.test(set1$no3.ic.re ~ set1$treatment) kruskal.test(set1$po4.hl.re ~ set1$treatment) wilcox.test(set1$po4.hl.re[which(set1$sterile == 0)], set1$po4.hl.re[which(set1$sterile == 1)]) mean(set1$po4.hl.re[which(set1$sterile == 1)]) mean(set1$po4.hl.re[which(set1$sterile == 0)]) #kruskal.test(set1$po4.ic.re ~ set1$treatment) set2 <- nutrients[which(nutrients$algae == "Ao" & nutrients$treatment != "MM"),] kruskal.test(set2$no3.hl.re ~ set1$treatment) wilcox.test(set2$no3.hl.re[which(set2$sterile == 0)], set2$no3.hl.re[which(set2$sterile == 1)]) mean(set2$no3.hl.re[which(set2$sterile == 1)]) mean(set2$no3.hl.re[which(set2$sterile == 0)]) #kruskal.test(set2$no3.ic.re ~ set2$treatment) kruskal.test(set2$po4.hl.re ~ set2$treatment) wilcox.test(set2$po4.hl.re[which(set2$sterile == 0)], set2$po4.hl.re[which(set2$sterile == 1)]) mean(set2$po4.hl.re[which(set2$sterile == 1)]) mean(set2$po4.hl.re[which(set2$sterile == 0)]) #kruskal.test(set2$po4.ic.re ~ set2$treatment) wilcox.test(set1$no3.hl.re[which(set1$sterile == 1)], set2$no3.hl.re[which(set2$sterile == 1)]) wilcox.test(set1$no3.hl.re[which(set1$sterile == 0)], set2$no3.hl.re[which(set2$sterile == 0)]) wilcox.test(set1$po4.hl.re[which(set1$sterile == 1)], set2$po4.hl.re[which(set2$sterile == 1)]) wilcox.test(set1$po4.hl.re[which(set1$sterile == 0)], set2$po4.hl.re[which(set2$sterile == 0)]) answer <- readline("Print more figures? (y/n) ") continue <- ifelse(answer == "y",T,F) } #Nutrient removal and growth if(choice == "e"){ pdf(file="figures/figure5.pdf", width = 16*0.394, height = 16*0.394, colormodel = "rgb", pointsize = 10) par(mfrow = c(4,4), mar = c(0,0,0,0), oma = c(4.5,4.5,0.5,0.5)) logax <- c("x", "x", "x", "x", "x") for(i in 1:nlevels(nutrients$algae)){ temp1 <- nutrients[which(nutrients$algae == levels(nutrients$algae)[i]),] for(k in 19:20){ for(j in 5:9){ if(j == 6) next # skip fluorescence data plot(temp1[,j],temp1[,k], xlab = "", ylab = "", ylim = c(c(60,80)[k-18],c(104,102)[k-18]), log = logax[j-4], las = 1, yaxt = "n", xaxt = "n", xlim = c(c(0.05,1,1e6,0.5,1e3)[j-4],c(2.5,1,2e9,12,2e6)[j-4])) if(j == 5){ par(xpd = NA) title(ylab = paste(c("nitrate","phosphate")[k-18],"removal efficiency (%)")) par(xpd = F) axis(2, at = (0:4)*c(10,5)[k-18]+c(60,80)[k-18], labels = (0:4)*c(10,5)[k-18]+c(60,80)[k-18], las = 1) }else{ axis(2, at = (0:4)*c(10,5)[k-18]+c(60,80)[k-18], labels = FALSE) } if(i == 2 & k == 20){ par(xpd = NA) title(xlab = c("Absorbance (750 nm)", "Fluorescence", "Culture density (cells/ml)", "Dry weight (g/l)", "Protozoa (per ml)")[j-4]) par(xpd = F) if(j == 5){ axis(1, at = c(0.05,0.1,0.2,0.5,1,2), labels = c(0.05,0.1,0.2,0.5,1,2)) } if(j == 7){ axis(1, at = c(1e6,1e7,1e8,1e9), labels = c(expression(10^6),expression(10^7),expression(10^8),expression(10^9))) } if(j == 8){ axis(1, at = c(0.5,1,2,5,10), labels = c(0.5,1,2,5,10)) } if(j == 9){ axis(1, at = c(1e3,1e4,1e5,1e6), labels = c(expression(10^3),expression(10^4),expression(10^5),expression(10^6))) } }else{ if(j == 5){ axis(1, at = c(0.05,0.1,0.2,0.5,1,2), labels = F) } if(j == 7){ axis(1, at = c(1e6,1e7,1e8,1e9), labels = F) } if(j == 8){ axis(1, at = c(0.5,1,2,5,10), labels = F) } if(j == 9){ axis(1, at = c(1e3,1e4,1e5,1e6), labels = F) } } pvalue <- round(cor.test(temp1[,j],temp1[,k], method = "spearman")$p.value,digits = 3) if(pvalue < 0.001){ mtext(paste(" p < 0.001"),side = 1, line = -2, adj = 0, cex = 0.7) }else{ if(pvalue < 0.05){ mtext(paste(" p =",round(pvalue, digits = 3)),side = 1, line = -2, adj = 0, cex = 0.7) }else{ mtext(paste(" p =",round(pvalue, digits = 2)),side = 1, line = -2, adj = 0, cex = 0.7) } } } } } dev.off() answer <- readline("Print more figures? (y/n) ") continue <- ifelse(answer == "y",T,F) } #For presentation if(choice == "g"){ pdf(file = "figures/figure_for_presentation_1.pdf", width = 20*0.394, height = 8*0.394, colormodel = "rgb", pointsize = 10) par(mfcol = c(1,2), mar = c(3.5,5,2,0.5), bty = "n", cex = 1.2) lineweight <- c(1,2,2,2,2,2) mycolors2 <- c("black", rep(colorschemes$Categorical.12[6], 4), colorschemes$Categorical.12[10]) for(i in 1:nlevels(daily.means$algae)){ temp1 <- daily.means[which(daily.means$algae == levels(daily.means$algae)[i]),] plot(1, type = "n", xlab = "", ylab = "", xlim = c(0,20), ylim = c(0,2) , las = 1) title(xlab = "time (days)", line = 2) title(ylab = "Absorbance (750 nm)", line = 3) title(main = paste(temp1$algae.long[1])) mycolors2 <- add.alpha(mycolors2, alpha=0.1) temp1$treatment = factor(temp1$treatment,levels(temp1$treatment)[c(1,3:6,2)]) for(k in 1:nlevels(temp1$treatment)){ temp2 <- temp1[which(temp1$treatment == levels(temp1$treatment)[k]),] temp3 <- temp2[which(!is.na(temp2$abs.mean)),] polygon(c(temp3$day,rev(temp3$day)),c(temp3$abs.mean-temp3$abs.sd,rev(temp3$abs.mean+temp3$abs.sd)), border = NA, col = mycolors2[k]) } mycolors2 <- add.alpha(mycolors2, alpha=1) for(k in 1:nlevels(temp1$treatment)){ temp2 <- temp1[which(temp1$treatment == levels(temp1$treatment)[k]),] temp3 <- temp2[which(!is.na(temp2$abs.mean)),] lines(temp3$day,temp3$abs.mean, lwd = lineweight[k], col = mycolors2[k]) } } dev.off() days <- c(13,18) temp <- df[which(df$treatment != "MM"),] temp$algae <- ifelse(temp$algae == "Cv", "C. vulgaris", "A. obliquus") temp$algae <- factor(temp$algae) temp$algae <- factor(temp$algae, levels(temp$algae)[c(2,1)]) temp$sterile <- ifelse(temp$sterile == 1,"sterile","not sterile") temp$sterile <- factor(temp$sterile) temp$sterile <- factor(temp$sterile, levels(temp$sterile)[c(2,1)]) temp$Wastewater <- temp$sterile temp1 <- temp[which(temp$day == 13 & temp$algae == "C. vulgaris"),] temp <- rbind(temp1, temp[which(temp$day == 18 & temp$algae == "A. obliquus"),]) p1 <- ggplot(data=temp, aes(x=algae, y=dw, fill=Wastewater)) + geom_boxplot() + theme(axis.title.x = element_blank(), axis.text.x = element_text(size=12), axis.text.y = element_text(size = 12), legend.text = element_text(size = 12)) + ggtitle("Biomass dry weight") + ylab("Biomass dry weight (g/l)") + scale_fill_manual(values=c(colorschemes$Categorical.12[10], colorschemes$Categorical.12[6])) p2 <- ggplot(data=temp, aes(x=algae, y=cult.dens, fill=Wastewater)) + geom_boxplot() + theme(axis.title.x = element_blank(), axis.text.x = element_text(size=12), axis.text.y = element_text(size = 12), legend.text = element_text(size = 12)) + scale_y_log10() + ggtitle("Cell density") + ylab("Cells per ml") + scale_fill_manual(values=c(colorschemes$Categorical.12[10], colorschemes$Categorical.12[6])) pdf(file = "figures/figure_for_presentation_2.pdf", width = 24*0.394, height = 8*0.394, colormodel = "rgb", pointsize = 10) multiplot(p1,p2, cols = 2) dev.off() answer <- readline("Print more figures? (y/n) ") continue <- ifelse(answer == "y",T,F) } #Perform some statistics on abs, fluo, cult.dens, etc. if(choice == "s"){ temp1 <- data.frame(day = df$day, algae = df$algae, treatment = df$treatment, sterile = df$sterile, abs = df$abs, cult.dens = df$cult.dens, dw = df$dw, bacteria = df$bacteria, protozoa = df$protozoa) #temp1 <- temp1[which(temp1$treatment != "MM"),] #temp1$treatment <- factor(temp1$treatment,levels(temp1$treatment)[c(2,3,4,5,6)]) #bonf.corr <- nlevels(temp1$treatment)*(nlevels(temp1$treatment)-1)/2 days <- list(c(1:13),c(1:18),c(6,10,13),c(6,10,14,18),c(6,10,13),c(6,10,14,18)) counter <- 0 for(i in 5:7){ measurement <- names(temp1)[i] cat(paste("\nMeasured variable:",measurement)) for(j in 1:nlevels(temp1$algae)){ counter <- counter + 1 currentalgae <- levels(temp1$algae)[j] cat(paste("\n\tAlgae species:",currentalgae)) temp2 <- temp1[which(temp1$algae == currentalgae),] max.day <- max(temp2$day[which(!is.na(temp2[,i]))]) temp3 <- temp2[which(temp2$day == max.day),] ttest <- t.test(temp3[which(temp3$treatment == "MM"),i],temp3[which(temp3$treatment == "FT 0.2"),i], var.equal = F) ratio <- mean(temp3[which(temp3$treatment == "MM"),i])/mean(temp3[which(temp3$treatment == "FT 0.2"),i]) cat(paste("\n\t\tComparison MM/FT0.2:\n\t\t\tRatio at day ",max.day," = ",round(ratio,digits = 2),", df = ",round(ttest$parameter, digits = 2),", t = ",round(ttest$statistic, digits = 2),", p = ",signif(ttest$p.value, digits = 3), sep = "")) temp3 <- temp2[which(temp2$treatment != "MM"),] temp3$treatment <- factor(temp3$treatment,levels(temp3$treatment)[c(2,3,4,5,6)]) cat(paste("\n\t\tComparison sterile/non-sterile:", sep = "")) for(k in 1:lengths(days)[counter]){ temp4 <- temp3[which(temp3$day == days[[counter]][k]),] ttest <- t.test(temp4[,i] ~ temp4$sterile, var.equal = F) cat(paste("\n\t\t\tDay ",days[[counter]][k],": df = ",round(ttest$parameter, digits = 2),", t = ",round(ttest$statistic,digits = 2),", p = ",signif(ttest$p.value, digits = 3), sep = "")) } temp3 <- temp3[which(temp3$treatment != "FT 0.2"),] temp3$treatment <- factor(temp3$treatment,levels(temp3$treatment)[c(2,3,4,5)]) cat(paste("\n\t\tComparison UV versus (FT4-7, FT, DF):", sep = "")) temp3$uv <- ifelse(temp3$treatment == "UV",1,0) for(k in 1:lengths(days)[counter]){ temp4 <- temp3[which(temp3$day == days[[counter]][k]),] ttest <- t.test(temp4[,i] ~ temp4$uv, var.equal = F) mean.1 <- round(mean(temp4[which(temp4$uv == 1),i]), digits = 2) mean.2 <- round(mean(temp4[which(temp4$uv == 0),i]), digits = 2) cat(paste("\n\t\t\tDay ",days[[counter]][k],": uv = ",mean.1,", others: ",mean.2,", df = ",round(ttest$parameter, digits = 2),", t = ",round(ttest$statistic,digits = 2),", p = ",signif(ttest$p.value, digits = 3), sep = "")) } } } cat("\n\nFinal dry weight:") for(i in 1:nlevels(df$algae)){ cat(paste("\n\t",levels(df$algae)[i])) temp1 <- daily.means[which(daily.means$algae == levels(daily.means$algae)[i]),] max.day <- max(temp1$day[which(!is.na(temp1$dw.mean))]) temp1 <- temp1[which(temp1$day == max.day),] for(j in 1:nlevels(temp1$treatment)){ temp2 <- temp1[which(temp1$treatment == levels(temp1$treatment)[j]),] cat(paste("\n\t\t",levels(temp2$treatment)[j],":\t",round(temp2$dw.mean, digits = 2)," +/- ",round(temp2$dw.sd/2, digits = 2)," (mean and SEM)",sep = "")) } } temp1 <- df[which(df$treatment == "FT 0.2"),] set1 <- temp1$dw[which((temp1$day == 13) & (temp1$algae == "Cv"))] set2 <- temp1$dw[which((temp1$day == 18) & (temp1$algae == "Ao"))] cat("\n\tin sterile water, Cv is ",round(mean(set1)/mean(set2),digits = 2)," higher than To:") ttest <- t.test(set1, set2, var.equal = F) cat(paste(" df = ",round(ttest$parameter, digits = 2),", t = ",round(ttest$statistic,digits = 2),", p = ",signif(ttest$p.value, digits = 3), sep = "")) temp1 <- df[which(df$sterile == 0),] set1 <- temp1$dw[which((temp1$day == 13) & (temp1$algae == "Cv"))] set2 <- temp1$dw[which((temp1$day == 18) & (temp1$algae == "Ao"))] cat("\n\tin non-sterile water, To is ",round(mean(set2)/mean(set1),digits = 2)," higher than Cv:") ttest <- t.test(set2, set1, var.equal = F) cat(paste(" df = ",round(ttest$parameter, digits = 2),", t = ",round(ttest$statistic,digits = 2),", p = ",signif(ttest$p.value, digits = 3), sep = "")) answer <- readline("\nPrint more figures? (y/n) ") continue <- ifelse(answer == "y",T,F) } if(choice == "t"){ pdf(file = "figures/figure_with_growth_curves.pdf", width = 28*0.394, height = 10*0.394, colormodel = "rgb", pointsize = 10) par(mfrow = c(2,6), mar = c(4.5,5,2,0.2), oma = c(0,3,0,0), xpd = F) gr <- data.frame(algae = factor(levels = levels(df$algae)), treatment = factor(levels = levels(df$treatment)), replicate = factor(levels = levels(df$replicate)), max.rate = numeric(), max.rate.day = numeric(), group = factor(levels = 1:12)) group.no <- 1 for(i in 1:nlevels(df$algae)){ temp1 <- df[which(df$algae == levels(df$algae)[i]),] for(j in 1:nlevels(temp1$treatment)){ temp2 <- temp1[which(temp1$treatment == levels(temp1$treatment)[j]),] plot(temp2$abs ~ temp2$day, ylab = "Absorbance (750 nm)", xlab = "Time (days)", ylim = c(0, 2), las = 1) title(main = paste(levels(temp2$treatment)[j])) for(k in 1:nlevels(temp2$replicate)){ temp3 <- temp2[which(temp2$replicate == levels(temp2$replicate)[k]),] temp4 <- temp3[which(temp3$day <= c(8,10)[i]),] spl <- smooth.spline(temp4$day, temp4$abs, df = 6) lines(spl) maxy <- max(predict(spl, 0:(c(8,10)[i]*1000)/1000,deriv = 1)$y) maxx <- which.max(predict(spl, 0:(c(8,10)[i]*1000)/1000,deriv = 1)$y)/1000 arrows(maxx,0,maxx,0.2, length = 0.05, lwd = 0.5) newrow <- data.frame(algae = levels(df$algae)[i], treatment = levels(temp1$treatment)[j], replicate = levels(temp2$replicate)[k], max.rate = maxy, max.rate.day = maxx, group = group.no) gr <- rbind(gr, newrow) } group.no <- group.no + 1 } mtext(text = paste(c("C. vulgaris","T. obliquus")[i]), side = 2, outer = F, line = 87, cex = 1) } dev.off() gr$group <- factor(gr$group) avg.rates <- aggregate(gr$max.rate, by = list(gr$treatment, gr$algae), FUN = mean) avg.rates$sem <- aggregate(gr$max.rate, by = list(gr$treatment, gr$algae), FUN = sd)[,3]/sqrt(4) names(avg.rates) <- c("treatment","algae","mean","sem") for(i in 1:nlevels(gr$algae)){ pdf(file = paste("figures/figure_growth_rates_",levels(gr$algae)[i],".pdf",sep = ""), width = 12*0.394, height = 8*0.394, colormodel = "rgb", pointsize = 10) par(mar = c(2,4,0.2,0.2)) temp1 <- gr[which(gr$algae == levels(gr$algae)[i]),] boxplot(temp1$max.rate ~ temp1$treatment, las = 1, ylab = "Maximum growth rate per day", ylim = c(0.12,0.25)) dev.off() kruskal.test(temp1$max.rate[5:24] ~ temp1$treatment[5:24]) } model1 <- aov(gr$max.rate ~ gr$group) summary.aov(model1) TukeyHSD(model1) answer <- readline("\nPrint more figures? (y/n) ") continue <- ifelse(answer == "y",T,F) } if(choice == "x"){ continue <- FALSE } } cat("\n\nScript end.\n\n")