ggplot script for VizBi 2016

I have updated the ggplot script for VizBi.

It's based on Dean Hammond's previous script and includes a loop - often seen as bad practice in R but it seems to work here.

It makes this:

# START
# The plot with the six lines is a bit busy. 
# Let's make a ggplot equivalent to the base R multiplot.
# Here we do it using faceting.
# To make the script free standing the reshaping code is included. 

setwd("/Users/paulbrennan/Documents/RforBiochemistsScripts/VizBiMar2016")

library(reshape2)
library(ggplot2)
library(ggthemes)

# some data:
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)
no.Exp <- c("Exp.1", "Exp.2", "Exp.3", "Exp.4", "Exp.5", "Exp.6")
Sub <- c(0, 1, 2, 4, 8, 12, 16, 20, 30, 40)

enzdata <- as.data.frame(enzdata) # converted to a data.frame for plotting
colnames(enzdata) <- rep("Exp", times=6)  # add a column name
enzdata <- cbind(Sub, enzdata) # add the Substrate data

# so we have a data.frame and we want to draw some graphs. 
# as is always the case with R, there are many ways to do this. 

# using our dataframe we can draw the graphs one by one...

p <-  ggplot(data = enzdata, aes(x = Sub, y = Exp)) +
  geom_point() +
  theme_few() +
  ylab("velocity (nmol/s)") + xlab("substrate (mM)") +
  ggtitle("Experiment 1")
p

# one of the useful feature of ggplot is that you can force new data into an old plot
# we have the data from Exp 6 in an object. 
# if we make the same shaped dataframe with another set of data we can force that in.
new.data <- data.frame(enzdata$Sub, enzdata[,3])
colnames(new.data) <- c("Sub", "Exp")  # col names must be the same
p <- p %+% new.data   # this %+% is what forces in the new data...
p <- p + ggtitle("Experiment 2")
p


# we can use this feature and put making and saving the graphs inside a loop
for(i in 1:6) {  
  new.data <- data.frame(enzdata$Sub, enzdata[,i+1])
  colnames(new.data) <- c("Sub", "Exp")  # col names must be the same
  p <- p %+% new.data
  p <- p + ggtitle(paste("Experiment", i))
  filename <- paste0("Experiment_", i, ".pdf")
  p + ggsave(filename)
}
# it would probably be better style to write a function and vectorise it!


## there is also an 'even better' way to do it using facets - a useful feature of ggplot
# requires us to change the structure of the dataframe
# melt the data from wide to long:
colnames(enzdata) <- c("Sub", "Exp.1", "Exp.2", "Exp.3", "Exp.4", "Exp.5", "Exp.6")

melted_data <- melt(enzdata, id.vars = "Sub", value.name = "v", variable.name = "Exp")

# Check out the differences between these two files:
view(enzdata)
view(melted_data)

# the way ggplot draws colours depends on the required number of colours. The way it does it is just with equally spaced hues around the color wheel, starting from 15:
gg_color_hue <- function(n){
  hues = seq(15, 375, length = n + 1)
  hcl(h = hues, l = 65, c = 100)[1:n]
  }
  
# we have 6 experiments (so we want to create our own 6 colour "mini-palette")
cols <- gg_color_hue(6)  
# this gets the default ggplot colours when we want to plot 6 different variables.


# 1st, create an object with the ggplot 
# It has all the data simply as points, shaped and colour-coded by expt. of origin (as above). 
# We store it in an object called 'x':
x <- ggplot(data = melted_data, aes(x = Sub, y = v)) +
  geom_point(aes(colour = Exp, shape = Exp)) +
  theme_few() +
  ylab("velocity (nmol/s)") + xlab("substrate (mM)") +
  theme(axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        plot.title = element_blank(),            # we don't want individual plot titles as the facet "strip" will give us this
        legend.position = "none",                # we don't want a legend either
        panel.border = element_rect(fill = NA, color = "darkgrey", size = 1.25, linetype = "solid"),
        axis.ticks = element_line(colour = 'darkgrey', size = 1.25, linetype = 'solid'))     # here, I just alter to colour and thickness of the plot outline and tick marks. You generally have to do this when faceting, as well as alter the text sizes (= element_text() in theme also)

