############################################################################### ##load packages################################################################ ############################################################################### library(data.table) library(cowplot) library(lubridate) library(MASS) library(foreach) library(lme4) #set working directory: setwd("~/Dropbox/Helen-Priscilla/") ############################################################################### ##Figure 1##################################################################### ############################################################################### #read in csv for freezer temperature data as a data table temp<-fread("temp_data_usable.csv") #format time to start at 0 temp[,time:=dmy_hm(Time)] temp[,time.adjust:=as.numeric(time-min(time))/60] #graph and save as pdf pdf("C:/Users/Helen/Dropbox/Helen-Priscilla/almost_final/fig_1.pdf", height=3.15, width=3.15) ggplot(temp, aes(x=time.adjust, y=temperature, color=location))+geom_line(size=1.2)+labs(x="Minutes frozen", y="Temperature (°C)", color="Thermometer\nlocation")+lims(y=c(-6,0))+theme(legend.position=c(.5, .75)) dev.off() ############################################################################### ##Figure 2##################################################################### ############################################################################### #read in csv for Morven data as a data table dat<-fread("morven_master_sheet_V9.csv") #exclude outdoor cornmeal cages dat<-dat[cage%in%c(1:6) | cage%in%c("A", "B")] dat<-dat[collection_date != "10/15/18"] dat<-dat[collection_date !="10/30/18"] #drop F0 7/24 indoor point because outdoor data were insufficient to produce LT50 dat<-dat[!(collection_date == "7/24/18" & generation == "F0")] #drop F0 12/10 indoor point because outdoor data were insufficient to produce LT50 dat<-dat[collection_date != "12/10/18"] #calculate percent alive dat[,percent_alive:=ifelse(is.na(n_dead_after), n_alive_after/(n_start-n_dead_before), ((n_start-n_dead_before)-n_dead_after)/(n_start-n_dead_before))] #calculate starting number dat[,n:=n_start-n_dead_before] #do some renaming dat[cage==1| cage==3|cage==5| cage=="A", founder_cage:="A"] dat[cage==2| cage==4|cage==6| cage=="B", founder_cage:="B"] dat[cage%in%c(1:6), environment:="outdoor"] dat[cage%in%c("A", "B"), environment:="indoor"] dat[cage%in%c(1:6), location:="outdoor"] dat[cage%in%c("A", "B"), location:="indoor"] dat[cage%in%c(1:6), food:="fruit"] dat[cage%in%c("A", "B"), food:="trays"] #reorder collection dates into an ordered factor by date dat[,collection_date:=factor(dat$collection_date, levels=c("6/26/2018", "7/24/2018", "8/21/2018", "9/18/2018", "10/16/2018", "10/30/2018", "11/7/2018", "11/21/2018", "11/30/2018", "12/10/2018" ))] #loop through all data subsets and calculate an LT50 for each; return table o<-foreach(gen=c("F0", "F1", "F2"), .errorhandling="remove")%do%{ print(gen) dat1<-dat[generation==gen] #print(dat1) d<-foreach(cd=unique(dat1$collection_date), .errorhandling="remove")%do%{ print(cd) dat2<-dat1[collection_date==cd] #print(dat2) l<-foreach(loc=unique(dat2$environment), .errorhandling="remove")%do%{ print(loc) dat3<-dat2[environment==loc] #print(dat3) mod<-glm(percent_alive~minutes_frozen, data=dat3, weights=n, family="binomial") return(data.table(generation=gen, collection_date=cd, location=loc, lt50=dose.p(mod, cf=1:2, p=0.5)[1], lt50se=attr(dose.p(mod, cf=1:2, p=0.5),"SE")[1,1])) } return(rbindlist(l)) } return(rbindlist(d)) } o<-rbindlist(o) #change date to mdy o[,date:=mdy(collection_date)] #get rid of too-large error bar o[lt50se>100,lt50se:=0] #plot LT50s for each generation a<-ggplot(o[generation == "F0"], aes(x=date, y=lt50, ymin=45, ymax=300, color=location))+geom_errorbar(data=o[generation == "F0"], width=4, aes(x=date, ymin=(lt50-lt50se), ymax=(lt50+lt50se)))+geom_line(size=.3) + labs(x="", y="LT50 (minutes)")+theme(axis.text=element_text(size=12), axis.text.x=element_blank(), axis.title.x=element_blank(), axis.title=element_text(size=14))+scale_x_date(limits=as.Date(c('2018-06-22', '2018-12-14'))) li<-ggplot(o[generation == "F1"], aes(x=date, y=lt50, ymin=45, ymax=300, color=location))+geom_errorbar(data=o[generation == "F1"], width=4, aes(x=date, ymin=(lt50-lt50se), ymax=(lt50+lt50se)))+geom_line(size=.3) + labs(x="", y="LT50 (minutes)")+theme(axis.text.x=element_blank(), axis.title.x=element_blank(), axis.title=element_text(size=14), axis.text=element_text(size=12))+scale_x_date(limits=as.Date(c('2018-06-22', '2018-12-14'))) pi<-ggplot(o[generation == "F2"], aes(x=date, y=lt50, ymin=45, ymax=300, color=location))+geom_errorbar(data=o[generation == "F2"], width=4, aes(x=date, ymin=(lt50-lt50se), ymax=(lt50+lt50se)))+geom_line(size=.3) + labs(x="Collection date", y="LT50 (minutes)")+theme(axis.title=element_text(size=14), axis.text=element_text(size=12))+scale_x_date(limits=as.Date(c('2018-06-22', '2018-12-14'))) #graph and save as pdf pdf("C:/Users/Helen/Dropbox/Helen-Priscilla/almost_final/fig2.pdf", height=6, width=3.15) plot_grid(a+theme(legend.position="none"), li+theme(legend.position="none"), pi+theme(legend.position="none"), nrow=3, ncol=1, align="v", labels=c("A", "B", "C"), rel_heights=c(.31,.31,.38)) dev.off() ############################################################################### ##Figure 2 and Table 2 stats################################################### ############################################################################### sumtable<-foreach(gen=c("F0", "F1", "F2"), .combine="rbind", .errorhandling="remove") %do% { po<-foreach(d=unique(dat$collection_date), .combine="rbind", .errorhandling="remove") %do% { e<-summary(glmer(percent_alive~location+minutes_frozen+(location|cage), weights=n, family="binomial", data=dat[generation==gen & collection_date==d])) return(data.table(generation=gen, collection.date=d, estimate=e$coefficients[2,1], p=e$coefficients[2,4] )) } return(po) } sumtable[, sig:=p<(.05/nrow(sumtable))] sumtable[, dir:=ifelse(sign(estimate)==1,"outdoor greater", "indoor greater")] ############################################################################### ##Figure 3##################################################################### ############################################################################### #use code from Figure 2 and start with data table "o" #create wide table to calculate changes in lt50 p<-dcast(o, generation+collection_date+date~location, value.var=c("lt50", "lt50se")) p[,diff:=lt50_outdoor-lt50_indoor] #read in weather data w<-fread("weather_data.csv") #convert to Celsius w[,tempC:=(Avg_temp-32)*(5/9)] #format dates w[,date:=mdy(date)] w[,pretty_date:=paste(month(date, label=T), day(date), sep=" ")] w[,pretty_date:=factor(pretty_date, levels=c("Jun 26", "Aug 21", "Sep 18", "Oct 16", "Nov 7", "Nov 21"))] #merge seasonal data and weather data p<-merge(p, w, by="date") #make a graph wea<-ggplot(p[generation=="F0"&!is.na(diff)], aes(x=tempC, y=diff, color=as.factor(pretty_date)))+geom_point(size=4)+ scale_color_manual(values=c("orange1", "orangered", "red3", "mediumvioletred", "darkmagenta", "purple4"))+ geom_abline(slope=-7.658, intercept=62.514, linetype="dashed", size = 1)+labs(x="Avg. field temperature (°C)", y="Outdoor - indoor LT50 (min.)", color = "Date")+theme(axis.text=element_text(size=12), axis.title=element_text(size=14), legend.position=c(.55,.72)) #save the graph as a pdf pdf("C:/Users/Helen/Dropbox/Helen-Priscilla/almost_final/fig3.pdf", width=3.15, height=3.15) wea dev.off() ############################################################################### ##Figure 3 stats############################################################### ############################################################################### summary(lm(diff~tempC, data=p[generation=="F0"&collection_date!="11/21/2018"])) ############################################################################### ##Figure 4##################################################################### ############################################################################### #read in csv of cornmeal-molasses fed cages as data table new<-fread("Morv_F0_7to11_102918.csv") #calculate percent alive new[,percent_alive:=ifelse(is.na(n_additional_dead_after), n_alive_after/(n_start-n_dead_before), ((n_start-n_dead_before)-n_additional_dead_after)/(n_start-n_dead_before))] #calculate starting number new[,n:=n_start-n_dead_before] #do some renaming new[cage%in%c(7:11), environment:="Outdoor cornmeal-molasses-fed"] new[cage%in%c("A", "B"), environment:="Indoor cornmeal-molasses-fed"] #pool cages of same environment & summarize weighted averages across cages new.sum<-new[,.(mean.percent=weighted.mean(x=percent_alive, w=n)), .(environment, minutes_frozen)] #make a graph new_oct<-ggplot(data=new.sum, aes(x=minutes_frozen, y=mean.percent*100, color=environment))+geom_point(size=1.8)+geom_line(size=1)+labs(x="Minutes frozen", y="% survival", color="Cage Location and Food")+theme(axis.title=element_text(size=14), axis.text=element_text(size=12))+scale_color_manual(values=c("#F8766D", "#7CAE00")) #read in csv of fruit fed cages as data table old<-fread("morven_F0_1to6_103018.csv") #calculate percent alive old[,percent_alive:=ifelse(is.na(n_additional_dead_after), n_alive_after/(n_start-n_dead_before), ((n_start-n_dead_before)-n_additional_dead_after)/(n_start-n_dead_before))] #calculate starting number old[,n:=n_start-n_dead_before] #do some renaming old[cage%in%c(1:6), environment:="Outdoor fruit-fed"] old[cage%in%c("A", "B"), environment:="Indoor cornmeal-molasses-fed"] #pool cages of same environment & summarize weighted averages across cages old.sum<-old[,.(mean.percent=weighted.mean(x=percent_alive, w=n)), .(environment, minutes_frozen)] #make a graph old_oct<-ggplot(data=old.sum, aes(x=minutes_frozen, y=mean.percent*100, color=environment))+geom_point(size=1.8)+geom_line(size=1)+labs(x="Minutes frozen", y="% survival", color="Cage Location and Food")+theme(axis.text=element_text(size=12), axis.title=element_text(size=14)) #save combined graph as pdf pdf("C:/Users/Helen/Dropbox/Helen-Priscilla/almost_final/fig4.pdf", height=4, width=3.15) plot_grid(old_oct+theme(legend.position="none"), new_oct+theme(legend.position="none"), align="v", ncol=1, labels = c("A", "B")) dev.off() ############################################################################### ##Figure 4 and Table 3 stats################################################### ############################################################################### #stats for cages 1-6, A-B summary(aov(glm(percent_alive~minutes_frozen+environment+cage, weights=n, family="binomial", data=old))) #stats for cages 7-11, A-B summary(aov(glm(percent_alive~minutes_frozen+environment+cage, weights=n, family="binomial", data=new))) ############################################################################### ##Figure 5##################################################################### ############################################################################### #read in csv for nutrition/antimicrobial assay as a data table fd<-fread("food_freeze_040119.csv") #calculate percent alive fd[,percent_alive:=ifelse(is.na(n_dead_after), n_alive_after/(n_start-n_dead_before), ((n_start-n_dead_before)-n_dead_after)/(n_start-n_dead_before))] #calculate starting number fd[,n:=n_start-n_dead_before] #pool cages & summarize weighted averages across cages fd.sum<-fd[,.(mean.percent=weighted.mean(x=percent_alive, w=n)), .(substrate, antimicrobials, minutes_frozen)] #do some renaming fd.sum[substrate=="ban", sub:="banana"] fd.sum[substrate=="cm", sub:="cornmeal-molasses"] fd.sum[antimicrobials=="yes", ant:="+antimicrobials"] fd.sum[antimicrobials=="no", ant:="-antimicrobials"] #make a graph food<-ggplot(data=fd.sum, aes(x=minutes_frozen, y=mean.percent*100, color=sub, linetype=ant))+geom_point(size=1.8) + geom_line(size=1)+labs(x="Minutes frozen", y="% survival", color="Developmental food", linetype=" ")+theme(axis.text=element_text(size=12), axis.title=element_text(size=14))+scale_color_manual(values=c("deepskyblue3", "deeppink3")) #read in csv for yeast assay as a data table ye<-fread("yeast_freeze_32419.csv") #calculate percent alive ye[,percent_alive:=ifelse(is.na(n_dead_after), n_alive_after/(n_start-n_dead_before), ((n_start-n_dead_before)-n_dead_after)/(n_start-n_dead_before))] #calculate starting number ye[,n:=n_start-n_dead_before] #pool cages & summarize weighted averages across cages ye.sum<-ye[,.(mean.percent=weighted.mean(x=percent_alive, w=n)), .(yeast_added, minutes_frozen)] #do some renaming ye.sum[yeast_added=="yes", yeast:="+yeast"] ye.sum[yeast_added=="no", yeast:="-yeast"] #make a graph yeast<-ggplot(data=ye.sum, aes(x=minutes_frozen, y=mean.percent*100, color=yeast))+geom_point(size=1.8) + geom_line(size=1)+labs(x="Minutes frozen", y="% survival", color=" ")+theme(axis.text=element_text(size=12), axis.title=element_text(size=14))+scale_color_manual(values=c("chocolate4", "orange3")) #read in csv for aging assay as a data table ag<-fread("aging_flies_PC_freeze_111318.csv") #calculate percent alive ag[,percent_alive:=ifelse(is.na(n_dead_after), n_alive_after/(n_start-n_dead_before), ((n_start-n_dead_before)-n_dead_after)/(n_start-n_dead_before))] #calculate starting number ag[,n:=(n_start-n_dead_before)] #do some renaming ag[lay_date=="16-Oct", time:="4 days"] ag[lay_date=="18-Sep", time:="32 days"] ag[lay_date=="2-Oct", time:="18 days"] #make time a factor ag[,time:=factor(ag$time, levels=c("4 days", "18 days", "32 days"))] #pool cages & summarize weighted averages across cages ag.sum<-ag[,.(mean.alive=weighted.mean(x=percent_alive, w=n)), .(time, minutes_frozen)] #make a graph age<-ggplot(data=ag.sum, aes(x=minutes_frozen, y=mean.alive*100, color=time, group=time))+geom_point(size=1.8)+geom_line(size=1)+labs(x="Minutes frozen", y="% survival", color="Age at precool")+theme(axis.text=element_text(size=12), axis.title=element_text(size=14))+scale_color_manual(values=c("#C77CFF", "dodgerblue4", "violetred4")) #read in csv for density assay as a data table de<-fread("C:/Users/Helen/Dropbox/Helen-Priscilla/density_flies.csv") #calculate percent alive de[,percent_alive:=ifelse(is.na(n_dead_after), n_alive_after/(n_start-n_dead_before), ((n_start-n_dead_before)-n_dead_after)/(n_start-n_dead_before))] #calculate starting number de[,n:=n_start-n_dead_before] #exclude density with insufficient data de<-de[Density!="300 eggs/mL"] #make density a factor de[, Density:=factor(de$Density, levels=c("120 eggs/mL", "40 eggs/mL", "5 eggs/mL"))] #pool cages & summarize weighted averages across cages de.sum<-de[,.(mean.alive=weighted.mean(x=percent_alive, w=n)), .(minutes_frozen, Density)] #make a graph density<-ggplot(de.sum, aes(x=minutes_frozen, y=mean.alive*100, color=Density, group=interaction(Density)))+geom_line(size=1)+geom_point(size=1.8)+labs(x="Minutes frozen", y="% survival")+theme(axis.text=element_text(size=12), axis.title=element_text(size=14))+scale_color_manual(values=c("orangered1", "seagreen3", "royalblue")) #compile graphs fig_five <-plot_grid(food+theme(legend.position="none"), yeast+theme(legend.position="none"), age+theme(legend.position="none"), density+theme(legend.position="none"), ncol=1, nrow=4, labels=c("A", "B", "C", "D"), align="v") #save graph as pdf pdf("C:/Users/Helen/Dropbox/Helen-Priscilla/almost_final/fig5.pdf", height=8, width=3.5) fig_five dev.off() ############################################################################### ##Figure 5 and Table 3 stats################################################### ############################################################################### ##container experiment #read in csv for bottles as a data table b=fread("Morven_F1_bottles_101618.csv") #read in csv for vials as a data table v=fread("Morven_F1_vials_101518.csv") #do some renaming b[,container:="bottle"] v[,container:="vial"] #combine a<-rbind(b,v) #more renaming a[cage%in%c(1:11), location:="outdoor"] a[cage%in%c("A", "B"), location:="indoor"] #calculate percent alive a[,percent_alive:=ifelse(is.na(n_additional_dead_after), n_alive_after/(n_start-n_dead_before), ((n_start-n_dead_before)-n_additional_dead_after)/(n_start-n_dead_before))] #calculate starting number a[,n:=n_start-n_dead_before] #container stats summary(aov(glm(percent_alive~minutes_frozen+container+location+cage, data=a, weights=n, family = "binomial"))) ##nutrition and antimicrobials summary(aov(glm(percent_alive~minutes_frozen+substrate+antimicrobials+cage, weights=n, family="binomial", data=fd))) ##yeast summary(aov(glm(percent_alive~minutes_frozen+yeast_added+cage, weights=n, family="binomial", data=ye[minutes_frozen>250]))) ##aging summary(aov(glm(percent_alive~minutes_frozen+time+cage, weights=n, family="binomial", data=ag))) ##density summary(aov(glm(percent_alive~minutes_frozen+Density+cage, weights=n, family="binomial", data=de[minutes_frozen>200])))