## FIGURES library(ggplot2) library(plyr) library(plotrix) library(RColorBrewer) library(gplots) # Import the raw data # These data do not include correlated response assays CellCountsRaw1 <- read.csv("~/Desktop/Figures & Analyses/CellCountsRaw1.csv") View(CellCountsRaw1) ## CELL DIVISIONS # Take the mean of the technical replcates of the flow cytometer (Canto Reps) CellCounts<-ddply(CellCountsRaw1, .(SelecEnv, AssayEnv, Xaxis, Line, ResType, Biorep, AsRep, Block, PopID), summarize, meanT0=mean(CellsPerMlT0, na.rm=TRUE), meanT1=mean(CellsPerMlT1, na.rm=TRUE), MeanMicrons=mean(Microns, na.rm=TRUE), MeanSscT1=mean(SSCT1, na.rm=TRUE), MeanChlT1=mean(RelChl, na.rm=TRUE)) # Calculate cell divisions per day CellDivisions<-ddply(CellCounts, .(SelecEnv, AssayEnv, Xaxis, Line, ResType, Biorep, AsRep, Block, PopID, meanT0, meanT1), summarise, DivPerDay=log((meanT1/meanT0), base=2)/7) # Mean cell divisions per population (biorep) CellDivisionsMean<-ddply(CellDivisions, .(SelecEnv, AssayEnv, Xaxis, Line, ResType, Biorep), summarise, mean=mean(DivPerDay), std.error=std.error(DivPerDay)) # Provide titles for the factes x_names<-c('A'='Control\nControl', 'B'= 'Low light\nLow light','C'= 'Low light\nControl','D'= 'Control\nLow light','E'= 'Low salt\nLow salt','F'= 'Low salt\nControl','G'= 'Control\nLow salt','H'= 'Low PO4\nLow PO4','I'= 'Low PO4\nControl','J'= 'Control\nLow PO4','K'='Random\nControl') # BoxPlot (PlotGrowth<- ggplot()+ geom_boxplot(data = CellDivisionsMean, aes(x=ResType, y = mean, col=ResType, group=ResType, panel.background="white"), size=0.5, width=0.9, position=pd) + geom_point(data=CellDivisionsMean, aes(x=ResType, y=mean, group=ResType, shape=Line, col=ResType), size = 2, position=pd) + scale_shape_manual(values=c(1,2,3,4,5,6,7,8,9), breaks=c("NG5", "NG'13", "NG26", "NG'10", "NG'16", "NG27", "NG'2", "NG'3", "NG'4")) + scale_colour_discrete(name="ResType", breaks=c("NG5", "NG'13", "NG26", "NG'10", "NG'16", "NG27", "NG'2", "NG'3", "NG'4")) + theme_bw() + theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank()) + xlab("Resistance Type") + ylab("Cell divisions (per day)") +facet_grid(.~Xaxis, labeller=as_labeller(x_names)) +theme(strip.background=element_rect(fill="white",colour="black",size=0.35)) + ylim(0.45,1.5) + theme(axis.title=element_text(size=16), axis.text=element_text(size=14), strip.text.x=element_text(size=13))+theme(axis.title.y=element_text(vjust=1.5))+ geom_hline(aes(yintercept=1), data=CellDivisionsMean, linetype="dashed", size=0.6, col="black")+theme(axis.title.x=element_text(vjust=-0.5))) ## CELL SIZE CellSize<-ddply(CellCounts, .(SelecEnv, AssayEnv, Xaxis, Line, PopID, ResType, Biorep), summarise, Size=mean(MeanMicrons), std.error=std.error(MeanMicrons)) (PlotSize<- ggplot()+ geom_boxplot(data = CellSize, aes(x=ResType, y = Size, col=ResType, group=ResType, panel.background="white"), size=0.5, width=0.8, position=pd) + geom_point(data=CellSize, aes(x=ResType, y=Size, group=ResType, shape=Line, col=ResType), size = 2, position=pd) + scale_shape_manual(values=c(1,2,3,4,5,6,7,8,9)) + theme_bw() + theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank()) + xlab("Resistance Type") + ylab("Cell diameter (μm)") +facet_grid(.~Xaxis, labeller=as_labeller(x_names)) + theme(axis.title=element_text(size=15), axis.text=element_text(size=12), strip.text.x=element_text(size=13))+theme(axis.title.y=element_text(vjust=1.5))+theme(axis.title.x=element_text(vjust=-0.5)) + theme(strip.background=element_rect(fill="white",colour="black",size=0.35)) + ylim(0.83,1.01) +theme(axis.title.x=element_text(margin=margin(12)))) ## CHLOROPHYLL Chl<-ddply(CellCounts, .(SelecEnv, AssayEnv, Xaxis, Line, PopID, ResType, Biorep), summarise, ChlMean=mean(MeanChlT1), std.error=std.error(MeanChlT1)) (PlotChl<- ggplot()+ geom_boxplot(data = Chl, aes(x=ResType, y = ChlMean, col=ResType, group=ResType, panel.background="white"), size=0.5, width=0.8, position=pd) + geom_point(data=Chl, aes(x=ResType, y=ChlMean, group=ResType, shape=Line, col=ResType), size = 2, position=pd) + scale_shape_manual(values=c(1,2,3,4,5,6,7,8,9)) + theme_bw() + theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank()) + xlab("Resistance Type") + ylab("Relative Chlorophyll") +facet_grid(.~Xaxis, labeller=as_labeller(x_names)) + theme(axis.title=element_text(size=15), axis.text=element_text(size=12), strip.text.x=element_text(size=13))+theme(axis.title.y=element_text(vjust=1.5))+theme(axis.title.x=element_text(vjust=-0.5)) + theme(strip.background=element_rect(fill="white",colour="black",size=0.35)) + geom_hline(aes(yintercept=1), data=CellDivisionsMean, linetype="dashed", size=0.6, col="grey") + ylim(0.7,1.5) +theme(axis.title.x=element_text(margin=margin(12)))) ## STATISTICS library(lme4) library(plyr) library("lmerTest", lib.loc="~/Library/R/3.2/library") # Import the raw data # These data include correlated response assays CellCountsRaw <- read.csv("~/Desktop/Figures & Analyses/CellCountsRaw.csv") View(CellCountsRaw) # Take the means of the technical replicates (cantorep) and then calculate cell divisions CellCounts<-ddply(CellCountsRaw, .(SelecEnv, AssayEnv, Assay, Line, Xaxis, ResType, Biorep, AsRep, PopID, Block), summarize, meanT0=mean(CellsPerMlT0, na.rm=TRUE), meanT1=mean(CellsPerMlT1, na.rm=TRUE), MeanFscT1=mean(FSCT1, na.rm=TRUE), MeanSscT1=mean(SSCT1, na.rm=TRUE), MeanChlT1=mean(RelChl, na.rm=TRUE)) CellDivisions<-ddply(CellCounts, .(SelecEnv, AssayEnv, Assay, Line, Xaxis, ResType, Biorep, AsRep, PopID, Block, meanT0, meanT1), summarise, DivPerDay=log((meanT1/meanT0), base=2)/7) View(CellDivisions) ## Have a look at the data by eye str(CellDivisions) summary(CellDivisions) plot(CellDivisions$ResType, CellDivisions$DivPerDay) plot(CellDivisions$SelecEnv, CellDivisions$DivPerDay) plot(CellDivisions$AssayEnv, CellDivisions$DivPerDay) ## Check whether the data are normally distributed hist(CellDivisions$DivPerDay, 300) hist(CellDivisions$DivPerDay) qqnorm(CellDivisions$DivPerDay):qqline(CellDivisions$DivPerDay) shapiro.test(CellDivisions$DivPerDay) # Data are not normally distributed, but this is ok because, as we see after the model is run, the residuals are normally dirstributed # Run the model growth<-lmer(DivPerDay~SelecEnv * AssayEnv * ResType + (1|PopID) + (1|Block), data=CellDivisions) # When this model is run, there is a warning message # fixed-effect model matrix is rank deficient so dropping 27 columns / coefficients # Such shrinkage is a good sign because it ensures numerical stability # Check the residuals plot(growth) hist(resid(growth), 100) # Residuals are normally distributed anova(growth) summary(growth) # Fitting a model to compare growth when populations were assayed in the selection environment vs assayed in a different environment from the selection environment Assay<-lmer(DivPerDay~Assay * ResType + (1|PopID) + (1|Block), data=CellDivisions) anova(Assay) # Direct response to seletion direct<-lmer(DirectResponse~SelecEnv * ResType + (1|PopID) + (1|Block), data=DirectResponse) anova(direct) summary(direct) ## CELL SIZE AND CHLOROPHYLL CONTENT CellCounts<-ddply(CellCountsRaw, .(SelecEnv, AssayEnv, Xaxis, Line, ResType, Biorep, AsRep, PopID, Block), summarize, meanT0=mean(CellsPerMlT0, na.rm=TRUE), meanT1=mean(CellsPerMlT1, na.rm=TRUE), MeanSize=mean(Microns, na.rm=TRUE), MeanChl=mean(RelChl, na.rm=TRUE)) View(CellCounts) # Cell size # Have a look at the data by eye and check normality plot(CellCounts$SelecEnv, CellCounts$MeanSize) plot(CellCounts$AssayEnv, CellCounts$MeanSize) plot(CellCounts$ResType, CellCounts$MeanSize) hist(CellSize$Size, 100) qqnorm(CellSize$Size):qqline(CellSize$Size) # Model for cell size size<-lmer(MeanSize~SelecEnv * AssayEnv * ResType + (1|PopID) + (1|Block), data=CellCounts) anova(size) plot(size) hist(resid(size), 100) summary(size) # Chlorophyll # Have a look at the data by eye and check normality plot(CellCounts$ResType, CellCounts$MeanChl) plot(CellCounts$AssayEnv, CellCounts$MeanChl) plot(CellCounts$SelecEnv, CellCounts$MeanChl) hist(Chlorophyll$Chl, 100) qqnorm(Chlorophyll$Chl) # Model for chlorophyll chl<-lmer(MeanChl~SelecEnv * AssayEnv * ResType + (1|PopID) + (1|Block), data=CellCounts) anova(chl) summary(chl) plot(chl) hist(resid(chl), 100)