# have a quick look at x
x
# see all the data in different colours

# Next, let's modify 'x' by faceting based on the Experiment Name, 
# We also specify that we want a panel 3 plots wide x 2 plots high.
# You can change this, obviously. Can also 'facet_grid' too:
x <- x + facet_wrap( ~ Exp, ncol = 3)

# Finally, let's apply the Michaelis Menten fitting to our faceted data, adding the best fit line in each case:
x <- x + geom_smooth(method = "nls",                 
                     method.args = list(formula = y ~ Vmax * x / (Km + x), 
                                        start = list(Vmax = 50, Km = 2)),
                     se = F, colour = 'black', size = 0.5)

# Finally, just call 'x' to show the plots:
x

# if you want to save it as a .pdf, just add
# ggsave("GGPLOT_FACET_EnzKin.pdf"), when you call 'x'
x + ggsave("GGPLOT_FACET_EnzKin.pdf")

Are my fitted enzyme kinetics lines significantly different?

Updated: 14th March 2017 because of changes in the geom_smooth() requirements.

R is good at statistics. I've been learning how to calculate F ratios in R to allow me to compare the lines fitted using glucose hydrogenase wild type and H297F mutant data that we graphed up in the previous blog post.  I want to add the p-value to the plot:

Note the addition of the F ratio and the p-value (very small)

The question from the statistical point of view: is fitting two lines better than fitting just one.

• Our null hypothesis says that one line fits all the data well. 
• The alternative hypothesis says that two lines fits the data well.

The question from an biochemistry experimental point of view: does the mutant glucose dehydrogenase H297F have different properties to wild type. If the mutation has no effect on the enzyme kinetics, then one line should fit all the data well. If not then two lines should fit the line better. This test won't tell us which or how the lines are different, just that they are.

Key steps to answer this question:

1. Fit one line to all the data  (fitTotal)
2. Fit the two lines to the two separate sets of data (fitWT and fitH297F). 
3. Extract the differences (deviance) between the data and the lines (aka sum of squares.)  
1. SS.null <- fitTotal$m$deviance() 
2. SS.alt <- fitWT$m$deviance() + fitH297F$m$deviance()
4. Extract the degrees of freedom for null and alternative hypothesis
1. df.null <- df.residual(fitTotal)
2. df.alt <- df.residual(fitWT) + df.residual(fitH297F)
5. Divide the deviance by the degrees of freedom which is this case is the number of samples less the number of lines we have drawn. This gives us the mean square - the average sum of squares.
1. mean.squares.null <- SS.null/df.null
2. mean.squares.alt <- SS.alt/df.alt
6. The F-ratio is the mean square for the null hypothesis (one line fits all the data) divided by mean squares for the alternative hypothesis.
1. Fratio <- mean.squares.null/mean.squares.alt
7. I calculated the p value if you know the F-ratio and the degrees of freedom using the pf() function.
1. p.val <- pf(Fratio, df.alt, df.null, log.p = TRUE)
8. I added the F-ratio and the p-value to the ggplot using the annotate() function.

If the variation is the basically the same in both cases - no extra value from plotting two lines, then the F-ratio will be approximately one. The higher the F-ratio the better. The p-value corresponding to the F-ratio depends on the degrees of freedom. R can calculate this.

I've learned most of this from Andy Field's book who has, in my opinion, written the easiest to read statistics books. One or two parts even made me laugh and I didn't expect that. Lots of the material is available through his website. You can download his information about ANOVA and F-ratios here.

In case you hadn't guessed from the plots, the F-ratio is significant and the null hypothesis is unlikely to be correct. Two lines describe the data better than one line. Our enzyme kinetic values are different.

Here is the script that does the calculations:

# START OF SCRIPT
# compare two models. 
# need to do an ANOVA...

library(ggplot2)

library(ggthemes)

# data from Dr C Bennett, University of Bath
# Link to published data: http://www.jbc.org/content/285/44/33701.full
# and her thesis: http://opus.bath.ac.uk/27220/

Enz <- c("WT","WT","WT","WT","WT",
         "WT","WT","WT","WT","WT",
         "WT","WT","WT",
         "H297F","H297F","H297F",
         "H297F","H297F","H297F",
         "H297F","H297F")
