Measuring active ghrelin (calibration curve)

by Bartek Rajwa

Computing active ghrelin concentration

For an explanation of difference between acyl- (active) and des-acyl ghrelin (DAG) see AG-vs-DAG.

Similarly to a procedure for computation of a total gherkin concentration we will use R.

The data for the calibration curve should be read into a variable. For the purpose of demonstration we will use the calibration data below:

Sample name Dilution factor Ghrelin replicate A Ghrelin replicate B Concentration (pg/ml)
Blank 1 0.041 0.044 0
Standard 6 64 0.052 0.051 27.5
Standard 5 32 0.057 0.054 55
Standard 4 16 0.078 0.07 110
Standard 3 8 0.111 0.109 220
Standard 2 4 0.242 0.244 440
Standard 1 2 0.619 0.634 880
Standard 0 1 1.642 1.674 1760

We will use library “drc” to perform computations.


Loading the data… <- read.csv("Millipore_active_ghrelin_2017-02-12_standard_curve.csv") <-[,-c(1,2]
colnames( <- c("Concentration","Response_1", "Response_2")

Subtracting the blank

tmp <- apply([,2:3], 1, function(x){[1,2:3]} )[,2:3] <- matrix(data=unlist(tmp), ncol=2, byrow=T) <-[-1,]

Reorganize the data – from wide to the long format, suitable for curve fitting <- reshape(, varying=c("Response_1","Response_2"),  direction="long",  v.names=c("Response")) <-[,c("Concentration", "Response")]

Fitting the model (4-parameter log-logistic function). We use the following model:

y=c+\frac{d-c}{1+\exp \left( b\left( \log (x)-\log (e) \right) \right)} <- drm(Response ~ Concentration, data =, fct = LL.4())

Plotting the calibration curve

par(pty="s", mar=c(5,6,1,1), mgp = c(4, 1, 0))
plot(, type="confidence", cex.lab=2, axes=T, cex.axis=1.2, log="xy")
plot(, type="all", add=T, pch=21, col="red", lwd=1, cex=2, bg="green")


At this point, we can read in the actual measured data.

ghrelin <- read.csv("Millipore_active_ghrelin_2017-02-12_samples.csv")
ghrelin <- ghrelin[1:12,c(1,4,5)]
ghrelin[,1] <- unclass(ghrelin[,1])
ghrelin[,1] <- c(rep("Fasted",6), rep("Re-fed",6)) 
colnames(ghrelin) <- c("Feeding", "OD1", "OD2")
tmp <- apply(ghrelin[,-1], 1, function(x){mean(x)})
ghrelin <- data.frame(ghrelin[,1], tmp)
colnames(ghrelin) <- c("Feeding", "OD")

The computation of the concentration is performed by the “ED” function of the “drc” library

hormone.calc <- ED(, respLev=ghrelin[,2], interval="delta", type="absolute", level=0.95)
colnames(hormone.calc) <- c("yi", "sei", "lb", "ub")

After the computation is done, we may want to find the best estimate of active ghrelin both investigated rat group. We use random effect model calculation available in R-library “metafor”

fasted <- rma(yi=yi, sei=sei, data=hormone.calc[1:6,], method="DL")
forest(fasted, slab=c("Fasted rat 1","Fasted rat 2", "Fasted rat 3", "Fasted rat 4", "Fasted rat 5", "Fasted rat 6"), efac=2, width=2)
re.fed <- rma(yi=yi, sei=sei, data=hormone.calc[7:12,], method="DL")
forest(re.fed, slab=c("Re-fed rat 1","Re-fed rat 2", "Re-fed rat 3", "Re-fed rat 4", "Re-fed rat 5", "Re-fed rat 6"), efac=2, width=2)

The results can be plotted. We will plot the measured values, and the established confidence intervals.

par(mar=c(5,5,1,1),mgp = c(3, 1, 0))
plot(1, hormone.calc[1,1], xlim=c(1,14), ylim=c(100, 4000), xlab="", ylab="Concentration", pch=15, cex=2, xaxt="n", cex.lab=1.5)
axis(side=1, at=c(1:14), labels=c(paste("Fasted rat", seq(1:6)), paste("Re-red rat", seq(1:6)), "Fasted average", "Re-fed average"), las=3)
points(1:6, hormone.calc[1:6,1], cex=2, pch=15, col="red")
points(7:12, hormone.calc[7:12,1], cex=2, pch=15, col="green")
arrows(x0=c(1:12), y0=hormone.calc[,3], x1=c(1:12), y1=hormone.calc[,4],  length=0.15, angle=90, code=3)
points(13, fasted$b, cex=2, pch=19, col="red")
points(14, re.fed$b, cex=2, pch=19, col="green")
arrows(x0=13, y0=fasted$, x1=13, y1=fasted$ci.ub, length=0.15, angle=90, code=3)
arrows(x0=14, y0=re.fed$, x1=14,$ci.ub, length=0.15, angle=90, code=3)


Created on , Last modified on