#DATASET= RAW READS FROM A TECAN d <- read.delim("RawData.txt", sep="\t") #transpose the excel file (so that data is organized in columns, now rows) #change the column names names(d) <- c(paste("C", seq(length=ncol(d)), sep="")) #add in a column for time t <-seq(0,length(d[,1])-1, by=1)*0.25 #bind the time column to the main data frame d <- cbind(t, d) head(d) ########################################################################################## #PULL OUT THE OPTICAL DENSITY AT 60 HOURS which(t==60) OD <- as.numeric(d[241,2:97]) #Find the first data point where OD is greater than 0.3, using the built-in help menu first3 <- apply(d>0.3,2,which.max) #TANGENT - Find help help(which.max) ?which.max help.search("max") #find examples of functions example(which.max) #make a single data frame with time (ancestral/evolved), OD (at 60 hours), and first (the first time point that OD > 0.5) df <- data.frame(line=c(1:48, 1:48),time=c(rep("anc", 48), rep("evol", 48)), OD, first3=first3[2:97]) #subset the data anc <- subset(df, time=="anc") evol <- subset(df, time=="evol") ############################################################################################### ############################################################################################### #1) Pull out the optical density at 48 hours. Store it to a vector called OD48 in the dataframe df. #2) Make the vector first3 refer to the hour at which od 0.3 was passed (right now it simply points you to the row number ############################################################################################### ############################################################################################### dev.off() #Many different ways to plot the data: #Look at raw curves for one pair of wells plot(d$t, d[,2], ylim=c(0, 1.2)) points(d$t, d[,53]) #Many graphic options (look at par) # main="growth curves" # add a title above the graph # pch=16 # change plot symbol to a filled circle # col="red" # change the fill color # xlim=c(-10,10) # change limits of the x-axis (horizontal axis) # ylim=c(0,100) # change limits of the y-axis (vertical axis) # lty=2 # change line type to dashed # cex=1.5 # magnify the plotting symbols 1.5-fold # xlab="Body size" # label for the x-axis # ylab="Frequency" # label for the y-axis # yaxt="n", xaxt="n" # suppress the y, x axes # las=2 # label axis perpendicular plot(d$t, d[,2], ylim=c(0,1.5), xlab="Time (hours)", ylab="OD", yaxt="n", xaxt="n", col="red") axis(1, at=c(0, 12, 24, 36, 48, 60)) axis(2, las=2) points(d$t, d[,53], col="blue") legend("topleft", legend=c("anc", "evol"), pch=19, col=c("red", "blue"), inset=0.05) #plot all the OD data plot(df$line, df$OD, ylim=c(0, 1.4)) plot(df$line, df$OD, ylim=c(0, 1.4), col=df$time) #boxplot: boxplot(df$OD~df$time) plot(df$time, df$OD) #histogram: hist(df$OD) #plot anc and evol separately par(mfrow=c(2,1)) #two plots one on top of each other hist(anc$OD) hist(evol$OD) #plot anc and evol separately, on the same histogram dev.off() #clear the plot workspace hist(anc$OD, ylim=c(0, 20), xlim=c(0,1.5), col=rgb(1,0,0, alpha=0.65), main="", ylab="Frequency", xlab="OD at 60 hours") par(new=T) #plot a second histogram in the same area hist(evol$OD, ylim=c(0,20), xlim=c(0,1.5), breaks=10, col=rgb(0,0,1, alpha=0.5), main="", ylab="", xlab="") #note the colour designation - this allows you to use transparent colours #stripchart to show different pairs stripchart(df$OD~df$time, vertical=TRUE, pch=21, xlim=c(0.75, 2.25), xaxt="n", yaxt="n") arrows(rep(1,48), anc$OD, rep(2,48), evol$OD, length=0) #add lines between the points axis(1, at=c(1,2), labels=c("anc", "evol")) axis(2, las=2) #scatterplot plot(df$OD, df$first3) plot(df$OD, df$first3, col=df$time) ################################################################## #Summary statistics: #Look at the means and standard deviation of the two fitness assays d.ag <- aggregate(df[c("OD","first3")], df[c("time")], mean) d.ag d.sd <- aggregate(df[c("OD", "first3")], df[c("time")], sd) d.sd #join the se information onto the main dataframe d.ag <- cbind(d.ag, sd.OD=d.sd$OD, sd.first3=d.sd$first3) d.ag plot(d.ag$time, d.ag$OD) #use "seq_along" to for the graph to be points instead of horizontal lines (becuse R wants to plot a box plot) plot(seq_along(d.ag$time), d.ag$OD, ylim=c(0, 1.3), xlim=c(0.75, 2.25)) arrows(seq_along(d.ag$time), d.ag$OD-d.ag$sd.OD, seq_along(d.ag$time), d.ag$OD+d.ag$sd.OD, length=0) #simple t-test treating all 48 strains as equal (ignore different starting colonies) boxplot(df$first3~df$time) od.t<-t.test(anc$first3, evol$first3) od.t #paired t-test stripchart(df$first3~df$time, vertical=TRUE, pch=21, xlim=c(0.75, 2.25), xaxt="n", yaxt="n", ylim=c(0, 60)) arrows(rep(1,48), anc$first3, rep(2,48), evol$first3, length=0) #add lines between the points axis(1, at=c(1,2), labels=c("anc", "evol")) axis(2, las=2) od.tp<-t.test(anc$first3, evol$first3, paired=TRUE) od.tp #correlation plot(df$OD, df$first3) cor.test(df$OD, df$first3) cor.test(df$first3, df$OD) #note, it doesn't matter which comes first palette(c("red", "blue")) plot(df$OD, df$first3, col=df$time) reg1<-lm(df$first3~df$OD) summary(reg1) reg2<-lm(df$first3~df$OD, subset = df$time=="anc") summary(reg2) reg3<-lm(df$first3~df$OD, subset = df$time=="evol") summary(reg3) abline(reg1, col="black") abline(reg2, col="red") abline(reg3, col="blue") dfsub <- subset(df, time=="evol" & first3>2) reg4<-lm(dfsub$first3~dfsub$OD) summary(reg4) abline(reg4, col="blue", lty=2) reg5<-lm(dfsub$OD~dfsub$first3) #anova head(john) str(john) plot(as.factor(john$timeEnv), john$DT) ano <- aov(john$DT ~ john$timeEnv) summary(ano) ano2 <- aov(john$DT ~ john$Time*john$Env*john$PloidyFac) summary(ano2) TukeyHSD(ano) ################################################################################## ################################################################################## #1) Are ploidy and doubling time correlated in this data set? Does time (ancestral versus evolved) affect the relationship between ploidy and doubling time? #Make a plot to visualize the relationship(s). #2) Make a plot to visualize the relationship between doubling time and PloidyFac in evolved strains #Is there a significant relationship between them? Does it depend on the assay environment? ################################################################################## ##################################################################################