S <- c(2.00, 1.00, 0.60, 0.50, 0.40, 
         0.30, 0.20, 0.10, 0.09, 0.08, 
         0.06, 0.04, 0.02, 
         0.05, 0.10, 0.20, 
         0.30, 0.40, 0.50, 
         1.00, 2.00)
v <- c(59.01, 58.29, 54.17, 51.82, 49.76, 
         45.15, 36.88, 26.10, 23.50, 22.26, 
         16.45, 13.67, 6.14, 
         11.8, 19.9, 30.3, 
         36.6, 40.2, 42.1, 
         47.8, 50.0)

# assemble the data into a data.frame
enzdata <- as.data.frame(Enz)
enzdata$S <- S
enzdata$v <- v

# fit all the data to one NLS...
MMcurve<-formula(v~Vmax*S/(Km+S))
fitTotal <- nls(MMcurve, enzdata, start=list(Vmax=50,Km=0.2), subset=)
# residual sum-of-square: 328   fitTotal$m$deviance()
# from here: https://stat.ethz.ch/pipermail/r-help/2010-August/249065.html
SS.null <- fitTotal$m$deviance()
# from here: http://stackoverflow.com/questions/21734248/how-to-return-only-the-degrees-of-freedom-from-a-summary-of-a-regression-in-r
df.null <- df.residual(fitTotal)

# split the data and fit WT and H297F data
WT <- subset(enzdata, Enz=="WT")
fitWT <- nls(MMcurve, WT, start=list(Vmax=50,Km=0.2))
# residual sum-of-square: 23.36  

H297F <- subset(enzdata, Enz=="H297F")
fitH297F <- nls(MMcurve, H297F, start=list(Vmax=50,Km=0.2))
# residual sum-of-square: 5.523
# add these together because they are one (alternative) model
SS.alt <- fitWT$m$deviance() + fitH297F$m$deviance()
df.alt <- df.residual(fitWT) + df.residual(fitH297F)

mean.squares.null <- SS.null/df.null
mean.squares.alt <- SS.alt/df.alt

# calculate the F-statistic
# http://stats.stackexchange.com/questions/12398/how-to-interpret-f-and-p-value-in-anova
Fratio <- mean.squares.null/mean.squares.alt

# use pf(q, df1, df2, ncp, lower.tail = TRUE, log.p = FALSE)
p.val <- pf(Fratio, df.alt, df.null, log.p = TRUE)
p.val <- format(p.val, dig=3)
           
# Now, draw the plot AND add the p-value using annotate()

ggplot(data=enzdata,         
       aes(x=S,            
           y=v,            
           colour = Enz)) +  
  geom_point() +            
  xlab("Substrate (mM)") +  
  ylab("Velocity (uM/min/mg.enzyme)") +    
  ggtitle("Glucose Dehydrogenase \n wild type and mutant") +  

  geom_smooth(method = "nls", 
              method.args = list(formula = y ~ Vmax * x / (Km + x), 
                start = list(Vmax = 50, Km = 0.2)),
             se = F, size = 0.5,
             data = subset(enzdata, Enz=="WT"))  +

  geom_smooth(method = "nls", 
                method.args = list(formula = y ~ Vmax * x / (Km + x), 
                                   start = list(Vmax = 50, Km = 0.2)),
                se = F, size = 0.5, 
              data = subset(enzdata, Enz=="H297F")) +
  theme_few() +

            annotate("text", x = 1.15, y = 25, label = "F ratio: ") +
            annotate("text", x = 1.4, y = 25, label = format(Fratio, dig=3)) +
            annotate("text", x = 1.15, y = 20, label = "P-value: ") +
            annotate("text", x = 1.5, y = 20, label = p.val)

# END OF SCRIPT

Plotting two enzyme plots with ggplot...

UPDATE: As of ggplot 2.0.0, released in Dec 2015, to use the geom_smooth() ggplot function, there is a need to put the method arguments (method.args = list()) into a list as detailed below. I have corrected this plot to include the new information and it works now.

ggplot2 is a powerful graphing package in R. It was created by Hadley Wickham, an expert in R. Books and websites are dedicated to ggplot2.

Most of my work has used base graphics but I'm trying to learn how to use ggplot2 and the syntax involved in controlling it.

