在 R 中如何将来自 dr4pl 模型的剂量反应 IC50 值合并到 ggplot2 曲线中？

问题描述

我在此站点上得到了帮助，将 dr4pl（替代 '''drc'''）4 参数逻辑剂量反应模型与 ggplot 结合使用。它工作得很好，但我希望将每条曲线的 IC50 合并到图中，并输出一个带有值的表格。 示例数据：

    #load packages
    library(ggplot2)
    library(data.table)
    library(dr4pl)
    library(car)
    
    #example data
curve = c("C1","C1","C2","C3","C3")
POC =c(1.07129314,0.91126280,0.97914297,0.95904437,0.88509670,0.84338263,0.75843762,0.61319681,0.52635571,0.84563087,1.24435113,1.11757648,0.82383523,0.82763447,0.72585483,0.31953609,0.15056989,0.10057988,0.57384256,0.65984339,0.81439758,0.84572057,0.62797088,0.30800934,0.08957274,0.06360764,0.04451161)
dose = c(0.078125,0.156250,0.312500,0.625000,1.250000,2.500000,5.000000,10.000000,20.000000,0.078125,20.000000)
example2<-data.frame(POC,dose,curve)

#this code will write  model that can be incorporated into 
predict.dr4pl <- function (object,newdata=NULL,se.fit=FALSE,level,interval) {
  xseq <- if (is.null(newdata)) object$data$Dose else newdata$x
  pred <- MeanResponse(xseq,object$parameters)
  if (!se.fit) {
    return(pred)
  }
  qq <- qnorm((1+level)/2)
  se <- sapply(xseq,function(x) car::deltaMethod(object,"UpperLimit + (LowerLimit - UpperLimit)/(1 + (x/IC50)^Slope)")[["Estimate"]])
  return(list(fit=data.frame(fit=pred,lwr=pred-qq*se,upr=pred+qq*se),se.fit=se))
}

#ggplot the example
ggplot(example2,aes(dose,POC,col=curve))+ 
  geom_point(size=4,shape=1) +
  geom_smooth(method="dr4pl",se=F)+ 
  coord_trans(x="log10")+
  theme_bw()+
  scale_x_continuous(breaks = c(0.01,0.1,1,10,100))+
  theme(plot.title = element_text(lineheight = 0.9,face="bold",size=20,hjust=0.5))+
  ggtitle("Dose Response")+
  theme(axis.title = element_text(face="bold",size = 14))+
  theme(axis.text = element_text(face="bold",size = 12,colour="black"))

以上将生成一个包含三条曲线的图，每条曲线都有一个拟合的 4pl 模型。我需要 IC50 显示在图中，或其他地方的表格中，以便我可以自己将其放入图中。

提前感谢您的帮助。

2021 年 6 月编辑：

有助于调试所提供脚本的新数据集：

new_data<-structure(list(cell_line = c("MCF7","MCF7","MCF7-A","T47D","T47D-A","T47D-pR","T47D-pR"),compound = c("Drug A","Drug A","Drug A"),dose = c(0.01,0.15625,0.3125,0.625,1.25,2.5,5,20,0.01,20),parental = c("MCF7","T47D"),medPOC = c(1.00934409313115,0.956036980819649,0.954030667407294,0.711853180626466,0.241538799146477,0.0567749571032392,0.0155968806238752,0.00581957552880206,0.00378717257100532,0.0059174440467226,0.945368038063252,0.929888217437452,0.835779167562444,0.586348465064098,0.262177612337144,0.0824082662867902,0.0308960382876894,0.00894423842501589,0.0061419835150226,0.00374687836102447,0.964364310824993,0.95414565025418,1.00205847809834,0.973855651871399,0.876721635476499,0.782657842062487,0.485817405545171,0.240963289831941,0.0807368351333095,0.0218556352405241,0.988828883403082,1.12687653176236,1.14174538911381,1.09793468012842,0.9371711107187,0.683580746738641,0.349521000688784,0.153807969259144,0.0557784035696989,0.00951109986306135,0.945736502721647,1.03384512851939,1.13359282933971,1.02798863227664,0.731271070765377,0.394548651675076,0.139037002094163,0.0594879837491901,0.0233289720871743,0.00163350004591345),sd = c(0.0487946611701357,0.0652579737057601,0.0447000542436252,0.0761742909565082,0.0271364685090581,0.0157757699801731,0.00457291598057361,0.00618036276893572,0.0031691999231949,0.00189564811623345,0.0286788195412834,0.0726775067251048,0.0598425497920862,0.11240260462323,0.0876844650331781,0.016269001522037,0.00768323265792931,0.0027283359541071,0.00530943514487439,0.00340553154613564,0.0513524691501586,0.0984457469912835,0.0451143466008112,0.0818218693028968,0.0934768958498176,0.0540191453047165,0.0608084812187634,0.0238245870875823,0.011089245968621,0.00995992003975256,0.0717550694254788,0.0510188619416171,0.0901696938469758,0.100897280181198,0.0426100049631484,0.0403803080203038,0.0383100464040449,0.0383778124586436,0.0153293502919647,0.0108310224482204,0.110583411728824,0.0714281429520549,0.127292213170872,0.121050868088146,0.15773243775593,0.144314150863785,0.0584119879286731,0.0179382293021653,0.0142390344682626,0.0148403765582003),dose1000 = c(10,78.125,156.25,312.5,625,1250,2500,5000,10000,20000,20000)),row.names = c(NA,-50L),groups = structure(list(cell_line = c("MCF7",.rows = structure(list(1L,2L,3L,4L,5L,6L,7L,8L,9L,10L,11L,12L,13L,14L,15L,16L,17L,18L,19L,20L,21L,22L,23L,24L,25L,26L,27L,28L,29L,30L,31L,32L,33L,34L,35L,36L,37L,38L,39L,40L,41L,42L,43L,44L,45L,46L,47L,48L,49L,50L),ptype = integer(0),class = c("vctrs_list_of","vctrs_vctr","list"))),class = c("tbl_df","tbl","data.frame"),.drop = TRUE),class = c("grouped_df","tbl_df","data.frame"))

解决方法

似乎 dr4pl 不像 drc 包通过 curveid 参数那样处理多剂量响应。

我的解决方案是使用 dr4pl::IC() 提取 IC50，然后将它们作为垂直线添加到图中。

这是一个提取 IC50 并将其保存为 data.frame 的函数：

multiIC <- function(data,colDose,colResp,colID,inhib.percent,...) {
  
  # Get curve IDs
  locID <- unique(data[[colID]])
  
  # Prepare a vector to store IC50s
  locIC <- rep(NA,length(locID))
  
  # Calculate IC50 separately for every curve
  for (ii in seq_along(locID)) {
    # Subset a single dose response
    locSub <- data[get(colID) == locID[[ii]],]
    
    # Calculate IC50
    locIC[[ii]] <- dr4pl::IC(
      dr4pl::dr4pl(dose = locSub[[colDose]],response = locSub[[colResp]],...),inhib.percent)
  }
  
  return(data.frame(id = locID,x = locIC))
}

如下使用：

dfIC50 <- multiIC(data = example2,colDose = "dose",colResp = "POC",colID = "curve",inhib.percent = 50)

获得：

> dfIC50
  id         x
1 C1 29.065593
2 C2  3.269516
3 C3  2.186479

然后将此行添加到您的 ggplot：

geom_vline(data = dfIC50,aes(xintercept = x,color = id))