### replace with appropriate path setwd("/Users/dominik/Google Drive/ETH/Projects/lambda Rex/dryad") library(dichromat) reads <- read.csv("120524_1.txt", skip = 5, header = FALSE, nrows = 201, sep = "\t", quote="\"", dec=".", fill = TRUE, comment.char="", na.strings = "NA", stringsAsFactors = FALSE)[,3:98] reads <- cbind("time" = 2*(1:nrow(reads))-2, reads) treatm <- c("0R999","0G999","0R99","0G99","0R95","0G95","0R90","0G90","0R50","0G50","0R10","0G10","TR999","TG999","TR99","TG99","TR95","TG95","TR90","TG90","TR50","TG50","TR10","TG10") names(reads) <- c("time",rep(treatm,4)) for(i in 1:24){ cols <- c(i+1,i+25,i+49,i+73) n <- length(cols) temp1 <- rowMeans(reads[,cols]) temp2 <- apply(reads[,cols],1,var) reads <- cbind(reads,temp1,sqrt(temp2/n)) names(reads)[ncol(reads)-1] <- paste(treatm[i]) names(reads)[ncol(reads)] <- paste(treatm[i],"_sem",sep="") } # get correct initial frequencies fd <- read.csv("120524_97-108.csv", header = TRUE, sep = ";", quote="\"", dec=".", fill = TRUE, comment.char="", na.strings = "NA") fd$lambda.color <- as.factor(fd$lambda.color) correction <- T for(i in 1:nrow(fd)){ corrected <- F if(fd$lambda.color[i] == "green"){ fd$lambda.cells[i] <- fd$green.cells[i] fd$HK97.cells[i] <- fd$red.cells[i] }else{ fd$lambda.cells[i] <- fd$red.cells[i] fd$HK97.cells[i] <- fd$green.cells[i] } if(!is.na(fd$HK97.cells[i]) & !is.na(fd$lambda.cells[i])){ if(correction == T){ if(fd$HK97.cells[i] == 0){ fd$HK97.cells[i] <- 1 corrected <- T } if(fd$lambda.cells[i] == 0){ fd$lambda.cells[i] <- 1 corrected <- T } } } fd$corrected[i] <- corrected } fd$freq.lambda <- fd$lambda.cells/(fd$lambda.cells + fd$HK97.cells) fd <- fd[,c(1,5,7:9,11)] before <- fd[order(fd$lambda.color,fd$prop.lambda),] init.freq <- c(round(100*before$freq.lambda[1:6],1),round(100*before$freq.lambda[7:12],1)) # final frequencies fd <- read.csv("120525_exp391.csv", header = TRUE, sep = ";", quote="\"", dec=".", fill = TRUE, comment.char="", na.strings = "NA") fd$lambda.color <- as.factor(fd$lambda.color) correction <- T for(i in 1:nrow(fd)){ corrected <- F if(fd$lambda.color[i] == "green"){ fd$lambda.cells[i] <- fd$green.cells[i] fd$HK97.cells[i] <- fd$red.cells[i] }else{ fd$lambda.cells[i] <- fd$red.cells[i] fd$HK97.cells[i] <- fd$green.cells[i] } if(!is.na(fd$HK97.cells[i]) & !is.na(fd$lambda.cells[i])){ if(correction == T){ if(fd$HK97.cells[i] == 0){ fd$HK97.cells[i] <- 1 corrected <- T } if(fd$lambda.cells[i] == 0){ fd$lambda.cells[i] <- 1 corrected <- T } } } fd$corrected[i] <- corrected fd$extinct[i] <- ifelse(fd$green.cells[i] + fd$red.cells[i] < 1000,TRUE,FALSE) } fd$freq.lambda <- fd$lambda.cells/(fd$lambda.cells + fd$HK97.cells) fd$logit <- log(fd$freq.lambda/(1-fd$freq.lambda)) fd <- fd[order(fd$lambda.color,fd$prop.lambda,fd$T4rII),] summary <- NULL for(i in 1:12){ set.ctrl <- fd$logit[(i*8-7):(i*8-4)] set.treat <- fd$logit[(i*8-3):(i*8)] temp <- t.test(set.ctrl,set.treat,var.equal = T) p.star <- ifelse(temp$p.value <= 0.05,"*","") p.star <- ifelse(temp$p.value <= 0.01,"**",p.star) p.star <- ifelse(temp$p.value <= 0.001,"***",p.star) summary <- rbind(summary,data.frame(lambda.color = fd$lambda.color[i*8],prop.lambda = fd$prop.lambda[i*8],df = round(temp$parameter,2), t = round(temp$statistic,2),p = round(temp$p.value,6), signif = p.star)) } summary <- cbind(summary,data.frame(init.freq = init.freq)) ctrl <- aggregate(fd$freq.lambda[which(fd$T4rII == 0)],list(fd$prop.lambda[fd$T4rII == 0],fd$lambda.color[fd$T4rII == 0],fd$T4rII[fd$T4rII == 0]),mean)[,c(1,2,4)] ctrl <- cbind(ctrl,aggregate(fd$logit[fd$T4rII == 0],list(fd$prop.lambda[fd$T4rII == 0],fd$lambda.color[fd$T4rII == 0],fd$T4rII[fd$T4rII == 0]),mean)[,4]) names(ctrl) <- c("prop.lambda","lambda.color","freq.lambda","logit") summary <- cbind(summary,ctrl[,3:4]) treat <- fd[which(fd$T4rII == 10000),c(3,6,9:14)] for(i in 1:nrow(summary)){ color <- summary$lambda.color[i] prop.lambda <- summary$prop.lambda[i] freq.lambda <- summary$freq.lambda[i] rows <- which(treat$lambda.color == color & treat$prop.lambda == prop.lambda) treat$fitness[rows] <- exp((treat$logit[rows]-summary$logit[i])/10) treat$fitness_alt[rows] <- exp(treat$logit[rows]-summary$logit[i]) treat$fitness_len[rows] <- log(treat$freq.lambda[rows]/summary$freq.lambda[i])/log((1-treat$freq.lambda[rows])/(1-summary$freq.lambda[i])) treat$fitness_gard[rows] <- (treat$freq.lambda[rows]*(1-summary$freq.lambda[i]))/((1-treat$freq.lambda[rows])*(summary$freq.lambda[i])) treat$fitness_m <- log(treat$fitness_gard) } mreads <- reads[,98:145] mreads <- mreads[,c(47,48,43,44,39,40,35,36,31,32,27,28,45,46,41,42,37,38,33,34,29,30,25,26)] mreads <- cbind(mreads,data.frame(time = reads$time)) before <- before[order(rev(before$lambda.color),rev(before$prop.lambda)),] # plot figure width.inch <- 8.4 * 0.393700787402 height.inch <- 13 * 0.393700787402 pdf(file = "figure2.pdf", width = width.inch, height = height.inch, family = "Times", pointsize = 9, colormodel = "rgb") layout(matrix(c(1,2), 2, 1, byrow = F), heights = c(2,1.2)) par(bty = "n", xpd = NA, mar=c(4,4,1,2), oma=c(0,0,0,0)) plot(1,type = "n", xlim = c(0,nrow(mreads)*2), ylim = c(0.09,0.45), xlab = "Time (min)", ylab = "Optical density (600 nm)", las=1, yaxt = "n") mtext("(a)",side = 2, at = 0.47, line = 3, las = 2) axis(2, at=c(0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45), label=c(0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45),las=1) col.pGFP <- colorschemes$Categorical.12[8] col.pDsRed <- colorschemes$Categorical.12[2] color <- c(col.pGFP,col.pDsRed) for(j in 1:2){ for(i in 1:6){ column <- 2*i - 1 + (j-1)*12 polygon(c(mreads$time[1:nrow(mreads)],rev(mreads$time[1:nrow(mreads)])),c(mreads[1:nrow(mreads),column]+mreads[1:nrow(mreads),column+1],rev(mreads[1:nrow(mreads),column])-rev(mreads[1:nrow(mreads),column+1])), border = NA, col = color[j]) } } for(j in 1:2){ for(i in 1:6){ column <- 2*i - 1 + (j-1)*12 lines(mreads$time[1:nrow(mreads)], mreads[1:nrow(mreads),column]) summary$y[i + (j-1)*6] <- mreads[201,column] } } text(405, summary$y[6], ">99%", adj = 0) text(405, summary$y[11], paste(summary$init.freq[11],"%", sep = ""), adj = 0) text(405, summary$y[4], paste(summary$init.freq[4],"%", sep = ""), adj = 0) text(405, summary$y[3]+0.005, paste(summary$init.freq[3],"%", sep = ""), adj = 0) text(405, summary$y[10]-0.005, paste(summary$init.freq[10],"%", sep = ""), adj = 0) text(405, summary$y[2], "<91%", adj = 0) treat <- treat[order(treat$lambda.color,treat$prop.lambda),] plot(1, type = "n", ylim = c(-4,6), xlim = c(0.5,7.5), ylab = "Malthusian fitness", xaxt = "n", xlab = expression(paste("Initial frequency of ",italic(E.)," ",italic(coli)," ",lambda)), las = 1) mtext("(b)",side = 2, at = 6.9, line = 3, las = 2) for(i in 1:7){ ycoord <- treat$fitness_m[(order(init.freq)[i+5]*4-3):(order(init.freq)[i+5]*4)] points(rep(i,4),ycoord,col = ifelse(as.character(treat$lambda.color[order(init.freq)[i+5]*4]) == "green", color[1], color[2]), pch = ifelse(as.character(treat$lambda.color[order(init.freq)[i+5]*4]) == "green", 1, 2)) } axis(1, at = c(1:7), labels = paste(init.freq[rev(order(init.freq,decreasing = T))][6:12],"%", sep = ""), las = 2, tick = F, line = -1) text(c(1:7),5,summary$signif[order(init.freq)[6:12]]) lines(c(0.5,7.5),c(0,0), lty="dotted") dev.off()