Here is a script and some of the graphs made with ggplot2. I have used enzyme data supplied by Charlie Bennett from University of Bath who is attending the R for Biochemists Training Day.

This script is show with various graphs intermingled. The idea is show how the various parts of the script change the graphs.

I hope it's useful. I have learned a little more about ggplot preparing it.

# SCRIPT STARTS

library(ggplot2)
library(ggthemes)

# data from Dr C Bennett, University of Bath
# Link to published data: http://www.jbc.org/content/285/44/33701.full
# and her thesis: http://opus.bath.ac.uk/27220/

# This is the data
Enz <- c("WT","WT","WT","WT","WT",
         "WT","WT","WT","WT","WT",
         "WT","WT","WT",
         "H297F","H297F","H297F",
         "H297F","H297F","H297F",
         "H297F","H297F")
sub <- c(2.00, 1.00, 0.60, 0.50, 0.40, 
         0.30, 0.20, 0.10, 0.09, 0.08, 
         0.06, 0.04, 0.02, 
         0.05, 0.10, 0.20, 
         0.30, 0.40, 0.50, 
         1.00, 2.00)
vel <- c(59.01, 58.29, 54.17, 51.82, 49.76, 
         45.15, 36.88, 26.10, 23.50, 22.26, 
         16.45, 13.67, 6.14, 
         11.8, 19.9, 30.3, 
         36.6, 40.2, 42.1, 
         47.8, 50.0)

# assemble the data into a data.frame
enzdata <- as.data.frame(Enz)
enzdata$sub <- sub
enzdata$vel <- vel

# plot the data with ggplot
# most of the syntax seems relatively easy to understand...
ggplot(data=enzdata,         # give the ggplot() function the data
       aes(x=sub,            # data.frame col with values for x-axis
           y=vel,            # and the y-axis
           colour = Enz)) +  # colour by WT or H297F **NOTE the "+"
  geom_point() +             # key function that show points
  xlab("Substrate (mM)") +   # label x-axis
  ylab("Velocity (uM/min/mg.enzyme)") +    # label y-axis
  ggtitle("Glucose Dehydrogenase \n wild type and mutant") 

# so far, so good. 

# we can create an object with the plot in it... 
enz.plot <- ggplot(data=enzdata,         
                   aes(x=sub,            
                       y=vel,            
                       colour = Enz)) +  
            geom_point() +            
            xlab("Substrate (mM)") +  
            ylab("Velocity (uM/min/mg.enzyme)") +    
            ggtitle("Glucose Dehydrogenase \n wild type and mutant")  


# then we can apply a different theme to change the style of the plot. 
enz.plot + theme_bw()

# or 
enz.plot + theme_few()

# or change arguments in the geom_point() function to make the points larger
enz.plot + geom_point(size=5)

# finally, we can add the enzyme kinetic lines using the geom_smooth() function
enz.plot  + geom_smooth(method = "nls", 
              method.args = list(formula = y ~ Vmax * x / (Km + x), 
                                 start = list(Vmax = 50, Km = 0.2)),
                        se = F, size = 0.5, 
                        data = subset(enzdata, Enz=="WT")) +
            geom_smooth(method = "nls", 
              method.args = list(formula = y ~ Vmax * x / (Km + x), 
                                 start = list(Vmax = 50, Km = 0.2)),
                        se = F, size = 0.5, 
                        data = subset(enzdata, Enz=="H297F")) +
            theme_few()

# SCRIPT ENDS

The following sources were also useful:

• Dean's script for plotting six enzyme plots with ggplot
• http://www.cookbook-r.com/Graphs/Bar_and_line_graphs_(ggplot2)/

Draw six enzymology graphs with ggplot...

Added March 2016:

The latest version of ggplot2 required a change in this script which has been added below. 

The plot with six separate enzymology lines on it looked a bit busy and difficult to appreciate. We made previous used base graphics to make six individual plots. Here Dean has used ggplot to make a nicer plot. 

First, we make an object containing the combined graph. It's a list. To learn more about data structures, please check out our earlier blog article. 

Then we use faceting to split the object into six separate graphs. 

Finally, we use the geom_smooth()function with the argument method = "nls" to add the lines. 

Here is the graph: 

