#Before running any of the code below, download all the text files #from http://depts.washington.edu/kerrpost/Public/RifRampCode without #changing any of their names, and place them in a new folder. #Write the full path to this folder here: wd<-"C:/[Write in the relevant path here]" #After executing the above assignment, in your new folder, create #the following subfolders: # Evolved_PDFS # Single_PDFS # ANC_PDFS # G13_PDFS # G18_PDFS # M358_PDFS # M367_PDFS # S497_PDFS # S1236_PDFS #****************************************************************** #Figure 1a: Rifampicin concentrations for the three main treatments #------------------------------------------------------------------ #User-specified parameters: Scol<-"tomato4" #color for Sudden treatment Mcol<-"tomato3" #color for Moderate treatment Gcol<-"tomato1" #color for Gradual treatment Spt<-1.6 #point size for Sudden treatment Mpt<-1.2 #point size for Moderate treatment Gpt<-.7 #point size for Gradual treatment ltck<-2 #line thickness legsize<-1.5 #size of text in legend axissize<-1.2 #size of text on axes labsize<-1.5 #size of text on axis labels #Read in data setwd(wd) RifCon.df<-read.table("RifConc.txt", header=TRUE) transfer<-RifCon.df[RifCon.df$treatment=="Sudden",]$transfer Sconc<-log10(RifCon.df[RifCon.df$treatment=="Sudden",] $rif.conc[2:length(transfer)]) Mconc<-log10(RifCon.df[RifCon.df$treatment=="Moderate",] $rif.conc[2:length(transfer)]) Gconc<-log10(RifCon.df[RifCon.df$treatment=="Gradual",] $rif.conc[2:length(transfer)]) val0<-(-8) Sconc<-c(val0,Sconc) Mconc<-c(val0,Mconc) Gconc<-c(val0,Gconc) #plot the rif concentrations of the three main treatments par(mar=c(5,5,4,2)) plot(transfer,Sconc,col="white",cex.lab=labsize, ylab=expression(paste("rifampicin concentration (", mu, "g/mL)",sep="")),axes=FALSE) box() axis(1,cex.axis=axissize) axis(2,at=c(val0,-6,-4,-2,0,2),labels=c(0,expression(10^-6), expression(10^-4),expression(10^-2),1,100),las=1,cex.axis=axissize) lines(transfer,Sconc,col=Scol,lwd=ltck) points(transfer,Sconc,col=Scol,cex=Spt,pch=19) lines(transfer,Mconc,col=Mcol,lwd=ltck) points(transfer,Mconc,col="white",cex=Mpt,pch=19) points(transfer,Mconc,col=Mcol,cex=Mpt,pch=1) lines(transfer,Gconc,col=Gcol,lwd=ltck) points(transfer,Gconc,col=Gcol,cex=Gpt,pch=19) #place break in the y axis bt<-0.1 rect(-3,-7-bt,3,-7+bt,col="white",border=FALSE) lines(c(-3,-1.1),c(-7-bt,-7-bt)) lines(c(-3,-1.1),c(-7+bt,-7+bt)) lines(c(-0.25,.75),c(-7-bt,-7-bt)) lines(c(-0.25,.75),c(-7+bt,-7+bt)) #place legend on graph lowx<-18 highx<-30 lowy<-(-6) highy<-(-3) rect(lowx-1,lowy,highx+1,highy,col="white") lines(c(lowx,lowx+2),c(highy-(.25*(highy-lowy)),highy-(.25*(highy-lowy))), col=Scol,lwd=ltck) points(lowx+1,highy-(.25*(highy-lowy)),col=Scol,cex=Spt,pch=19) lines(c(lowx,lowx+2),c(highy-(.5*(highy-lowy)),highy-(.5*(highy-lowy))), col=Mcol,lwd=ltck) points(lowx+1,highy-(.5*(highy-lowy)),col="white",cex=Mpt,pch=19) points(lowx+1,highy-(.5*(highy-lowy)),col=Mcol,cex=Mpt,pch=1) lines(c(lowx,lowx+2),c(highy-(.75*(highy-lowy)),highy-(.75*(highy-lowy))), col=Gcol,lwd=ltck) points(lowx+1,highy-(.75*(highy-lowy)),col=Gcol,cex=Gpt,pch=19) text(lowx+3,highy-(.25*(highy-lowy)),"Sudden",col=Scol, pos=4, cex=legsize) text(lowx+3,highy-(.5*(highy-lowy)),"Moderate",col=Mcol, pos=4, cex=legsize) text(lowx+3,highy-(.75*(highy-lowy)),"Gradual",col=Gcol, pos=4, cex=legsize) #End of Fig 1a----------------------------------------------------- #****************************************************************** #****************************************************************** #****************************************************************** #Figure 1b: Survivorship fraction in each treatment #------------------------------------------------------------------ #User-specified parameters: Scol<-"tomato4" #color for Rapid treatment Mcol<-"tomato3" #color for Moderate treatment Gcol<-"tomato1" #color for Gradual treatment axissize<-1.3 #size of text on axes labsize<-1.5 #size of text on axis labels #read and process the data: setwd(wd) Survive.df<-read.table("Survivorship.txt", header=TRUE) Sinit<-Survive.df[Survive.df$treatment=="Sudden",]$initial.N.pops Sfin<-Survive.df[Survive.df$treatment=="Sudden",]$final.N.pops Slost<-Survive.df[Survive.df$treatment=="Sudden",]$lost.N.pops Minit<-Survive.df[Survive.df$treatment=="Moderate",]$initial.N.pops Mfin<-Survive.df[Survive.df$treatment=="Moderate",]$final.N.pops Mlost<-Survive.df[Survive.df$treatment=="Moderate",]$lost.N.pops Ginit<-Survive.df[Survive.df$treatment=="Gradual",]$initial.N.pops Gfin<-Survive.df[Survive.df$treatment=="Gradual",]$final.N.pops Glost<-Survive.df[Survive.df$treatment=="Gradual",]$lost.N.pops Ssur<-((Sfin+Slost)/Sinit)*100 Msur<-((Mfin+Mlost)/Minit)*100 Gsur<-((Gfin+Glost)/Ginit)*100 #plot the survivorships par(mar=c(5,5,4,2)) lab<-c("Sudden","Moderate","Gradual") xcoor<-barplot(c(Ssur,Msur,Gsur),names.arg=lab, col=c(Scol,Mcol,Gcol), ylab="percent of populations surviving",xlab="treatment", cex.names=axissize,cex.lab=labsize,cex.axis=axissize,ylim=c(0,100),las=1) del<-5 text(xcoor[1,1],Ssur+del,paste("(",Sinit,")",sep=""),col=Scol, cex=axissize) text(xcoor[2,1],Msur+del,paste("(",Minit,")",sep=""),col=Mcol, cex=axissize) text(xcoor[3,1],Gsur+del,paste("(",Ginit,")",sep=""),col=Gcol, cex=axissize) #End of Fig 1b----------------------------------------------------- #****************************************************************** #****************************************************************** #Data Analysis for Fig 1b~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Set working directory setwd(wd) #Read in survivorship data Survive.df<-read.table("Survivorship.txt", header=TRUE) Sinit<-Survive.df[Survive.df$treatment=="Sudden",]$initial.N.pops Sfin<-Survive.df[Survive.df$treatment=="Sudden",]$final.N.pops Slost<-Survive.df[Survive.df$treatment=="Sudden",]$lost.N.pops Minit<-Survive.df[Survive.df$treatment=="Moderate",]$initial.N.pops Mfin<-Survive.df[Survive.df$treatment=="Moderate",]$final.N.pops Mlost<-Survive.df[Survive.df$treatment=="Moderate",]$lost.N.pops Ginit<-Survive.df[Survive.df$treatment=="Gradual",]$initial.N.pops Gfin<-Survive.df[Survive.df$treatment=="Gradual",]$final.N.pops Glost<-Survive.df[Survive.df$treatment=="Gradual",]$lost.N.pops #Create survivorship by treatment matrix Pop.data<-c(Sinit-Sfin-Slost,Minit-Mfin-Mlost,Ginit-Gfin-Glost,Sfin+Slost,Mfin+Mlost,Gfin+Glost) dim(Pop.data)<-c(3,2) #Pearson's chi-squared test chisq.test(Pop.data) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #****************************************************************** #Figure 1c: Growth rates in each treatment (RUN FIG 1B FIRST) #------------------------------------------------------------------ #INSTRUCTIONS: #The data file that this program uses has time points in the rows #(times listed in the first column) and each column lists the #absorbances for each time point for a single replicate of a given #evolved population (the pop/rep name heads the column). All growth #assays were performed at the highest concentration of rif. A separate #file contains all the names of the strains. #SET AND RESET WORKING DIRECTORY (AFTER READING IN SOME DATA) setwd(wd) spec.df<-read.table("Evolved_RifRamp.txt", header=TRUE) names.df<-read.table("PopulationNames.txt", header=TRUE) wd2<-paste(wd,"/Evolved_PDFS",sep="") setwd(wd2) names<-names.df$strain types<-names.df$type gens<-length(names) #PARAMETERS TO CHANGE reps<-4 #number of replicates (must be even) window<-5 #number of time points in sliding window (must be odd) period<-30 #number of minutes between OD readings lb<-FALSE #whether the log of the OD values are used #FUNCTION: posify #This function changes all non-positive values in a #vector to the smallest positive value in the vector #or to a set positive value if all entries are #non-positive. #--------------------------------------------------- posify<-function(x,emer.value=0.0001) { y<-numeric(length(x)) smallest.strict.pos<-max(x) if(smallest.strict.pos<=0) { y<-rep(emer.value,length(x)) } if(smallest.strict.pos>0) { for(i in 1:length(x)) { if(x[i]>0 & x[i](maxval-meanval), which.min(x),which.max(x)) if(ind==1) { y<-x[2:length(x)] return(y) } if(ind==(length(x))) { y<-x[1:(length(x)-1)] return(y) } y<-c(x[1:(ind-1)],x[(ind+1):length(x)]) return(y) } #--------------------------------------------------- #FUNCTION: logornot #This function either gives the log of the argument #or just returns the argument, depending on the second #(boolean) argument. #--------------------------------------------------- logornot<-function(x,logbool=TRUE) { ifelse(logbool,y<-log(x),y<-x) y } #--------------------------------------------------- #Create the fitness arrays fit1<-array(data = numeric(gens*reps), dim = c(gens,reps)) fit2<-array(data = numeric(gens*reps), dim = c(gens,reps)) fit3<-array(data = numeric(gens*reps), dim = c(gens,reps)) time1<-array(data = numeric(gens*reps), dim = c(gens,reps)) time2<-array(data = numeric(gens*reps), dim = c(gens,reps)) time3<-array(data = numeric(gens*reps), dim = c(gens,reps)) FIT<-array(data = numeric(gens*reps), dim = c(gens,reps)) NAMES<-array(data = numeric(gens*reps), dim = c(gens,reps)) #Calculate maximal slopes and plot graphs to check slopes for(x in 1:gens) { for(z in 1:reps) { #Get column name (for strain/replicate) ifelse(z==1, geno<-paste(names[x],"a",sep=""), ifelse(z==2, geno<-paste(names[x],"b",sep=""), ifelse(z==3, geno<-paste(names[x],"c",sep=""), geno<-paste(names[x],"d",sep="")))) NAMES[x,z]=geno #Get times and absorbance data for focal strain/replicate t<-spec.df$time od<-spec.df[,geno] od<-posify(od) #Calculate slopes across the OD trace times<-numeric(length(od)-window+1) slopes<-numeric(length(od)-window+1) for(i in (1+((window-1)/2)):(length(od)-((window-1)/2))){ slopes[i-(window-1)/2]<-lsfit(t[(i-(window-1)/2):(i+(window-1)/2)],logornot(od[(i-(window-1)/2):(i+(window-1)/2)],lb))$coefficients[2] times[i-(window-1)/2]<-t[i] } #Find the three highest slopes across the OD trace and average fit1[x,z]<-max(slopes) tind<-which.max(slopes) time1[x,z]<-times[tind] slopes[tind]<-(min(slopes)-1) fit2[x,z]<-max(slopes) tind<-which.max(slopes) time2[x,z]<-times[tind] slopes[tind]<-(min(slopes)-1) fit3[x,z]<-max(slopes) tind<-which.max(slopes) time3[x,z]<-times[tind] slopes[tind]<-(min(slopes)-1) FIT[x,z]=(fit1[x,z]+fit2[x,z]+fit3[x,z])/3.0 #Plot the OD traces if(z==1) { dev.set(3) pdf(paste(names[x],".pdf", sep="")) par(mfrow=c(2,reps/2)) } plot(times,logornot(od[(1+((window-1)/2)): (length(od)-((window-1)/2))],lb),xlab="time (min)", ylab="OD",main=paste(names[x],":",z,"\n","slope=", FIT[x,z]),col="white") points(times,logornot(od[(1+((window-1)/2)): (length(od)-((window-1)/2))],lb),pch=19) points(time1[x,z],logornot(od[(time1[x,z]/period)+1],lb), col="red",cex=2) points(time1[x,z],logornot(od[(time1[x,z]/period)+1],lb), col="red",cex=3) lines(times,fit1[x,z]*(times-time1[x,z])+ logornot(od[(time1[x,z]/period)+1],lb),lwd=2,col="red") points(time2[x,z],logornot(od[(time2[x,z]/period)+1],lb), col="red",cex=2) points(time2[x,z],logornot(od[(time2[x,z]/period)+1],lb), col="red",cex=3) lines(times,fit2[x,z]*(times-time2[x,z])+ logornot(od[(time2[x,z]/period)+1],lb),lwd=2,col="red") points(time3[x,z],logornot(od[(time3[x,z]/period)+1],lb), col="red",cex=2) points(time3[x,z],logornot(od[(time3[x,z]/period)+1],lb), col="red",cex=3) lines(times,fit3[x,z]*(times-time3[x,z])+ logornot(od[(time3[x,z]/period)+1],lb),lwd=2,col="red") if(z==reps) { dev.off() } } } #Create vectors for the mean and sd fitnesses fit.mean<-numeric(gens) fit.sde<-numeric(gens) names.short<-numeric(gens) #Throw out the most deviant slope and calculate mean and sd for(x in 1:gens) { fit.mean[x]=mean(toss.outlier(FIT[x,1:reps])) fit.sde[x]=(sd(toss.outlier(FIT[x,1:reps]))/sqrt(reps-1)) names.short[x]=NAMES[x,1] } #Export names of strains in order of decreasing mean fitness write(names.short[order(fit.mean,decreasing=TRUE)],file = "out.txt",ncolumns = 1,sep = "\t") write(sort(fit.mean,decreasing=TRUE),file = "out2.txt",ncolumns = 1,sep = "\t") #Collect the mean fitnesses of each treatment into vectors G.gro<-fit.mean[which(types=="Gradual")] M.gro<-fit.mean[which(types=="Moderate")] S.gro<-fit.mean[which(types=="Sudden")] Sh.gro<-fit.mean[which(types=="Shock")] #User-specified parameters: setwd(wd) Scol<-"tomato4" #color for Rapid treatment Mcol<-"tomato3" #color for Moderate treatment Gcol<-"tomato1" #color for Gradual treatment pt.cex<-1.2 #point size ltck<-2 #line thickness legsize<-1.5 #size of text in legend axissize<-1.1 #size of text on axes labsize<-1.5 #size of text on axis labels N.treatments<-3 Hist.range<-.25 Plot.spacer<-.1 xmin<-0 xmax.plot<-0.0012 ymin<-0 ymax<-(N.treatments*Hist.range+(N.treatments-1)*Plot.spacer) Nbars<-25 #Put in the zero growth strains (note this requires running Fig 1B above) S.gro.all<-c(S.gro,rep(0,Sinit-length(S.gro))) M.gro.all<-c(M.gro,rep(0,Minit-length(M.gro))) G.gro.all<-c(G.gro,rep(0,Ginit-length(G.gro))) #Calculate bar length and reserve histogram vectors xmax<-max(c(S.gro.all,M.gro.all,G.gro.all)) BarLength<-(xmax/Nbars) SHist<-rep(0,Nbars) MHist<-rep(0,Nbars) GHist<-rep(0,Nbars) #Fill the histograms for(i in 1:length(S.gro.all)) { ifelse(S.gro.all[i]==0, SHist[1]<-(SHist[1]+(1.0/(length(S.gro.all)))), SHist[ceiling(S.gro.all[i]/BarLength)]<-(SHist[ceiling(S.gro.all[i]/BarLength)]+(1.0/(length(S.gro.all))))) } for(i in 1:length(M.gro.all)) { ifelse(M.gro.all[i]==0, MHist[1]<-(MHist[1]+(1.0/(length(M.gro.all)))), MHist[ceiling(M.gro.all[i]/BarLength)]<-(MHist[ceiling(M.gro.all[i]/BarLength)]+(1.0/(length(M.gro.all))))) } for(i in 1:length(G.gro.all)) { ifelse(G.gro.all[i]==0, GHist[1]<-(GHist[1]+(1.0/(length(G.gro.all)))), GHist[ceiling(G.gro.all[i]/BarLength)]<-(GHist[ceiling(G.gro.all[i]/BarLength)]+(1.0/(length(G.gro.all))))) } #Plot the axes par(mar=c(5,5,4,2)) plot(c(xmin,xmax.plot),c(ymin,ymax),col="white",cex.lab=labsize,xlab="growth rate (AU/min)",ylab="frequency",axes=FALSE) axis(1,at=c(0,.0004,.0008,.0012),cex.axis=axissize) axis(2,at= c(0,(1/3)*Hist.range,(2/3)*Hist.range,Hist.range, Plot.spacer+Hist.range, Plot.spacer+Hist.range+(1/3)*Hist.range, Plot.spacer+Hist.range+(2/3)*Hist.range, Plot.spacer+2*Hist.range, 2*Plot.spacer+2*Hist.range, 2*Plot.spacer+2*Hist.range+(1/3)*Hist.range, 2*Plot.spacer+2*Hist.range+(2/3)*Hist.range, 2*Plot.spacer+3*Hist.range), labels=c(c(0,.05,.1,.15), c(0,.05,.1,.6), c(0,.05,.1,1)), las=1, cex.axis=axissize) #breaks in the y-axis: epi1<-0.01 epi2<-0.00175 rect(xmin-epi1,Hist.range+epi2,xmin+epi1,Hist.range+Plot.spacer-epi2,col="white",border=FALSE) rect(xmin-epi1,2*Hist.range+Plot.spacer+epi2,xmin+epi1,2*Plot.spacer+2*Hist.range-epi2,col="white",border=FALSE) #add lines and boxes: lines(c(xmin,xmax.plot),c(0,0)) lines(c(xmin,xmax.plot),c(Plot.spacer+Hist.range,Plot.spacer+Hist.range)) lines(c(xmin,xmax.plot),c(2*Plot.spacer+2*Hist.range,2*Plot.spacer+2*Hist.range)) #reassign zero bars: SHist[1]=0.15*(SHist[1]) MHist[1]=0.15*(MHist[1]/0.6) for(i in 1:Nbars) { rect(xmin+BarLength*(i-1),2*Plot.spacer+2*Hist.range,xmin+BarLength*i, 2*Plot.spacer+2*Hist.range+(SHist[i]*Hist.range)/.15,col=Scol) rect(xmin+BarLength*(i-1),Plot.spacer+Hist.range,xmin+BarLength*i, Plot.spacer+Hist.range+(MHist[i]*Hist.range)/.15,col=Mcol) rect(xmin+BarLength*(i-1),0,xmin+BarLength*i, (GHist[i]*Hist.range)/.15,col=Gcol) } #place breaks in the y axis epi1<-0.000025 epi2<-0.0075 epi3<-0.00001 rect(-2*BarLength,2*Plot.spacer+2*Hist.range+(5/6)*Hist.range-epi2,BarLength+epi1, 2*Plot.spacer+2*Hist.range+(5/6)*Hist.range+epi2,col="white",border=FALSE) lines(c(-2*BarLength,-epi1-epi3),c(2*Plot.spacer+2*Hist.range+(5/6)*Hist.range-epi2, 2*Plot.spacer+2*Hist.range+(5/6)*Hist.range-epi2)) lines(c(-2*BarLength,-epi1-epi3),c(2*Plot.spacer+2*Hist.range+(5/6)*Hist.range+epi2, 2*Plot.spacer+2*Hist.range+(5/6)*Hist.range+epi2)) lines(c(-epi3,BarLength+epi3),c(2*Plot.spacer+2*Hist.range+(5/6)*Hist.range-epi2, 2*Plot.spacer+2*Hist.range+(5/6)*Hist.range-epi2)) lines(c(-epi3,BarLength+epi3),c(2*Plot.spacer+2*Hist.range+(5/6)*Hist.range+epi2, 2*Plot.spacer+2*Hist.range+(5/6)*Hist.range+epi2)) rect(-2*BarLength,Plot.spacer+Hist.range+(5/6)*Hist.range-epi2,BarLength+epi1, Plot.spacer+Hist.range+(5/6)*Hist.range+epi2,col="white",border=FALSE) lines(c(-2*BarLength,-epi1-epi3),c(Plot.spacer+Hist.range+(5/6)*Hist.range-epi2, Plot.spacer+Hist.range+(5/6)*Hist.range-epi2)) lines(c(-2*BarLength,-epi1-epi3),c(Plot.spacer+Hist.range+(5/6)*Hist.range+epi2, Plot.spacer+Hist.range+(5/6)*Hist.range+epi2)) lines(c(-epi3,BarLength+epi3),c(Plot.spacer+Hist.range+(5/6)*Hist.range-epi2, Plot.spacer+Hist.range+(5/6)*Hist.range-epi2)) lines(c(-epi3,BarLength+epi3),c(Plot.spacer+Hist.range+(5/6)*Hist.range+epi2, Plot.spacer+Hist.range+(5/6)*Hist.range+epi2)) #Draw triangles giving the mean growth rates ignoring the zeroes (filled) #and incorporating the zeroes (open) del<-.02 points(mean(S.gro.all),2*Plot.spacer+2*Hist.range-del,col=Scol,cex=pt.cex,pch=2) points(mean(S.gro),2*Plot.spacer+2*Hist.range-del,col=Scol,cex=pt.cex+.2,pch=17) points(mean(M.gro.all),Plot.spacer+Hist.range-del,col=Mcol,cex=pt.cex,pch=2) points(mean(M.gro),Plot.spacer+Hist.range-del,col=Mcol,cex=pt.cex+.2,pch=17) points(mean(G.gro.all),-del,,col=Gcol,cex=pt.cex,pch=2) points(mean(G.gro),-del,col=Gcol,cex=pt.cex+.2,pch=17) #Plot labels delt<-0.05 text(0.0011,2*Plot.spacer+3*Hist.range-delt,"Sudden",col=Scol, cex=legsize) text(0.0011,Plot.spacer+2*Hist.range-delt,"Moderate",col=Mcol, cex=legsize) text(0.0011,Hist.range-delt,"Gradual",col=Gcol, cex=legsize) #FUNCTION: diracZero #This function gives a value of 1 if x is zero and #a value of zero otherwise #--------------------------------------------------- diracZero.f<- function(x) { ifelse(x==0,1,0) } #--------------------------------------------------- #Here we define a negative log likelihood function which assumes that #the data is distributed according to a zero-inflated normal distribution NLL.f <- function(data, param) { -sum(log(param[1]*diracZero.f(data)+(1-param[1])*dnorm(data,mean=param[2],sd=param[3]))) } #Find the maximum likelihood estimates for distribution parameters #Here we do the analysis independently for each treatment (but see below) #Note that we multiple all data by 1000 (but then rescale back when plotting) opt.S<-optim(c(.98,.8,.1),NLL.f,data=1000*S.gro.all) opt.M<-optim(c(.55,.7,.1),NLL.f,data=1000*M.gro.all) opt.G<-optim(c(.1,.7,.1),NLL.f,data=1000*G.gro.all) #Name the ML parameters S.ml.p<-opt.S$par[1] S.ml.m<-opt.S$par[2] S.ml.s<-opt.S$par[3] M.ml.p<-opt.M$par[1] M.ml.m<-opt.M$par[2] M.ml.s<-opt.M$par[3] G.ml.p<-opt.G$par[1] G.ml.m<-opt.G$par[2] G.ml.s<-opt.G$par[3] #Draw the normal curves over our histograms lines(seq(0,1.2,.01)/1000, 2*Plot.spacer+2*Hist.range+(Hist.range/.15)*BarLength*(1-opt.S$par[1])* dnorm(seq(0,1.2,.01)/1000,mean=(opt.S$par[2]/1000),sd=(opt.S$par[3]/1000)), lwd=.5,col="black",ylim=c(0,.15)) lines(seq(0,1.2,.01)/1000, Plot.spacer+Hist.range+(Hist.range/.15)*BarLength*(1-opt.M$par[1])* dnorm(seq(0,1.2,.01)/1000,mean=(opt.M$par[2]/1000),sd=(opt.M$par[3]/1000)), lwd=.5,col="black",ylim=c(0,.15)) lines(seq(0,1.2,.01)/1000, (Hist.range/.15)*BarLength*(1-opt.G$par[1])* dnorm(seq(0,1.2,.01)/1000,mean=(opt.G$par[2]/1000),sd=(opt.G$par[3]/1000)), lwd=.5,col="black",ylim=c(0,.15)) #End of Fig 1c----------------------------------------------------- #****************************************************************** #****************************************************************** #Data Analysis for Fig 1c (RUN FIG 1C FIRST)~~~~~~~~~~~~~~~~~~~~~~~ #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Here we find the ML parameters for a null model assuming that #all treatments have the same parameters (zero-inflation probability, #mean of normal, sd of normal) opt.SMG<-optim(c(.5,3,.5),NLL.f,data=1000*c(S.gro.all,M.gro.all,G.gro.all)) SMG.ml.p<-opt.SMG$par[1] SMG.ml.m<-opt.SMG$par[2] SMG.ml.s<-opt.SMG$par[3] #The growth data is put into a vector growth<-1000*c(S.gro.all,M.gro.all,G.gro.all) #Each treatment is given a unique code codep<-c(rep(1,length(S.gro.all)),rep(2,length(M.gro.all)),rep(3,length(G.gro.all))) #The starting parameters for each of the alternative models we consider #Model "SM.G" for instance is the model where the Sudden and Moderate #treatments have the same parameters for the mean of the normal, but the #Gradual treatment has a unique mean. We always assume each treatment has #a unique zero-inflation probability and an identical sd. par.SM.G=c(S.ml.p,M.ml.p,G.ml.p,SMG.ml.m,SMG.ml.m,M.ml.s) par.MG.S=c(S.ml.p,M.ml.p,G.ml.p,SMG.ml.m,SMG.ml.m,M.ml.s) par.SG.M=c(S.ml.p,M.ml.p,G.ml.p,SMG.ml.m,SMG.ml.m,M.ml.s) par.S.M.G=c(S.ml.p,M.ml.p,G.ml.p,S.ml.m,M.ml.m,G.ml.m,M.ml.s) #Each model has a relevant code codem.SM.G<-c(rep(1,length(S.gro.all)),rep(1,length(M.gro.all)),rep(2,length(G.gro.all))) codem.MG.S<-c(rep(1,length(S.gro.all)),rep(2,length(M.gro.all)),rep(2,length(G.gro.all))) codem.SG.M<-c(rep(1,length(S.gro.all)),rep(2,length(M.gro.all)),rep(1,length(G.gro.all))) codem.S.M.G<-c(rep(1,length(S.gro.all)),rep(2,length(M.gro.all)),rep(3,length(G.gro.all))) #The negative log likelihood for model where there are 2 unique means NLL2.f <- function(data, code.p, code.m, param) { p <- param[1:3][code.p] m <- param[4:5][code.m] s <- param[6] -sum(log(p*diracZero.f(data)+(1-p)*dnorm(data,mean=m,sd=s))) } #The negative log likelihood for model where there are 3 unique means NLL3.f <- function(data, code.p, code.m, param) { p <- param[1:3][code.p] m <- param[4:6][code.m] s <- param[7] -sum(log(p*diracZero.f(data)+(1-p)*dnorm(data,mean=m,sd=s))) } #We maximize the likelihood (or minimize the negative log likelihood) opt.SM.G<-optim(par.SM.G,NLL2.f,data=growth,code.p=codep,code.m=codem.SM.G,method="CG") opt.MG.S<-optim(par.MG.S,NLL2.f,data=growth,code.p=codep,code.m=codem.MG.S,method="CG") opt.SG.M<-optim(par.SG.M,NLL2.f,data=growth,code.p=codep,code.m=codem.SG.M,method="CG") opt.S.M.G<-optim(par.S.M.G,NLL3.f,data=growth,code.p=codep,code.m=codem.S.M.G,method="CG") #The log likelihoods: LL.SMG<-(-opt.SMG$value) LL.SM.G<-(-opt.SM.G$value) LL.MG.S<-(-opt.MG.S$value) LL.SG.M<-(-opt.SG.M$value) LL.S.M.G<-(-opt.S.M.G$value) #The relevant chi-square values for the log-likelihood ratio tests 1-pchisq(2*(LL.S.M.G-LL.MG.S),1) 1-pchisq(2*(LL.S.M.G-LL.SM.G),1) 1-pchisq(2*(LL.S.M.G-LL.SG.M),1) 1-pchisq(2*(LL.MG.S-LL.SMG),1) 1-pchisq(2*(LL.SM.G-LL.SMG),1) 1-pchisq(2*(LL.SG.M-LL.SMG),1) 1-pchisq(2*(LL.S.M.G-LL.SMG),2) #ANALYSIS OF SUDDEN AND GRADUAL######################## #The growth data is put into a vector growth<-1000*c(S.gro.all,G.gro.all) #Codes code.diff<-c(rep(1,length(S.gro.all)),rep(2,length(G.gro.all))) code.same<-c(rep(1,length(S.gro.all)),rep(1,length(G.gro.all))) #The negative log likelihood for models with unique or identical means and sd's NLL.null.f <- function(data, code.diff, code.same, param) { p <- param[1:2][code.diff] m <- param[3:4][code.same] s <- param[5:6][code.same] -sum(log(p*diracZero.f(data)+(1-p)*dnorm(data,mean=m,sd=s))) } NLL.diffm.f <- function(data, code.diff, code.same, param) { p <- param[1:2][code.diff] m <- param[3:4][code.diff] s <- param[5:6][code.same] -sum(log(p*diracZero.f(data)+(1-p)*dnorm(data,mean=m,sd=s))) } NLL.diffs.f <- function(data, code.diff, code.same, param) { p <- param[1:2][code.diff] m <- param[3:4][code.same] s <- param[5:6][code.diff] -sum(log(p*diracZero.f(data)+(1-p)*dnorm(data,mean=m,sd=s))) } NLL.diffms.f <- function(data, code.diff, code.same, param) { p <- param[1:2][code.diff] m <- param[3:4][code.diff] s <- param[5:6][code.diff] -sum(log(p*diracZero.f(data)+(1-p)*dnorm(data,mean=m,sd=s))) } par=c(S.ml.p,G.ml.p,S.ml.m,G.ml.m,S.ml.s,G.ml.s) #We maximize the likelihood (or minimize the negative log likelihood) opt.null<-optim(par,NLL.null.f,data=growth,code.diff=code.diff,code.same=code.same,method="CG") opt.mdif<-optim(par,NLL.diffm.f,data=growth,code.diff=code.diff,code.same=code.same,method="CG") opt.sdif<-optim(par,NLL.diffs.f,data=growth,code.diff=code.diff,code.same=code.same,method="CG") opt.msdif<-optim(par,NLL.diffms.f,data=growth,code.diff=code.diff,code.same=code.same,method="CG") opt.msdif$par #The log likelihoods: LL.null<-(-opt.null$value) LL.mdif<-(-opt.mdif$value) LL.sdif<-(-opt.sdif$value) LL.msdif<-(-opt.msdif$value) 1-pchisq(2*(LL.mdif-LL.null),1) 1-pchisq(2*(LL.sdif-LL.null),1) 1-pchisq(2*(LL.msdif-LL.mdif),1) 1-pchisq(2*(LL.msdif-LL.sdif),1) 1-pchisq(2*(LL.msdif-LL.null),2) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #****************************************************************** #Figure 2a: Mutation locations #------------------------------------------------------------------ #This function will plot the mutations in a set of strains from a set of treatments #The arguments are described below PlotMutations<-function(raw.seq.df, N.bases, section.colors, treatments, treatment.colors, treatment.labels, gene.color="ivory",gene.width=0.025,min.gap.length=30,section.buffer=1,mag.spacer=0.075,tab.spacer=0.01, treatment.label.cex=1.0,treatment.label.spacer=0.02, first.mutation.color="green", following.mutation.color="white", hist.spacer=0.015, hist.unit.height=0.015, hist.axis.spacer=0.02, hist.ticks.length=0.005, hist.numbers.cex=0.6, hist.numbers.spacer=0.01, hist.label="number of\nmutations", hist.label.spacer=0.05, hist.label.cex=0.6,xmin=-0.05,xmax=1.1,ymin=-.1,ymax=1) { #Function arguments: #------------------- # raw.seq.df is the data frame with all the mutations (must have the following columns: treatment, strain, base, and first, a bool indicating whether the mutation was first to arise in the strain) (note: each mutation would be listed as its own row) # N.bases is number of bases in the gene # section.colors is a vector giving the colors of the highlighted sections (make this long enough to cover all possible sections) # treatments is a vector of all the treatments as they appear in the data frame # treatment.colors is a vector of the colors for each treatment # treatment.labels is a vector with the labels to place on the mutation table # gene.color is the color of the gene in the figure (default="ivory") # gene.width is the width of the gene in the figure (default=0.05) # min.gap.length is minimum number of bases consecutive mutations need to spaced in order to be in different sections (default=30) # section.buffer is number of bases subtracted and added to the ends of each section to buffer the sections in the figure (default=1) # mag.spacer is width of the spacing between the gene and the section header of the table of mutations (default=0.1) # tab.spacer is width of the spacing between the section header and the table of mutations (default=0.01) # treatment.label.cex is the size of the text of the treatment labels (default=1.0) # treatment.label.spacer is the spacing between the tabel and the treatment labels (default=0.02) # mutation.color is the color of mutations in the table (default="white") # hist.spacer is the spacing between the table and the histogram (default=0.02) # hist.unit.height is the height of "one mutation" in the histogram (default=0.015) # hist.axis.spacer is the spacing between the histogram and its axis (default=0.02) # hist.ticks.length is the length of the ticks on the histogram axis (default=0.005) # hist.numbers.cex is the size of the text of the numbers on the histogram axis (default=0.6) # hist.numbers.spacer is the spacing between the ticks and the numbers (default=0.01) # hist.label is the label on the histogram axis (default="number of\nmutations") # hist.label.spacer is the spacing between the numbers and label on the histogram (default=0.05) # hist.label.cex is the size of the text of the label on the histogram axis (default=0.6) # xmin is the minimum x value of the plot (default=-0.05) # xmax is the maximum x value of the plot (default=1.1) # ymin is the minimum y value of the plot (default=-0.1) # ymax is the maximum y value of the plot (default=1) #------------------- #Reserve space for the plot plot(c(xmin,xmax),c(ymin,ymax),col="white",axes=FALSE,ylab="",xlab="") #Draw full gene rect(0,1-gene.width,1,1,col=gene.color) #Pick the specified treatments out of the data frame focal.raw.seq.df<-raw.seq.df[is.element(raw.seq.df$treatment,treatments),] #Construct site vectors/parameters sites<-sort(unique(focal.raw.seq.df$base)) Nsit<-length(sites) #Determine the number and boundaries of genetic sections to plot sec.init<-(sites[1]-section.buffer) if(Nsit==1) { sec.fin<-(sites[1]+section.buffer) N.sec<-1 } if(Nsit>1) { N.sec<-1 j<-1 for(i in 2:Nsit) { if((sites[i]-sites[i-1])>min.gap.length) { sec.init<-c(sec.init,sites[i]-section.buffer) ifelse(j==1,sec.fin<-sites[i-1]+section.buffer,sec.fin<-c(sec.fin,sites[i-1]+section.buffer)) j<-j+1 N.sec<-N.sec+1 } } ifelse(j==1,sec.fin<-sites[Nsit]+section.buffer,sec.fin<-c(sec.fin,sites[Nsit]+section.buffer)) } #Draw sections on gene for(i in 1:N.sec) { rect(sec.init[i]/N.bases,1-gene.width,sec.fin[i]/N.bases,1,col=section.colors[i],border=FALSE) } lines(c(0,1),c(1,1)) #Compute the lengths of each section and total length of all sections new.run<-numeric(N.sec) tot<-0 for(i in 1:N.sec) { new.run[i]<-sec.fin[i]-sec.init[i] tot<-tot+new.run[i] } #Put together "new" sections that are confluent on the interval [0,1], preserving relative lengths new.sec.init<-numeric(N.sec) new.sec.fin<-numeric(N.sec) tot2<-0 for(i in 1:N.sec) { new.sec.init[i]<-tot2 new.sec.fin[i]<-tot2+(new.run[i]/tot) tot2<-tot2+(new.run[i]/tot) } #"Magnify" the sections for(i in 1:N.sec) { polygon(c(sec.init[i]/N.bases,new.sec.init[i],new.sec.fin[i], sec.fin[i]/N.bases), c(1-gene.width,1-gene.width-mag.spacer,1-gene.width-mag.spacer,1-gene.width), col=section.colors[i],border=FALSE) } lines(c(0,1),c(1-gene.width,1-gene.width)) #Function to organize strains in a list first by number of mutations #then by location of the first mutation (this makes the mutation table #a bit prettier) OrderStrains<-function(X.df) { Xstrains<-sort(unique(X.df$strain)) len<-length(Xstrains) x<-rep(0,len) y<-rep(0,len) z<-rep(0,len) SO.df<-data.frame(strain=x,base=y,N.mut=z) for(i in 1:len) { SO.df$strain[i]<-X.df[X.df$strain==Xstrains[i],]$strain[1] SO.df$base[i]<-min(X.df[X.df$strain==Xstrains[i],]$base) SO.df$N.mut[i]<-length(X.df[X.df$strain==Xstrains[i],]$base) } SO.df<-SO.df[order(SO.df$N.mut,SO.df$base),] SO.df } #Organize strains in all treatments and collect information on the number #of mutations in each. for(i in 1:length(treatments)) { T.df<-OrderStrains(raw.seq.df[raw.seq.df$treatment==treatments[i],]) ifelse(i==1,strains<-T.df$strain,strains<-c(strains,T.df$strain)) ifelse(i==1,mutHist<-T.df$N.mut,mutHist<-c(mutHist,T.df$N.mut)) } Nstr<-length(strains) #the number of strains to be plotted row.width<-(1-gene.width-mag.spacer-tab.spacer)/(Nstr+1) #the width of each row (strain) #Draw confluent sections at row height for(i in 1:N.sec) { rect(new.sec.init[i],1-gene.width-mag.spacer-row.width,new.sec.fin[i], 1-gene.width-mag.spacer,col=section.colors[i]) } #Put together a list of all the bases in the table focal.bases<-c((sec.init[1]):(sec.fin[1])) if(N.sec>1) { for(i in 2:N.sec) { nex<-c((sec.init[i]):(sec.fin[i])) focal.bases<-c(focal.bases,nex) } } #Draw the proper background colors for each strain current.treatment<-raw.seq.df[raw.seq.df$strain==strains[1],]$treatment[1] for(j in 1:length(treatments)) { if(current.treatment==treatments[j]) { current.color<-treatment.colors[j] } } row.of.last.color<-0 for(i in 1:Nstr) { if(raw.seq.df[raw.seq.df$strain==strains[i],]$treatment[1] != current.treatment) { rect(0,1-gene.width-mag.spacer-row.width-tab.spacer-(row.width*(i-1)),1, 1-gene.width-mag.spacer-row.width-tab.spacer-(row.width*(row.of.last.color)),col=current.color,border=FALSE) for(j in 1:length(treatments)) { if(raw.seq.df[raw.seq.df$strain==strains[i],]$treatment[1]==treatments[j]) { current.color<-treatment.colors[j] } } row.of.last.color<-(i-1) current.treatment<-raw.seq.df[raw.seq.df$strain==strains[i],]$treatment[1] } if(i==Nstr) { rect(0,1-gene.width-mag.spacer-row.width-tab.spacer-(row.width*(i)),1, 1-gene.width-mag.spacer-row.width-tab.spacer-(row.width*(row.of.last.color)),col=current.color,border=FALSE) } } #Draw the mutations for each strain for(i in 1:Nstr) { mutations<-raw.seq.df[raw.seq.df$strain==strains[i],]$base firstbools<-raw.seq.df[raw.seq.df$strain==strains[i],]$first lfb<-length(focal.bases) for(j in 1:length(mutations)) { pos<-which(focal.bases==mutations[j]) if(firstbools[j]) { rect((pos-1)/lfb,1-gene.width-mag.spacer-row.width-tab.spacer-(row.width*(i)),pos/lfb, 1-gene.width-mag.spacer-row.width-tab.spacer-(row.width*(i-1)),col=first.mutation.color,border=FALSE) } if(!firstbools[j]) { rect((pos-1)/lfb,1-gene.width-mag.spacer-row.width-tab.spacer-(row.width*(i)),pos/lfb, 1-gene.width-mag.spacer-row.width-tab.spacer-(row.width*(i-1)),col=following.mutation.color,border=FALSE) } } } #Draw boundaries around (and in) mutation table for(i in 1:N.sec) { rect(new.sec.init[i],1-gene.width-mag.spacer-row.width-tab.spacer-(row.width*Nstr),new.sec.fin[i], 1-gene.width-mag.spacer-row.width-tab.spacer) } #Draw histogram to the right of the table for(i in 1:Nstr) { for(j in 1:length(treatments)) { if(raw.seq.df[raw.seq.df$strain==strains[i],]$treatment[1]==treatments[j]) { color<-treatment.colors[j] } } rect(1+hist.spacer,1-gene.width-mag.spacer-row.width-tab.spacer-(row.width*(i)), 1+hist.spacer+mutHist[i]*hist.unit.height, 1-gene.width-mag.spacer-row.width-tab.spacer-(row.width*(i-1)), col=color) } #Plot the treatment labels to the left of the table Tlen<-numeric(length(treatments)) Totlen<-0 for(i in 1:length(treatments)) { Tlen[i]<-length(unique(raw.seq.df[raw.seq.df$treatment==treatments[i],]$strain)) Totlen<-(Totlen+Tlen[i]) } Ty<-numeric(length(treatments)) RunTot<-0 for(i in 1:length(treatments)) { Ty[i]<-(1-gene.width-mag.spacer-row.width-tab.spacer)*(1-RunTot-0.5*(Tlen[i]/Totlen)) text(-treatment.label.spacer, Ty[i], treatment.labels[i], col=treatment.colors[i],cex=treatment.label.cex,srt=90) RunTot<-(RunTot+(Tlen[i]/Totlen)) } #Plot the histogram axis, numbers, and label arrows(1+hist.spacer,-hist.axis.spacer,1+hist.spacer+(max(mutHist)+1)*hist.unit.height,-hist.axis.spacer,length=2*hist.unit.height) for(i in 1:max(mutHist)) { lines(c(1+hist.spacer+(i)*hist.unit.height,1+hist.spacer+(i)*hist.unit.height), c(-hist.axis.spacer,-hist.axis.spacer-hist.ticks.length)) text(1+hist.spacer+(i)*hist.unit.height,-hist.numbers.spacer,i,cex=hist.numbers.cex,pos=1) } text(1+hist.spacer+(max(mutHist)+1)*.5*hist.unit.height,-hist.label.spacer,hist.label,cex=hist.label.cex,pos=1) } #Read in sequence data and set arguments setwd(wd) raw.seq.df<-read.table("MutationData.txt", header=TRUE) N.bases<-4029 treatments<-c("Sudden","Moderate","Gradual") treatment.colors<-c("tomato4","tomato3","tomato1") treatment.labels<-c("Sudden","Moderate","Gradual") section.colors<-rep(rev(c("cornsilk3","cornsilk4")),25) #Plot mutation table PlotMutations(raw.seq.df, N.bases, section.colors, treatments, treatment.colors, treatment.labels,min.gap.length=15,section.buffer=5,hist.unit.height=0.015, following.mutation.color="white",first.mutation.color="skyblue") #------------------------------------------------------------------ #End Fig 2a #****************************************************************** #****************************************************************** #Figure 2b: Nucleotide Diversity (Run Figure 2a first) #------------------------------------------------------------------ #Construct site and strain vectors/parameters sites<-sort(unique(raw.seq.df$base)) strains<-sort(unique(raw.seq.df$strain)) Nsit<-length(sites) Nstr<-length(strains) #Build the processed sequence data frame treatment<-rep("blank",Nstr) strain<-strains for(j in 1:Nstr) { treatment[j]<-paste(raw.seq.df[raw.seq.df$strain==strain[j],]$treatment[1]) } seq.df<-data.frame(treatment,strain) for(i in 1:Nsit) { x<-rep("blank",Nstr) for(j in 1:Nstr) { ifelse(any(raw.seq.df[raw.seq.df$strain==strain[j],]$base==sites[i]), x[j]<-paste(raw.seq.df[raw.seq.df$strain==strain[j] & raw.seq.df$base==sites[i],]$mut.base[1]), x[j]<-paste(raw.seq.df[raw.seq.df$base==sites[i],]$anc.base[1])) } seq.df<-data.frame(seq.df,x) } seqnames<-c("treatment","strain") for(i in 1:Nsit) { seqnames<-c(seqnames,paste("b",sites[i],sep="")) } names(seq.df)<-seqnames G.df<-seq.df[seq.df$treatment=="Gradual",] M.df<-seq.df[seq.df$treatment=="Moderate",] S.df<-seq.df[seq.df$treatment=="Sudden",] Calculate.ND<-function(X.df,treat,Nsites) { X.NucDiv<-0.0 NX<-length(X.df$treat) for(i in 1:(NX-1)) { for(j in (i+1):NX) { hap1<-rep("blank",Nsites) hap2<-rep("blank",Nsites) for(k in 1:Nsites) { hap1[k]<-paste(X.df[i,k+2]) hap2[k]<-paste(X.df[j,k+2]) } X.NucDiv<-X.NucDiv+(1.0/NX)*(1.0/NX)*(sum(hap1!=hap2)) } } X.NucDiv } #Calculate the ND for each treatment G.NucDiv<-Calculate.ND(G.df,treatment,Nsit) M.NucDiv<-Calculate.ND(M.df,treatment,Nsit) S.NucDiv<-Calculate.ND(S.df,treatment,Nsit) #plot the NDs par(mar=c(5,4,4,2)) axissize<-1.3 #size of text on axes labsize<-1.5 #size of text on axis labels Scol<-"tomato4" #color for Rapid treatment Mcol<-"tomato3" #color for Moderate treatment Gcol<-"tomato1" #color for Gradual treatment lab<-c("Sudden","Moderate","Gradual") xcoor<-barplot(c(S.NucDiv,M.NucDiv,G.NucDiv),names.arg=lab, col=treatment.colors, ylab="nucleotide diversity",xlab="treatment", cex.names=axissize,cex.lab=labsize,cex.axis=axissize,ylim=c(0,2),las=1) del<-.1 text(xcoor[1,1],S.NucDiv+del, "a", cex=1.5) text(xcoor[2,1],M.NucDiv+del, "b", cex=1.5) text(xcoor[3,1],G.NucDiv+del, "b", cex=1.5) #------------------------------------------------------------------ #End Fig 2b #****************************************************************** #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #Analysis for Figure 2b (requires running Figure 2b first) #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ N.perms<-1000 #Permutation test to compare Gradual and Sudden GS.seq.df<-seq.df[seq.df$treatment=="Gradual" | seq.df$treatment=="Sudden",] N.ext<-0 for(i in 1:N.perms) { shuffle.treatment<-sample(GS.seq.df$treatment) new.GS.seq.df<-data.frame(GS.seq.df,shuffle.treatment) new.G.df<-new.GS.seq.df[new.GS.seq.df$shuffle.treatment=="Gradual",] new.S.df<-new.GS.seq.df[new.GS.seq.df$shuffle.treatment=="Sudden",] if(abs(Calculate.ND(new.G.df,shuffle.treatment,Nsit)-Calculate.ND(new.S.df,shuffle.treatment,Nsit))>abs(G.NucDiv-S.NucDiv)) { N.ext<-N.ext+1 } } GSP<-N.ext/(N.perms) GSP #Permutation test to compare Gradual and Moderate GM.seq.df<-seq.df[seq.df$treatment=="Gradual" | seq.df$treatment=="Moderate",] N.ext<-0 for(i in 1:N.perms) { shuffle.treatment<-sample(GM.seq.df$treatment) new.GM.seq.df<-data.frame(GM.seq.df,shuffle.treatment) new.G.df<-new.GM.seq.df[new.GM.seq.df$shuffle.treatment=="Gradual",] new.M.df<-new.GM.seq.df[new.GM.seq.df$shuffle.treatment=="Moderate",] if(abs(Calculate.ND(new.G.df,shuffle.treatment,Nsit)-Calculate.ND(new.M.df,shuffle.treatment,Nsit))>abs(G.NucDiv-M.NucDiv)) { N.ext<-N.ext+1 } } GMP<-N.ext/(N.perms) GMP #Permutation test to compare Sudden and Moderate SM.seq.df<-seq.df[seq.df$treatment=="Sudden" | seq.df$treatment=="Moderate",] N.ext<-0 for(i in 1:N.perms) { shuffle.treatment<-sample(SM.seq.df$treatment) new.SM.seq.df<-data.frame(SM.seq.df,shuffle.treatment) new.S.df<-new.SM.seq.df[new.SM.seq.df$shuffle.treatment=="Sudden",] new.M.df<-new.SM.seq.df[new.SM.seq.df$shuffle.treatment=="Moderate",] if(abs(Calculate.ND(new.S.df,shuffle.treatment,Nsit)-Calculate.ND(new.M.df,shuffle.treatment,Nsit))>abs(S.NucDiv-M.NucDiv)) { N.ext<-N.ext+1 } } SMP<-N.ext/(N.perms) SMP Nboots<-1000 R.NucDiv<-Calculate.ND(R.df,treatment,Nsit) All.seq.df<-seq.df[seq.df$treatment=="Rapid" | seq.df$treatment=="Moderate" | seq.df$treatment=="Gradual",] len<-length(All.seq.df$strain) NDBoots<-numeric(Nboots) run<-0 for(i in 1:Nboots) { Resam.seq.df<-seq.df[seq.df$treatment=="Rapid",] for(j in 1:length(seq.df[seq.df$treatment=="Rapid",]$strain)) { Resam.seq.df[j,]<-All.seq.df[ceiling(len*runif(1)),] } NDBoots[i]<-Calculate.ND(Resam.seq.df,treatment,Nsit) if(NDBoots[i]abs(G.NucDiv-S.NucDiv)) { N.ext<-N.ext+1 } } GSP<-N.ext/(N.perms) GSP #Permutation test to compare Gradual and Moderate GM.seq.df<-seq.df[seq.df$treatment=="Gradual" | seq.df$treatment=="Moderate",] N.ext<-0 for(i in 1:N.perms) { shuffle.treatment<-sample(GM.seq.df$treatment) new.GM.seq.df<-data.frame(GM.seq.df,shuffle.treatment) new.G.df<-new.GM.seq.df[new.GM.seq.df$shuffle.treatment=="Gradual",] new.M.df<-new.GM.seq.df[new.GM.seq.df$shuffle.treatment=="Moderate",] if(abs(Calculate.ND(new.G.df,shuffle.treatment,Nsit)-Calculate.ND(new.M.df,shuffle.treatment,Nsit))>abs(G.NucDiv-M.NucDiv)) { N.ext<-N.ext+1 } } GMP<-N.ext/(N.perms) GMP #Permutation test to compare Sudden and Moderate SM.seq.df<-seq.df[seq.df$treatment=="Sudden" | seq.df$treatment=="Moderate",] N.ext<-0 for(i in 1:N.perms) { shuffle.treatment<-sample(SM.seq.df$treatment) new.SM.seq.df<-data.frame(SM.seq.df,shuffle.treatment) new.S.df<-new.SM.seq.df[new.SM.seq.df$shuffle.treatment=="Sudden",] new.M.df<-new.SM.seq.df[new.SM.seq.df$shuffle.treatment=="Moderate",] if(abs(Calculate.ND(new.S.df,shuffle.treatment,Nsit)-Calculate.ND(new.M.df,shuffle.treatment,Nsit))>abs(S.NucDiv-M.NucDiv)) { N.ext<-N.ext+1 } } SMP<-N.ext/(N.perms) SMP #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ #****************************************************************** #Figures 2c: Growth rates for single mutants #------------------------------------------------------------------ setwd(wd) spec.df<-read.table("Single_24h.txt", header=TRUE) gen<-c("g1546t.anc","g436t.anc","t1532g.anc","t437a.anc","c1527a.anc","c1721t.anc","g1546t","g436t","t1532g","t437a","c1527a","c1721t") wd2<-paste(wd,"/Single_PDFS",sep="") setwd(wd2) par(mar=c(5,5,4,2)) axissize<-1.3 #size of text on axes labsize<-1.5 #size of text on axis labels Scol<-"tomato4" #color for Rapid treatment Mcol<-"tomato3" #color for Moderate treatment Gcol<-"tomato1" #color for Gradual treatment con<-c("c1","c2","c4","c6","c8","c10","c12","c14","c16","c18") gencol<-c("black","black","black","black","black","black",Scol,Scol,Mcol,Mcol,Gcol,Gcol) gensdcol<-rep("lightgray",12) reps<-4 #must be even window<-5 #must be odd period<-30 cons<-length(con) gens<-length(gen) #FUNCTION: posify #This function changes all non-positive values in a #vector to the smallest positive value in the vector #or to a set positive value if all entries are #non-positive. #--------------------------------------------------- posify<-function(x,emer.value=0.0001) { y<-numeric(length(x)) smallest.strict.pos<-max(x) if(smallest.strict.pos<=0) { y<-rep(emer.value,length(x)) } if(smallest.strict.pos>0) { for(i in 1:length(x)) { if(x[i]>0 & x[i](maxval-meanval),which.min(x),which.max(x)) if(ind==1) { y<-x[2:length(x)] return(y) } if(ind==(length(x))) { y<-x[1:(length(x)-1)] return(y) } y<-c(x[1:(ind-1)],x[(ind+1):length(x)]) return(y) } #--------------------------------------------------- #Create the fitness arrays fit1<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) fit2<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) fit3<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time1<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time2<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time3<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) FIT<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) for(x in 1:gens) { for(y in 1:cons) { for(z in 1:reps) { geno<-gen[x] conc<-con[y] t<-spec.df[spec.df$rep==z & spec.df$strain==geno,]$time od<-spec.df[spec.df$rep==z & spec.df$strain==geno,4+y] od<-posify(od) times<-numeric(length(od)-window+1) slopes<-numeric(length(od)-window+1) for(i in (1+((window-1)/2)):(length(od)-((window-1)/2))){ slopes[i-(window-1)/2]<-lsfit(t[(i-(window-1)/2):(i+(window-1)/2)],od[(i-(window-1)/2):(i+(window-1)/2)])$coefficients[2] times[i-(window-1)/2]<-t[i] } fit1[x,y,z]<-max(slopes) tind<-which.max(slopes) time1[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) fit2[x,y,z]<-max(slopes) tind<-which.max(slopes) time2[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) fit3[x,y,z]<-max(slopes) tind<-which.max(slopes) time3[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) FIT[x,y,z]=(fit1[x,y,z]+fit2[x,y,z]+fit3[x,y,z])/3.0 if(z==1) { dev.set(3) pdf(paste(gen[x],"___",con[y],".pdf", sep="")) par(mfrow=c(2,reps/2)) } plot(times,(od[(1+((window-1)/2)):(length(od)-((window-1)/2))]),xlab="time (min)",ylab="OD",main=paste(gen[x],":",con[y],":",z,"\n","slope=",FIT[x,y,z]),col="white") points(times,(od[(1+((window-1)/2)):(length(od)-((window-1)/2))]),pch=19) points(time1[x,y,z],(od[(time1[x,y,z]/period)+1]),col="red",cex=2) points(time1[x,y,z],(od[(time1[x,y,z]/period)+1]),col="red",cex=3) lines(times,fit1[x,y,z]*(times-time1[x,y,z])+(od[(time1[x,y,z]/period)+1]),lwd=2,col="red") points(time2[x,y,z],(od[(time2[x,y,z]/period)+1]),col="red",cex=2) points(time2[x,y,z],(od[(time2[x,y,z]/period)+1]),col="red",cex=3) lines(times,fit2[x,y,z]*(times-time2[x,y,z])+(od[(time2[x,y,z]/period)+1]),lwd=2,col="red") points(time3[x,y,z],(od[(time3[x,y,z]/period)+1]),col="red",cex=2) points(time3[x,y,z],(od[(time3[x,y,z]/period)+1]),col="red",cex=3) lines(times,fit3[x,y,z]*(times-time3[x,y,z])+(od[(time3[x,y,z]/period)+1]),lwd=2,col="red") if(z==reps) { dev.off() } } } } fit.mean<-array(data = numeric(gens*cons), dim = c(gens,cons)) fit.sde<-array(data = numeric(gens*cons), dim = c(gens,cons)) for(x in 1:gens) { for(y in 1:cons) { fit.mean[x,y]=mean(toss.outlier(FIT[x,y,1:reps])) fit.sde[x,y]=(sd(toss.outlier(FIT[x,y,1:reps]))/sqrt(reps-1)) } } windows() setwd(wd) RifCon.df<-read.table("RifConc.txt", header=TRUE) con.val<-RifCon.df[RifCon.df$treatment=="Moderate",]$rif.conc[c(1,2,4,6,8,10,12,14,16,18)] axissize<-1.2 #size of text on axes labsize<-1.5 #size of text on axis labels Spt<-1.6 #point size for Sudden treatment Mpt<-1.2 #point size for Moderate treatment Gpt<-.7 #point size for Gradual treatment par(mar=c(5,5,4,2)) plot(c(1,length(con.val)),c(min(fit.mean)-max(fit.sde),max(fit.mean)+max(fit.sde)),col="white",xaxt="n",yaxt="n", xlab=expression(paste("rifampicin concentration (",mu,"g/mL)",sep="")),ylab=expression(paste("growth rate (",10^-4," AU/min)",sep="")),cex.lab=labsize, cex.axis=2) axis(1, at=1:length(con), labels=con.val,cex.axis=axissize) axis(2, at=c(0,0.0002,0.0004,0.0006,0.0008,0.001), labels=c(0,2,4,6,8,10), las=1, cex.axis=axissize) lines(1:cons,fit.mean[1,1:cons],col=gencol[1],lwd=2,lty=2) points(1:cons,fit.mean[1,1:cons],col=gencol[1],pch=15,cex=1.6) for(x in c(7,8,9,10,11,12)){ #polygon(c(1:cons,cons:1),c(fit.mean[x,1:cons]-fit.sde[x,1:cons],rev(fit.mean[x,1:cons]+fit.sde[x,1:cons])),col=gensdcol[x],border=FALSE) lines(1:cons,fit.mean[x,1:cons],col=gencol[x],lwd=2) ifelse(spec.df[spec.df$strain==gen[x],]$treatment[1]=="Sudden", ptsize<-Spt, ifelse(spec.df[spec.df$strain==gen[x],]$treatment[1]=="Moderate", ptsize<-Mpt, ptsize<-Gpt) ) points(1:cons,fit.mean[x,1:cons],col=gencol[x],pch=19,cex=ptsize) if(ptsize==Mpt) { points(1:cons,fit.mean[x,1:cons],col="white",pch=19,cex=ptsize) points(1:cons,fit.mean[x,1:cons],col=Mcol,pch=1,cex=ptsize) } } #End of Fig 2c----------------------------------------------------- #------------------------------------------------------------------ #****************************************************************** #****************************************************************** #Figures 4a: Growth rates for G13 #------------------------------------------------------------------ setwd(wd) spec.df<-read.table("G13_24h.txt", header=TRUE) gen<-c("g428a.anc","c1721t.anc","g428a.c1721t.anc","g428a","c1721t","g428a.c1721t") wd2<-paste(wd,"/G13_PDFS",sep="") setwd(wd2) con<-c("c1","c2","c4","c6","c8","c10","c12","c14","c16","c18") gencol<-c("black","black","black","green4","blue","red") gensdcol<-c("lightgray","lightgray","lightgray","lightgreen","lightblue","pink") reps<-4 #must be even window<-5 #must be odd period<-30 cons<-length(con) gens<-length(gen) #FUNCTION: posify #This function changes all non-positive values in a #vector to the smallest positive value in the vector #or to a set positive value if all entries are #non-positive. #--------------------------------------------------- posify<-function(x,emer.value=0.0001) { y<-numeric(length(x)) smallest.strict.pos<-max(x) if(smallest.strict.pos<=0) { y<-rep(emer.value,length(x)) } if(smallest.strict.pos>0) { for(i in 1:length(x)) { if(x[i]>0 & x[i](maxval-meanval),which.min(x),which.max(x)) if(ind==1) { y<-x[2:length(x)] return(y) } if(ind==(length(x))) { y<-x[1:(length(x)-1)] return(y) } y<-c(x[1:(ind-1)],x[(ind+1):length(x)]) return(y) } #--------------------------------------------------- #Create the fitness arrays fit1<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) fit2<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) fit3<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time1<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time2<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time3<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) FIT<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) for(x in 1:gens) { for(y in 1:cons) { for(z in 1:reps) { geno<-gen[x] conc<-con[y] t<-spec.df[spec.df$rep==z & spec.df$strain==geno,]$time od<-spec.df[spec.df$rep==z & spec.df$strain==geno,3+y] od<-posify(od) times<-numeric(length(od)-window+1) slopes<-numeric(length(od)-window+1) for(i in (1+((window-1)/2)):(length(od)-((window-1)/2))){ slopes[i-(window-1)/2]<-lsfit(t[(i-(window-1)/2):(i+(window-1)/2)],od[(i-(window-1)/2):(i+(window-1)/2)])$coefficients[2] times[i-(window-1)/2]<-t[i] } fit1[x,y,z]<-max(slopes) tind<-which.max(slopes) time1[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) fit2[x,y,z]<-max(slopes) tind<-which.max(slopes) time2[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) fit3[x,y,z]<-max(slopes) tind<-which.max(slopes) time3[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) FIT[x,y,z]=(fit1[x,y,z]+fit2[x,y,z]+fit3[x,y,z])/3.0 if(z==1) { dev.set(3) pdf(paste(gen[x],"___",con[y],".pdf", sep="")) par(mfrow=c(2,reps/2)) } plot(times,(od[(1+((window-1)/2)):(length(od)-((window-1)/2))]),xlab="time (min)",ylab="OD",main=paste(gen[x],":",con[y],":",z,"\n","slope=",FIT[x,y,z]),col="white") points(times,(od[(1+((window-1)/2)):(length(od)-((window-1)/2))]),pch=19) points(time1[x,y,z],(od[(time1[x,y,z]/period)+1]),col="red",cex=2) points(time1[x,y,z],(od[(time1[x,y,z]/period)+1]),col="red",cex=3) lines(times,fit1[x,y,z]*(times-time1[x,y,z])+(od[(time1[x,y,z]/period)+1]),lwd=2,col="red") points(time2[x,y,z],(od[(time2[x,y,z]/period)+1]),col="red",cex=2) points(time2[x,y,z],(od[(time2[x,y,z]/period)+1]),col="red",cex=3) lines(times,fit2[x,y,z]*(times-time2[x,y,z])+(od[(time2[x,y,z]/period)+1]),lwd=2,col="red") points(time3[x,y,z],(od[(time3[x,y,z]/period)+1]),col="red",cex=2) points(time3[x,y,z],(od[(time3[x,y,z]/period)+1]),col="red",cex=3) lines(times,fit3[x,y,z]*(times-time3[x,y,z])+(od[(time3[x,y,z]/period)+1]),lwd=2,col="red") if(z==reps) { dev.off() } } } } fit.mean<-array(data = numeric(gens*cons), dim = c(gens,cons)) fit.sde<-array(data = numeric(gens*cons), dim = c(gens,cons)) for(x in 1:gens) { for(y in 1:cons) { fit.mean[x,y]=mean(toss.outlier(FIT[x,y,1:reps])) fit.sde[x,y]=(sd(toss.outlier(FIT[x,y,1:reps]))/sqrt(reps-1)) } } windows() setwd(wd) RifCon.df<-read.table("RifConc.txt", header=TRUE) con.val<-RifCon.df[RifCon.df$treatment=="Moderate",]$rif.conc[c(1,2,4,6,8,10,12,14,16,18)] axissize<-1.2 #size of text on axes labsize<-1.5 #size of text on axis labels par(mar=c(5,5,4,2)) plot(c(1,length(con.val)),c(min(fit.mean)-max(fit.sde),max(fit.mean)+max(fit.sde)),col="white",xaxt="n",yaxt="n", xlab=expression(paste("rifampicin concentration (",mu,"g/mL)",sep="")),ylab=expression(paste("growth rate (",10^-4," AU/min)",sep="")),cex.lab=labsize, cex.axis=2) axis(1, at=1:length(con), labels=con.val,cex.axis=axissize) axis(2, at=c(0,0.0002,0.0004,0.0006,0.0008,0.001), labels=c(0,2,4,6,8,10), las=1, cex.axis=axissize) for(x in c(2,4,5,6)){ polygon(c(1:cons,cons:1),c(fit.mean[x,1:cons]-fit.sde[x,1:cons],rev(fit.mean[x,1:cons]+fit.sde[x,1:cons])),col=gensdcol[x],border=FALSE) lines(1:cons,fit.mean[x,1:cons],col=gencol[x],lwd=2) points(1:cons,fit.mean[x,1:cons],col=gencol[x],pch=19) } text(2,0.0002,"ancestors",col=gencol[1]) text(2,0.00015,gen[4],col=gencol[4]) text(2,0.0001,gen[5],col=gencol[5]) text(2,0.00005,gen[6],col=gencol[6]) #End of Fig 4a----------------------------------------------------- #------------------------------------------------------------------ #****************************************************************** #****************************************************************** #Figures 4c: Growth rates for G18 #------------------------------------------------------------------ setwd(wd) spec.df<-read.table("G18_24h.txt", header=TRUE) gen<-c("a443t.anc","c1527a.anc","a443t.c1527a.anc","a443t","c1527a","a443t.c1527a") wd2<-paste(wd,"/G18_PDFS",sep="") setwd(wd2) #********************************************************************************** con<-c("c1","c2","c4","c6","c8","c10","c12","c14","c16","c18") gencol<-c("black","black","black","green4","blue","red") gensdcol<-c("lightgray","lightgray","lightgray","lightgreen","lightblue","pink") reps<-4 #must be even window<-5 #must be odd period<-30 cons<-length(con) gens<-length(gen) #FUNCTION: posify #This function changes all non-positive values in a #vector to the smallest positive value in the vector #or to a set positive value if all entries are #non-positive. #--------------------------------------------------- posify<-function(x,emer.value=0.0001) { y<-numeric(length(x)) smallest.strict.pos<-max(x) if(smallest.strict.pos<=0) { y<-rep(emer.value,length(x)) } if(smallest.strict.pos>0) { for(i in 1:length(x)) { if(x[i]>0 & x[i](maxval-meanval),which.min(x),which.max(x)) if(ind==1) { y<-x[2:length(x)] return(y) } if(ind==(length(x))) { y<-x[1:(length(x)-1)] return(y) } y<-c(x[1:(ind-1)],x[(ind+1):length(x)]) return(y) } #--------------------------------------------------- #Create the fitness arrays fit1<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) fit2<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) fit3<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time1<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time2<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time3<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) FIT<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) for(x in 1:gens) { for(y in 1:cons) { for(z in 1:reps) { geno<-gen[x] conc<-con[y] t<-spec.df[spec.df$rep==z & spec.df$strain==geno,]$time od<-spec.df[spec.df$rep==z & spec.df$strain==geno,3+y] od<-posify(od) times<-numeric(length(od)-window+1) slopes<-numeric(length(od)-window+1) for(i in (1+((window-1)/2)):(length(od)-((window-1)/2))){ slopes[i-(window-1)/2]<-lsfit(t[(i-(window-1)/2):(i+(window-1)/2)],od[(i-(window-1)/2):(i+(window-1)/2)])$coefficients[2] times[i-(window-1)/2]<-t[i] } fit1[x,y,z]<-max(slopes) tind<-which.max(slopes) time1[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) fit2[x,y,z]<-max(slopes) tind<-which.max(slopes) time2[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) fit3[x,y,z]<-max(slopes) tind<-which.max(slopes) time3[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) FIT[x,y,z]=(fit1[x,y,z]+fit2[x,y,z]+fit3[x,y,z])/3.0 if(z==1) { dev.set(3) pdf(paste(gen[x],"___",con[y],".pdf", sep="")) par(mfrow=c(2,reps/2)) } plot(times,(od[(1+((window-1)/2)):(length(od)-((window-1)/2))]),xlab="time (min)",ylab="OD",main=paste(gen[x],":",con[y],":",z,"\n","slope=",FIT[x,y,z]),col="white") points(times,(od[(1+((window-1)/2)):(length(od)-((window-1)/2))]),pch=19) points(time1[x,y,z],(od[(time1[x,y,z]/period)+1]),col="red",cex=2) points(time1[x,y,z],(od[(time1[x,y,z]/period)+1]),col="red",cex=3) lines(times,fit1[x,y,z]*(times-time1[x,y,z])+(od[(time1[x,y,z]/period)+1]),lwd=2,col="red") points(time2[x,y,z],(od[(time2[x,y,z]/period)+1]),col="red",cex=2) points(time2[x,y,z],(od[(time2[x,y,z]/period)+1]),col="red",cex=3) lines(times,fit2[x,y,z]*(times-time2[x,y,z])+(od[(time2[x,y,z]/period)+1]),lwd=2,col="red") points(time3[x,y,z],(od[(time3[x,y,z]/period)+1]),col="red",cex=2) points(time3[x,y,z],(od[(time3[x,y,z]/period)+1]),col="red",cex=3) lines(times,fit3[x,y,z]*(times-time3[x,y,z])+(od[(time3[x,y,z]/period)+1]),lwd=2,col="red") if(z==reps) { dev.off() } } } } fit.mean<-array(data = numeric(gens*cons), dim = c(gens,cons)) fit.sde<-array(data = numeric(gens*cons), dim = c(gens,cons)) for(x in 1:gens) { for(y in 1:cons) { fit.mean[x,y]=mean(toss.outlier(FIT[x,y,1:reps])) fit.sde[x,y]=(sd(toss.outlier(FIT[x,y,1:reps]))/sqrt(reps-1)) } } windows() setwd(wd) RifCon.df<-read.table("RifConc.txt", header=TRUE) con.val<-RifCon.df[RifCon.df$treatment=="Moderate",]$rif.conc[c(1,2,4,6,8,10,12,14,16,18)] axissize<-1.2 #size of text on axes labsize<-1.5 #size of text on axis labels par(mar=c(5,5,4,2)) plot(c(1,length(con.val)),c(min(fit.mean)-max(fit.sde),max(fit.mean)+max(fit.sde)),col="white",xaxt="n",yaxt="n", xlab=expression(paste("rifampicin concentration (",mu,"g/mL)",sep="")),ylab=expression(paste("growth rate (",10^-4," AU/min)",sep="")),cex.lab=labsize, cex.axis=2) axis(1, at=1:length(con), labels=con.val,cex.axis=axissize) axis(2, at=c(0,0.0002,0.0004,0.0006,0.0008,0.001), labels=c(0,2,4,6,8,10), las=1, cex.axis=axissize) for(x in c(5,1,4,6)){ polygon(c(1:cons,cons:1),c(fit.mean[x,1:cons]-fit.sde[x,1:cons],rev(fit.mean[x,1:cons]+fit.sde[x,1:cons])),col=gensdcol[x],border=FALSE) lines(1:cons,fit.mean[x,1:cons],col=gencol[x],lwd=2) points(1:cons,fit.mean[x,1:cons],col=gencol[x],pch=19) } text(2,0.0002,"ancestors",col=gencol[1]) text(2,0.00015,gen[4],col=gencol[4]) text(2,0.0001,gen[5],col=gencol[5]) text(2,0.00005,gen[6],col=gencol[6]) #End of Fig 4c----------------------------------------------------- #------------------------------------------------------------------ #****************************************************************** #****************************************************************** #Figures 4b,d: Generating Mathematica text files #------------------------------------------------------------------ setwd(wd) viab.df<-read.table("Viability.txt", header=TRUE) reps<-6 #must be even crit.frac<-(1/6) con<-c("c1","c2","c4","c6","c8","c10","c12","c14","c16","c18") cons<-length(con) gen18<-c("a443t.anc","c1527a.anc","a443t.c1527a.anc","a443t","c1527a","a443t.c1527a") gen13<-c("g428a.anc","c1721t.anc","g428a.c1721t.anc","g428a","c1721t","g428a.c1721t") PV<-array(data = numeric(4*cons), dim = c(4,cons)) EV<-array(data = numeric(4*cons), dim = c(4,cons)) rownames(PV)<-c("ab","Ab","aB","AB") rownames(EV)<-c("ab","Ab","aB","AB") colnames(PV)<-con colnames(EV)<-con gen<-gen13 lin<-"G13" gens<-length(gen) PERSIST<-array(data = numeric(gens*cons), dim = c(gens,cons)) EMERGE<-array(data = numeric(gens*cons), dim = c(gens,cons)) for(x in 1:gens) { for(y in 1:cons) { geno<-gen[x] conc<-con[y] PERSIST[x,y]<-sum(viab.df[viab.df$lineage==lin & viab.df$strain==geno & viab.df$concentration==conc & viab.df$dilution==0.025,]$growth)/reps EMERGE[x,y]<-sum(viab.df[viab.df$lineage==lin & viab.df$strain==geno & viab.df$concentration==conc & viab.df$dilution==0.00001,]$growth)/reps } } for(x in 1:4) { for(y in 1:cons) { if(x==1) { PV[x,y]<-(PERSIST[1,y]+PERSIST[2,y]+PERSIST[3,y])/3 EV[x,y]<-(EMERGE[1,y]+EMERGE[2,y]+EMERGE[3,y])/3 } if(x>1) { PV[x,y]<-PERSIST[x+2,y] EV[x,y]<-EMERGE[x+2,y] } if(PV[x,y]<=crit.frac) { PV[x,y]<-0 } if(EV[x,y]<=crit.frac) { EV[x,y]<-0 } } } PV13<-PV EV13<-EV write.table(PV13,"PV13.txt",row.names=FALSE,col.names=FALSE) gen<-gen18 lin<-"G18" gens<-length(gen) PERSIST<-array(data = numeric(gens*cons), dim = c(gens,cons)) EMERGE<-array(data = numeric(gens*cons), dim = c(gens,cons)) for(x in 1:gens) { for(y in 1:cons) { geno<-gen[x] conc<-con[y] PERSIST[x,y]<-sum(viab.df[viab.df$lineage==lin & viab.df$strain==geno & viab.df$concentration==conc & viab.df$dilution==0.025,]$growth)/reps EMERGE[x,y]<-sum(viab.df[viab.df$lineage==lin & viab.df$strain==geno & viab.df$concentration==conc & viab.df$dilution==0.00001,]$growth)/reps } } for(x in 1:4) { for(y in 1:cons) { if(x==1) { PV[x,y]<-(PERSIST[1,y]+PERSIST[2,y]+PERSIST[3,y])/3 EV[x,y]<-(EMERGE[1,y]+EMERGE[2,y]+EMERGE[3,y])/3 } if(x>1) { PV[x,y]<-PERSIST[x+2,y] EV[x,y]<-EMERGE[x+2,y] } if(PV[x,y]<=crit.frac) { PV[x,y]<-0 } if(EV[x,y]<=crit.frac) { EV[x,y]<-0 } } } PV18<-PV EV18<-EV write.table(PV18,"PV18.txt",row.names=FALSE,col.names=FALSE) edge<-c("ab.Ab","ab.aB","Ab.AB","aB.AB") edges<-length(edge) EdgeV<-array(data = numeric(edges*cons), dim = c(edges,cons)) rownames(EdgeV)<-edge colnames(EdgeV)<-con PV<-PV13 EV<-EV13 for(x in edge) { for(y in 1:cons) { if(x=="ab.Ab") { ifelse(PV["ab",y]>crit.frac & PV["Ab",y]>crit.frac, EdgeV[x,y]<-1, ifelse( (EV["ab",y]>crit.frac | EV["Ab",y]>crit.frac), EdgeV[x,y]<-1, EdgeV[x,y]<-0)) } if(x=="ab.aB") { ifelse(PV["ab",y]>crit.frac & PV["aB",y]>crit.frac, EdgeV[x,y]<-1, ifelse( (EV["ab",y]>crit.frac | EV["aB",y]>crit.frac), EdgeV[x,y]<-1, EdgeV[x,y]<-0)) } if(x=="Ab.AB") { ifelse(PV["Ab",y]>crit.frac & PV["AB",y]>crit.frac, EdgeV[x,y]<-1, ifelse( (EV["Ab",y]>crit.frac | EV["AB",y]>crit.frac), EdgeV[x,y]<-1, EdgeV[x,y]<-0)) } if(x=="aB.AB") { ifelse(PV["aB",y]>crit.frac & PV["AB",y]>crit.frac, EdgeV[x,y]<-1, ifelse( (EV["aB",y]>crit.frac | EV["AB",y]>crit.frac), EdgeV[x,y]<-1, EdgeV[x,y]<-0)) } } } EdgeV13<-EdgeV write.table(EdgeV13,"EdgeV13.txt",row.names=FALSE,col.names=FALSE) PV<-PV18 EV<-EV18 for(x in edge) { for(y in 1:cons) { if(x=="ab.Ab") { ifelse(PV["ab",y]>crit.frac & PV["Ab",y]>crit.frac, EdgeV[x,y]<-1, ifelse( (EV["ab",y]>crit.frac | EV["Ab",y]>crit.frac), EdgeV[x,y]<-1, EdgeV[x,y]<-0)) } if(x=="ab.aB") { ifelse(PV["ab",y]>crit.frac & PV["aB",y]>crit.frac, EdgeV[x,y]<-1, ifelse( (EV["ab",y]>crit.frac | EV["aB",y]>crit.frac), EdgeV[x,y]<-1, EdgeV[x,y]<-0)) } if(x=="Ab.AB") { ifelse(PV["Ab",y]>crit.frac & PV["AB",y]>crit.frac, EdgeV[x,y]<-1, ifelse( (EV["Ab",y]>crit.frac | EV["AB",y]>crit.frac), EdgeV[x,y]<-1, EdgeV[x,y]<-0)) } if(x=="aB.AB") { ifelse(PV["aB",y]>crit.frac & PV["AB",y]>crit.frac, EdgeV[x,y]<-1, ifelse( (EV["aB",y]>crit.frac | EV["AB",y]>crit.frac), EdgeV[x,y]<-1, EdgeV[x,y]<-0)) } } } EdgeV18<-EdgeV write.table(EdgeV18,"EdgeV18.txt",row.names=FALSE,col.names=FALSE) setwd(wd) comp18.df<-read.table("G18EdgeCompetitions.txt", header=TRUE) comp13.df<-read.table("G13EdgeCompetitions.txt", header=TRUE) con<-c(1,2,4,6,8,10,12,14,16,18) cons<-length(con) comp.df<-comp13.df gen<-gen13 EdgeV<-EdgeV13 reps<-4 #Create the fitness arrays FIT<-array(data = numeric(edges*cons*reps), dim = c(edges,cons,reps)) for(x in 1:edges) { for(y in 1:cons) { edg<-edge[x] conc<-con[y] if(edg=="ab.Ab") { str1=gen[1] str2=gen[4] } if(edg=="ab.aB") { str1=gen[2] str2=gen[5] } if(edg=="Ab.AB") { str1=gen[4] str2=gen[6] } if(edg=="aB.AB") { str1=gen[5] str2=gen[6] } cur.comp0<-comp.df[comp.df$strain1==str1 & comp.df$strain2==str2 & comp.df$concentration==conc & comp.df$time==0,] cur.comp48<-comp.df[comp.df$strain1==str1 & comp.df$strain2==str2 & comp.df$concentration==conc & comp.df$time==48,] cfu1i=(cur.comp0$count1)*(10^(-cur.comp0$dilution1+1))*(1/80) cfu2i=(cur.comp0$count2)*(10^(-cur.comp0$dilution2+1))*(1/80) cfu1f=(cur.comp48$count1)*(10^(-cur.comp48$dilution1+1)) cfu2f=(cur.comp48$count2)*(10^(-cur.comp48$dilution2+1)) fit<-(cfu1f/cfu1i)/(cfu2f/cfu2i) FIT[x,y,1]<-fit[1] FIT[x,y,2]<-fit[2] FIT[x,y,3]<-fit[3] FIT[x,y,4]<-fit[4] } } EdgeDir<-array(data = numeric(edges*cons), dim = c(edges,cons)) rownames(EdgeDir)<-edge colnames(EdgeDir)<-con corrected.alpha<-0.05/sum(EdgeV) for(x in 1:edges) { for(y in 1:cons) { ifelse(EdgeV[x,y]==1, ifelse(t.test(FIT[x,y,1:reps],mu=1,alternative="two.sided")[3]$p.value1,EdgeDir[x,y]<-(-1),EdgeDir[x,y]<-1), EdgeDir[x,y]<-0), EdgeDir[x,y]<-0) } } EdgeDir13<-EdgeDir write.table(EdgeDir13,"EdgeDir13.txt",row.names=FALSE,col.names=FALSE) comp.df<-comp18.df gen<-gen18 EdgeV<-EdgeV18 reps<-4 #Create the fitness arrays FIT<-array(data = numeric(edges*cons*reps), dim = c(edges,cons,reps)) for(x in 1:edges) { for(y in 1:cons) { edg<-edge[x] conc<-con[y] if(edg=="ab.Ab") { str1=gen[1] str2=gen[4] } if(edg=="ab.aB") { str1=gen[2] str2=gen[5] } if(edg=="Ab.AB") { str1=gen[4] str2=gen[6] } if(edg=="aB.AB") { str1=gen[5] str2=gen[6] } cur.comp0<-comp.df[comp.df$strain1==str1 & comp.df$strain2==str2 & comp.df$concentration==conc & comp.df$time==0,] cur.comp48<-comp.df[comp.df$strain1==str1 & comp.df$strain2==str2 & comp.df$concentration==conc & comp.df$time==48,] cfu1i=(cur.comp0$count1)*(10^(-cur.comp0$dilution1+1))*(1/80) cfu2i=(cur.comp0$count2)*(10^(-cur.comp0$dilution2+1))*(1/80) cfu1f=(cur.comp48$count1)*(10^(-cur.comp48$dilution1+1)) cfu2f=(cur.comp48$count2)*(10^(-cur.comp48$dilution2+1)) fit<-(cfu1f/cfu1i)/(cfu2f/cfu2i) FIT[x,y,1]<-fit[1] FIT[x,y,2]<-fit[2] FIT[x,y,3]<-fit[3] FIT[x,y,4]<-fit[4] } } EdgeDir<-array(data = numeric(edges*cons), dim = c(edges,cons)) rownames(EdgeDir)<-edge colnames(EdgeDir)<-con corrected.alpha<-0.05/sum(EdgeV) for(x in 1:edges) { for(y in 1:cons) { ifelse(EdgeV[x,y]==1, ifelse(t.test(FIT[x,y,1:reps],mu=1,alternative="two.sided")[3]$p.value1,EdgeDir[x,y]<-(-1),EdgeDir[x,y]<-1), EdgeDir[x,y]<-0), EdgeDir[x,y]<-0) } } EdgeDir18<-EdgeDir write.table(EdgeDir18,"EdgeDir18.txt",row.names=FALSE,col.names=FALSE) #End of Fig 4b,d--------------------------------------------------- #------------------------------------------------------------------ #****************************************************************** #****************************************************************** #Supplementary Figure 1: Ara Marker Check #------------------------------------------------------------------ par(mfrow=c(1,1),mar=c(5,5,4,2),oma=c(1,1,1,1)) setwd(wd) ara.df<-read.table("AraCompetitionResults.txt", header=TRUE) strains<-paste(unique(ara.df$strain)) pvals<-rep(-1,length(unique(ara.df$strain))) for(i in 1:length(unique(ara.df$strain))) { x<-ara.df[ara.df$strain==unique(ara.df$strain)[i],]$fitness pvals[i]<-t.test(x,mu=1,alternative="two.sided")[3]$p.value } ara.stats.df<-data.frame(strains,pvals) legsize<-1.5 #size of text in legend axissize<-1.2 #size of text on axes labsize<-1.5 #size of text on axis labels ab13<-ara.df[ara.df$strain==unique(ara.df$strain)[1],]$fitness Ab13<-ara.df[ara.df$strain==unique(ara.df$strain)[2],]$fitness aB13<-ara.df[ara.df$strain==unique(ara.df$strain)[4],]$fitness AB13<-ara.df[ara.df$strain==unique(ara.df$strain)[5],]$fitness ab18<-ara.df[ara.df$strain==unique(ara.df$strain)[6],]$fitness Ab18<-ara.df[ara.df$strain==unique(ara.df$strain)[7],]$fitness aB18<-ara.df[ara.df$strain==unique(ara.df$strain)[9],]$fitness AB18<-ara.df[ara.df$strain==unique(ara.df$strain)[10],]$fitness boxplot(ab13,Ab13,aB13,AB13,ab18,Ab18,aB18,AB18,names=rep(c("ab","Ab","aB","AB"),2), col=c(rep("red",4),rep("blue",4)),ylab="fitness",xlab="genotype", cex.names=axissize,cex.lab=labsize,cex.axis=axissize,ylim=c(0,2),las=1) for(i in 1:8) { text(i,1.5,"ns") } #End of Supp Fig 1------------------------------------------------- #------------------------------------------------------------------ #****************************************************************** #****************************************************************** #Supplementary Figure 2: Ancestral reconstruction #------------------------------------------------------------------ setwd(wd) spec.df<-read.table("AncestorReconstruction.txt", header=TRUE) spec.df[spec.df$strain=="BK26",] gen<-c("BK26") #gen<-c("g428a.anc","c1721t.anc","g428a.c1721t.anc","g428a","c1721t","g428a.c1721t") gens<-length(gen) wd2<-paste(wd,"/ANC_PDFS",sep="") setwd(wd2) con<-c("c1","c2","c4","c6","c8","c10","c12","c14","c16","c18") gencol<-c("black") gensdcol<-c("lightgray") reps<-20 #must be even window<-5 #must be odd period<-30 cons<-length(con) #FUNCTION: posify #This function changes all non-positive values in a #vector to the smallest positive value in the vector #or to a set positive value if all entries are #non-positive. #--------------------------------------------------- posify<-function(x,emer.value=0.0001) { y<-numeric(length(x)) smallest.strict.pos<-max(x) if(smallest.strict.pos<=0) { y<-rep(emer.value,length(x)) } if(smallest.strict.pos>0) { for(i in 1:length(x)) { if(x[i]>0 & x[i](maxval-meanval),which.min(x),which.max(x)) if(ind==1) { y<-x[2:length(x)] return(y) } if(ind==(length(x))) { y<-x[1:(length(x)-1)] return(y) } y<-c(x[1:(ind-1)],x[(ind+1):length(x)]) return(y) } #--------------------------------------------------- #Create the fitness arrays fit1<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) fit2<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) fit3<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time1<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time2<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time3<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) FIT<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) for(x in 1:gens) { for(y in 1:cons) { for(z in 1:reps) { geno<-gen[x] conc<-con[y] t<-spec.df[spec.df$rep==z & spec.df$strain==geno,]$time od<-spec.df[spec.df$rep==z & spec.df$strain==geno,3+y] od<-posify(od) times<-numeric(length(od)-window+1) slopes<-numeric(length(od)-window+1) for(i in (1+((window-1)/2)):(length(od)-((window-1)/2))){ slopes[i-(window-1)/2]<-lsfit(t[(i-(window-1)/2):(i+(window-1)/2)],od[(i-(window-1)/2):(i+(window-1)/2)])$coefficients[2] times[i-(window-1)/2]<-t[i] } fit1[x,y,z]<-max(slopes) tind<-which.max(slopes) time1[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) fit2[x,y,z]<-max(slopes) tind<-which.max(slopes) time2[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) fit3[x,y,z]<-max(slopes) tind<-which.max(slopes) time3[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) FIT[x,y,z]=(fit1[x,y,z]+fit2[x,y,z]+fit3[x,y,z])/3.0 } } } fit.max<-array(data = numeric(gens*cons), dim = c(gens,cons)) fit.min<-array(data = numeric(gens*cons), dim = c(gens,cons)) fit.mean<-array(data = numeric(gens*cons), dim = c(gens,cons)) fit.sde<-array(data = numeric(gens*cons), dim = c(gens,cons)) for(x in 1:gens) { for(y in 1:cons) { v<-toss.outlier(FIT[x,y,1:reps]) fit.max[x,y]=max(v) fit.min[x,y]=min(v) fit.mean[x,y]=mean(v) fit.sde[x,y]=(sd(v)/sqrt(reps-1)) } } setwd(wd) RifCon.df<-read.table("RifConc.txt", header=TRUE) con.val<-RifCon.df[RifCon.df$treatment=="Moderate",]$rif.conc[c(1,2,4,6,8,10,12,14,16,18)] axissize<-1.2 #size of text on axes labsize<-1.5 #size of text on axis labels par(mar=c(5,5,4,2)) plot(c(1,length(con.val)),c(min(fit.mean)-max(fit.sde),max(fit.max[x,1:cons])),col="white",xaxt="n",yaxt="n", xlab=expression(paste("rifampicin concentration (",mu,"g/mL)",sep="")),ylab=expression(paste("growth rate (",10^-4," AU/min)",sep="")),cex.lab=labsize, cex.axis=2) axis(1, at=1:length(con), labels=con.val,cex.axis=axissize) axis(2, at=c(0,0.0002,0.0004,0.0006,0.0008,0.001), labels=c(0,2,4,6,8,10), las=1, cex.axis=axissize) lines(1:cons,fit.max[x,1:cons],col="gray",lwd=2, lty=2) lines(1:cons,fit.min[x,1:cons],col="gray",lwd=2, lty=2) for(x in c(1)){ polygon(c(1:cons,cons:1),c(fit.mean[x,1:cons]-fit.sde[x,1:cons],rev(fit.mean[x,1:cons]+fit.sde[x,1:cons])),col=gensdcol[x],border=FALSE) lines(1:cons,fit.mean[x,1:cons],col=gencol[x],lwd=2) points(1:cons,fit.mean[x,1:cons],col=gencol[x],pch=19) } labs<-c("a","b","c","d","e","f","g","h","i","j","k","l","m","n","o") j<-1 par(mfrow=c(5,3),mar=c(1,1,1,1),oma=c(1,1,1,1)) for(i in c(15,13,14,1,2,3,4,5,6,10,11,12,7,8,9)) { gen<-paste(unique(spec.df$strain)[i]) gens<-length(gen) reps<-4 #Create the fitness arrays fit1<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) fit2<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) fit3<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time1<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time2<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time3<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) FIT<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) for(x in 1:gens) { for(y in 1:cons) { for(z in 1:reps) { geno<-gen[x] conc<-con[y] t<-spec.df[spec.df$rep==z & spec.df$strain==geno,]$time od<-spec.df[spec.df$rep==z & spec.df$strain==geno,3+y] od<-posify(od) times<-numeric(length(od)-window+1) slopes<-numeric(length(od)-window+1) for(i in (1+((window-1)/2)):(length(od)-((window-1)/2))){ slopes[i-(window-1)/2]<-lsfit(t[(i-(window-1)/2):(i+(window-1)/2)],od[(i-(window-1)/2):(i+(window-1)/2)])$coefficients[2] times[i-(window-1)/2]<-t[i] } fit1[x,y,z]<-max(slopes) tind<-which.max(slopes) time1[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) fit2[x,y,z]<-max(slopes) tind<-which.max(slopes) time2[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) fit3[x,y,z]<-max(slopes) tind<-which.max(slopes) time3[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) FIT[x,y,z]=(fit1[x,y,z]+fit2[x,y,z]+fit3[x,y,z])/3.0 } } } fit.mean<-array(data = numeric(gens*cons), dim = c(gens,cons)) fit.sde<-array(data = numeric(gens*cons), dim = c(gens,cons)) for(x in 1:gens) { for(y in 1:cons) { v<-toss.outlier(FIT[x,y,1:reps]) fit.mean[x,y]=mean(v) fit.sde[x,y]=(sd(v)/sqrt(reps-1)) } } plot(c(1,length(con.val)),c(min(fit.mean)-max(fit.sde),max(fit.max[x,1:cons])),col="white",xaxt="n",yaxt="n", xlab="",ylab="",cex.lab=labsize, cex.axis=2) axis(1, at=1:length(con), labels=rep("",length(con)),cex.axis=axissize) axis(2, at=c(0,0.0002,0.0004,0.0006,0.0008,0.001), labels=c("","","","","",""), las=1, cex.axis=axissize) lines(1:cons,fit.max[x,1:cons],col="red",lwd=1, lty=2) lines(1:cons,fit.min[x,1:cons],col="red",lwd=1, lty=2) for(x in c(1)){ polygon(c(1:cons,cons:1),c(fit.mean[x,1:cons]-fit.sde[x,1:cons],rev(fit.mean[x,1:cons]+fit.sde[x,1:cons])),col=gensdcol[x],border=FALSE) lines(1:cons,fit.mean[x,1:cons],col=gencol[x],lwd=1.5) points(1:cons,fit.mean[x,1:cons],col=gencol[x],pch=19, cex=.75) } text(9.5,0.001,labs[j],cex=1.25) j<-j+1 } #End of Supp Fig 2------------------------------------------------- #------------------------------------------------------------------ #****************************************************************** #****************************************************************** #Supplemental Figure 3a: Growth rates for M358 #------------------------------------------------------------------ setwd(wd) spec.df<-read.table("M358_24h.txt", header=TRUE) gen<-c("t437a.anc","a1685c.anc","t437a.a1685c.anc","t437a","a1685c","t437a.a1685c") wd2<-paste(wd,"/M358_PDFS",sep="") setwd(wd2) con<-c("c1","c2","c4","c6","c8","c10","c12","c14","c16","c18") gencol<-c("black","black","black","green4","blue","red") gensdcol<-c("lightgray","lightgray","lightgray","lightgreen","lightblue","pink") reps<-4 #must be even window<-5 #must be odd period<-30 cons<-length(con) gens<-length(gen) #FUNCTION: posify #This function changes all non-positive values in a #vector to the smallest positive value in the vector #or to a set positive value if all entries are #non-positive. #--------------------------------------------------- posify<-function(x,emer.value=0.0001) { y<-numeric(length(x)) smallest.strict.pos<-max(x) if(smallest.strict.pos<=0) { y<-rep(emer.value,length(x)) } if(smallest.strict.pos>0) { for(i in 1:length(x)) { if(x[i]>0 & x[i](maxval-meanval),which.min(x),which.max(x)) if(ind==1) { y<-x[2:length(x)] return(y) } if(ind==(length(x))) { y<-x[1:(length(x)-1)] return(y) } y<-c(x[1:(ind-1)],x[(ind+1):length(x)]) return(y) } #--------------------------------------------------- #Create the fitness arrays fit1<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) fit2<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) fit3<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time1<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time2<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time3<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) FIT<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) for(x in 1:gens) { for(y in 1:cons) { for(z in 1:reps) { geno<-gen[x] conc<-con[y] t<-spec.df[spec.df$rep==z & spec.df$strain==geno,]$time od<-spec.df[spec.df$rep==z & spec.df$strain==geno,3+y] od<-posify(od) times<-numeric(length(od)-window+1) slopes<-numeric(length(od)-window+1) for(i in (1+((window-1)/2)):(length(od)-((window-1)/2))){ slopes[i-(window-1)/2]<-lsfit(t[(i-(window-1)/2):(i+(window-1)/2)],od[(i-(window-1)/2):(i+(window-1)/2)])$coefficients[2] times[i-(window-1)/2]<-t[i] } fit1[x,y,z]<-max(slopes) tind<-which.max(slopes) time1[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) fit2[x,y,z]<-max(slopes) tind<-which.max(slopes) time2[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) fit3[x,y,z]<-max(slopes) tind<-which.max(slopes) time3[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) FIT[x,y,z]=(fit1[x,y,z]+fit2[x,y,z]+fit3[x,y,z])/3.0 if(z==1) { dev.set(3) pdf(paste(gen[x],"___",con[y],".pdf", sep="")) par(mfrow=c(2,reps/2)) } plot(times,(od[(1+((window-1)/2)):(length(od)-((window-1)/2))]),xlab="time (min)",ylab="OD",main=paste(gen[x],":",con[y],":",z,"\n","slope=",FIT[x,y,z]),col="white") points(times,(od[(1+((window-1)/2)):(length(od)-((window-1)/2))]),pch=19) points(time1[x,y,z],(od[(time1[x,y,z]/period)+1]),col="red",cex=2) points(time1[x,y,z],(od[(time1[x,y,z]/period)+1]),col="red",cex=3) lines(times,fit1[x,y,z]*(times-time1[x,y,z])+(od[(time1[x,y,z]/period)+1]),lwd=2,col="red") points(time2[x,y,z],(od[(time2[x,y,z]/period)+1]),col="red",cex=2) points(time2[x,y,z],(od[(time2[x,y,z]/period)+1]),col="red",cex=3) lines(times,fit2[x,y,z]*(times-time2[x,y,z])+(od[(time2[x,y,z]/period)+1]),lwd=2,col="red") points(time3[x,y,z],(od[(time3[x,y,z]/period)+1]),col="red",cex=2) points(time3[x,y,z],(od[(time3[x,y,z]/period)+1]),col="red",cex=3) lines(times,fit3[x,y,z]*(times-time3[x,y,z])+(od[(time3[x,y,z]/period)+1]),lwd=2,col="red") if(z==reps) { dev.off() } } } } fit.mean<-array(data = numeric(gens*cons), dim = c(gens,cons)) fit.sde<-array(data = numeric(gens*cons), dim = c(gens,cons)) for(x in 1:gens) { for(y in 1:cons) { fit.mean[x,y]=mean(toss.outlier(FIT[x,y,1:reps])) fit.sde[x,y]=(sd(toss.outlier(FIT[x,y,1:reps]))/sqrt(reps-1)) } } windows() setwd(wd) RifCon.df<-read.table("RifConc.txt", header=TRUE) con.val<-RifCon.df[RifCon.df$treatment=="Moderate",]$rif.conc[c(1,2,4,6,8,10,12,14,16,18)] axissize<-1.2 #size of text on axes labsize<-1.5 #size of text on axis labels par(mar=c(5,5,4,2)) plot(c(1,length(con.val)),c(min(fit.mean)-max(fit.sde),max(fit.mean)+max(fit.sde)),col="white",xaxt="n",yaxt="n", xlab=expression(paste("rifampicin concentration (",mu,"g/mL)",sep="")),ylab=expression(paste("growth rate (",10^-4," AU/min)",sep="")),cex.lab=labsize, cex.axis=2) axis(1, at=1:length(con), labels=con.val,cex.axis=axissize) axis(2, at=c(0,0.0002,0.0004,0.0006,0.0008,0.001), labels=c(0,2,4,6,8,10), las=1, cex.axis=axissize) for(x in c(2,4,5,6)){ polygon(c(1:cons,cons:1),c(fit.mean[x,1:cons]-fit.sde[x,1:cons],rev(fit.mean[x,1:cons]+fit.sde[x,1:cons])),col=gensdcol[x],border=FALSE) lines(1:cons,fit.mean[x,1:cons],col=gencol[x],lwd=2) points(1:cons,fit.mean[x,1:cons],col=gencol[x],pch=19) } text(2,0.0002,"ancestors",col=gencol[1]) text(2,0.00015,gen[4],col=gencol[4]) text(2,0.0001,gen[5],col=gencol[5]) text(2,0.00005,gen[6],col=gencol[6]) #End of Supp Fig 3a------------------------------------------------ #------------------------------------------------------------------ #****************************************************************** #****************************************************************** #Supplemental Figure 3c: Growth rates for M367 #------------------------------------------------------------------ #PARAMETERS TO CHANGE************************************************************** setwd(wd) spec.df<-read.table("M367_24h.txt", header=TRUE) gen<-c("t1532g.anc","g1546a.anc","t1532g.g1546a.anc","t1532g","g1546a","t1532g.g1546a") wd2<-paste(wd,"/M367_PDFS" ,sep="") setwd(wd2) #********************************************************************************** con<-c("c1","c2","c4","c6","c8","c10","c12","c14","c16","c18") gencol<-c("black","black","black","green4","blue","red") gensdcol<-c("lightgray","lightgray","lightgray","lightgreen","lightblue","pink") reps<-4 #must be even window<-5 #must be odd period<-30 cons<-length(con) gens<-length(gen) #FUNCTION: posify #This function changes all non-positive values in a #vector to the smallest positive value in the vector #or to a set positive value if all entries are #non-positive. #--------------------------------------------------- posify<-function(x,emer.value=0.0001) { y<-numeric(length(x)) smallest.strict.pos<-max(x) if(smallest.strict.pos<=0) { y<-rep(emer.value,length(x)) } if(smallest.strict.pos>0) { for(i in 1:length(x)) { if(x[i]>0 & x[i](maxval-meanval),which.min(x),which.max(x)) if(ind==1) { y<-x[2:length(x)] return(y) } if(ind==(length(x))) { y<-x[1:(length(x)-1)] return(y) } y<-c(x[1:(ind-1)],x[(ind+1):length(x)]) return(y) } #--------------------------------------------------- #Create the fitness arrays fit1<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) fit2<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) fit3<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time1<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time2<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time3<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) FIT<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) for(x in 1:gens) { for(y in 1:cons) { for(z in 1:reps) { geno<-gen[x] conc<-con[y] t<-spec.df[spec.df$rep==z & spec.df$strain==geno,]$time od<-spec.df[spec.df$rep==z & spec.df$strain==geno,3+y] od<-posify(od) times<-numeric(length(od)-window+1) slopes<-numeric(length(od)-window+1) for(i in (1+((window-1)/2)):(length(od)-((window-1)/2))){ slopes[i-(window-1)/2]<-lsfit(t[(i-(window-1)/2):(i+(window-1)/2)],od[(i-(window-1)/2):(i+(window-1)/2)])$coefficients[2] times[i-(window-1)/2]<-t[i] } fit1[x,y,z]<-max(slopes) tind<-which.max(slopes) time1[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) fit2[x,y,z]<-max(slopes) tind<-which.max(slopes) time2[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) fit3[x,y,z]<-max(slopes) tind<-which.max(slopes) time3[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) FIT[x,y,z]=(fit1[x,y,z]+fit2[x,y,z]+fit3[x,y,z])/3.0 if(z==1) { dev.set(3) pdf(paste(gen[x],"___",con[y],".pdf", sep="")) par(mfrow=c(2,reps/2)) } plot(times,(od[(1+((window-1)/2)):(length(od)-((window-1)/2))]),xlab="time (min)",ylab="OD",main=paste(gen[x],":",con[y],":",z,"\n","slope=",FIT[x,y,z]),col="white") points(times,(od[(1+((window-1)/2)):(length(od)-((window-1)/2))]),pch=19) points(time1[x,y,z],(od[(time1[x,y,z]/period)+1]),col="red",cex=2) points(time1[x,y,z],(od[(time1[x,y,z]/period)+1]),col="red",cex=3) lines(times,fit1[x,y,z]*(times-time1[x,y,z])+(od[(time1[x,y,z]/period)+1]),lwd=2,col="red") points(time2[x,y,z],(od[(time2[x,y,z]/period)+1]),col="red",cex=2) points(time2[x,y,z],(od[(time2[x,y,z]/period)+1]),col="red",cex=3) lines(times,fit2[x,y,z]*(times-time2[x,y,z])+(od[(time2[x,y,z]/period)+1]),lwd=2,col="red") points(time3[x,y,z],(od[(time3[x,y,z]/period)+1]),col="red",cex=2) points(time3[x,y,z],(od[(time3[x,y,z]/period)+1]),col="red",cex=3) lines(times,fit3[x,y,z]*(times-time3[x,y,z])+(od[(time3[x,y,z]/period)+1]),lwd=2,col="red") if(z==reps) { dev.off() } } } } fit.mean<-array(data = numeric(gens*cons), dim = c(gens,cons)) fit.sde<-array(data = numeric(gens*cons), dim = c(gens,cons)) for(x in 1:gens) { for(y in 1:cons) { fit.mean[x,y]=mean(toss.outlier(FIT[x,y,1:reps])) fit.sde[x,y]=(sd(toss.outlier(FIT[x,y,1:reps]))/sqrt(reps-1)) } } windows() setwd(wd) RifCon.df<-read.table("RifConc.txt", header=TRUE) con.val<-RifCon.df[RifCon.df$treatment=="Moderate",]$rif.conc[c(1,2,4,6,8,10,12,14,16,18)] axissize<-1.2 #size of text on axes labsize<-1.5 #size of text on axis labels par(mar=c(5,5,4,2)) plot(c(1,length(con.val)),c(min(fit.mean)-max(fit.sde),max(fit.mean)+max(fit.sde)),col="white",xaxt="n",yaxt="n", xlab=expression(paste("rifampicin concentration (",mu,"g/mL)",sep="")),ylab=expression(paste("growth rate (",10^-4," AU/min)",sep="")),cex.lab=labsize, cex.axis=2) axis(1, at=1:length(con), labels=con.val,cex.axis=axissize) axis(2, at=c(0,0.0002,0.0004,0.0006,0.0008,0.001), labels=c(0,2,4,6,8,10), las=1, cex.axis=axissize) for(x in c(2,4,5,6)){ polygon(c(1:cons,cons:1),c(fit.mean[x,1:cons]-fit.sde[x,1:cons],rev(fit.mean[x,1:cons]+fit.sde[x,1:cons])),col=gensdcol[x],border=FALSE) lines(1:cons,fit.mean[x,1:cons],col=gencol[x],lwd=2) points(1:cons,fit.mean[x,1:cons],col=gencol[x],pch=19) } text(2,0.0002,"ancestors",col=gencol[1]) text(2,0.00015,gen[4],col=gencol[4]) text(2,0.0001,gen[5],col=gencol[5]) text(2,0.00005,gen[6],col=gencol[6]) #End of Supp Fig 3c------------------------------------------------ #------------------------------------------------------------------ #****************************************************************** #****************************************************************** #Supplementary Figure 3b,d: Generating Mathematica text files #------------------------------------------------------------------ setwd(wd) viab.df<-read.table("Viability.txt", header=TRUE) reps<-6 #must be even crit.frac<-(1/6) con<-c("c1","c2","c4","c6","c8","c10","c12","c14","c16","c18") cons<-length(con) gen358<-c("t437a.anc","a1685c.anc","t437a.a1685c.anc","t437a","a1685c","t437a.a1685c") gen367<-c("t1532g.anc","g1546a.anc","t1532g.g1546a.anc","t1532g","g1546a","t1532g.g1546a") PV<-array(data = numeric(4*cons), dim = c(4,cons)) EV<-array(data = numeric(4*cons), dim = c(4,cons)) rownames(PV)<-c("ab","Ab","aB","AB") rownames(EV)<-c("ab","Ab","aB","AB") colnames(PV)<-con colnames(EV)<-con gen<-gen358 lin<-"M358" gens<-length(gen) PERSIST<-array(data = numeric(gens*cons), dim = c(gens,cons)) EMERGE<-array(data = numeric(gens*cons), dim = c(gens,cons)) for(x in 1:gens) { for(y in 1:cons) { geno<-gen[x] conc<-con[y] PERSIST[x,y]<-sum(viab.df[viab.df$lineage==lin & viab.df$strain==geno & viab.df$concentration==conc & viab.df$dilution==0.025,]$growth)/reps EMERGE[x,y]<-sum(viab.df[viab.df$lineage==lin & viab.df$strain==geno & viab.df$concentration==conc & viab.df$dilution==0.00001,]$growth)/reps } } for(x in 1:4) { for(y in 1:cons) { if(x==1) { PV[x,y]<-(PERSIST[1,y]+PERSIST[2,y]+PERSIST[3,y])/3 EV[x,y]<-(EMERGE[1,y]+EMERGE[2,y]+EMERGE[3,y])/3 } if(x>1) { PV[x,y]<-PERSIST[x+2,y] EV[x,y]<-EMERGE[x+2,y] } if(PV[x,y]<=crit.frac) { PV[x,y]<-0 } if(EV[x,y]<=crit.frac) { EV[x,y]<-0 } } } PV358<-PV EV358<-EV write.table(PV358,"PV358.txt",row.names=FALSE,col.names=FALSE) gen<-gen367 lin<-"M367" gens<-length(gen) PERSIST<-array(data = numeric(gens*cons), dim = c(gens,cons)) EMERGE<-array(data = numeric(gens*cons), dim = c(gens,cons)) for(x in 1:gens) { for(y in 1:cons) { geno<-gen[x] conc<-con[y] PERSIST[x,y]<-sum(viab.df[viab.df$lineage==lin & viab.df$strain==geno & viab.df$concentration==conc & viab.df$dilution==0.025,]$growth)/reps EMERGE[x,y]<-sum(viab.df[viab.df$lineage==lin & viab.df$strain==geno & viab.df$concentration==conc & viab.df$dilution==0.00001,]$growth)/reps } } for(x in 1:4) { for(y in 1:cons) { if(x==1) { PV[x,y]<-(PERSIST[1,y]+PERSIST[2,y]+PERSIST[3,y])/3 EV[x,y]<-(EMERGE[1,y]+EMERGE[2,y]+EMERGE[3,y])/3 } if(x>1) { PV[x,y]<-PERSIST[x+2,y] EV[x,y]<-EMERGE[x+2,y] } if(PV[x,y]<=crit.frac) { PV[x,y]<-0 } if(EV[x,y]<=crit.frac) { EV[x,y]<-0 } } } PV367<-PV EV367<-EV write.table(PV367,"PV367.txt",row.names=FALSE,col.names=FALSE) edge<-c("ab.Ab","ab.aB","Ab.AB","aB.AB") edges<-length(edge) EdgeV<-array(data = numeric(edges*cons), dim = c(edges,cons)) rownames(EdgeV)<-edge colnames(EdgeV)<-con PV<-PV358 EV<-EV358 for(x in edge) { for(y in 1:cons) { if(x=="ab.Ab") { ifelse(PV["ab",y]>crit.frac & PV["Ab",y]>crit.frac, EdgeV[x,y]<-1, ifelse( (EV["ab",y]>crit.frac | EV["Ab",y]>crit.frac), EdgeV[x,y]<-1, EdgeV[x,y]<-0)) } if(x=="ab.aB") { ifelse(PV["ab",y]>crit.frac & PV["aB",y]>crit.frac, EdgeV[x,y]<-1, ifelse( (EV["ab",y]>crit.frac | EV["aB",y]>crit.frac), EdgeV[x,y]<-1, EdgeV[x,y]<-0)) } if(x=="Ab.AB") { ifelse(PV["Ab",y]>crit.frac & PV["AB",y]>crit.frac, EdgeV[x,y]<-1, ifelse( (EV["Ab",y]>crit.frac | EV["AB",y]>crit.frac), EdgeV[x,y]<-1, EdgeV[x,y]<-0)) } if(x=="aB.AB") { ifelse(PV["aB",y]>crit.frac & PV["AB",y]>crit.frac, EdgeV[x,y]<-1, ifelse( (EV["aB",y]>crit.frac | EV["AB",y]>crit.frac), EdgeV[x,y]<-1, EdgeV[x,y]<-0)) } } } EdgeV358<-EdgeV write.table(EdgeV358,"EdgeV358.txt",row.names=FALSE,col.names=FALSE) PV<-PV367 EV<-EV367 for(x in edge) { for(y in 1:cons) { if(x=="ab.Ab") { ifelse(PV["ab",y]>crit.frac & PV["Ab",y]>crit.frac, EdgeV[x,y]<-1, ifelse( (EV["ab",y]>crit.frac | EV["Ab",y]>crit.frac), EdgeV[x,y]<-1, EdgeV[x,y]<-0)) } if(x=="ab.aB") { ifelse(PV["ab",y]>crit.frac & PV["aB",y]>crit.frac, EdgeV[x,y]<-1, ifelse( (EV["ab",y]>crit.frac | EV["aB",y]>crit.frac), EdgeV[x,y]<-1, EdgeV[x,y]<-0)) } if(x=="Ab.AB") { ifelse(PV["Ab",y]>crit.frac & PV["AB",y]>crit.frac, EdgeV[x,y]<-1, ifelse( (EV["Ab",y]>crit.frac | EV["AB",y]>crit.frac), EdgeV[x,y]<-1, EdgeV[x,y]<-0)) } if(x=="aB.AB") { ifelse(PV["aB",y]>crit.frac & PV["AB",y]>crit.frac, EdgeV[x,y]<-1, ifelse( (EV["aB",y]>crit.frac | EV["AB",y]>crit.frac), EdgeV[x,y]<-1, EdgeV[x,y]<-0)) } } } EdgeV367<-EdgeV write.table(EdgeV367,"EdgeV367.txt",row.names=FALSE,col.names=FALSE) setwd(wd) comp358.df<-read.table("M358EdgeCompetitions.txt", header=TRUE) comp367.df<-read.table("M367EdgeCompetitions.txt", header=TRUE) con<-c(1,2,4,6,8,10,12,14,16,18) cons<-length(con) comp.df<-comp358.df gen<-gen358 EdgeV<-EdgeV358 reps<-4 #Create the fitness arrays FIT<-array(data = numeric(edges*cons*reps), dim = c(edges,cons,reps)) for(x in 1:edges) { for(y in 1:cons) { edg<-edge[x] conc<-con[y] if(edg=="ab.Ab") { str1=gen[1] str2=gen[4] } if(edg=="ab.aB") { str1=gen[2] str2=gen[5] } if(edg=="Ab.AB") { str1=gen[4] str2=gen[6] } if(edg=="aB.AB") { str1=gen[5] str2=gen[6] } cur.comp0<-comp.df[comp.df$strain1==str1 & comp.df$strain2==str2 & comp.df$concentration==conc & comp.df$time==0,] cur.comp48<-comp.df[comp.df$strain1==str1 & comp.df$strain2==str2 & comp.df$concentration==conc & comp.df$time==48,] cfu1i=(cur.comp0$count1)*(10^(-cur.comp0$dilution1+1))*(1/80) cfu2i=(cur.comp0$count2)*(10^(-cur.comp0$dilution2+1))*(1/80) cfu1f=(cur.comp48$count1)*(10^(-cur.comp48$dilution1+1)) cfu2f=(cur.comp48$count2)*(10^(-cur.comp48$dilution2+1)) fit<-(cfu1f/cfu1i)/(cfu2f/cfu2i) FIT[x,y,1]<-fit[1] FIT[x,y,2]<-fit[2] FIT[x,y,3]<-fit[3] FIT[x,y,4]<-fit[4] } } EdgeDir<-array(data = numeric(edges*cons), dim = c(edges,cons)) rownames(EdgeDir)<-edge colnames(EdgeDir)<-con corrected.alpha<-0.05/sum(EdgeV) for(x in 1:edges) { for(y in 1:cons) { ifelse(EdgeV[x,y]==1, ifelse(t.test(FIT[x,y,1:reps],mu=1,alternative="two.sided")[3]$p.value1,EdgeDir[x,y]<-(-1),EdgeDir[x,y]<-1), EdgeDir[x,y]<-0), EdgeDir[x,y]<-0) } } EdgeDir358<-EdgeDir write.table(EdgeDir358,"EdgeDir358.txt",row.names=FALSE,col.names=FALSE) comp.df<-comp367.df gen<-gen367 EdgeV<-EdgeV367 reps<-4 #Create the fitness arrays FIT<-array(data = numeric(edges*cons*reps), dim = c(edges,cons,reps)) for(x in 1:edges) { for(y in 1:cons) { edg<-edge[x] conc<-con[y] if(edg=="ab.Ab") { str1=gen[1] str2=gen[4] } if(edg=="ab.aB") { str1=gen[2] str2=gen[5] } if(edg=="Ab.AB") { str1=gen[4] str2=gen[6] } if(edg=="aB.AB") { str1=gen[5] str2=gen[6] } cur.comp0<-comp.df[comp.df$strain1==str1 & comp.df$strain2==str2 & comp.df$concentration==conc & comp.df$time==0,] cur.comp48<-comp.df[comp.df$strain1==str1 & comp.df$strain2==str2 & comp.df$concentration==conc & comp.df$time==48,] cfu1i=(cur.comp0$count1)*(10^(-cur.comp0$dilution1+1))*(1/80) cfu2i=(cur.comp0$count2)*(10^(-cur.comp0$dilution2+1))*(1/80) cfu1f=(cur.comp48$count1)*(10^(-cur.comp48$dilution1+1)) cfu2f=(cur.comp48$count2)*(10^(-cur.comp48$dilution2+1)) fit<-(cfu1f/cfu1i)/(cfu2f/cfu2i) FIT[x,y,1]<-fit[1] FIT[x,y,2]<-fit[2] FIT[x,y,3]<-fit[3] FIT[x,y,4]<-fit[4] } } EdgeDir<-array(data = numeric(edges*cons), dim = c(edges,cons)) rownames(EdgeDir)<-edge colnames(EdgeDir)<-con corrected.alpha<-0.05/sum(EdgeV) for(x in 1:edges) { for(y in 1:cons) { ifelse(EdgeV[x,y]==1, ifelse(t.test(FIT[x,y,1:reps],mu=1,alternative="two.sided")[3]$p.value1,EdgeDir[x,y]<-(-1),EdgeDir[x,y]<-1), EdgeDir[x,y]<-0), EdgeDir[x,y]<-0) } } EdgeDir367<-EdgeDir write.table(EdgeDir367,"EdgeDir367.txt",row.names=FALSE,col.names=FALSE) #End of Sup Fig 4b,d----------------------------------------------- #------------------------------------------------------------------ #****************************************************************** #****************************************************************** #Supplemental Figure 4a: Growth rates for S497 #------------------------------------------------------------------ setwd(wd) spec.df<-read.table("Sudden_24h.txt", header=TRUE) gen<-c("g1546t.anc","g1546t") wd2<-paste(wd,"/S497_PDFS" ,sep="") setwd(wd2) con<-c("c1","c2","c4","c6","c8","c10","c12","c14","c16","c18") gencol<-c("black","green4") gensdcol<-c("lightgray","lightgreen") reps<-4 #must be even window<-5 #must be odd period<-30 cons<-length(con) gens<-length(gen) #FUNCTION: posify #This function changes all non-positive values in a #vector to the smallest positive value in the vector #or to a set positive value if all entries are #non-positive. #--------------------------------------------------- posify<-function(x,emer.value=0.0001) { y<-numeric(length(x)) smallest.strict.pos<-max(x) if(smallest.strict.pos<=0) { y<-rep(emer.value,length(x)) } if(smallest.strict.pos>0) { for(i in 1:length(x)) { if(x[i]>0 & x[i](maxval-meanval),which.min(x),which.max(x)) if(ind==1) { y<-x[2:length(x)] return(y) } if(ind==(length(x))) { y<-x[1:(length(x)-1)] return(y) } y<-c(x[1:(ind-1)],x[(ind+1):length(x)]) return(y) } #--------------------------------------------------- #Create the fitness arrays fit1<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) fit2<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) fit3<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time1<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time2<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time3<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) FIT<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) for(x in 1:gens) { for(y in 1:cons) { for(z in 1:reps) { geno<-gen[x] conc<-con[y] t<-spec.df[spec.df$rep==z & spec.df$strain==geno,]$time od<-spec.df[spec.df$rep==z & spec.df$strain==geno,3+y] od<-posify(od) times<-numeric(length(od)-window+1) slopes<-numeric(length(od)-window+1) for(i in (1+((window-1)/2)):(length(od)-((window-1)/2))){ slopes[i-(window-1)/2]<-lsfit(t[(i-(window-1)/2):(i+(window-1)/2)],od[(i-(window-1)/2):(i+(window-1)/2)])$coefficients[2] times[i-(window-1)/2]<-t[i] } fit1[x,y,z]<-max(slopes) tind<-which.max(slopes) time1[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) fit2[x,y,z]<-max(slopes) tind<-which.max(slopes) time2[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) fit3[x,y,z]<-max(slopes) tind<-which.max(slopes) time3[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) FIT[x,y,z]=(fit1[x,y,z]+fit2[x,y,z]+fit3[x,y,z])/3.0 if(z==1) { dev.set(3) pdf(paste(gen[x],"___",con[y],".pdf", sep="")) par(mfrow=c(2,reps/2)) } plot(times,(od[(1+((window-1)/2)):(length(od)-((window-1)/2))]),xlab="time (min)",ylab="OD",main=paste(gen[x],":",con[y],":",z,"\n","slope=",FIT[x,y,z]),col="white") points(times,(od[(1+((window-1)/2)):(length(od)-((window-1)/2))]),pch=19) points(time1[x,y,z],(od[(time1[x,y,z]/period)+1]),col="red",cex=2) points(time1[x,y,z],(od[(time1[x,y,z]/period)+1]),col="red",cex=3) lines(times,fit1[x,y,z]*(times-time1[x,y,z])+(od[(time1[x,y,z]/period)+1]),lwd=2,col="red") points(time2[x,y,z],(od[(time2[x,y,z]/period)+1]),col="red",cex=2) points(time2[x,y,z],(od[(time2[x,y,z]/period)+1]),col="red",cex=3) lines(times,fit2[x,y,z]*(times-time2[x,y,z])+(od[(time2[x,y,z]/period)+1]),lwd=2,col="red") points(time3[x,y,z],(od[(time3[x,y,z]/period)+1]),col="red",cex=2) points(time3[x,y,z],(od[(time3[x,y,z]/period)+1]),col="red",cex=3) lines(times,fit3[x,y,z]*(times-time3[x,y,z])+(od[(time3[x,y,z]/period)+1]),lwd=2,col="red") if(z==reps) { dev.off() } } } } fit.mean<-array(data = numeric(gens*cons), dim = c(gens,cons)) fit.sde<-array(data = numeric(gens*cons), dim = c(gens,cons)) for(x in 1:gens) { for(y in 1:cons) { fit.mean[x,y]=mean(toss.outlier(FIT[x,y,1:reps])) fit.sde[x,y]=(sd(toss.outlier(FIT[x,y,1:reps]))/sqrt(reps-1)) } } windows() setwd(wd) RifCon.df<-read.table("RifConc.txt", header=TRUE) con.val<-RifCon.df[RifCon.df$treatment=="Moderate",]$rif.conc[c(1,2,4,6,8,10,12,14,16,18)] axissize<-1.2 #size of text on axes labsize<-1.5 #size of text on axis labels par(mar=c(5,5,4,2)) plot(c(1,length(con.val)),c(min(fit.mean)-max(fit.sde),max(fit.mean)+max(fit.sde)),col="white",xaxt="n",yaxt="n", xlab=expression(paste("rifampicin concentration (",mu,"g/mL)",sep="")),ylab=expression(paste("growth rate (",10^-4," AU/min)",sep="")),cex.lab=labsize, cex.axis=2) axis(1, at=1:length(con), labels=con.val,cex.axis=axissize) axis(2, at=c(0,0.0002,0.0004,0.0006,0.0008,0.001), labels=c(0,2,4,6,8,10), las=1, cex.axis=axissize) for(x in c(1,2)){ polygon(c(1:cons,cons:1),c(fit.mean[x,1:cons]-fit.sde[x,1:cons],rev(fit.mean[x,1:cons]+fit.sde[x,1:cons])),col=gensdcol[x],border=FALSE) lines(1:cons,fit.mean[x,1:cons],col=gencol[x],lwd=2) points(1:cons,fit.mean[x,1:cons],col=gencol[x],pch=19) } text(2,0.0002,"ancestor",col=gencol[1]) text(2,0.00015,gen[2],col=gencol[2]) #End of Supp Fig 4a------------------------------------------------ #------------------------------------------------------------------ #****************************************************************** #****************************************************************** #Supplemental Figure 4b: Growth rates for S1236 #------------------------------------------------------------------ setwd(wd) spec.df<-read.table("Sudden_24h.txt", header=TRUE) gen<-c("g436t.anc","g436t") wd2<-paste(wd,"/S1236_PDFS",sep="") setwd(wd2) con<-c("c1","c2","c4","c6","c8","c10","c12","c14","c16","c18") gencol<-c("black","green4") gensdcol<-c("lightgray","lightgreen") reps<-4 #must be even window<-5 #must be odd period<-30 cons<-length(con) gens<-length(gen) #FUNCTION: posify #This function changes all non-positive values in a #vector to the smallest positive value in the vector #or to a set positive value if all entries are #non-positive. #--------------------------------------------------- posify<-function(x,emer.value=0.0001) { y<-numeric(length(x)) smallest.strict.pos<-max(x) if(smallest.strict.pos<=0) { y<-rep(emer.value,length(x)) } if(smallest.strict.pos>0) { for(i in 1:length(x)) { if(x[i]>0 & x[i](maxval-meanval),which.min(x),which.max(x)) if(ind==1) { y<-x[2:length(x)] return(y) } if(ind==(length(x))) { y<-x[1:(length(x)-1)] return(y) } y<-c(x[1:(ind-1)],x[(ind+1):length(x)]) return(y) } #--------------------------------------------------- #Create the fitness arrays fit1<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) fit2<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) fit3<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time1<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time2<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) time3<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) FIT<-array(data = numeric(gens*cons*reps), dim = c(gens,cons,reps)) for(x in 1:gens) { for(y in 1:cons) { for(z in 1:reps) { geno<-gen[x] conc<-con[y] t<-spec.df[spec.df$rep==z & spec.df$strain==geno,]$time od<-spec.df[spec.df$rep==z & spec.df$strain==geno,3+y] od<-posify(od) times<-numeric(length(od)-window+1) slopes<-numeric(length(od)-window+1) for(i in (1+((window-1)/2)):(length(od)-((window-1)/2))){ slopes[i-(window-1)/2]<-lsfit(t[(i-(window-1)/2):(i+(window-1)/2)],od[(i-(window-1)/2):(i+(window-1)/2)])$coefficients[2] times[i-(window-1)/2]<-t[i] } fit1[x,y,z]<-max(slopes) tind<-which.max(slopes) time1[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) fit2[x,y,z]<-max(slopes) tind<-which.max(slopes) time2[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) fit3[x,y,z]<-max(slopes) tind<-which.max(slopes) time3[x,y,z]<-times[tind] slopes[tind]<-(min(slopes)-1) FIT[x,y,z]=(fit1[x,y,z]+fit2[x,y,z]+fit3[x,y,z])/3.0 if(z==1) { dev.set(3) pdf(paste(gen[x],"___",con[y],".pdf", sep="")) par(mfrow=c(2,reps/2)) } plot(times,(od[(1+((window-1)/2)):(length(od)-((window-1)/2))]),xlab="time (min)",ylab="OD",main=paste(gen[x],":",con[y],":",z,"\n","slope=",FIT[x,y,z]),col="white") points(times,(od[(1+((window-1)/2)):(length(od)-((window-1)/2))]),pch=19) points(time1[x,y,z],(od[(time1[x,y,z]/period)+1]),col="red",cex=2) points(time1[x,y,z],(od[(time1[x,y,z]/period)+1]),col="red",cex=3) lines(times,fit1[x,y,z]*(times-time1[x,y,z])+(od[(time1[x,y,z]/period)+1]),lwd=2,col="red") points(time2[x,y,z],(od[(time2[x,y,z]/period)+1]),col="red",cex=2) points(time2[x,y,z],(od[(time2[x,y,z]/period)+1]),col="red",cex=3) lines(times,fit2[x,y,z]*(times-time2[x,y,z])+(od[(time2[x,y,z]/period)+1]),lwd=2,col="red") points(time3[x,y,z],(od[(time3[x,y,z]/period)+1]),col="red",cex=2) points(time3[x,y,z],(od[(time3[x,y,z]/period)+1]),col="red",cex=3) lines(times,fit3[x,y,z]*(times-time3[x,y,z])+(od[(time3[x,y,z]/period)+1]),lwd=2,col="red") if(z==reps) { dev.off() } } } } fit.mean<-array(data = numeric(gens*cons), dim = c(gens,cons)) fit.sde<-array(data = numeric(gens*cons), dim = c(gens,cons)) for(x in 1:gens) { for(y in 1:cons) { fit.mean[x,y]=mean(toss.outlier(FIT[x,y,1:reps])) fit.sde[x,y]=(sd(toss.outlier(FIT[x,y,1:reps]))/sqrt(reps-1)) } } windows() setwd(wd) RifCon.df<-read.table("RifConc.txt", header=TRUE) con.val<-RifCon.df[RifCon.df$treatment=="Moderate",]$rif.conc[c(1,2,4,6,8,10,12,14,16,18)] axissize<-1.2 #size of text on axes labsize<-1.5 #size of text on axis labels par(mar=c(5,5,4,2)) plot(c(1,length(con.val)),c(min(fit.mean)-max(fit.sde),max(fit.mean)+max(fit.sde)),col="white",xaxt="n",yaxt="n", xlab=expression(paste("rifampicin concentration (",mu,"g/mL)",sep="")),ylab=expression(paste("growth rate (",10^-4," AU/min)",sep="")),cex.lab=labsize, cex.axis=2) axis(1, at=1:length(con), labels=con.val,cex.axis=axissize) axis(2, at=c(0,0.0002,0.0004,0.0006,0.0008,0.001), labels=c(0,2,4,6,8,10), las=1, cex.axis=axissize) for(x in c(1,2)){ polygon(c(1:cons,cons:1),c(fit.mean[x,1:cons]-fit.sde[x,1:cons],rev(fit.mean[x,1:cons]+fit.sde[x,1:cons])),col=gensdcol[x],border=FALSE) lines(1:cons,fit.mean[x,1:cons],col=gencol[x],lwd=2) points(1:cons,fit.mean[x,1:cons],col=gencol[x],pch=19) } text(2,0.0002,"ancestor",col=gencol[1]) text(2,0.00015,gen[2],col=gencol[2]) #End of Supp Fig 4b------------------------------------------------ #------------------------------------------------------------------ #****************************************************************** #****************************************************************** #Supplementary Figure 4c: Generating Mathematica text files #------------------------------------------------------------------ setwd(wd) viab.df<-read.table("Viability.txt", header=TRUE) reps<-6 #must be even crit.frac<-(1/6) con<-c("c1","c2","c4","c6","c8","c10","c12","c14","c16","c18") cons<-length(con) gen497<-c("g1546t.anc","g1546t") gen1236<-c("g436t.anc","g436t") PV<-array(data = numeric(2*cons), dim = c(2,cons)) EV<-array(data = numeric(2*cons), dim = c(2,cons)) rownames(PV)<-c("a","A") rownames(EV)<-c("a","A") colnames(PV)<-con colnames(EV)<-con gen<-gen497 lin<-"S497" gens<-length(gen) PERSIST<-array(data = numeric(gens*cons), dim = c(gens,cons)) EMERGE<-array(data = numeric(gens*cons), dim = c(gens,cons)) for(x in 1:gens) { for(y in 1:cons) { geno<-gen[x] conc<-con[y] PERSIST[x,y]<-sum(viab.df[viab.df$lineage==lin & viab.df$strain==geno & viab.df$concentration==conc & viab.df$dilution==0.025,]$growth)/reps EMERGE[x,y]<-sum(viab.df[viab.df$lineage==lin & viab.df$strain==geno & viab.df$concentration==conc & viab.df$dilution==0.00001,]$growth)/reps } } for(x in 1:2) { for(y in 1:cons) { PV[x,y]<-PERSIST[x,y] EV[x,y]<-EMERGE[x,y] if(PV[x,y]<=crit.frac) { PV[x,y]<-0 } if(EV[x,y]<=crit.frac) { EV[x,y]<-0 } } } PV497<-PV EV497<-EV write.table(PV497,"PV497.txt",row.names=FALSE,col.names=FALSE) gen<-gen1236 lin<-"S1236" gens<-length(gen) PERSIST<-array(data = numeric(gens*cons), dim = c(gens,cons)) EMERGE<-array(data = numeric(gens*cons), dim = c(gens,cons)) for(x in 1:gens) { for(y in 1:cons) { geno<-gen[x] conc<-con[y] PERSIST[x,y]<-sum(viab.df[viab.df$lineage==lin & viab.df$strain==geno & viab.df$concentration==conc & viab.df$dilution==0.025,]$growth)/reps EMERGE[x,y]<-sum(viab.df[viab.df$lineage==lin & viab.df$strain==geno & viab.df$concentration==conc & viab.df$dilution==0.00001,]$growth)/reps } } for(x in 1:2) { for(y in 1:cons) { PV[x,y]<-PERSIST[x,y] EV[x,y]<-EMERGE[x,y] if(PV[x,y]<=crit.frac) { PV[x,y]<-0 } if(EV[x,y]<=crit.frac) { EV[x,y]<-0 } } } PV1236<-PV EV1236<-EV write.table(PV1236,"PV1236.txt",row.names=FALSE,col.names=FALSE) edge<-c("a.A") edges<-length(edge) EdgeV<-array(data = numeric(edges*cons), dim = c(edges,cons)) rownames(EdgeV)<-edge colnames(EdgeV)<-con PV<-PV497 EV<-EV497 for(x in edge) { for(y in 1:cons) { ifelse(PV["a",y]>crit.frac & PV["A",y]>crit.frac, EdgeV[x,y]<-1, ifelse( (EV["a",y]>crit.frac | EV["A",y]>crit.frac), EdgeV[x,y]<-1, EdgeV[x,y]<-0 ) ) } } EdgeV497<-EdgeV write.table(EdgeV497,"EdgeV497.txt",row.names=FALSE,col.names=FALSE) PV<-PV1236 EV<-EV1236 for(x in edge) { for(y in 1:cons) { ifelse(PV["a",y]>crit.frac & PV["A",y]>crit.frac, EdgeV[x,y]<-1, ifelse( (EV["a",y]>crit.frac | EV["A",y]>crit.frac), EdgeV[x,y]<-1, EdgeV[x,y]<-0 ) ) } } EdgeV1236<-EdgeV write.table(EdgeV1236,"EdgeV1236.txt",row.names=FALSE,col.names=FALSE) setwd(wd) compS.df<-read.table("SEdgeCompetitions.txt", header=TRUE) con<-c(1,2,4,6,8,10,12,14,16,18) cons<-length(con) comp.df<-compS.df gen<-gen497 EdgeV<-EdgeV497 reps<-4 #Create the fitness arrays FIT<-array(data = numeric(edges*cons*reps), dim = c(edges,cons,reps)) for(x in 1:edges) { for(y in 1:cons) { edg<-edge[x] conc<-con[y] if(edg=="a.A") { str1=gen[1] str2=gen[2] } cur.comp0<-comp.df[comp.df$strain1==str1 & comp.df$strain2==str2 & comp.df$concentration==conc & comp.df$time==0,] cur.comp48<-comp.df[comp.df$strain1==str1 & comp.df$strain2==str2 & comp.df$concentration==conc & comp.df$time==48,] cfu1i=(cur.comp0$count1)*(10^(-cur.comp0$dilution1+1))*(1/80) cfu2i=(cur.comp0$count2)*(10^(-cur.comp0$dilution2+1))*(1/80) cfu1f=(cur.comp48$count1)*(10^(-cur.comp48$dilution1+1)) cfu2f=(cur.comp48$count2)*(10^(-cur.comp48$dilution2+1)) fit<-(cfu1f/cfu1i)/(cfu2f/cfu2i) FIT[x,y,1]<-fit[1] FIT[x,y,2]<-fit[2] FIT[x,y,3]<-fit[3] FIT[x,y,4]<-fit[4] } } EdgeDir<-array(data = numeric(edges*cons), dim = c(edges,cons)) rownames(EdgeDir)<-edge colnames(EdgeDir)<-con corrected.alpha<-0.05/sum(EdgeV) for(x in 1:edges) { for(y in 1:cons) { ifelse(EdgeV[x,y]==1, ifelse(t.test(FIT[x,y,1:reps],mu=1,alternative="two.sided")[3]$p.value1,EdgeDir[x,y]<-(-1),EdgeDir[x,y]<-1), EdgeDir[x,y]<-0), EdgeDir[x,y]<-0) } } EdgeDir497<-EdgeDir write.table(EdgeDir497,"EdgeDir497.txt",row.names=FALSE,col.names=FALSE) comp.df<-compS.df gen<-gen1236 EdgeV<-EdgeV1236 reps<-4 #Create the fitness arrays FIT<-array(data = numeric(edges*cons*reps), dim = c(edges,cons,reps)) for(x in 1:edges) { for(y in 1:cons) { edg<-edge[x] conc<-con[y] if(edg=="a.A") { str1=gen[1] str2=gen[2] } cur.comp0<-comp.df[comp.df$strain1==str1 & comp.df$strain2==str2 & comp.df$concentration==conc & comp.df$time==0,] cur.comp48<-comp.df[comp.df$strain1==str1 & comp.df$strain2==str2 & comp.df$concentration==conc & comp.df$time==48,] cfu1i=(cur.comp0$count1)*(10^(-cur.comp0$dilution1+1))*(1/80) cfu2i=(cur.comp0$count2)*(10^(-cur.comp0$dilution2+1))*(1/80) cfu1f=(cur.comp48$count1)*(10^(-cur.comp48$dilution1+1)) cfu2f=(cur.comp48$count2)*(10^(-cur.comp48$dilution2+1)) fit<-(cfu1f/cfu1i)/(cfu2f/cfu2i) FIT[x,y,1]<-fit[1] FIT[x,y,2]<-fit[2] FIT[x,y,3]<-fit[3] FIT[x,y,4]<-fit[4] } } EdgeDir<-array(data = numeric(edges*cons), dim = c(edges,cons)) rownames(EdgeDir)<-edge colnames(EdgeDir)<-con corrected.alpha<-0.05/sum(EdgeV) for(x in 1:edges) { for(y in 1:cons) { ifelse(EdgeV[x,y]==1, ifelse(t.test(FIT[x,y,1:reps],mu=1,alternative="two.sided")[3]$p.value1,EdgeDir[x,y]<-(-1),EdgeDir[x,y]<-1), EdgeDir[x,y]<-0), EdgeDir[x,y]<-0) } } EdgeDir1236<-EdgeDir write.table(EdgeDir1236,"EdgeDir1236.txt",row.names=FALSE,col.names=FALSE) #End of Sup Fig 4c------------------------------------------------- #------------------------------------------------------------------ #******************************************************************