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.

install.packages("drc")
require(drc}

hormone.data <- hormone.data[,-c(1,2]
colnames(hormone.data) <- c("Concentration","Response_1", "Response_2")

Subtracting the blank

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

Reorganize the data – from wide to the long format, suitable for curve fitting

hormone.data <- reshape(hormone.data, varying=c("Response_1","Response_2"),  direction="long",  v.names=c("Response"))
hormone.data <- hormone.data[,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)}$

hormone.data.model <- drm(Response ~ Concentration, data = hormone.data, fct = LL.4())
summary(hormone.data.model)

Plotting the calibration curve

par(pty="s", mar=c(5,6,1,1), mgp = c(4, 1, 0))
plot(hormone.data.model, type="confidence", cex.lab=2, axes=T, cex.axis=1.2, log="xy")
plot(hormone.data.model, 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 <- 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(hormone.data.model, 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$ci.lb, x1=13, y1=fasted$ci.ub, length=0.15, angle=90, code=3)
arrows(x0=14, y0=re.fed$ci.lb, x1=14, y1=re.red$ci.ub, length=0.15, angle=90, code=3)