Drawing a cell death curve and calculating an LD50

I'm going to show you the this analysis before the script.
We do a lot of drug testing in our laboratory at Cardiff and it's a valuable thing for biochemists and pharmacologists. To do the analysis, most of the lab usually use a different (expensive) statistical package but I like R.
The experimental set up, done by a student in the lab, is as follows:

• using some cells and add various concentrations of a novel drug
• leave for 48 hours
• then measure how many cells are dead
• the experiment was done four times with the doses improved each time

Here is the graph and calculated LD50 - a measure of how good the drug is:

Here is the script I used to generate this graph:

### START of SCRIPT 

###  this is the data   ###

# data from four experiments

conc <- c(5.00E-07, 1.00E-06, 1.00E-05, 

          5.00E-07, 1.00E-06, 5.00E-06, 1.00E-05, 2.00E-05, 

          5.00E-07, 1.00E-06, 2.50E-06, 5.00E-06, 1.00E-05, 

          5.00E-07, 1.00E-06, 2.50E-06, 5.00E-06, 1.00E-05)

dead.cells <- c(34.6, 47.7, 81.7, 

                37.6, 55.7, 89.1, 84.3, 85.2, 

                34.4, 46.1, 76.2, 84.3, 84.1, 

                24.5, 26.1, 60.6, 82.7, 87)

# transform the data to make it postive and put into a data frame for fitting 

data <- as.data.frame(conc)   # create the data frame

data$dead.cells <- dead.cells

data$nM <- data$conc * 1000000000

data$log.nM <- log10(data$nM) 

###  fit the data  ###

# make sure logconc remains positive, otherwise multiply to keep positive values

# (such as in this example where the conc is multiplied by 1000000

fit <- nls(dead.cells ~ bot+(top-bot)/(1+(log.nM/LD50)^slope),

             data = data,

             start=list(bot=20, top=95, LD50=3, slope=-12))

###  Plot the results  ###

#this lets you graph your calculated equations nice and pretty

x <- seq(0,6, length=100)

y <- (coef(fit)["bot"]+(coef(fit)["top"]-coef(fit)["bot"])/(1+(x/coef(fit)["LD50"])^coef(fit)["slope"]))

m <- coef(fit)

# plot the points first

plot(data$nM, data$dead.cells, 

     log="x", 

     main="Drug Dose Response and LD50", 

     xlab="Drug concentration (microM)", 

     ylab="Dead cells (% of cells)",

     xlim= c(500, 10000),

     ylim= c(20,100),

     xaxt = "n")   # suppresses the labels on the x-axis

# adds the x-axis labels they way I want it to look

axis(1, at=c(500, 1000, 2500, 5000, 10000),  lab=c("0.5","1","2.5","5","10"))

# adds the fitted non-linear line 

lines(10^x,y, lty="dotted", col="red", lwd =2)

# add the LD50 in the legend which allows nice positioning. 

rp = vector('expression',1)

rp[1] = substitute(expression(LD50(microM) == MYVALUE), 

                   list(MYVALUE = format((10^m[3])/1000,dig=3)))[2]

legend('topleft', legend = rp, bty = 'n')

# END OF SCRIPT

To prepare this script, I used the following sources:

• http://www.carlyhuitema.com/r.html
• http://www.harding.edu/fmccown/r/
• http://lukemiller.org/index.php/2012/10/adding-p-values-and-r-squared-values-to-a-plot-using-expression/

Thanks Dr Carly Huitema, Dr Frank McCown and Dr Luke Miller for the resources. 

Thanks to Mel for the data.