剂量反应曲线的参数和 AUC 和 IC50答案

【问题标题】：Parameters and AUC and IC50 of a dose response curve剂量反应曲线的参数和 AUC 和 IC50
【发布时间】：2016-11-02 00:12:13
【问题描述】：

我有这些剂量反应数据：

df <- data.frame(viability=c(14,81,58,78,71,83,64,16,32,100,100,81,86,83,100,90,15,100,38,100,91,84,92,100),
                 dose=c(10,0.62,2.5,0.16,0.039,0.0024,0.0098,0.00061,10,0.62,2.5,0.16,0.039,0.0024,0.0098,0.00061,10,0.62,2.5,0.16,0.039,0.0024,0.0098,0.00061),
                 stringsAsFactors=F)

然后我使用 drc 包的 drm 函数来拟合这些数据的对数逻辑曲线：

library(drc)
fit <- drm(viability~dose,data=df,fct=LL.4(names=c("slope","low","high","ED50")),type="continuous")

> summary(fit)

Model fitted: Log-logistic (ED50 as parameter) (4 parms)

Parameter estimates:

                  Estimate Std. Error  t-value p-value
slope:(Intercept)  5.15328   18.07742  0.28507  0.7785
low:(Intercept)   20.19430   12.61122  1.60130  0.1250
high:(Intercept)  83.33181    4.96736 16.77586  0.0000
ED50:(Intercept)   2.98733    1.99685  1.49602  0.1503

Residual standard error:

 21.0743 (20 degrees of freedom)

然后我生成预测，以便绘制曲线：

pred.df <- expand.grid(dose=exp(seq(log(max(df$dose)),log(min(df$dose)),length=100))) 
pred <- predict(fit,newdata=pred.df,interval="confidence") 
pred.df$viability <- pmax(pred[,1],0)
pred.df$viability <- pmin(pred.df$viability,100)
pred.df$viability.low <- pmax(pred[,2],0)
pred.df$viability.low <- pmin(pred.df$viability.low,100)
pred.df$viability.high <- pmax(pred[,3],0)
pred.df$viability.high <- pmin(pred.df$viability.high,100)

我还使用 PharmacoGx Bioconductor 包来计算曲线及其上限和下限的 AUC 和 IC50：

library(PharmacoGx)
auc.mid <- computeAUC(rev(pred.df$dose),rev(pred.df$viability))/((max(pred.df$viability)-min(pred.df$viability))*(max(pred.df$dose)-min(pred.df$dose)))
auc.low <- computeAUC(rev(pred.df$dose),rev(pred.df$viability.low))/((max(pred.df$viability.low)-min(pred.df$viability.low))*(max(pred.df)-min(pred.df$dose)))
auc.high <- computeAUC(rev(pred.df$dose),rev(pred.df$viability.high))/((max(pred.df$viability.high)-min(pred.df$viability.high))*(max(pred.df$dose)-min(pred.df$dose)))
ic50.mid <- computeIC50(rev(pred.df$dose),rev(pred.df$viability))
ic50.low <- computeIC50(rev(pred.df$dose),rev(pred.df$viability.low))
ic50.high <- computeIC50(rev(pred.df$dose),rev(pred.df$viability.high))

创建一个包含所有参数的表格，以便我可以将所有内容绘制在一起：

ann.df <- data.frame(param=c("slope","low","high","ED50","auc.mid","auc.high","auc.low","ic50.mid","ic50.high","ic50.low"),value=signif(c(summary(fit)$coefficient[,1],auc.mid,auc.high,auc.low,ic50.mid,ic50.high,ic50.low),2),stringsAsFactors=F)

最后绘制出来：

library(ggplot2)
library(grid)
library(gridExtra)
pl <- ggplot(df,aes(x=dose,y=viability))+geom_point()+geom_ribbon(data=pred.df,aes(x=dose,y=viability,ymin=viability.low,ymax=viability.high),alpha=0.2)+labs(y="viability")+
  geom_line(data=pred.df,aes(x=dose,y=viability))+coord_trans(x="log")+theme_bw()+scale_x_continuous(name="dose",breaks=sort(unique(df$dose)),labels=format(signif(sort(unique(df$dose)),3),scientific=T))
ggdraw(pl)+draw_grob(tableGrob(ann.df,rows=NULL),x=0.1,y=0.175,width=0.3,height=0.4)

这给出了：

我的问题是：

我认为斜率应该是负数。怎么是5.2？
auc.mid、auc.high 和 auc.low累积为：

auc.mid

分别给出 21.47818、37.52389 和 2.678228。

由于这些不在 [0,1] 范围内，我认为将它们除以最高相应生存能力下的区域将给出我正在寻找的结果，即相对 AUC，但这些值相对于如图所示。那么这些 AUC 是什么？

还有，auc.mid > auc.low > auc.high 怎么来的？我认为应该是auc.high > auc.mid > auc.low

IC50 值似乎也有点低。它们有意义吗？

额外问题：如何避免图中slope、low、high、ED50、ic50.mid 和 ic50.high 中的尾随零？

【问题讨论】：

作为你的斜坡summary(fit)$coefficient[,1]。在大多数模型中，这将是截距。

标签： r ggplot2 auc drc

【解决方案1】：

你拉出来的参数是山坡参数，或者指数中浓度变量前面的系数，不是曲线的实际斜率。
对于曲线上方的区域，提供的 AUC 在 [0-100] 范围内。我运行代码并得到 auc.low>auc.mid>auc.high 的顺序。传统上报告响应曲线下的面积，或 1-活力。
需要注意的是，PharmacoGx 包使用 3 参数山坡模型，类似于drc 中的 LL.3。因此，该图将不对应于由PharmacoGx 拟合的函数来计算 IC50 或 AUC。

来源：PharmacoGx 开发。

【讨论】：