library(gplots) # get averages and standard deviations for phage timecourse data dataTablePhage <- read.table( paste('experimentalData','phageTimecourse20111110-11_M9succinate_M9glucose.txt', sep= .Platform$file.sep), sep='\t', header = TRUE) dataTablePhageAcetate <- read.table( paste('experimentalData','phageTimecourse20120223_Acetate.txt', sep= .Platform$file.sep), sep='\t', header = TRUE) dataTablePhageRich <- read.table( paste('experimentalData','phageTimecourse_20110826_TB.txt', sep= .Platform$file.sep), sep='\t', header = TRUE) dataTablePhageRich <- subset( dataTablePhageRich, time < 0.28 ) dataTablePhage$G1[2] sdPhageGlucose <- 0 avePhageGlucose <- 0 for (i in 1:length(dataTablePhage$G1)) { sdPhageGlucose[i] <- sd( c(dataTablePhage$G1[i], dataTablePhage$G2[i], dataTablePhage$G3[i]) ) avePhageGlucose[i] <- mean( c(dataTablePhage$G1[i], dataTablePhage$G2[i], dataTablePhage$G3[i]) ) } sdPhageSuccinate <- 0 avePhageSuccinate <- 0 for (i in 1:length(dataTablePhage$S1)) { sdPhageSuccinate[i] <- sd( c(dataTablePhage$S1[i], dataTablePhage$S2[i], dataTablePhage$S3[i]) ) avePhageSuccinate[i] <- mean( c(dataTablePhage$S1[i], dataTablePhage$S2[i], dataTablePhage$S3[i]) ) } sdPhageTryptone <- 0 avePhageTryptone <- 0 for (i in 1:length(dataTablePhageRich$phage1)) { sdPhageTryptone[i] <- sd( c(dataTablePhageRich$phage1[i], dataTablePhageRich$phage2[i], dataTablePhageRich$phage3[i]) ) avePhageTryptone[i] <- mean( c(dataTablePhageRich$phage1[i], dataTablePhageRich$phage2[i], dataTablePhageRich$phage3[i]) ) } sdPhageAcetate <- 0 avePhageAcetate <- 0 for (i in 1:length(dataTablePhageAcetate$t)) { sdPhageAcetate[i] <- sd( c(dataTablePhageAcetate$A1[i], dataTablePhageAcetate$A2[i], dataTablePhageAcetate$A3[i]) ) avePhageAcetate[i] <- mean( c(dataTablePhageAcetate$A1[i], dataTablePhageAcetate$A2[i], dataTablePhageAcetate$A3[i]) ) } tComp <- 17 muData <- c(0.27,0.45,0.66,1.5) phageData <- c(0,0,0,0) phageSD <- c(0,0,0,0) # Acetate for bar plot idxMedia <- 1 idxAcetateData <- which.min( abs(dataTablePhageAcetate$t-tComp) ) phageData[idxMedia] <- avePhageAcetate[idxAcetateData] phageSD[idxMedia] <- sdPhageAcetate[idxAcetateData] # Succinate for bar plot idxMedia <- 2 idxSuccinateData <- which.min( abs(dataTablePhage$time-tComp) ) phageData[idxMedia] <- avePhageSuccinate[idxSuccinateData] phageSD[idxMedia] <- sdPhageSuccinate[idxSuccinateData] # GLucose for bar plot idxMedia <- 3 idxGlucoseData <- which.min( abs(dataTablePhage$time-tComp) ) phageData[idxMedia] <- avePhageGlucose[idxGlucoseData] phageSD[idxMedia] <- sdPhageGlucose[idxGlucoseData] # Tryptone for bar plot idxMedia <- 4 idxTryptoneData <- which.min( abs(dataTablePhageRich$time-tComp) ) phageData[idxMedia] <- avePhageTryptone[idxTryptoneData] phageSD[idxMedia] <- sdPhageTryptone[idxTryptoneData] # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - folderString <- paste( 'simulationsText', 'fluxBoundsSet_metabolicLimit' , sep= .Platform$file.sep) files <- list.files( path = folderString, pattern = ".txt" ) nFile <- length( files ) metLimitMus <- rep(0,nFile) metLimitPhage <- rep(0,nFile) for( iFile in 1:nFile ){ dataTableTemp <- read.table(paste(folderString,files[iFile],sep= .Platform$file.sep), sep='\t', header = TRUE) dataTableTemp <- subset( dataTableTemp, hostFlag == 1 ) idxDataTemp <- which.min( abs(dataTableTemp$t*60-tComp) ) metLimitMus[iFile] <- dataTableTemp$bBIOMASS[2] metLimitPhage[iFile] <- dataTableTemp$CONCvT7[idxDataTemp] } # corners of region # # if convex # metLimitMus <- c( metLimitMus, min(metLimitMus), max(metLimitMus) ) # metLimitPhage <- c( metLimitPhage, 0, 0 ) # metLimitRegion <- chull( x = metLimitMus, y = metLimitPhage ) # xMetLimit <- metLimitMus[metLimitRegion] # yMetLimit <- metLimitPhage[metLimitRegion] # smoothing fSmooth <- 0.1 smoothed <- lowess(x=metLimitMus,y=metLimitPhage,f=fSmooth) xMetLimit <- c(min(metLimitMus), smoothed$x ,max(metLimitMus)) yMetLimit <- c(0, smoothed$y ,0) # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - folderString <- paste( 'simulationsText', 'fluxBoundsSet_mechineryLimit' , sep= .Platform$file.sep) files <- list.files( path = folderString, pattern = ".txt" ) nFile <- length( files ) intLimitMus <- rep(0,nFile) intLimitPhage <- rep(0,nFile) mechLimitPhage <- rep(0,nFile) for( iFile in 1:nFile ){ dataTableTemp <- read.table(paste(folderString,files[iFile],sep= .Platform$file.sep), sep='\t', header = TRUE) dataTableTemp <- subset( dataTableTemp, hostFlag == 1 ) idxDataTemp <- which.min( abs(dataTableTemp$t*60-tComp) ) intLimitMus[iFile] <- dataTableTemp$bBIOMASS[2] intLimitPhage[iFile] <- dataTableTemp$CONCvT7[idxDataTemp] mechLimitPhage[iFile] <- dataTableTemp$CONCvT7Ref[idxDataTemp] } # corners of region mechLimitSortIdx <- order( intLimitMus ) xMechLimitRegion <- c( min(intLimitMus), intLimitMus[mechLimitSortIdx] ,max(intLimitMus) ) yMechLimitRegion <- c( 0, mechLimitPhage[mechLimitSortIdx] ,0) xInfeasRegion <- c( min(intLimitMus), intLimitMus[mechLimitSortIdx] ) yInfeasRegion <- c( max(mechLimitPhage), mechLimitPhage[mechLimitSortIdx]) # sort integrated data intSortIdx <- order( intLimitMus ) intLimitMus <- intLimitMus[intSortIdx] intLimitPhage <- intLimitPhage[intSortIdx] # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - figureName <- paste('figure7_floatLegend_8o7cm_width_',tComp,'min.eps', sep = '' ) # postscript(figureName, width = 7.0, height = 6.0, horizontal = FALSE, onefile = FALSE, paper = "special") postscript(paste('figureOutput',figureName, sep= .Platform$file.sep), width = 3.27, height = 3.27, horizontal = FALSE, onefile = FALSE, paper = "special", pointsize = 8) par( mar = (c(5, 5, 0, 0) + 0.1)) xAxis <- c(min(intLimitMus), 1.65) yAxis <- c(0, 240) xScale <- seq(from = 0.25, to = 1.5, by = 0.25) yScale <- seq(from = 0, to = 200, by = 50) yString <- paste('Phage T7 Produced (phage/infected host) \n at ', tComp, ' Minutes Post Infection', sep='') plot(x=NULL, y=NULL, xlab = 'Growth Rate (/hour)', xlim = xAxis, ylab = yString, ylim = yAxis, axes = FALSE) # print(xMechLimitRegion) # print(xMetLimit) polygon( x=xInfeasRegion, y =yInfeasRegion, col = 'grey95' , border = NA) polygon( x=xMechLimitRegion, y =yMechLimitRegion, col = 'grey70' , border = NA) # polygon( x=xMetLimit, y = yMetLimit , col = 'olivedrab1', border = NA ) # polygon( x=xMechLimitRegion, y =yMechLimitRegion, col = 'black', density = 10, lwd = 1, angle = 45 ) polygon( x=xMetLimit, y = yMetLimit , col = 'black' , density = 10, lwd = 1, angle = 135 ) # points(x=xMechLimitRegion, y =yMechLimitRegion,pch = 4, col = 'grey30') # points( x=xMetLimit, y = yMetLimit,pch = 4, col = 'olivedrab3') lines(x=intLimitMus,y=intLimitPhage, lwd = 4 , col = 'darkolivegreen3' ) # 'orangered1' plotCI( x = muData, y = phageData, uiw = phageSD, col = 'cornflowerblue' , barcol = 'cornflowerblue' , add = TRUE, cex = 1.5, lwd = 1.5) axis(1, at = xScale ) axis(2, at = yScale) # legend(x=1.75, y=250, legend=c('Mechanistically Infeasible','Mechanistically Feasible','Metabolicly Feasible ','Integrated Simulation','Expirimental Data'), # pch = c(NA, NA, NA, NA, 1), lty = c(NA,NA,NA,1, NA), lwd = c(1,1,1,4,2),col=c(NA, NA ,"black",'grey20','black'), # fill=c('grey90','grey50',NA,NA,NA), bty="n", border = c(NA,NA, 'black', NA, NA), pt.cex = c(4,4,4,NA,1), seg.len = 1, # density = c(NA, NA, 10, NA, NA), angle = c(NA, NA, 135, NA, NA) ) wLeg <- 0.12 hLeg <- 15 yLeg <- max(yInfeasRegion) - 10 - hLeg#235 xLeg <- min(xInfeasRegion)+0.1 #0.25 sLeg <- 1.75 rect(xleft = (xLeg-0.05), ybottom = (yLeg-4*sLeg*hLeg-0.35*hLeg ) , xright= (xLeg+(0.8-0.2)), ytop = (yLeg+1.35*hLeg), col = 'white', border = NA) yLegBot <- yLeg-(0:2)*(sLeg*hLeg) -2 yLegTop <- yLegBot + hLeg + 2 rect(xleft = xLeg, ybottom = yLegBot , xright= (xLeg+wLeg), ytop = yLegTop, density = c(NA, NA, 10), col=c('grey95','grey70','black'),border = c(NA,NA, 'black'), angle = c(NA, NA, 135)) yLegBot <- yLeg-(0:4)*(sLeg*hLeg) yLegTop <- yLegBot + hLeg text(c('Machinery \nInfeasible','Machinery \nFeasible','Metabolically \nFeasible','Integrated \nSimulation','Experimental \nData'), x= (xLeg+wLeg*1.2), y = (yLegBot+ yLegTop)/2,adj = c(0,0.5)) lines(x=c(xLeg,(xLeg+wLeg)),y=rep(((yLegBot[4]+ yLegTop[4])/2),2), lwd = 4 , col = 'darkolivegreen3' ) points(x=c((xLeg+wLeg/2)),y=(yLegBot[5]+ yLegTop[5])/2, col = 'cornflowerblue', cex = 1.5, lwd =1.5 ) dev.off()