2

我有这些剂量反应数据：

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)

这使： 在此处输入图像描述

我的问题是：

  1. 我认为斜率应该是负数。怎么是5.2？

  2. ,auc.mid和累积为auc.high：auc.low

    auc.mid <- computeAUC(rev(pred.df$dose),rev(pred.df$viability)) auc.low <- computeAUC(rev(pred.df$dose),rev(pred.df$viability.low )) auc.high <- computeAUC(rev(pred.df$dose),rev(pred.df$viability.high))

分别给出 21.47818、37.52389 和 2.678228。

由于这些不在 [0,1] 范围内，我认为将它们除以最高相应生存能力下的区域将给出我正在寻找的东西，即相对 AUC，但这些值相对于数字而言似乎太低了显示。那么这些 AUC 是什么？

还有，怎么来的auc.mid> auc.low> auc.high？我觉得应该是auc.high>>auc.midauc.low

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

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

4

1 回答 1

3

  1. 您要提取的参数是山坡参数，或者指数中浓度变量前面的系数，而不是曲线的实际斜率。

  2. 对于曲线上方的区域，提供的 AUC 在 [0-100] 范围内。我运行了代码并得到了 auc.low>auc.mid>auc.high 的顺序。传统上报告响应曲线下的面积，或 1-生存力。

  3. 需要注意的是，该PharmacoGx包使用 3 参数山坡模型，类似于 LL.3 in drc。因此，该图将不对应于PharmacoGx计算 IC50 或 AUC 的函数拟合。

资料来源：PharmacoGx 开发。

于 2017-04-19T20:41:32.217 回答