Temperature time course

Tracking changes over time is a very useful way of understanding and analysing your system. This applies whether your system is biological or otherwise. As part of my hobby of baking bread, I have built a brick oven in my garden. I have installed thermocouples that allow me to measure the temperature at the top and the bottom of my oven. To investigate the performance of my oven, I have recorded temperatures when I fired the oven in May, June and August of this year.

The data for this is available on github and is downloaded as part of the script.

I have written a script to download, manipulate and graph the data. 

The manipulations required were:

• turning my date and time into a date an time object that R could understand
• subtracting the start time to calculate the elapsed time
• the elapsed time was seconds which was converted to hours
• the time was then converted to a number to allow for graphing by ggplot
• the data was melted into a format for ggplot

Within the script needs to be applied separately to each of the three months of data. This can be done easily by changing the read.table() function to the appropriate file. 

By doing this, I graph and analyse the data in the same way making the graphs easier to compare. Here are the three graphs:

Here is the script for generating these graphs:

library(reshape2)

library(ggplot2)

library(RCurl) # allows us to download data through urls 

# URL for May data:

May <- getURL("https://raw.githubusercontent.com/brennanpincardiff/RforBiochemists/master/data/ovenTempMay.tsv")

# URL for Jun data:

Jun <- getURL("https://raw.githubusercontent.com/brennanpincardiff/RforBiochemists/master/data/ovenTempJun.tsv")

# URL for Aug data: 

Aug <- getURL("https://raw.githubusercontent.com/brennanpincardiff/RforBiochemists/master/data/ovenTempAug.tsv")

# put the URLs together in a vector

urls <- c(May, Jun, Aug)

months <- c("May", "Jun", "Aug")

# do it in a loop to apply to each of the three data sets

for(i in 1:3){

# read in the data for the relevant month - it's a tab separated file 

    data <- read.table(text = urls[i], stringsAsFactors=FALSE, sep = "\t", header = TRUE)

    

    

# convert the data and time into a format that R will understand. make it into a POSIXct object. 

# so as not to mess with the original data put it in a new column called time.P 

# because my dates are separted by slashes "/", I need to tell R about the format. 

    data$time.P <- as.POSIXct(data$time, format = "%d/%m/%Y %H:%M")

    

# the first value is the start.time for this temperature profile. 

    start.time <- data$time.P[1]

    

# calculate the elapsed time 

# substract the start time from each value. 

# This returns the time in seconds

    data$e.time.sec <- data$time.P - start.time

# convert this into an hour time and change back into a number for graphing purposes. 

    data$e.time.hour.num <- as.numeric(data$e.time.sec/3600)

    

    

 # take out the data we need for the plot

    data.subset <- as.data.frame(data$e.time.hour.num)

    colnames(data.subset)[1] <- "e.time.hour.num"

    data.subset$top <- data$top

    data.subset$bot <- data$bot

    

# melt it into a format for ggplot using melt() function

   data.subset.melt <- melt(data.subset, 
                           id.vars = "e.time.hour.num")

# put in nice column names
  colnames(data.subset.melt) <- c("elapsed.time", "place", "temp")
  
# make the graph object
  p <- ggplot(data.subset.melt, aes(x=elapsed.time, 
                                    y= temp, 
                                    colour = factor(place, labels = c("Top", "Bottom")))) + 

# colour = factor and the labels allows us to customize the legend

        geom_line(size=1) +

        geom_point() +

        labs(color = "Place") + # customizes the legend title

        scale_colour_manual(values=c("black","red")) +

        ylab("Temperature") + # y-label

        xlab("Elapsed time (hours)") + # x-label

        ylim(0,400) +

        scale_x_continuous(limits=c(0, 48), 

                           breaks=c(0,1,4,8,12, 24, 48)) +

        theme_bw()

# position the legend  

    p <- p + theme(legend.position=c(1,1), # move to the top right

                   legend.justification=c(1,1), # move it in a bit

                   legend.text=element_text(size = 12), # increase size of text

                   legend.title=element_text(size = 12)) # and the title

    

    p <- p + theme(axis.title.y = element_text(size = 14 )) + 

      theme(axis.text = element_text(size = 12))

    

# add the appropriate title

    p <- p + ggtitle(paste0("Oven Temp (", months[i], " Firing)")) 

# print the object - you have to do this because of the loop

    print(p)

}

Helpful resources, I used to prepare this script:

