（在 R 中）如何将剂量反应 IC50 值从 dr4pl 模型合并到 ggplot2 曲线中？
(in R) How to incorporate dose response IC50 values, from dr4pl models into a ggplot2 curve?
我在这个网站上得到帮助，将 dr4pl（'''drc''' 的替代品）4 参数逻辑剂量反应模型与 ggplot 相结合。它工作得很好，但我希望将每条曲线的 IC50 合并到图中，并输出 table 和这些值。
示例数据：
#load packages
library(ggplot2)
library(data.table)
library(dr4pl)
library(car)
#example data
curve = c("C1","C1","C1","C1","C1","C1","C1","C1","C1","C2","C2","C2","C2","C2","C2","C2","C2","C2","C3","C3","C3","C3","C3","C3","C3","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,0.156250,0.312500,0.625000,
1.250000,2.500000,5.000000,10.000000,20.000000,0.078125,
0.156250,0.312500,0.625000,1.250000,2.500000,5.000000,10.000000,
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 出现在情节中，或者出现在其他地方的 table 中，这样我就可以自己将它放入情节中。
在此先感谢您的帮助。
编辑 2021 年 6 月：
帮助调试脚本的新数据集：
new_data<-structure(list(cell_line = c("MCF7", "MCF7", "MCF7", "MCF7",
"MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7-A", "MCF7-A",
"MCF7-A", "MCF7-A", "MCF7-A", "MCF7-A", "MCF7-A", "MCF7-A", "MCF7-A",
"MCF7-A", "T47D", "T47D", "T47D", "T47D", "T47D", "T47D", "T47D",
"T47D", "T47D", "T47D", "T47D-A", "T47D-A", "T47D-A", "T47D-A",
"T47D-A", "T47D-A", "T47D-A", "T47D-A", "T47D-A", "T47D-A", "T47D-pR",
"T47D-pR", "T47D-pR", "T47D-pR", "T47D-pR", "T47D-pR", "T47D-pR",
"T47D-pR", "T47D-pR", "T47D-pR"), compound = c("Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A"), dose = c(0.01, 0.078125, 0.15625, 0.3125, 0.625,
1.25, 2.5, 5, 10, 20, 0.01, 0.078125, 0.15625, 0.3125, 0.625,
1.25, 2.5, 5, 10, 20, 0.01, 0.078125, 0.15625, 0.3125, 0.625,
1.25, 2.5, 5, 10, 20, 0.01, 0.078125, 0.15625, 0.3125, 0.625,
1.25, 2.5, 5, 10, 20, 0.01, 0.078125, 0.15625, 0.3125, 0.625,
1.25, 2.5, 5, 10, 20), parental = c("MCF7", "MCF7", "MCF7", "MCF7",
"MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7",
"MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7",
"T47D", "T47D", "T47D", "T47D", "T47D", "T47D", "T47D", "T47D",
"T47D", "T47D", "T47D", "T47D", "T47D", "T47D", "T47D", "T47D",
"T47D", "T47D", "T47D", "T47D", "T47D", "T47D", "T47D", "T47D",
"T47D", "T47D", "T47D", "T47D", "T47D", "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, 10, 78.125, 156.25,
312.5, 625, 1250, 2500, 5000, 10000, 20000, 10, 78.125, 156.25,
312.5, 625, 1250, 2500, 5000, 10000, 20000, 10, 78.125, 156.25,
312.5, 625, 1250, 2500, 5000, 10000, 20000, 10, 78.125, 156.25,
312.5, 625, 1250, 2500, 5000, 10000, 20000)), row.names = c(NA,
-50L), groups = structure(list(cell_line = c("MCF7", "MCF7",
"MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7",
"MCF7-A", "MCF7-A", "MCF7-A", "MCF7-A", "MCF7-A", "MCF7-A", "MCF7-A",
"MCF7-A", "MCF7-A", "MCF7-A", "T47D", "T47D", "T47D", "T47D",
"T47D", "T47D", "T47D", "T47D", "T47D", "T47D", "T47D-A", "T47D-A",
"T47D-A", "T47D-A", "T47D-A", "T47D-A", "T47D-A", "T47D-A", "T47D-A",
"T47D-A", "T47D-pR", "T47D-pR", "T47D-pR", "T47D-pR", "T47D-pR",
"T47D-pR", "T47D-pR", "T47D-pR", "T47D-pR", "T47D-pR"), compound = c("Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A"), dose = c(0.01, 0.078125, 0.15625, 0.3125, 0.625,
1.25, 2.5, 5, 10, 20, 0.01, 0.078125, 0.15625, 0.3125, 0.625,
1.25, 2.5, 5, 10, 20, 0.01, 0.078125, 0.15625, 0.3125, 0.625,
1.25, 2.5, 5, 10, 20, 0.01, 0.078125, 0.15625, 0.3125, 0.625,
1.25, 2.5, 5, 10, 20, 0.01, 0.078125, 0.15625, 0.3125, 0.625,
1.25, 2.5, 5, 10, 20), .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"))), row.names = c(NA, -50L), class = c("tbl_df",
"tbl", "data.frame"), .drop = TRUE), class = c("grouped_df",
"tbl_df", "tbl", "data.frame"))
似乎 dr4pl 不像 drc 包通过 curveid 参数处理多剂量反应。
我的解决方案是使用 dr4pl::IC() 提取 IC50，然后将它们作为垂直线添加到图中。
这是一个提取 IC50 并将其保存为 data.frame 的函数：
library(data.table)
library(dr4pl)
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))
}
输入数据应该是data.table。使用代码如下：
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))
我在这个网站上得到帮助，将 dr4pl（'''drc''' 的替代品）4 参数逻辑剂量反应模型与 ggplot 相结合。它工作得很好，但我希望将每条曲线的 IC50 合并到图中，并输出 table 和这些值。
示例数据：
#load packages
library(ggplot2)
library(data.table)
library(dr4pl)
library(car)
#example data
curve = c("C1","C1","C1","C1","C1","C1","C1","C1","C1","C2","C2","C2","C2","C2","C2","C2","C2","C2","C3","C3","C3","C3","C3","C3","C3","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,0.156250,0.312500,0.625000,
1.250000,2.500000,5.000000,10.000000,20.000000,0.078125,
0.156250,0.312500,0.625000,1.250000,2.500000,5.000000,10.000000,
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 出现在情节中，或者出现在其他地方的 table 中，这样我就可以自己将它放入情节中。
在此先感谢您的帮助。
编辑 2021 年 6 月：
帮助调试脚本的新数据集：
new_data<-structure(list(cell_line = c("MCF7", "MCF7", "MCF7", "MCF7",
"MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7-A", "MCF7-A",
"MCF7-A", "MCF7-A", "MCF7-A", "MCF7-A", "MCF7-A", "MCF7-A", "MCF7-A",
"MCF7-A", "T47D", "T47D", "T47D", "T47D", "T47D", "T47D", "T47D",
"T47D", "T47D", "T47D", "T47D-A", "T47D-A", "T47D-A", "T47D-A",
"T47D-A", "T47D-A", "T47D-A", "T47D-A", "T47D-A", "T47D-A", "T47D-pR",
"T47D-pR", "T47D-pR", "T47D-pR", "T47D-pR", "T47D-pR", "T47D-pR",
"T47D-pR", "T47D-pR", "T47D-pR"), compound = c("Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A"), dose = c(0.01, 0.078125, 0.15625, 0.3125, 0.625,
1.25, 2.5, 5, 10, 20, 0.01, 0.078125, 0.15625, 0.3125, 0.625,
1.25, 2.5, 5, 10, 20, 0.01, 0.078125, 0.15625, 0.3125, 0.625,
1.25, 2.5, 5, 10, 20, 0.01, 0.078125, 0.15625, 0.3125, 0.625,
1.25, 2.5, 5, 10, 20, 0.01, 0.078125, 0.15625, 0.3125, 0.625,
1.25, 2.5, 5, 10, 20), parental = c("MCF7", "MCF7", "MCF7", "MCF7",
"MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7",
"MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7",
"T47D", "T47D", "T47D", "T47D", "T47D", "T47D", "T47D", "T47D",
"T47D", "T47D", "T47D", "T47D", "T47D", "T47D", "T47D", "T47D",
"T47D", "T47D", "T47D", "T47D", "T47D", "T47D", "T47D", "T47D",
"T47D", "T47D", "T47D", "T47D", "T47D", "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, 10, 78.125, 156.25,
312.5, 625, 1250, 2500, 5000, 10000, 20000, 10, 78.125, 156.25,
312.5, 625, 1250, 2500, 5000, 10000, 20000, 10, 78.125, 156.25,
312.5, 625, 1250, 2500, 5000, 10000, 20000, 10, 78.125, 156.25,
312.5, 625, 1250, 2500, 5000, 10000, 20000)), row.names = c(NA,
-50L), groups = structure(list(cell_line = c("MCF7", "MCF7",
"MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7", "MCF7",
"MCF7-A", "MCF7-A", "MCF7-A", "MCF7-A", "MCF7-A", "MCF7-A", "MCF7-A",
"MCF7-A", "MCF7-A", "MCF7-A", "T47D", "T47D", "T47D", "T47D",
"T47D", "T47D", "T47D", "T47D", "T47D", "T47D", "T47D-A", "T47D-A",
"T47D-A", "T47D-A", "T47D-A", "T47D-A", "T47D-A", "T47D-A", "T47D-A",
"T47D-A", "T47D-pR", "T47D-pR", "T47D-pR", "T47D-pR", "T47D-pR",
"T47D-pR", "T47D-pR", "T47D-pR", "T47D-pR", "T47D-pR"), compound = c("Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A", "Drug A", "Drug A", "Drug A",
"Drug A"), dose = c(0.01, 0.078125, 0.15625, 0.3125, 0.625,
1.25, 2.5, 5, 10, 20, 0.01, 0.078125, 0.15625, 0.3125, 0.625,
1.25, 2.5, 5, 10, 20, 0.01, 0.078125, 0.15625, 0.3125, 0.625,
1.25, 2.5, 5, 10, 20, 0.01, 0.078125, 0.15625, 0.3125, 0.625,
1.25, 2.5, 5, 10, 20, 0.01, 0.078125, 0.15625, 0.3125, 0.625,
1.25, 2.5, 5, 10, 20), .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"))), row.names = c(NA, -50L), class = c("tbl_df",
"tbl", "data.frame"), .drop = TRUE), class = c("grouped_df",
"tbl_df", "tbl", "data.frame"))
似乎 dr4pl 不像 drc 包通过 curveid 参数处理多剂量反应。
我的解决方案是使用 dr4pl::IC() 提取 IC50，然后将它们作为垂直线添加到图中。
这是一个提取 IC50 并将其保存为 data.frame 的函数：
library(data.table)
library(dr4pl)
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))
}
输入数据应该是data.table。使用代码如下：
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))