Multi-panel enzymology data plotted with ggplot 

Is the ggplot graphic nicer than the base graphics plot? It's personal choice in large part.

Multi-panel enzymology data plotted with base graphics

Here is the script to make the multi panel plot with ggplot.

# The plot with the six lines is a bit busy. 
# Let's make a ggplot equivalent to the base R multiplot.
# Here we do it using faceting.
# To make the script free standing the reshaping code is included. 

library(reshape2)
library(ggplot2)
library (ggthemes)

# the data:
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)
no.Exp <- c("Exp 1", "Exp 2", "Exp 3", "Exp 4", "Exp 5", "Exp 6")
Sub <- c(0, 1, 2, 4, 8, 12, 16, 20, 30, 40)

enzdata <- as.data.frame(enzdata) # converted to a data.frame
colnames(enzdata) <- no.Exp # add column names
enzdata <- cbind(Sub, enzdata) # add the Substrate data
# melt the data from wide to long:
melted_data <- melt(enzdata, id.vars = "Sub", value.name = "v", variable.name = "Exp")

# the way ggplot draws colours depends on the required number of colours. The way it does it is just with equally spaced hues around the color wheel, starting from 15:
gg_color_hue <- function(n){
  hues = seq(15, 375, length = n + 1)
  hcl(h = hues, l = 65, c = 100)[1:n]
  }
  
# we have 6 experiments (so we want to create our own 6 colour "mini-palette")
cols <- gg_color_hue(6)  
# this gets the default ggplot colours when we want to plot 6 different variables.


# 1st, create an object with the ggplot 
# It has all the data simply as points, shaped and colour-coded by expt. of origin (as above). 
# We store it in an object called 'x':
x <- ggplot(data = melted_data, aes(x = Sub, y = v)) +
  geom_point(aes(colour = Exp, shape = Exp)) +
  theme_few() +
  ylab("velocity (nmol/s)") + xlab("substrate (mM)") +
  theme(axis.title.x = element_text(size = 15),
        axis.title.y = element_text(size = 15),
        axis.text.x = element_text(size = 15),
        axis.text.y = element_text(size = 15),
        plot.title = element_blank(),            # we don't want individual plot titles as the facet "strip" will give us this
        legend.position = "none",                # we don't want a legend either
        panel.border = element_rect(fill = NA, color = "darkgrey", size = 1.25, linetype = "solid"),
        axis.ticks = element_line(colour = 'darkgrey', size = 1.25, linetype = 'solid'))     # here, I just alter to colour and thickness of the plot outline and tick marks. You generally have to do this when faceting, as well as alter the text sizes (= element_text() in theme also)

This is what we have made so far

# Next, let's modify 'x' by faceting based on the Experiment Name, 
# We also specify that we want a panel 3 plots wide x 2 plots high.
# You can change this, obviously. Can also 'facet_grid' too:
x <- x + facet_wrap( ~ Exp, ncol = 3)

The facet_wrap() function turns it into this

# Finally, let's apply the Michaelis Menten fitting to our faceted data, adding the best fit line in each case:
x <- x + geom_smooth(method = "nls",                 
         method.args = list(formula = y ~ Vmax * x / (Km + x), 
                            start = list(Vmax = 50, Km = 2)),
                     se = F, colour = 'black', size = 0.5)  

# Finally, just call 'x' to show the plots:
x

# if you want to save it as a .pdf, just add + ggsave("GGPLOT_FACET_EnzKin.pdf"), when you call 'x'
x + ggsave("GGPLOT_FACET_EnzKin.pdf")

Plotting enzyme data with ggplot - Part I

ggplot is a very powerful graphics package.  The construction of the graphs is quite different to base graphics. "It uses the grammar of graphics" particularly layers. 

Dr Dean Hammond has written this script that illustrates how ggplot can be used to draw multiple enzymatic data plots with colours. Part II will draw six separate plots. 

Using ggplot sometimes requires some reshaping of the data (data munging or data wrangling) which is shown in the script below. This uses the melt() function.

The points are plotted alone first creating a layer with geom_point():
geom_point(aes(color = Exp, shape = Exp)).
aes() refers to the aesthetics of the plot.  

