Thursday, 28 May 2015

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