• http://www.noamross.net/blog/2014/2/10/using-times-and-dates-in-r---presentation-code.html
• http://www.markhneedham.com/blog/2014/07/20/r-ggplot-dont-know-how-to-automatically-pick-scale-for-object-of-type-difftime-discrete-value-supplied-to-continuous-scale/
• http://stackoverflow.com/questions/7531868/how-to-rename-a-single-column-in-a-data-frame-in-r
• http://docs.ggplot2.org/current/scale_continuous.html
• http://stackoverflow.com/questions/14441729/read-a-csv-from-github-into-r
• http://stackoverflow.com/questions/23897027/pass-multiple-variables-into-ggtitle-r
• https://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-do-lattice_002ftrellis-graphics-not-work_003f

Using a for loop to compare data from multiple experiments

One of the key benefits of computer programming is allowing the programme to repeat steps (so you don't have to).
Dr Dean Hammond has written this script that illustrates how to use a for loop to plot data from six separate enzyme kinetic experiments to allow a comparison. It builds on the previous example of plotting enzymatic data.

The data is all contained in a matrix, a type of two dimensional data structure where all the data is of the same type.

Here is the output that is produced:

Here is the script that makes it:

# Following on from Prof. Beynon's example with enzymatic data... 

# For multiplotting 6 enzymology data-sets, using base R 

# Data from six experiments

no.Exp <- c("Exp 1","Exp 2", "Exp 3", "Exp 4","Exp 5", "Exp 6")

# Substrate concentrations:

Sub <- c(0, 1, 2, 4, 8, 12, 16, 20, 30, 40)

# Data in a matrix - 2D object with data of the same class (numeric)

enzdata <- matrix(c(0, 17.36667, 31.97143, 52.68889, 61.95385, 74.2, 77.97143, 84.28, 99.91429, 93.66667, 

                    0, 15.7, 29.42286, 45.64, 62.60615, 75.78118, 69.88, 75.256, 89.59429, 86.84, 

                    0, 27.10667, 42.12, 63.48, 69.56, 74.26857, 79.44444, 83.29091, 87.1, 82.08571, 

                    0, 24.72, 39.07, 47.4, 57.928, 67.6, 71.35556, 67, 75.79375, 70.86667, 

                    0, 5.723636, 11.48, 17.697143, 28.813333, 37.567273, 42.483077, 40.68, 52.81, 56.92, 

                    0, 2.190476, 5.254545, 8.95, 15.628571, 20.8, 25.355556, 26.55, 32.44, 33.333333),

                  nrow=10,

                  ncol=6)

# specify plotting parameters for our multiplot page:

par(mfrow = c(3, 2),  # 6 plots in a 2 column x 3 row format

    oma = c(1,1,1,1), # oma = outer margin in lines, of each plots (bottom, left, top, right)

    mar = c(3,3,2,1),  # mar = no. of lines to be specified on the four sides of each plot

    cex.main = 0.9,    # main text size

    las = 1)           # all axis labels horizontal

# here's a for loop 

# to plot the data for each enzymatic reaction (Expt), 

# get the values of Km and Vmax from the theoretical formula. 

# Then, build a theoretical line defining the best fit curve.

for(i in 1:length(no.Exp)){    # for every experiment - one col

    v <- enzdata[, i]          # get the velocity ('v') data

    data <- cbind(Sub, v)      # create a data-set for each expt  

    fit <- nls(v ~ Vmax * (Sub / (Km + Sub)),

               start = list(Vmax = 50, Km = 2))

    

    # write a title for each peptide plot, based on colnames: 

    title = paste(no.Exp[i])    # use exp name as a plot title

    

    SconcRange <- seq(0, 50, 0.1)

    theorLine <- predict(fit, list(Sub = SconcRange))

    

    # draw each plot, with points, 

    #applying the correct title adjusting font sizes accordingly:

    plot(Sub, v, main = title,

         col = 'red', pch = 16,

         cex.lab = 0.6, cex.axis = 0.8,

         xlab = NA, ylab = NA)      # omit drawing axis labels

    

    # add the theoretical lines: 

    lines(SconcRange, theorLine, lwd = 1.5, col = 'blue')

  }  # this curly bracket is the end of the for loop. 

# use mtext to add single x- and y-axis labels for all plots, close to the margin:

mtext("Velocity (nmol/s)", side = 2, las = 0, outer = TRUE, line = -1)

mtext("Substrate (mM)", side = 1, outer = TRUE, line = -0.5)