The line fitting is done within ggplot using geom_smooth():
geom_smooth(method = "nls", formula = y ~ Vmax * x / (Km + x), start = list(Vmax = 50, Km = 2),

I particularly like the ease with which the plot can be saved using ggsave():
ggsave("All_points_plus_fits.pdf")

The final graph is here:

# Following on from my code to multiplot 6 enzymology data-sets using base R in a for loop,
# here's how to create a single plot with all fitted curves on it too
# I have used ggplot and faceting.

# Some manipulation of the data is required, essentially to melt it from 'wide' to 'long' (reshape2)

# I personally don't like the standard ggplot look, so
# I also load package ggthemes and think Stephen Few's theme is the nicest (hence, theme_few())

# You may need to do this once: install.packages("ggplot2", "reshape2", "ggthemes")

library(reshape2)  # allows us to use the packages
library(ggplot2)
library(ggthemes)

# lets read the data in:
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)
no.Exp <- c("Exp 1", "Exp 2", "Exp 3", "Exp 4", "Exp 5", "Exp 6")
Sub <- c(0, 1, 2, 4, 8, 12, 16, 20, 30, 40)
 
# convert it to a data.frame for melting:
enzdata <- as.data.frame(enzdata)

# add column names:
colnames(enzdata) <- no.Exp

# add the Sunstrate data as an additional column to the data.frame:
enzdata <- cbind(Sub, enzdata)

# melt the data from wide to long:
melted_data <- melt(enzdata, id.vars = "Sub", value.name = "v", variable.name = "Exp")

# an examination of the two data files shows what has happened. 
View(enzdata)
View(melted_data)

Data before "melting"

Some of the data after melting

# The way ggplot draws colours depends on the required number of colours. 
# The way it does it is just with equally spaced hues around the color wheel, starting from 15:
gg_color_hue <- function(n){
  hues = seq(15, 375, length = n + 1)
  hcl(h = hues, l = 65, c = 100)[1:n]
  }
  
# we have 6 experiments (so we want to create our own 6 colour "mini-palette")
cols <- gg_color_hue(6)  
# this gets the default ggplot colours when we want to plot 6 different variables.

# Let's plot:
ggplot(melted_data, aes(x = Sub, y = v)) +
  ylab("velocity (nmol/s)") + xlab("substrate (mM)") +
  theme_few() +
  # let's plot the data simply as points, shaped and colour-coded by expt. of origin:
  geom_point(aes(color = Exp, shape = Exp)) +

Just the points - no fitted lines yet. 

  # then add each best-fit line from our Michaelis Menten nls equation (according to Expt), colouring based-on our colour palette to match ggplots own colour-coding:
  geom_smooth(method = "nls", formula = y ~ Vmax * x / (Km + x), start = list(Vmax = 50, Km = 2),
              se = F, colour = cols[1], size = 0.5, data = filter(melted_data, Exp == "Exp 1")) +
  geom_smooth(method = "nls", formula = y ~ Vmax * x / (Km + x), start = list(Vmax = 50, Km = 2),
              se = F, colour = cols[2], size = 0.5, data = filter(melted_data, Exp == "Exp 2")) +
  geom_smooth(method = "nls", formula = y ~ Vmax * x / (Km + x), start = list(Vmax = 50, Km = 2),
              se = F, colour = cols[3], size = 0.5, data = filter(melted_data, Exp == "Exp 3")) +
  geom_smooth(method = "nls", formula = y ~ Vmax * x / (Km + x), start = list(Vmax = 50, Km = 2),
              se = F, colour = cols[4], size = 0.5, data = filter(melted_data, Exp == "Exp 4")) +
  geom_smooth(method = "nls", formula = y ~ Vmax * x / (Km + x), start = list(Vmax = 50, Km = 2),
              se = F, colour = cols[5], size = 0.5, data = filter(melted_data, Exp == "Exp 5")) +
  geom_smooth(method = "nls", formula = y ~ Vmax * x / (Km + x), start = list(Vmax = 50, Km = 2),
              se = F, colour = cols[6], size = 0.5, data = filter(melted_data, Exp == "Exp 6")) +
  # Let's add a title to the plot
  ggtitle("Enzyme Kinetics Data") +
  # and save it as a .pdf (optional, obviously)
  ggsave("All_points_plus_fits.pdf")

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)