0

我有一些代码可以将 dr4pl 模型拟合到数据集中的剂量反应曲线，并生成一个包含相关 IC50 的表格，或者将曲线拟合到 ggplot。

但是，当 dr4pl 模型无法将曲线拟合到特定数据位时，它会失败，并且代码不会运行。因此，如果我在包含 100 条曲线的数据上运行代码，其中一条包含错误数据，那么整个过程都会失败，而不会告诉我哪里出了问题。

我需要帮助弄清楚如何 i) 让它运行，并忽略它无法拟合的曲线，或者 ii) 标记数据集的有问题的部分。

library(dr4pl)
library(car)

这是可以添加到 ggplot 的模型：

    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(dat2, aes(dose, medPOC, col=cell_palbo))+
  geom_point(size=2.5) +
  geom_smooth(method="dr4pl",se=F)+ 
  coord_trans(x="log10")+
  scale_x_continuous(breaks = c(0.01, 0.1, 1, 10))+
  theme_bw()+
  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"))

示例图： 左=曲线拟合成功，右=曲线拟合失败

这产生了一个函数，允许我将模型传递给数据集并生成 IC50：#builds IC50 table

    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))
}

这将在数​​据集上运行该函数以生成包含 IC50 的表：

dfIC50 <- multiIC(data = dataset, 
                  colDose = "dose", 
                  colResp = "medPOC", 
                  colID = "ID", 
                  inhib.percent = 50)

加载两个不同的数据集，每个数据集有 2 条曲线。Dat1 将使函数失败并且不生成输出。Dat2 工作得很好。

dat1<-structure(list(cell_palbo = c("T47D-", "T47D-", "T47D-", "T47D-", 
                              "T47D-", "T47D-", "T47DpR-", "T47DpR-", "T47DpR-", "T47DpR-", 
                              "T47DpR-", "T47DpR-"), medPOC = c(1, 0.859642814335371, 0.826920771904591, 
                                                                0.611051180630469, 0.0391945343401654, 0.284310200167805, 1, 
                                                                0.905880254474107, 0.790891253938998, 0.624650692669005, 1.01212913966348, 
                                                                0.365955169748499), dose = c(0.01, 0.1, 0.3, 1, 10, 3, 0.01, 
                                                                                             0.1, 0.3, 1, 10, 3), ID = c("SN1051760333 - T47D-", "SN1051760333 - T47D-", 
                                                                                                                         "SN1051760333 - T47D-", "SN1051760333 - T47D-", "SN1051760333 - T47D-", 
                                                                                                                         "SN1051760333 - T47D-", "SN1051760333 - T47DpR-", "SN1051760333 - T47DpR-", 
                                                                                                                         "SN1051760333 - T47DpR-", "SN1051760333 - T47DpR-", "SN1051760333 - T47DpR-", 
                                                                                                                         "SN1051760333 - T47DpR-")), row.names = c(NA, -12L), class = c("data.table", 
                                                                                                                                                                                        "data.frame"))

dat2<-structure(list(cell_palbo = c("T47D-", "T47D-", "T47D-", "T47D-", 
                              "T47D-", "T47D-", "T47DpR-", "T47DpR-", "T47DpR-", "T47DpR-", 
                              "T47DpR-", "T47DpR-"), medPOC = c(1, 0.897544504013507, 0.0792630550042949, 
                                                                0.0194307040668227, 0.00746423387932822, 0.0135067089244987, 
                                                                1, 0.0396359365825015, 0.0404633534404527, 0.00371003042758768, 
                                                                0.00184166978060108, 0.0024555597074681), dose = c(0.01, 0.1, 
                                                                                                                   0.3, 1, 10, 3, 0.01, 0.1, 0.3, 1, 10, 3), ID = c("SN1023563430 - T47D-", 
                                                                                                                                                                    "SN1023563430 - T47D-", "SN1023563430 - T47D-", "SN1023563430 - T47D-", 
                                                                                                                                                                    "SN1023563430 - T47D-", "SN1023563430 - T47D-", "SN1023563430 - T47DpR-", 
                                                                                                                                                                    "SN1023563430 - T47DpR-", "SN1023563430 - T47DpR-", "SN1023563430 - T47DpR-", 
                                                                                                                                                                    "SN1023563430 - T47DpR-", "SN1023563430 - T47DpR-")), row.names = c(NA, 
                                                                                                                                                                                                                                        -12L), class = c("data.table", "data.frame"))

正如您将看到的，如果您运行这些并尝试生成 IC50 表：

dfIC50a <- multiIC(data = dat1, 
                   colDose = "dose", 
                   colResp = "medPOC", 
                   colID = "ID", 
                   inhib.percent = 50)

dfIC50b <- multiIC(data = dat2, 
                   colDose = "dose", 
                   colResp = "medPOC", 
                   colID = "ID", 
                   inhib.percent = 50)

我真的很感激一些帮助，在此先感谢。

4

0 回答 0