我有一些代码可以将 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)
我真的很感激一些帮助，在此先感谢。