-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathdrug_sensitivity_scoring.R
166 lines (151 loc) · 6.19 KB
/
drug_sensitivity_scoring.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
# Curve fitting function -----------------------------------------------------
#' Fit logistic curve and calculate IC50/EC50
#'
#' \code{dss} Fir logistic regression to dose response data
#'
#' @param viaobj named vector of viability scores
#' @return returns a dataframe containing IC50/EC50, the slope at the IC50
#' @export
fit_logistic <- function(via) {
viadf <- data.frame(viability = via,
inhibition = 100-via,
dose = as.numeric(names(via))) %>%
filter(!is.na(dose))
fit <- drc::drm(viability ~ dose, data = viadf,
fct = drc::LL.4(fixed = c(NA, NA, NA, NA),
names = c("slope","min","max","EC50")),
control = drc::drmc(errorm = F))
if (class(fit) != "drc") {
fit <- drc::drm(viability ~ dose, data = viadf,
fct = drc::LL.4(fixed = c(NA, NA, 100, NA),
names = c("slope","min","max","EC50")),
control = drc::drmc(errorm = F))
}
return(fit)
}
# DSS function ---------------------------------------------------------------
#' Calculate drug sensitivity score
#'
#' \code{dss} Takes as input: ec50, slope, max_response, min_concentration, max_concentration
#'
#' @param ic50 calculated IC50/EC50
#' @param slope slope
#' @param maxr maximum response
#' @param minc minimum concentration tested
#' @param maxc maximum concentration tested
#' @param th minimum response threshold (default = 10)
#' @param minr minimum response (default = 0)
#' @param log whether to log input values (default = T)
#' @return returns the EC50, slope, max_response, DSS1, DSS2 and DSS3 scores as a dataframe
#' @export
calc_dss <- function(ec50, slope, maxr, minc, maxc, th=10, minr=0, log=T) {
# log inputs if needed
if (log) {
minc <- log10(minc)
maxc <- log10(maxc)
ec50 <- log10(ec50)
}
# check response over threshold
if (ec50 > maxc | maxr < th | slope < 0) {
dss1=0; dss2=0; dss3=0
return(data.frame(EC50 = 10^ec50,
slope = slope,
minc = 10^minc,
maxc = 10^maxc,
maxr = maxr,
DSS1 = dss1,
DSS2 = dss2,
DSS3 = dss3))
}
# x1 = minimum concentration that exceeds minimum response threshold
if (th != 0) {
x1 <- (ec50 - ((log(maxr-th) - log(th-minr)) / (slope*log(10))))
if (x1 < minc) { x1 <- minc }
if (x1 > maxc) { x1 <- maxc }
} else {
x1 <- minc
}
# x2 = minimum concentration that exceeds maximum response threshold
maxth <- maxr*0.9
x2 <- (ec50 - ((log(maxr-maxth) - log(maxth-minr)) / (slope*log(10))))
if (x2 < minc) { x2 <- minc }
if (x2 > maxc) { x2 <- maxc }
# integration
int_y <- (((((maxr-minr) * log(1+10^(slope*(ec50-maxc)))) / (slope*log(10))) + maxr*maxc) -
((((maxr-minr) * log(1+10^(slope*(ec50-x1)))) / (slope*log(10))) + maxr*x1)) - (th*(maxc-x1))
if (int_y < 0) { int_y = 0 }
# total area
total_area = (maxc - minc) * (100-th)
max_area = ((maxc-x2) / (maxc-minc)) * 2
max_area <- ifelse(max_area > 1, 1, max_area)
response_area = ((x2-x1)/(maxc-minc)) * 2
response_area <- ifelse(response_area > 1, 1, response_area)
# dss
dss1 = round(((int_y/total_area)*100), 3)
# dss2 = round(dss1 / log10(maxr), 3)
# dss3 = round(dss1 * (log10(100)/log10(maxr)) * ((maxc-x1)/(maxc-minc)), 3)
dss2 = round(dss1 * max_area, 3)
dss3 = round(dss2 * response_area, 3)
# return
return(data.frame(row.names = NULL,
EC50 = round(10^ec50, 3),
slope = round(slope, 3),
minc = 10^minc,
maxc = 10^maxc,
maxr = round(maxr, 3),
DSS1 = dss1,
DSS2 = dss2,
DSS3 = dss3))
}
# DRC plotting ---------------------------------------------------------------
#' Dose-response curve plotting function
#'
#' \code{dss} Takes as input: drc model generated by drm
#'
#' @param fit drc model
#' @param xlab label for x axis (default = "Concentration")
#' @param ylab label for y axis (default = "Viability %")
#' @param conf whether to plot confidence intervals (default = TRUE)
#' @param points whether to plot observed data points (default = TRUE)
#' @return returns a ggplot object
#' @export
plot_fit <- function(fit, xlab = "Concentration", ylab = "Viability (%)", conf = T, points = T) {
obsdat <- fit$data[,1:2]
colnames(obsdat) <- c("conc","Prediction")
minc <- min(obsdat$conc, na.rm=T)
maxc <- max(obsdat$conc, na.rm=T)
preddat <- expand.grid(conc = exp(seq(log(maxc), log(minc), length=100)))
preddat <- cbind(preddat, predict(fit, newdata=preddat, interval="confidence"))
p <- ggplot() +
geom_line(data = preddat, aes(x = conc, y = Prediction), lwd = 0.5) +
theme_bw(base_size = 14) +
theme(panel.grid = element_blank()) +
xlab(xlab) + ylab(ylab) +
annotation_logticks(sides = "b", alpha = 0.5) +
scale_x_continuous(trans = "log10")
if (conf) p <- p + geom_ribbon(data = preddat, alpha = 0.25,
aes(x = conc, y = Prediction, ymin = Lower, ymax = Upper))
if (points) p <- p + geom_point(data = obsdat,
aes(x = conc, y = Prediction), alpha = 0.7)
return(p)
}
# ---------------------------------------------------------------
#' Get predicted dose-response curve from drc model object
#'
#' \code{dss} Takes as input: drc model generated by drm
#'
#' @param fit drc model
#' @param n number of points to predict across the range
#' @param minc minimum concentration for range - taken from fit if NA
#' @param maxc maximum concentration for range - taken from fit if NA
#' @return returns a named vector of predicted response values over the dose range
#' @export
pred_dr_curve <- function(fit, n = 100, minc = NA, maxc = NA) {
obsdat <- fit$data[,1:2]
colnames(obsdat) <- c("conc","Prediction")
if (is.na(minc)) minc <- min(obsdat$conc, na.rm=T)
if (is.na(maxc)) maxc <- max(obsdat$conc, na.rm=T)
preddat <- expand.grid(conc = exp(seq(log(maxc), log(minc), length=n)))
preddat <- cbind(preddat, predict(fit, newdata=preddat, interval="confidence"))
return(preddat)
}