#Analyses for identifying the Minimum Inhibitory Concentration (MIC) of 6gingerol #6 Concentrations tested = 100, 75, 50, 25, 10, 5 ug/mL 6gingerol Gingerol<- read.csv("~DataFile4_6gingerol.csv") #Install packages library(ggplot2) library(ggpubr) library(DescTools) library(nlme) library(plyr) library(psych) #Making "Treatment" (6gingerol concentration or negative control), "Plate_Date" (three trials on three separate days), and "Plate" (plate number within trial) factors Gingerol$Treatment<-as.factor(Gingerol$Treatment) Gingerol$Plate_Date<-as.factor(Gingerol$Plate_Date) Gingerol$Plate<-as.factor(Gingerol$Plate) str(Gingerol) #Boxplot for 6gingerol concentrations across all trials and plates Gingerol$Treatment <- factor(Gingerol$Treatment, levels = c("5", "10", "25","50","75","100","Negative")) Gingerol_Boxplot<-ggplot(Gingerol, aes(x=Treatment, y=Conc,color=Treatment,palette="jco")) +geom_boxplot() +geom_boxplot()+ labs(y="Zoospore Viability (% of Positive Control)", x="Treatment (ug/mL)",font="bold",color="black",size=16)+ theme(legend.position="none",text=element_text(face="bold",colour="black",size=16),axis.text.y=element_text(face="bold",colour="black"),axis.text.x=element_text(face="bold",colour="black",size = 13)) Gingerol_Boxplot #Comparing six generalized least-squares (GLS) models with normal errors #(Treatment+Plate_Date+Plate, Treatment*Plate_Date, Treatment*Plate, Treatment*Plate_Date+Plate, Treatment*Plate+Plate_Date, Treatment*Plate_Date*Plate) #Making the negative control the reference level Gingerol$Treatment<-relevel(Gingerol$Treatment,ref="Negative") #1 Treatment+Plate_Date+Plate G_gls_mod_1<-gls(Conc~Treatment+Plate+Plate_Date,weights=varIdent(form=~1|Treatment),method = "ML",data=Gingerol) anova(G_gls_mod_1) summary(G_gls_mod_1) plot(G_gls_mod_1) qqnorm(resid(G_gls_mod_1, type='pearson')) qqline(resid(G_gls_mod_1, type='pearson')) #2 Treatment*Plate_Date G_gls_mod_2<-gls(Conc~Treatment*Plate_Date,weights=varIdent(form=~1|Treatment),method = "ML",data=Gingerol) anova(G_gls_mod_2) summary(G_gls_mod_2) plot(G_gls_mod_2) qqnorm(resid(G_gls_mod_2, type='pearson')) qqline(resid(G_gls_mod_2, type='pearson')) #3 Treatment*Plate G_gls_mod_3<-gls(Conc~Treatment*Plate,weights=varIdent(form=~1|Treatment),method = "ML",data=Gingerol) anova(G_gls_mod_3) summary(G_gls_mod_3) plot(G_gls_mod_3) qqnorm(resid(G_gls_mod_3, type='pearson')) qqline(resid(G_gls_mod_3, type='pearson')) #4 Treatment*Plate_Date + Plate G_gls_mod_4<-gls(Conc~Treatment*Plate_Date+Plate,weights=varIdent(form=~1|Treatment),method = "ML",data=Gingerol) anova(G_gls_mod_4) summary(G_gls_mod_4) plot(G_gls_mod_4) qqnorm(resid(G_gls_mod_4, type='pearson')) qqline(resid(G_gls_mod_4, type='pearson')) #5 Treatment*Plate + Plate_Date G_gls_mod_5<-gls(Conc~Treatment*Plate+Plate_Date,weights=varIdent(form=~1|Treatment),method = "ML",data=Gingerol) anova(G_gls_mod_5) summary(G_gls_mod_5) plot(G_gls_mod_5) qqnorm(resid(G_gls_mod_5, type='pearson')) qqline(resid(G_gls_mod_5, type='pearson')) #6 Treatment*Plate*Plate_Date G_gls_mod_6<-gls(Conc~Treatment*Plate*Plate_Date,weights=varIdent(form=~1|Treatment),method = "ML",data=Gingerol) anova(G_gls_mod_6) summary(G_gls_mod_6) plot(G_gls_mod_6) qqnorm(resid(G_gls_mod_6, type='pearson')) qqline(resid(G_gls_mod_6, type='pearson')) shapiro.test(G_gls_mod_6$residuals) #Comparing AIC between the six models AIC(G_gls_mod_1, G_gls_mod_2, G_gls_mod_3, G_gls_mod_4, G_gls_mod_5,G_gls_mod_6) #Model 6 (Treatment*Plate_Date*Plate) is our best fit model (lowest AIC). #Subset data into nine different plates (Plates 1-9) GPlate1<-subset(Gingerol,Plate_Date=="7.1.21" & Plate=="1") GPlate2<-subset(Gingerol,Plate_Date=="7.1.21" & Plate=="2") GPlate3<-subset(Gingerol,Plate_Date=="7.1.21" & Plate=="3") GPlate4<-subset(Gingerol,Plate_Date=="7.2.21" & Plate=="1") GPlate5<-subset(Gingerol,Plate_Date=="7.2.21" & Plate=="2") GPlate6<-subset(Gingerol,Plate_Date=="7.2.21" & Plate=="3") GPlate7<-subset(Gingerol,Plate_Date=="7.3.21" & Plate=="1") GPlate8<-subset(Gingerol,Plate_Date=="7.3.21" & Plate=="2") GPlate9<-subset(Gingerol,Plate_Date=="7.3.21" & Plate=="3") #Subset data into three plate dates (Trials 1-3) GPlate_Date1<-subset(Gingerol,Plate_Date=="7.1.21") GPlate_Date2<-subset(Gingerol,Plate_Date=="7.2.21") GPlate_Date3<-subset(Gingerol,Plate_Date=="7.3.21") #GLS models & boxplots for plates 1-9 #Plate 1 Plate1_boxplot<-ggplot(GPlate1, aes(x=Treatment, y=Conc,color=Treatment,palette="jco")) +geom_boxplot()+theme(legend.position = "none")+ ggtitle("Plate 1") Plate1_boxplot #Model G_P1_gls<-gls(Conc~Treatment,weights=varIdent(form=~1|Treatment),data=GPlate1) summary(G_P1_gls) plot(G_P1_gls) qqnorm(resid(G_P1_gls, type='pearson')) qqline(resid(G_P1_gls, type='pearson')) G_pVal1<-summary(G_P1_gls)$tTable[,4] G_pVal1 #Plate 2 #Boxplot Plate2_boxplot<-ggplot(GPlate2, aes(x=Treatment, y=Conc,color=Treatment,palette="jco")) +geom_boxplot()+theme(legend.position = "none")+ ggtitle("Plate 2") Plate2_boxplot #Model G_P2_gls<-gls(Conc~Treatment,weights=varIdent(form=~1|Treatment),data=GPlate2) summary(G_P2_gls) plot(G_P2_gls) qqnorm(resid(G_P2_gls, type='pearson')) qqline(resid(G_P2_gls, type='pearson')) G_pVal2<-summary(G_P2_gls)$tTable[,4] G_pVal2 #Plate 3 Plate3_boxplot<-ggplot(GPlate3, aes(x=Treatment, y=Conc,color=Treatment,palette="jco")) +geom_boxplot()+theme(legend.position = "none")+ ggtitle("Plate 3") Plate3_boxplot #Model G_P3_gls<-gls(Conc~Treatment,weights=varIdent(form=~1|Treatment),data=GPlate3) summary(G_P3_gls) plot(G_P3_gls) qqnorm(resid(G_P3_gls, type='pearson')) qqline(resid(G_P3_gls, type='pearson')) G_pVal3<-summary(G_P3_gls)$tTable[,4] G_pVal3 #Plate 4 #Boxplot Plate4_boxplot<-ggplot(GPlate4, aes(x=Treatment, y=Conc,color=Treatment,palette="jco")) +geom_boxplot()+theme(legend.position = "none")+ ggtitle("Plate 4") Plate4_boxplot #Model G_P4_gls<-gls(Conc~Treatment,weights=varIdent(form=~1|Treatment),data=GPlate4) summary(G_P4_gls) plot(G_P4_gls) qqnorm(resid(G_P4_gls, type='pearson')) qqline(resid(G_P4_gls, type='pearson')) G_pVal4<-summary(G_P4_gls)$tTable[,4] G_pVal4 #Plate 5 #Boxplot Plate5_boxplot<-ggplot(GPlate5, aes(x=Treatment, y=Conc,color=Treatment,palette="jco")) +geom_boxplot()+theme(legend.position = "none")+ ggtitle("Plate 5") Plate5_boxplot #Model G_P5_gls<-gls(Conc~Treatment,weights=varIdent(form=~1|Treatment),data=GPlate5) summary(G_P5_gls) plot(G_P5_gls) qqnorm(resid(G_P5_gls, type='pearson')) qqline(resid(G_P5_gls, type='pearson')) G_pVal5<-summary(G_P5_gls)$tTable[,4] G_pVal5 #Plate 6 #Boxplot Plate6_boxplot<-ggplot(GPlate6, aes(x=Treatment, y=Conc,color=Treatment,palette="jco")) +geom_boxplot()+theme(legend.position = "none")+ ggtitle("Plate 6") Plate6_boxplot #Model G_P6_gls<-gls(Conc~Treatment,weights=varIdent(form=~1|Treatment),data=GPlate6) summary(G_P6_gls) plot(G_P6_gls) qqnorm(resid(G_P6_gls, type='pearson')) qqline(resid(G_P6_gls, type='pearson')) G_pVal6<-summary(G_P6_gls)$tTable[,4] G_pVal6 #Plate 7 #Boxplot Plate7_boxplot<-ggplot(GPlate7, aes(x=Treatment, y=Conc,color=Treatment,palette="jco")) +geom_boxplot()+theme(legend.position = "none")+ ggtitle("Plate 7") Plate7_boxplot #Model G_P7_gls<-gls(Conc~Treatment,weights=varIdent(form=~1|Treatment),data=GPlate7) summary(G_P7_gls) plot(G_P7_gls) qqnorm(resid(G_P7_gls, type='pearson')) qqline(resid(G_P7_gls, type='pearson')) G_pVal7<-summary(G_P7_gls)$tTable[,4] G_pVal7 #Plate 8 #Boxplot Plate8_boxplot<-ggplot(GPlate8, aes(x=Treatment, y=Conc,color=Treatment,palette="jco")) +geom_boxplot()+theme(legend.position = "none")+ ggtitle("Plate 8") Plate8_boxplot #Model G_P8_gls<-gls(Conc~Treatment,weights=varIdent(form=~1|Treatment),data=GPlate8) summary(G_P8_gls) plot(G_P8_gls) qqnorm(resid(G_P8_gls, type='pearson')) qqline(resid(G_P8_gls, type='pearson')) G_pVal8<-summary(G_P8_gls)$tTable[,4] G_pVal8 #Plate 9 #Boxplot Plate9_boxplot<-ggplot(GPlate9, aes(x=Treatment, y=Conc,color=Treatment,palette="jco")) +geom_boxplot()+theme(legend.position = "none")+ ggtitle("Plate 9") Plate9_boxplot #Model G_P9_gls<-gls(Conc~Treatment,weights=varIdent(form=~1|Treatment),data=GPlate9) summary(G_P9_gls) plot(G_P9_gls) qqnorm(resid(G_P9_gls, type='pearson')) qqline(resid(G_P9_gls, type='pearson')) G_pVal9<-summary(G_P9_gls)$tTable[,4] G_pVal9 #Renaming treatments for plate numbers (Plates 1-9) G_pVal1 <- rename(G_pVal1, replace = c("(Intercept)" = "InterceptP1", "Treatment100"="Treatment100P1", "Treatment75"="Treatment75P1", "Treatment50"="Treatment50P1", "Treatment25"="Treatment25P1", "Treatment10"="Treatment10P1", "Treatment5"="Treatment5P1")) G_pVal2 <- rename(G_pVal2, replace = c("(Intercept)" = "InterceptP2", "Treatment100"="Treatment100P2", "Treatment75"="Treatment75P2", "Treatment50"="Treatment50P2", "Treatment25"="Treatment25P2", "Treatment10"="Treatment10P2", "Treatment5"="Treatment5P2")) G_pVal3 <- rename(G_pVal3, replace = c("(Intercept)" = "InterceptP3", "Treatment100"="Treatment100P3", "Treatment75"="Treatment75P3", "Treatment50"="Treatment50P3", "Treatment25"="Treatment25P3", "Treatment10"="Treatment10P3", "Treatment5"="Treatment5P3")) G_pVal4 <- rename(G_pVal4, replace = c("(Intercept)" = "InterceptP4", "Treatment100"="Treatment100P4", "Treatment75"="Treatment75P4", "Treatment50"="Treatment50P4", "Treatment25"="Treatment25P4", "Treatment10"="Treatment10P4", "Treatment5"="Treatment5P4")) G_pVal5 <- rename(G_pVal5, replace = c("(Intercept)" = "InterceptP5", "Treatment100"="Treatment100P5", "Treatment75"="Treatment75P5", "Treatment50"="Treatment50P5", "Treatment25"="Treatment25P5", "Treatment10"="Treatment10P5", "Treatment5"="Treatment5P5")) G_pVal6 <- rename(G_pVal6, replace = c("(Intercept)" = "InterceptP6", "Treatment100"="Treatment100P6", "Treatment75"="Treatment75P6", "Treatment50"="Treatment50P6", "Treatment25"="Treatment25P6", "Treatment10"="Treatment10P6", "Treatment5"="Treatment5P6")) G_pVal7 <- rename(G_pVal7, replace = c("(Intercept)" = "InterceptP7", "Treatment100"="Treatment100P7", "Treatment75"="Treatment75P7", "Treatment50"="Treatment50P7", "Treatment25"="Treatment25P7", "Treatment10"="Treatment10P7", "Treatment5"="Treatment5P7")) G_pVal8 <- rename(G_pVal8, replace = c("(Intercept)" = "InterceptP8", "Treatment100"="Treatment100P8", "Treatment75"="Treatment75P8", "Treatment50"="Treatment50P8", "Treatment25"="Treatment25P8", "Treatment10"="Treatment10P8", "Treatment5"="Treatment5P8")) G_pVal9 <- rename(G_pVal9, replace = c("(Intercept)" = "InterceptP9", "Treatment100"="Treatment100P9", "Treatment75"="Treatment75P9", "Treatment50"="Treatment50P9", "Treatment25"="Treatment25P9", "Treatment10"="Treatment10P9", "Treatment5"="Treatment5P9")) #Combining p-values from each pairwise concentration comparison and the negative control All_Pvalues<-c(G_pVal1,G_pVal2,G_pVal3,G_pVal4,G_pVal5,G_pVal6,G_pVal7,G_pVal8,G_pVal9) All_Pvalues #Using a False Discovery Rate (FDR) correction to correct for multiple post-hoc comparisons #BH adjusted p-values for plates 1-9 BHAdjustedPValues<-p.adjust(All_Pvalues, method = "BH", n = length(All_Pvalues)) BHAdjustedPValues #Selecting p-values that are above 0.05 to identify lowest concentration (%) that is not significantly different from the negative control per plate nonsignifPvalue<-0.05 idx <- which(BHAdjustedPValues > nonsignifPvalue) # row numbers BHAdjustedPValues[idx] #Lowest 6gingerol concentration (MIC) per plate: #P1 = 25 ug/mL #P2 = 50 ug/mL #P3 = 50 ug/mL #P4 = 25 ug/mL #P5 = 25 ug/mL #P6 = 10 ug/mL (All concentrations exhibited lower Bsal zoospore viability than the negative control, indicated by a significant p-value) #P7 = 25 ug/mL #P8 = 25 ug/mL #P9 = 25 ug/mL #Overall lowest 6gingerol concentration (MIC) for plates 1-9 = 25 ug/mL #GLS models & boxplots for plate dates 1-3 (Trials 1-3) #Plate Date 1 #Boxplot G_PD1_boxplot<-ggplot(GPlate_Date1, aes(x=Treatment, y=Conc,color=Treatment,palette="jco")) +geom_boxplot()+theme(legend.position = "none")+ ggtitle("Plate Date 1") G_PD1_boxplot #Model G_PD1_gls<-gls(Conc~Treatment,weights=varIdent(form=~1|Treatment),data=GPlate_Date1) summary(G_PD1_gls) pVal_PD1<-summary(G_PD1_gls)$tTable[,4] pVal_PD1 plot(G_PD1_gls) qqnorm(resid(G_PD1_gls, type='pearson')) qqline(resid(G_PD1_gls, type='pearson')) #Plate Date 2 #Boxplot G_PD2_boxplot<-ggplot(GPlate_Date2, aes(x=Treatment, y=Conc,color=Treatment,palette="jco")) +geom_boxplot()+theme(legend.position = "none")+ ggtitle("Plate Date 2") G_PD2_boxplot #Model G_PD2_gls<-gls(Conc~Treatment,weights=varIdent(form=~1|Treatment),data=GPlate_Date2) summary(G_PD2_gls) pVal_PD2<-summary(G_PD2_gls)$tTable[,4] pVal_PD2 plot(G_PD2_gls) qqnorm(resid(G_PD2_gls, type='pearson')) qqline(resid(G_PD2_gls, type='pearson')) #Plate Date 3 #Boxplot G_PD3_boxplot<-ggplot(GPlate_Date2, aes(x=Treatment, y=Conc,color=Treatment,palette="jco")) +geom_boxplot()+theme(legend.position = "none")+ ggtitle("Plate Date 3") G_PD3_boxplot #Model G_PD3_gls<-gls(Conc~Treatment,weights=varIdent(form=~1|Treatment),data=GPlate_Date3) summary(G_PD3_gls) pVal_PD3<-summary(G_PD3_gls)$tTable[,4] pVal_PD3 plot(G_PD3_gls) qqnorm(resid(G_PD3_gls, type='pearson')) qqline(resid(G_PD3_gls, type='pearson')) #Renaming treatments for plate dates (Trials 1-3) pVal_PD1 <- rename(pVal_PD1, replace = c("(Intercept)" = "InterceptPD1", "Treatment100"="Treatment100PD1", "Treatment75"="Treatment75PD1", "Treatment50"="Treatment50PD1", "Treatment25"="Treatment25PD1", "Treatment10"="Treatment10PD1", "Treatment5"="Treatment5PD1")) pVal_PD2 <- rename(pVal_PD2, replace = c("(Intercept)" = "InterceptPD2", "Treatment100"="Treatment100PD2", "Treatment75"="Treatment75PD2", "Treatment50"="Treatment50PD2", "Treatment25"="Treatment25PD2", "Treatment10"="Treatment10PD2", "Treatment5"="Treatment5PD2")) pVal_PD3 <- rename(pVal_PD3, replace = c("(Intercept)" = "InterceptPD3", "Treatment100"="Treatment100PD3", "Treatment75"="Treatment75PD3", "Treatment50"="Treatment50PD3", "Treatment25"="Treatment25PD3", "Treatment10"="Treatment10PD3", "Treatment5"="Treatment5PD3")) #Combining p-values from each pairwise treatment comparison and the negative control All_Pvalues_PD<-c(pVal_PD1,pVal_PD2, pVal_PD3) All_Pvalues_PD #Using a False Discovery Rate (FDR) correction to correct for multiple post-hoc comparisons #BH adjusted p-values for plate dates 1-3 BHAdjustedValues_PD<-p.adjust(All_Pvalues_PD, method = "BH", n = length(All_Pvalues_PD)) BHAdjustedValues_PD #Selecting p-values that are above 0.05 to identify lowest concentration (%) that is not significantly different from the negative control per plate date nonsignifPvalue<-0.05 idx <- which(BHAdjustedValues_PD > nonsignifPvalue) # row numbers BHAdjustedValues_PD[idx] #Lowest 6gingerol concentration (MIC) per plate date: #PD1 = 25 ug/mL #PD2 = 25 ug/mL #PD3 = 25 ug/mL #Overall lowest 6gingerol concentration (MIC) for plate dates 1-3 = 25 ug/mL #Bootstrap analysis to test that our assumptions of normality did not bias our estimated lowest 6-gingerol concentration. #Generating 1000 re-sampled datasets, refitting the GLS model (Treatment+Plate_Date+Plate), and extracting the bootstrapped sampling distributions of our model coefficients. num_boot = 1000 coefs = array(NA, dim=c(num_boot, length(coef(G_gls_mod_1)))) n = nrow(Gingerol) # Start of bootstrap for(i in 1:num_boot){ # Resample data new_dat = Gingerol[sample(1:n, n, replace=T), ] # Refit the model from above tfit = gls(Conc~Treatment+Plate+Plate_Date,weights=varIdent(form=~1|Treatment),method = "ML",data=new_dat) # Extract coefficients coefs[i, ] = coef(tfit) } # Assigning coeficient names colnames(coefs) = names(coef(G_gls_mod_1)) # Check if sampling distributions are normal pairs.panels(coefs) # Compare 95% confidence intervals from the bootrap to the normal assumption model apply(coefs, 2, quantile, c(0.025, 0.5, 0.975)) summary(G_gls_mod_1) #Identified the lowest concentration where the 95% confidence interval of the bootstrapped sampling distribution overlapped zero (i.e., the lowest concentration that was not significantly different from the negative control) #Lowest 6gingerol concentration (MIC) from bootstrap analysis = 25 ug/mL