se <- function(x) sqrt(var(x)/(length(x) - 1)) library(plotrix) sediff <-function(y1,y2) sqrt(((((length(y1)-1)*var(y1)+(length(y2)-1)*var(y2))/(length(y1)+length(y2)-2)))/length(y1)+(((length(y1)-1)*var(y1)+(length(y2)-1)*var(y2))/(length(y1)+length(y2)-2)/length(y2))) datacount <- read.csv("datacount.csv", sep=",",header=T) #Ethanol Ne #ethanol ceth <- subset(datacount, Environment=="Ethanol") ceth$Num.Cell <- ceth$Num.Cell*10 eth.test <- aggregate(ceth$Num.Cell,list(Ploidy=ceth$Ploidy,Time=ceth$Time,Line=ceth$Line), mean) eth.test.1 <- subset(eth.test, Ploidy==1) eth.anc.1 <- subset(eth.test.1, Time==49) eth.anc.1.se <- se((eth.anc.1$x*0.0315)) eth.anc.1.mean <- mean((eth.anc.1$x*0.0315)) eth.evol.1 <- subset(eth.test.1, Time==189) eth.evol.1.se <- se(eth.evol.1$x*0.0315) eth.evol.1.mean <- mean(eth.evol.1$x*0.0315) eth.test.2 <- subset(eth.test, Ploidy==2) eth.anc.2 <- subset(eth.test.2, Time==49) eth.anc.2.se <- se(eth.anc.2$x*0.0315) eth.anc.2.mean <- mean(eth.anc.2$x*0.0315) eth.evol.2 <- subset(eth.test.2, Time==189) eth.evol.2.se <- se(eth.evol.2$x*0.0315) eth.evol.2.mean <- mean(eth.evol.2$x*0.0315) eth.graph <- c(eth.anc.1.mean, eth.anc.2.mean, eth.evol.1.mean, eth.evol.2.mean) eth.graph.se <-c(eth.anc.1.se, eth.anc.2.se, eth.evol.1.se, eth.evol.2.se) eth.testing <- rbind(eth.anc.1, eth.anc.2, eth.evol.1, eth.evol.2) ethtest <- aov(eth.testing$x~eth.testing$Ploidy*eth.testing$Time) #kohanol Ne ckoh <- subset(datacount, Environment=="KOH") ckoh$Num.Cell <- ckoh$Num.Cell*10 koh.test <- aggregate(ckoh$Num.Cell,list(Ploidy=ckoh$Ploidy,Time=ckoh$Time,Line=ckoh$Line), mean) koh.test.1 <- subset(koh.test, Ploidy==1) koh.anc.1 <- subset(koh.test.1, Time==49) koh.anc.1.se <- se((koh.anc.1$x*0.0315)) koh.anc.1.mean <- mean((koh.anc.1$x*0.0315)) koh.evol.1 <- subset(koh.test.1, Time==189) koh.evol.1.se <- se(koh.evol.1$x*0.0315) koh.evol.1.mean <- mean(koh.evol.1$x*0.0315) koh.test.2 <- subset(koh.test, Ploidy==2) koh.anc.2 <- subset(koh.test.2, Time==49) koh.anc.2.se <- se(koh.anc.2$x*0.0315) koh.anc.2.mean <- mean(koh.anc.2$x*0.0315) koh.evol.2 <- subset(koh.test.2, Time==189) koh.evol.2.se <- se(koh.evol.2$x*0.0315) koh.evol.2.mean <- mean(koh.evol.2$x*0.0315) koh.graph <- c(koh.anc.1.mean, koh.anc.2.mean, koh.evol.1.mean, koh.evol.2.mean) koh.graph.se <-c(koh.anc.1.se, koh.anc.2.se, koh.evol.1.se, koh.evol.2.se) koh.testing <- rbind(koh.anc.1, koh.anc.2, koh.evol.1, koh.evol.2) kohtest <- aov(koh.testing$x~koh.testing$Ploidy*koh.testing$Time) #NaCl cnac <- subset(datacount, Environment=="NaCl") cnac$Num.Cell <- cnac$Num.Cell*10 nac.test <- aggregate(cnac$Num.Cell,list(Ploidy=cnac$Ploidy,Time=cnac$Time,Line=cnac$Line), mean) nac.test.1 <- subset(nac.test, Ploidy==1) nac.anc.1 <- subset(nac.test.1, Time==49) nac.anc.1.se <- se((nac.anc.1$x*0.0315)) nac.anc.1.mean <- mean((nac.anc.1$x*0.0315)) nac.evol.1 <- subset(nac.test.1, Time==189) nac.evol.1.se <- se(nac.evol.1$x*0.0315) nac.evol.1.mean <- mean(nac.evol.1$x*0.0315) nac.test.2 <- subset(nac.test, Ploidy==2) nac.anc.2 <- subset(nac.test.2, Time==49) nac.anc.2.se <- se(nac.anc.2$x*0.0315) nac.anc.2.mean <- mean(nac.anc.2$x*0.0315) nac.evol.2 <- subset(nac.test.2, Time==189) nac.evol.2.se <- se(nac.evol.2$x*0.0315) nac.evol.2.mean <- mean(nac.evol.2$x*0.0315) nac.graph <- c(nac.anc.1.mean, nac.anc.2.mean, nac.evol.1.mean, nac.evol.2.mean) nac.graph.se <-c(nac.anc.1.se, nac.anc.2.se, nac.evol.1.se, nac.evol.2.se) nac.testing <- rbind(nac.anc.1, nac.anc.2, nac.evol.1, nac.evol.2) nactest <- aov(nac.testing$x~nac.testing$Ploidy*nac.testing$Time) #Caffeine ccaf <- subset(datacount, Environment=="Caffeine") ccaf$Num.Cell <- ccaf$Num.Cell*10 caf.test <- aggregate(ccaf$Num.Cell,list(Ploidy=ccaf$Ploidy,Time=ccaf$Time,Line=ccaf$Line), mean) caf.test.1 <- subset(caf.test, Ploidy==1) caf.anc.1 <- subset(caf.test.1, Time==49) caf.anc.1.se <- se((caf.anc.1$x*0.0315)) caf.anc.1.mean <- mean((caf.anc.1$x*0.0315)) caf.evol.1 <- subset(caf.test.1, Time==189) caf.evol.1.se <- se(caf.evol.1$x*0.0315) caf.evol.1.mean <- mean(caf.evol.1$x*0.0315) caf.test.2 <- subset(caf.test, Ploidy==2) caf.anc.2 <- subset(caf.test.2, Time==49) caf.anc.2.se <- se(caf.anc.2$x*0.0315) caf.anc.2.mean <- mean(caf.anc.2$x*0.0315) caf.evol.2 <- subset(caf.test.2, Time==189) caf.evol.2.se <- se(caf.evol.2$x*0.0315) caf.evol.2.mean <- mean(caf.evol.2$x*0.0315) caf.graph <- c(caf.anc.1.mean, caf.anc.2.mean, caf.evol.1.mean, caf.evol.2.mean) caf.graph.se <-c(caf.anc.1.se, caf.anc.2.se, caf.evol.1.se, caf.evol.2.se) caf.testing <- rbind(caf.anc.1, caf.anc.2, caf.evol.1, caf.evol.2) caftest <- aov(caf.testing$x~caf.testing$Ploidy*caf.testing$Time) #nys cnys <- subset(datacount, Environment=="Nystatin") cnys$Num.Cell <- cnys$Num.Cell*10 nys.test <- aggregate(cnys$Num.Cell,list(Ploidy=cnys$Ploidy,Time=cnys$Time,Line=cnys$Line), mean) nys.test.1 <- subset(nys.test, Ploidy==1) nys.anc.1 <- subset(nys.test.1, Time==49) nys.anc.1.se <- se((nys.anc.1$x*0.0315)) nys.anc.1.mean <- mean((nys.anc.1$x*0.0315)) nys.evol.1 <- subset(nys.test.1, Time==189) nys.evol.1.se <- se(nys.evol.1$x*0.0315) nys.evol.1.mean <- mean(nys.evol.1$x*0.0315) nys.test.2 <- subset(nys.test, Ploidy==2) nys.anc.2 <- subset(nys.test.2, Time==49) nys.anc.2.se <- se(nys.anc.2$x*0.0315) nys.anc.2.mean <- mean(nys.anc.2$x*0.0315) nys.evol.2 <- subset(nys.test.2, Time==189) nys.evol.2.se <- se(nys.evol.2$x*0.0315) nys.evol.2.mean <- mean(nys.evol.2$x*0.0315) nys.graph <- c(nys.anc.1.mean, nys.anc.2.mean, nys.evol.1.mean, nys.evol.2.mean) nys.graph.se <-c(nys.anc.1.se, nys.anc.2.se, nys.evol.1.se, nys.evol.2.se) nys.testing <- rbind(nys.anc.1, nys.anc.2, nys.evol.1, nys.evol.2) nystest <- aov(nys.testing$x~nys.testing$Ploidy*nys.testing$Time) #HCl chcl <- subset(datacount, Environment=="HCl") chcl$Num.Cell <- chcl$Num.Cell*10 hcl.test <- aggregate(chcl$Num.Cell,list(Ploidy=chcl$Ploidy,Time=chcl$Time,Line=chcl$Line), mean) hcl.test.1 <- subset(hcl.test, Ploidy==1) hcl.anc.1 <- subset(hcl.test.1, Time==49) hcl.anc.1.se <- se((hcl.anc.1$x*0.0315)) hcl.anc.1.mean <- mean((hcl.anc.1$x*0.0315)) hcl.evol.1 <- subset(hcl.test.1, Time==189) hcl.evol.1.se <- se(hcl.evol.1$x*0.0315) hcl.evol.1.mean <- mean(hcl.evol.1$x*0.0315) hcl.test.2 <- subset(hcl.test, Ploidy==2) hcl.anc.2 <- subset(hcl.test.2, Time==49) hcl.anc.2.se <- se(hcl.anc.2$x*0.0315) hcl.anc.2.mean <- mean(hcl.anc.2$x*0.0315) hcl.evol.2 <- subset(hcl.test.2, Time==189) hcl.evol.2.se <- se(hcl.evol.2$x*0.0315) hcl.evol.2.mean <- mean(hcl.evol.2$x*0.0315) hcl.graph <- c(hcl.anc.1.mean, hcl.anc.2.mean, hcl.evol.1.mean, hcl.evol.2.mean) hcl.graph.se <-c(hcl.anc.1.se, hcl.anc.2.se, hcl.evol.1.se, hcl.evol.2.se) hcl.testing <- rbind(hcl.anc.1, hcl.anc.2, hcl.evol.1, hcl.evol.2) hcltest <- aov(hcl.testing$x~hcl.testing$Ploidy*hcl.testing$Time) #YPD cypd <- subset(datacount, Environment=="YPD") cypd$Num.Cell <- cypd$Num.Cell*10 ypd.test <- aggregate(cypd$Num.Cell,list(Ploidy=cypd$Ploidy,Time=cypd$Time,Line=cypd$Line), mean) ypd.test.1 <- subset(ypd.test, Ploidy==1) ypd.anc.1 <- subset(ypd.test.1, Time==49) ypd.anc.1.se <- se((ypd.anc.1$x*0.0315)) ypd.anc.1.mean <- mean((ypd.anc.1$x*0.0315)) ypd.evol.1 <- subset(ypd.test.1, Time==189) ypd.evol.1.se <- se(ypd.evol.1$x*0.0315) ypd.evol.1.mean <- mean(ypd.evol.1$x*0.0315) ypd.test.2 <- subset(ypd.test, Ploidy==2) ypd.anc.2 <- subset(ypd.test.2, Time==49) ypd.anc.2.se <- se(ypd.anc.2$x*0.0315) ypd.anc.2.mean <- mean(ypd.anc.2$x*0.0315) ypd.evol.2 <- subset(ypd.test.2, Time==189) ypd.evol.2.se <- se(ypd.evol.2$x*0.0315) ypd.evol.2.mean <- mean(ypd.evol.2$x*0.0315) ypd.graph <- c(ypd.anc.1.mean, ypd.anc.2.mean, ypd.evol.1.mean, ypd.evol.2.mean) ypd.graph.se <-c(ypd.anc.1.se, ypd.anc.2.se, ypd.evol.1.se, ypd.evol.2.se) ypd.testing <- rbind(ypd.anc.1, ypd.anc.2, ypd.evol.1, ypd.evol.2) ypdtest <- aov(ypd.testing$x~ypd.testing$Ploidy*ypd.testing$Time) graphmeans.N.a <- c(ypd.graph[1], hcl.graph[1],eth.graph[1],koh.graph[1],nys.graph[1],nac.graph[1],caf.graph[1]) graphmeans.2N.a <- c(ypd.graph[2], hcl.graph[2],eth.graph[2],koh.graph[2],nys.graph[2],nac.graph[2],caf.graph[2]) graphmeans.N.e <- c(ypd.graph[3], hcl.graph[3],eth.graph[3],koh.graph[3],nys.graph[3],nac.graph[3],caf.graph[3]) graphmeans.2N.e <- c(ypd.graph[4], hcl.graph[4],eth.graph[4],koh.graph[4],nys.graph[4],nac.graph[4],caf.graph[4]) graphmeans.N.a.se <- c(ypd.graph.se[1], hcl.graph.se[1],eth.graph.se[1],koh.graph.se[1],nys.graph.se[1],nac.graph.se[1],caf.graph.se[1]) graphmeans.2N.a.se <- c(ypd.graph.se[2], hcl.graph.se[2],eth.graph.se[2],koh.graph.se[2],nys.graph.se[2],nac.graph.se[2],caf.graph.se[2]) graphmeans.N.e.se <- c(ypd.graph.se[3], hcl.graph.se[3],eth.graph.se[3],koh.graph.se[3],nys.graph.se[3],nac.graph.se[3],caf.graph.se[3]) graphmeans.2N.e.se <- c(ypd.graph.se[4], hcl.graph.se[4],eth.graph.se[4],koh.graph.se[4],nys.graph.se[4],nac.graph.se[4],caf.graph.se[4]) #black and white pdf("figure2", width=6.5, height=4, family="Times", pointsize=9) par(mar=c(5,7,1,1)) plot(1:7, graphmeans.N.a, col="black", pch=21, ylim=c(0,4e+6), xlim=c(0.9,7.3), xaxt="n", ylab="", xlab="Environment (stressor + YPD)", cex=2, yaxt="n", cex.lab=1.3, cex.axis=1.2) names<-c("YPD","HCl","Ethanol","KOH","Nystatin","NaCl","Caffeine") axis(1, at=1:7+0.125, labels=names, cex.axis=1.2, cex.lab=1.2) axis(2, at=c(0, 1e+06, 2e+06, 3e+06, 4e+06), labels=c("0",expression(1 %*% 10^6),expression(2 %*% 10^6),expression(3 %*% 10^6),expression(4 %*% 10^6)), cex.axis=1.2, las=2, cex.lab=1.3) title(ylab="Adaptive effective population size", line=5, cex.lab=1.3) points(1:7+0.05, graphmeans.N.e, col=grey(0.2), pch=24, cex=1.6) points(1:7+0.2, graphmeans.2N.a, col="black", cex=1.6, pch=19) points(1:7+0.25, graphmeans.2N.e, col=grey(0.2), cex=1.6, pch=24,bg=grey(0.2)) arrows(1:7, graphmeans.N.a-graphmeans.N.a.se, 1:7, graphmeans.N.a+graphmeans.N.a.se, length=0, col="black") arrows(1:7+0.05, graphmeans.N.e-graphmeans.N.e.se, 1:7+0.05, graphmeans.N.e+graphmeans.N.e.se, length=0, col=grey(0.2)) arrows(1:7+0.2, graphmeans.2N.a-graphmeans.2N.a.se, 1:7+0.2, graphmeans.2N.a+graphmeans.2N.a.se, length=0, col="black") arrows(1:7+0.25, graphmeans.2N.e-graphmeans.2N.e.se, 1:7+0.25, graphmeans.2N.e+graphmeans.2N.e.se, length=0, col=grey(0.2)) #axis.break(axis=2,breakpos=5.5e+06,bgcol="white",breakcol="black", # style="zigzag",brw=0.035) txt_l1 <- expression(paste(italic(N[h]), " ancestral")) txt_l2 <- expression(paste(italic(N[h]), " evolved")) txt_l3 <- expression(paste(italic(N[d]), " ancestral")) txt_l4 <- expression(paste(italic(N[d]), " evolved")) legend(0.75, 1.15e+06, legend=c(txt_l1, txt_l2, txt_l3, txt_l4), col=c("black","black","black",grey(0.2)),pt.bg=c("white","white","black",grey(0.2)), pch=c(21,24,19,24), cex=1.1, y.intersp=1.1) dev.off() system("open figure2")