I was helped on this site in combining the dr4pl (the replacement for '''drc''') 4 parameter logistic dose response model with ggplot. It works really well, but I'm looking to incorporate the the IC50 for each curve into the plots, AND to output a table with the values as well.
Example data:
#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"))
the above will generate a plot with three curves, each with a fitted 4pl model. I need the IC50 to show up in the plot, or in a table elsewhere so I can put it into the plot myself.
Grateful in advance for any help.
Edit June 2021:
new data set to help with debugging the script provided:
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"))
It seems that dr4pl doesn't handle multiple dose responses as the drc package does via the curveid parameter.
My solution would be to extract the IC50s using dr4pl::IC() and then add them as vertical lines to the plot.
Here's a function that extracts and saves IC50 as a 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))
}
The input data should be a data.table. Use the code as follows:
dfIC50 <- multiIC(data = example2,
colDose = "dose",
colResp = "POC",
colID = "curve",
inhib.percent = 50)
To obtain:
> dfIC50
id x
1 C1 29.065593
2 C2 3.269516
3 C3 2.186479
Then add this line to your ggplot:
geom_vline(data = dfIC50,
aes(xintercept = x,
color = id))