R/Plot_synergy.R
In struckma/synergyfinder: Calculate and Visualize Synergy Scores for Drug Combinations
#' f2si2
#'
#' Format numbers with SI units
#'
#' @param numbers vector of numbers
#' @param rounding round numbers
#' @param force_unit NA = auto detect suitable units per number, a unit: convert numbers to that unit
#' @param no_suffix omit the unit prefix in the output
#'
#' @return a vector of converted numbers as characters
#'
#' @export
#'
#' @author Stephan Struckmann \email{stephan.struckmann@expressence.de}
#' @author Roland \url{https://stackoverflow.com/users/1412059/roland)}
#' @author Jonas Stein \email{news@jonasstein.de} \url{https://github.com/jonasstein/sitools}
#'
#' @examples
#' f2si2(c(1000000.1, 10, 100, 1000, 0.001, .000001, .00000001), force_unit = "k", no_suffix = FALSE )
#' f2si2(c(1000000.1, 10, 100, 1000, 0.001, .000001, .00000001), force_unit = "k" )
#' f2si2(c(1000000.1, 10, 100, 1000, 0.001, .000001, .00000001) )
#'
#' @seealso https://stackoverflow.com/a/11341038
f2si2<-function (numbers, rounding=FALSE, force_unit = NA, no_suffix = !is.na(force_unit))
{
sistrings = c()
for (i in 1:length(numbers)) {
number = numbers[[i]]
if (number < 0) {
number <- - number
sgn = "-"
} else
sgn = ""
lut <- c(1e-24, 1e-21, 1e-18, 1e-15, 1e-12, 1e-09, 1e-06,
0.001, 1, 1000, 1e+06, 1e+09, 1e+12, 1e+15, 1e+18, 1e+21,
1e+24)
pre <- c("y", "z", "a", "f", "p", "n", "µ", "m", "", "k",
"M", "G", "T", "P", "E", "Z", "Y")
if (!is.na(force_unit)) {
ix <- which(pre == force_unit)
} else {
if (number == 0)
ix = 9
else
ix <- findInterval(number, lut)
}
if (ix == 0)
ix <- 1
if (lut[ix]!=1) {
if (no_suffix)
unit = ""
else
unit = pre[ix]
if (rounding==T) {
sistring <- paste0(sgn, paste(round(number/lut[ix]), unit, sep = ifelse(unit == "", "", " ")))
}
else {
sistring <- paste0(sgn, paste(number/lut[ix], unit, sep = ifelse(unit == "", "", " ")))
}
} else {
sistring <- paste0(sgn, as.character(number))
}
sistrings = c(sistrings, sistring)
}
return(sistrings)
}
#'Drug interaction landscape
#'
#'A function to visualize the synergy scores for drug combinations as 2D or 3D
#'interaction landscape over the dose-response matrix.
#'
#' @param data a list object generated by function \code{\link{CalculateSynergy}}.
#' @param knitr leave out the dev.off / dev.new calls to make output working in knitr files
#' @param type a parameter to specify the type of the interaction landscape, 2D, 3D
#' or both. By default, 2D interaction landscape is returned.
#' @param save.file a logical parameter to specify if the interaction landscape is saved as a pdf file
#' in the current working directory or returned as an R object. By default, it is FALSE.
#' @param pair.index a parameter to specify which drug combination if there are many drug combinations in the data.
#' By default, it is NULL so that the synergy score visualization of all the drug combinations in the data is returned.
#' @param legend.start a parameter to specify the starting point of the legend. By defualt, it is NULL so the legend
#' starting point is fixed by the data automatically.
#' @param legend.end a parameter to specify the ending point of the legend. By defualt, it is NULL so the legend
#' ending point is fixed by the data automatically.
#' @param x.range a parameter to specify the starting and ending concentration of the drug on x-axis. Use e.g., c(1, 3) to
#' specify that only from 1st to 3rd concentrations of the drug on x-axis are used. By default, it is NULl so all the
#' concentrations are used.
#' @param y.range a parameter to specify the starting and ending concentration of the drug on y-axis. Use e.g., c(1, 3) to
#' specify that only from 1st to 3rd concentrations of the drug on y-axis are used. By default, it is NULl so all the
#' concentrations are used.
#' @param x.unit unit for display
#' @param y.unit unit for display
#' @return a pdf file or the interaction landscapes are only displayed depending on
#' the save.file parameter.
#' @author Liye He \email{liye.he@helsinki.fi}
#' @examples
#' data("mathews_screening_data")
#' data <- ReshapeData(mathews_screening_data)
#' scores <- CalculateSynergy(data)
#' PlotSynergy(scores, "2D")
PlotSynergy <- function(data, knitr = FALSE, type = "2D", save.file = FALSE, pair.index = NULL, legend.start = NULL, legend.end = NULL, x.range = NULL, y.range = NULL,
x.unit = "µ", y.unit = "µ") {
if(!is.list(data)) {
stop("Input data is not a list format!")
}
scores <- data$scores
drug.pairs <- data$drug.pairs
num.pairs <- 1:length(scores)
plots <- list()
if(!is.null(pair.index)) {
num.pairs <- pair.index
}
for (i in num.pairs) {
scores.dose <- t(scores[[i]])
drug.col <- drug.pairs$drug.col[i]
drug.row <- drug.pairs$drug.row[i]
scores.tmp <- scores.dose
if(!is.null(x.range)) {
if(x.range[1] == 1) {
scores.tmp <- scores.tmp[(x.range[1] + 1):x.range[2], ]
} else {
scores.tmp <- scores.tmp[x.range[1]:x.range[2], ]
}
} else {
# delete the first row
scores.tmp <- scores.tmp[-1, ]
}
if(!is.null(y.range)) {
if(y.range[1] == 1) {
scores.tmp <- scores.tmp[, (y.range[1] + 1):y.range[2]]
} else {
scores.tmp <- scores.tmp[, y.range[1]:y.range[2]]
}
} else {
# delete the first column
scores.tmp <- scores.tmp[, -1]
}
#summary.score <- round(mean(scores.dose[-1, -1]), 3)
summary.score <- round(mean(scores.tmp, na.rm = TRUE), 3)
# drug.col: the col-wise drug
# drug.row: the row-wise drug
x.conc <- as.numeric(rownames(scores.dose)) ## col-wise drug concentration
y.conc <- as.numeric(colnames(scores.dose)) ## row-wise drug concentration
# drug.col.conc in index
x <- c(0:(length(x.conc) - 1))
# drug.row.conc in index
y <- c(0:(length(y.conc) - 1))
# smoothing by kriging
tmp <- cbind(expand.grid(x, y),c(as.matrix(scores.dose)))
x.vec <- tmp[, 1]
y.vec <- tmp[, 2]
scores.dose.vec = tmp[, 3]
#kriged = kriging(x.vec,y.vec,delta.dose.vec,lags=3,pixels=50)
pixels.num <- 5 * (length(x.conc) - 1) + 2
if (dim(scores.dose)[1] < 8) {
# kriged <- kriging(x.vec, y.vec, scores.dose.vec,lags = 2, pixels = pixels.num, model = "spherical")
kriged <- kriging(y.vec, x.vec, scores.dose.vec,lags = 2, pixels = pixels.num, model = "spherical")
} else {
# kriged <- kriging(x.vec, y.vec, scores.dose.vec,lags = 3, pixels = pixels.num, model = "spherical")
kriged <- kriging(y.vec, x.vec, scores.dose.vec,lags = 3, pixels = pixels.num, model = "spherical")
}
xseq <- round(kriged[["map"]]$x / kriged$pixel)
yseq <- round(kriged[["map"]]$y / kriged$pixel)
a <- min(xseq):max(xseq)
b <- min(yseq):max(yseq)
z.len <- length(kriged[["map"]]$pred) ## n
na <- length(a)
nb <- length(b)
res1 <- as.double(rep(0, na * nb)) ## c
res2 <- as.integer(rep(0,na * nb)) ## d
for(idx1 in 1:na) {
for(idx2 in 1:nb) {
for(idx3 in 1:z.len) {
if(xseq[idx3] == a[idx1] && yseq[idx3] == b[idx2]) {
res1[idx2+(idx1-1)*nb] <- kriged[["map"]]$pred[idx3]
res2[idx2+(idx1-1)*nb] <- 1
break
}
}
}
}
res1[res2 == 0] <- NA
c <- matrix(res1, na, nb, byrow = TRUE)
plot.title <- paste(data$method,"synergy score:", summary.score, sep = " ")
lut <- c(1e-24, 1e-21, 1e-18, 1e-15, 1e-12, 1e-09, 1e-06,
0.001, 1, 1000, 1e+06, 1e+09, 1e+12, 1e+15, 1e+18, 1e+21,
1e+24)
pre <- c("y", "z", "a", "f", "p", "n", "µ", "m", "", "k",
"M", "G", "T", "P", "E", "Z", "Y")
conc.unit <- drug.pairs$concUnit[i] ## concentration unit
conc.unit.prefix <- substr(conc.unit, 1, 1)
unit.text.x <- paste("(", x.unit, "M)", sep = "")
unit.text.y <- paste("(", y.unit, "M)", sep = "")
file.name <- paste(drug.row, drug.col, "synergy", drug.pairs$blockIDs[i], data$method, "pdf", sep = ".")
drug.row <- paste(drug.row, unit.text.x, sep = " ")
drug.col <- paste(drug.col, unit.text.y, sep = " ")
max.dose <- max(abs(max(scores.dose)), abs(min(scores.dose)))
color.range <- round(max.dose + 5, -1)
if(is.null(legend.start)){
start.point <- -color.range
} else {
start.point <- legend.start
}
if(is.null(legend.end)) {
end.point <- color.range
} else {
end.point <- legend.end
}
if (type == "3D") {
plots.3d <- c
if(!is.null(x.range)) {
row.start <- (x.range[1] - 1) * 5
row.end <- (x.range[2] - 1) * 5
plots.3d <- plots.3d[row.start:row.end, ]
}
if(!is.null(y.range)) {
col.start <- (y.range[1] - 1) * 5
col.end <- (y.range[2] - 1) * 5
plots.3d <- plots.3d[, col.start:col.end]
}
x.conc1 <- x.conc
y.conc1 <- y.conc
if(!is.null(x.range)) {
x.conc1 <- x.conc1[x.range[1]:x.range[2]]
}
if(!is.null(y.range)) {
y.conc1 <- y.conc1[y.range[1]:y.range[2]]
}
# conc1 -> column.values -- ncol(plots.3d)
# conc2 -> row.values -- nrow(plots.3d)
c1 = data.frame(point=seq(from = 1, to = ncol(plots.3d), length.out = length(x.conc1)), conc=x.conc1)
c2 = data.frame(point=seq(from = 1, to = nrow(plots.3d), length.out = length(y.conc1)), conc=y.conc1)
column.values = spline(c1$point, c1$conc, xout = (1:ncol(plots.3d)))
row.values = spline(c2$point, c2$conc, xout = (1:nrow(plots.3d)))
x.at = spline(c1$conc, c1$point, xout = x.conc1)$y
y.at = spline(c2$conc, c2$point, xout = y.conc1)$y
x_axtcks = axisTicks(c(min(y.conc1), max(y.conc1)), FALSE)
y_axtcks = axisTicks(c(min(x.conc1), max(x.conc1)), FALSE)
#axisTicks(c(min(b), max(b)), FALSE, nint=dev.size(units = "px")[[2]]/30)
fig <- wireframe(plots.3d,scales = list(arrows = FALSE,distance=c(0.8,0.8,0.8),col=1,cex=0.8,z = list(tick.number=7)
,
x=list(at = x_axtcks, labels = f2si2(x_axtcks*lut[pre == conc.unit.prefix], force_unit = x.unit)),
y=list(at = y_axtcks, labels = f2si2(y_axtcks*lut[pre == conc.unit.prefix], force_unit = y.unit))
),
drape = TRUE, colorkey = list(space="top",width=0.5),
screen = list(z = 30, x = -55),
zlab=list(expression("Synergy score"),rot=90,cex=1,axis.key.padding = 0),xlab=list(as.character(drug.row),cex=1, rot=20),ylab=list(as.character(drug.col),cex=1,rot=-50),
zlim=c(start.point, end.point),
col.regions=colorRampPalette(c("green","white","red"))(100),
main = plot.title,
at=do.breaks(c(start.point, end.point),100),
par.settings = list(axis.line=list(col="transparent")),
zoom = 1, aspect = 1,
row.values = row.values$y,
column.values = column.values$y
#par.settings=list(layout.widths=list(left.padding=0,right.padding=0), layout.heights=list(top.padding=0, bottom.padding=0)) # margin
)
print(fig)
fig <- recordPlot()
} else if (type == "2D") {
dev.new(noRStudioGD = TRUE)
layout(matrix(c(1, 2), nrow = 2L, ncol = 1L), heights = c(0.1, 1))
par(mar = c(0, 6.1, 2.1, 4.1))
suppressWarnings(par(mgp = c(3, -0.4, 0)))
levels <- seq(start.point, end.point, by = 2)
col <- colorpanel(end.point, low = "green", mid = "white", high = "red")
plot.new()
plot.window(ylim = c(0, 1), xlim = range(levels), xaxs = "i", yaxs = "i")
rect(levels[-length(levels)],0, levels[-1L],0.3, col = col, border = NA)
axis(3,tick = FALSE, at = do.breaks(c(start.point, end.point), end.point/10))
title(plot.title)
par(mar = c(5.1,4.1,2.1,2.1))
suppressWarnings(par(mgp = c(2,1,0)))
plot.new()
if(!is.null(x.range)) {
row.start <- (x.range[1] - 1) * 5
row.end <- (x.range[2] - 1) * 5
c <- c[row.start:row.end, ]
}
if(!is.null(y.range)) {
col.start <- (y.range[1] - 1) * 5
col.end <- (y.range[2] - 1) * 5
c <- c[, col.start:col.end]
}
x.2D <- (1:dim(c)[1] - 1) / (dim(c)[1] - 1)
y.2D <- (1:dim(c)[2] - 1) / (dim(c)[2] - 1)
plot.window(asp = NA, xlim = range(x.2D), ylim = range(y.2D), "", xaxs = "i", yaxs = "i")
.filled.contour(x.2D, y.2D, z = c, levels, col = col)
#grid(dim(c)[1] - 1, dim(c)[2] - 1, col = "gray", lty = "solid", lwd = 0.3)
box()
mtext(drug.col, 1, cex = 1, padj = 3)
mtext(drug.row, 2, cex = 1, las = 3, padj = -3)
if(!is.null(x.range)) {
x.conc <- x.conc[x.range[1]:x.range[2]]
}
if(!is.null(y.range)) {
y.conc <- y.conc[y.range[1]:y.range[2]]
}
axis(side = 1, at = seq(0, 1, by = 1/(length(x.conc) - 1)), labels = signif(x.conc, 1))
axis(side = 2, at = seq(0, 1, by = 1/(length(y.conc) - 1)), labels = signif(y.conc, 1))
fig <- recordPlot()
dev.off()
} else {
plots.3d <- c
if(!is.null(x.range)) {
row.start <- (x.range[1] - 1) * 5
row.end <- (x.range[2] - 1) * 5
plots.3d <- plots.3d[row.start:row.end, ]
}
if(!is.null(y.range)) {
col.start <- (y.range[1] - 1) * 5
col.end <- (y.range[2] - 1) * 5
plots.3d <- plots.3d[, col.start:col.end]
}
x.conc1 <- x.conc
y.conc1 <- y.conc
if(!is.null(x.range)) {
x.conc1 <- x.conc1[x.range[1]:x.range[2]]
}
if(!is.null(y.range)) {
y.conc1 <- y.conc1[y.range[1]:y.range[2]]
}
# conc1 -> column.values -- ncol(plots.3d)
# conc2 -> row.values -- nrow(plots.3d)
c1 = data.frame(point=seq(from = 1, to = ncol(plots.3d), length.out = length(x.conc1)), conc=x.conc1)
c2 = data.frame(point=seq(from = 1, to = nrow(plots.3d), length.out = length(y.conc1)), conc=y.conc1)
column.values = spline(c1$point, c1$conc, xout = (1:ncol(plots.3d)))
row.values = spline(c2$point, c2$conc, xout = (1:nrow(plots.3d)))
x.at = spline(c1$conc, c1$point, xout = x.conc1)$y
y.at = spline(c2$conc, c2$point, xout = y.conc1)$y
x_axtcks = axisTicks(c(min(y.conc1), max(y.conc1)), FALSE)
y_axtcks = axisTicks(c(min(x.conc1), max(x.conc1)), FALSE)
syn.3d.plot <- wireframe(plots.3d,
scales = list(arrows = FALSE,distance=c(0.8,0.8,0.8),col=1,cex=0.8,z = list(tick.number=7)
,
x=list(at = x_axtcks, labels = f2si2(x_axtcks*lut[pre == conc.unit.prefix], force_unit = x.unit)),
y=list(at = y_axtcks, labels = f2si2(y_axtcks*lut[pre == conc.unit.prefix], force_unit = y.unit))
),
drape = TRUE, colorkey = list(space="top",width=0.5),
screen = list(z = 30, x = -55),
zlab=list(expression("Synergy score"),rot=90,cex=1,axis.key.padding = 0),xlab=list(as.character(drug.row),cex=1, rot=20),ylab=list(as.character(drug.col),cex=1,rot=-50),
zlim=c(start.point, end.point),
col.regions=colorRampPalette(c("green","white","red"))(100),
main = plot.title,
at=do.breaks(c(start.point, end.point),100),
par.settings = list(axis.line=list(col="transparent")),
zoom = 1, aspect = 1,
row.values = row.values$y,
column.values = column.values$y
#par.settings=list(layout.widths=list(left.padding=0,right.padding=0), layout.heights=list(top.padding=0, bottom.padding=0)) # margin
)
layout(matrix(c(1, 2, 3, 3), nrow = 2L, ncol = 2L), heights = c(0.1, 1))
par(mar = c(0, 6.1, 2.1, 4.1))
suppressWarnings(par(mgp = c(3, -1, 0)))
levels <- seq(start.point, end.point, by = 2)
col <- colorpanel(end.point, low = "green", mid = "white", high = "red")
plot.new()
plot.window(ylim = c(0, 1), xlim = range(levels), xaxs = "i", yaxs = "i")
rect(levels[-length(levels)],0, levels[-1L],0.3, col = col, border = NA)
axis(3,tick = FALSE, at = do.breaks(c(start.point, end.point), end.point/10))
title(plot.title)
par(mar = c(5.1,4.1,2.1,2.1))
suppressWarnings(par(mgp = c(2,1,0)))
plot.new()
if(!is.null(x.range)) {
row.start <- (x.range[1] - 1) * 5
row.end <- (x.range[2] - 1) * 5
c <- c[row.start:row.end, ]
}
if(!is.null(y.range)) {
col.start <- (y.range[1] - 1) * 5
col.end <- (y.range[2] - 1) * 5
c <- c[, col.start:col.end]
}
x.2D <- (1:dim(c)[1] - 1) / (dim(c)[1] - 1)
y.2D <- (1:dim(c)[2] - 1) / (dim(c)[2] - 1)
plot.window(asp = NA, xlim = range(x.2D), ylim = range(y.2D), "", xaxs = "i", yaxs = "i")
.filled.contour(x.2D, y.2D, z = c, levels, col = col)
#grid(dim(c)[1] - 1, dim(c)[2] - 1, col = "gray", lty = "solid", lwd = 0.3)
box()
mtext(drug.col, 1, cex = 1, padj = 3)
mtext(drug.row, 2, cex = 1, las = 3, padj = -3)
if(!is.null(x.range)) {
x.conc <- x.conc[x.range[1]:x.range[2]]
}
if(!is.null(y.range)) {
y.conc <- y.conc[y.range[1]:y.range[2]]
}
axis(side = 1, at = seq(0, 1, by = 1/(length(x.conc) - 1)), labels = signif(x.conc, 1))
axis(side = 2, at = seq(0, 1, by = 1/(length(y.conc) - 1)), labels = signif(y.conc, 1))
print(syn.3d.plot, position = c(0.5,0, 1, 1), newpage = FALSE)
fig <- recordPlot()
}
plots[[i]] <- fig
if(save.file) {
if (type %in% c("2D", "3D")) {
pdf(file.name, width = 10, height = 10)
} else {
pdf(file.name, width = 12, height = 6)
}
print(fig)
dev.off()
}
}
if(!save.file) {
for(i in num.pairs) {
dev.new(noRStudioGD = TRUE)
replayPlot(plots[[i]])
}
}
}
PlotEffect <- function(data, knitr = FALSE, type = "2D", save.file = FALSE,
effect.label = "Effect",
pair.index = NULL, legend.start = NULL, legend.end = NULL, x.range = NULL, y.range = NULL,
x.unit = "µ", y.unit = "µ") {
if(!is.list(data)) {
stop("Input data is not a list format!")
}
scores <- data$dose.response.mats
drug.pairs <- data$drug.pairs
num.pairs <- 1:length(scores)
plots <- list()
if(!is.null(pair.index)) {
num.pairs <- pair.index
}
for (i in num.pairs) {
scores.dose <- t(scores[[i]])
drug.col <- drug.pairs$drug.col[i]
drug.row <- drug.pairs$drug.row[i]
scores.tmp <- scores.dose
if(!is.null(x.range)) {
if(x.range[1] == 1) {
scores.tmp <- scores.tmp[(x.range[1] + 1):x.range[2], ]
} else {
scores.tmp <- scores.tmp[x.range[1]:x.range[2], ]
}
} else {
# delete the first row
scores.tmp <- scores.tmp[-1, ]
}
if(!is.null(y.range)) {
if(y.range[1] == 1) {
scores.tmp <- scores.tmp[, (y.range[1] + 1):y.range[2]]
} else {
scores.tmp <- scores.tmp[, y.range[1]:y.range[2]]
}
} else {
# delete the first column
scores.tmp <- scores.tmp[, -1]
}
#summary.score <- round(mean(scores.dose[-1, -1]), 3)
summary.score <- round(mean(scores.tmp, na.rm = TRUE), 3)
# drug.col: the col-wise drug
# drug.row: the row-wise drug
x.conc <- as.numeric(rownames(scores.dose)) ## col-wise drug concentration
y.conc <- as.numeric(colnames(scores.dose)) ## row-wise drug concentration
# drug.col.conc in index
x <- c(0:(length(x.conc) - 1))
# drug.row.conc in index
y <- c(0:(length(y.conc) - 1))
# smoothing by kriging
tmp <- cbind(expand.grid(x, y),c(as.matrix(scores.dose)))
x.vec <- tmp[, 1]
y.vec <- tmp[, 2]
scores.dose.vec = tmp[, 3]
#kriged = kriging(x.vec,y.vec,delta.dose.vec,lags=3,pixels=50)
pixels.num <- 5 * (length(x.conc) - 1) + 2
if (dim(scores.dose)[1] < 8) {
kriged <- kriging(y.vec, x.vec, scores.dose.vec,lags = 2, pixels = pixels.num, model = "spherical")
# kriged <- kriging(x.vec, y.vec, scores.dose.vec,lags = 2, pixels = pixels.num, model = "spherical")
} else {
kriged <- kriging(y.vec, x.vec, scores.dose.vec,lags = 3, pixels = pixels.num, model = "spherical")
# kriged <- kriging(x.vec, y.vec, scores.dose.vec,lags = 3, pixels = pixels.num, model = "spherical")
}
xseq <- round(kriged[["map"]]$x / kriged$pixel)
yseq <- round(kriged[["map"]]$y / kriged$pixel)
a <- min(xseq):max(xseq)
b <- min(yseq):max(yseq)
z.len <- length(kriged[["map"]]$pred) ## n
na <- length(a)
nb <- length(b)
res1 <- as.double(rep(0, na * nb)) ## c
res2 <- as.integer(rep(0,na * nb)) ## d
for(idx1 in 1:na) {
for(idx2 in 1:nb) {
for(idx3 in 1:z.len) {
if(xseq[idx3] == a[idx1] && yseq[idx3] == b[idx2]) {
res1[idx2+(idx1-1)*nb] <- kriged[["map"]]$pred[idx3]
res2[idx2+(idx1-1)*nb] <- 1
break
}
}
}
}
res1[res2 == 0] <- NA
c <- matrix(res1, na, nb, byrow = TRUE)
# plot.title <- paste(data$method,"synergy score:", summary.score, "(Effect)", sep = " ")
plot.title <- paste("Dose-response plane", sep = " ")
lut <- c(1e-24, 1e-21, 1e-18, 1e-15, 1e-12, 1e-09, 1e-06,
0.001, 1, 1000, 1e+06, 1e+09, 1e+12, 1e+15, 1e+18, 1e+21,
1e+24)
pre <- c("y", "z", "a", "f", "p", "n", "µ", "m", "", "k",
"M", "G", "T", "P", "E", "Z", "Y")
conc.unit <- drug.pairs$concUnit[i] ## concentration unit
conc.unit.prefix <- substr(conc.unit, 1, 1)
unit.text.x <- paste("(", x.unit, "M)", sep = "")
unit.text.y <- paste("(", y.unit, "M)", sep = "")
file.name <- paste(drug.row, drug.col, "effect", drug.pairs$blockIDs[i], data$method, "pdf", sep = ".")
drug.row <- paste(drug.row, unit.text.x, sep = " ")
drug.col <- paste(drug.col, unit.text.y, sep = " ")
max.dose <- max(abs(max(scores.dose)), abs(min(scores.dose)))
color.range <- round(max.dose + 5, -1)
if(is.null(legend.start)){
start.point <- 0
} else {
start.point <- legend.start
}
if(is.null(legend.end)) {
end.point <- color.range
} else {
end.point <- legend.end
}
if (type == "3D") {
plots.3d <- c
if(!is.null(x.range)) {
row.start <- (x.range[1] - 1) * 5
row.end <- (x.range[2] - 1) * 5
plots.3d <- plots.3d[row.start:row.end, ]
}
if(!is.null(y.range)) {
col.start <- (y.range[1] - 1) * 5
col.end <- (y.range[2] - 1) * 5
plots.3d <- plots.3d[, col.start:col.end]
}
x.conc1 <- x.conc
y.conc1 <- y.conc
if(!is.null(x.range)) {
x.conc1 <- x.conc1[x.range[1]:x.range[2]]
}
if(!is.null(y.range)) {
y.conc1 <- y.conc1[y.range[1]:y.range[2]]
}
# conc1 -> column.values -- ncol(plots.3d)
# conc2 -> row.values -- nrow(plots.3d)
c1 = data.frame(point=seq(from = 1, to = ncol(plots.3d), length.out = length(x.conc1)), conc=x.conc1)
c2 = data.frame(point=seq(from = 1, to = nrow(plots.3d), length.out = length(y.conc1)), conc=y.conc1)
column.values = spline(c1$point, c1$conc, xout = (1:ncol(plots.3d)))
row.values = spline(c2$point, c2$conc, xout = (1:nrow(plots.3d)))
x.at = spline(c1$conc, c1$point, xout = x.conc1)$y
y.at = spline(c2$conc, c2$point, xout = y.conc1)$y
x_axtcks = axisTicks(c(min(y.conc1), max(y.conc1)), FALSE)
y_axtcks = axisTicks(c(min(x.conc1), max(x.conc1)), FALSE)
fig <- wireframe(plots.3d,scales = list(arrows = FALSE,distance=c(0.8,0.8,0.8),col=1,cex=0.8,z = list(tick.number=7)
,
x=list(at = x_axtcks, labels = f2si2(x_axtcks*lut[pre == conc.unit.prefix], force_unit = x.unit)),
y=list(at = y_axtcks, labels = f2si2(y_axtcks*lut[pre == conc.unit.prefix], force_unit = y.unit))
),
drape = TRUE, colorkey = list(space="top",width=0.5),
screen = list(z = 30, x = -55),
zlab=list(effect.label,rot=90,cex=1,axis.key.padding = 0),xlab=list(as.character(drug.row),cex=1, rot=20),ylab=list(as.character(drug.col),cex=1,rot=-50),
zlim=c(start.point, end.point),
col.regions=colorRampPalette(c("green","white","red"))(100),
main = plot.title,
at=do.breaks(c(start.point, end.point),100),
par.settings = list(axis.line=list(col="transparent")),
zoom = 1, aspect = 1,
row.values = row.values$y,
column.values = column.values$y
#par.settings=list(layout.widths=list(left.padding=0,right.padding=0), layout.heights=list(top.padding=0, bottom.padding=0)) # margin
)
print(fig)
fig <- recordPlot()
} else if (type == "2D") {
dev.new(noRStudioGD = TRUE)
layout(matrix(c(1, 2), nrow = 2L, ncol = 1L), heights = c(0.1, 1))
par(mar = c(0, 6.1, 2.1, 4.1))
suppressWarnings(par(mgp = c(3, -0.4, 0)))
levels <- seq(start.point, end.point, by = 2)
col <- colorpanel(end.point, low = "green", mid = "white", high = "red")
plot.new()
plot.window(ylim = c(0, 1), xlim = range(levels), xaxs = "i", yaxs = "i")
rect(levels[-length(levels)],0, levels[-1L],0.3, col = col, border = NA)
axis(3,tick = FALSE, at = do.breaks(c(start.point, end.point), end.point/10))
title(plot.title)
par(mar = c(5.1,4.1,2.1,2.1))
suppressWarnings(par(mgp = c(2,1,0)))
plot.new()
if(!is.null(x.range)) {
row.start <- (x.range[1] - 1) * 5
row.end <- (x.range[2] - 1) * 5
c <- c[row.start:row.end, ]
}
if(!is.null(y.range)) {
col.start <- (y.range[1] - 1) * 5
col.end <- (y.range[2] - 1) * 5
c <- c[, col.start:col.end]
}
x.2D <- (1:dim(c)[1] - 1) / (dim(c)[1] - 1)
y.2D <- (1:dim(c)[2] - 1) / (dim(c)[2] - 1)
plot.window(asp = NA, xlim = range(x.2D), ylim = range(y.2D), "", xaxs = "i", yaxs = "i")
.filled.contour(x.2D, y.2D, z = c, levels, col = col)
#grid(dim(c)[1] - 1, dim(c)[2] - 1, col = "gray", lty = "solid", lwd = 0.3)
box()
mtext(drug.col, 1, cex = 1, padj = 3)
mtext(drug.row, 2, cex = 1, las = 3, padj = -3)
if(!is.null(x.range)) {
x.conc <- x.conc[x.range[1]:x.range[2]]
}
if(!is.null(y.range)) {
y.conc <- y.conc[y.range[1]:y.range[2]]
}
axis(side = 1, at = seq(0, 1, by = 1/(length(x.conc) - 1)), labels = signif(x.conc, 1))
axis(side = 2, at = seq(0, 1, by = 1/(length(y.conc) - 1)), labels = signif(y.conc, 1))
fig <- recordPlot()
dev.off()
} else {
plots.3d <- c
if(!is.null(x.range)) {
row.start <- (x.range[1] - 1) * 5
row.end <- (x.range[2] - 1) * 5
plots.3d <- plots.3d[row.start:row.end, ]
}
if(!is.null(y.range)) {
col.start <- (y.range[1] - 1) * 5
col.end <- (y.range[2] - 1) * 5
plots.3d <- plots.3d[, col.start:col.end]
}
x.conc1 <- x.conc
y.conc1 <- y.conc
if(!is.null(x.range)) {
x.conc1 <- x.conc1[x.range[1]:x.range[2]]
}
if(!is.null(y.range)) {
y.conc1 <- y.conc1[y.range[1]:y.range[2]]
}
# conc1 -> column.values -- ncol(plots.3d)
# conc2 -> row.values -- nrow(plots.3d)
c1 = data.frame(point=seq(from = 1, to = ncol(plots.3d), length.out = length(x.conc1)), conc=x.conc1)
c2 = data.frame(point=seq(from = 1, to = nrow(plots.3d), length.out = length(y.conc1)), conc=y.conc1)
column.values = spline(c1$point, c1$conc, xout = (1:ncol(plots.3d)))
row.values = spline(c2$point, c2$conc, xout = (1:nrow(plots.3d)))
x.at = spline(c1$conc, c1$point, xout = x.conc1)$y
y.at = spline(c2$conc, c2$point, xout = y.conc1)$y
x_axtcks = axisTicks(c(min(y.conc1), max(y.conc1)), FALSE)
y_axtcks = axisTicks(c(min(x.conc1), max(x.conc1)), FALSE)
syn.3d.plot <- wireframe(plots.3d,
scales = list(arrows = FALSE,distance=c(0.8,0.8,0.8),col=1,cex=0.8,z = list(tick.number=7)
,
x=list(at = x_axtcks, labels = f2si2(x_axtcks*lut[pre == conc.unit.prefix], force_unit = x.unit)),
y=list(at = y_axtcks, labels = f2si2(y_axtcks*lut[pre == conc.unit.prefix], force_unit = y.unit))
),
drape = TRUE, colorkey = list(space="top",width=0.5),
screen = list(z = 30, x = -55),
zlab=list(expression("Synergy score"),rot=90,cex=1,axis.key.padding = 0),xlab=list(as.character(drug.row),cex=1, rot=20),ylab=list(as.character(drug.col),cex=1,rot=-50),
zlim=c(start.point, end.point),
col.regions=colorRampPalette(c("green","white","red"))(100),
main = plot.title,
at=do.breaks(c(start.point, end.point),100),
par.settings = list(axis.line=list(col="transparent")),
zoom = 1, aspect = 1,
row.values = row.values$y,
column.values = column.values$y
#par.settings=list(layout.widths=list(left.padding=0,right.padding=0), layout.heights=list(top.padding=0, bottom.padding=0)) # margin
)
layout(matrix(c(1, 2, 3, 3), nrow = 2L, ncol = 2L), heights = c(0.1, 1))
par(mar = c(0, 6.1, 2.1, 4.1))
suppressWarnings(par(mgp = c(3, -1, 0)))
levels <- seq(start.point, end.point, by = 2)
col <- colorpanel(end.point, low = "green", mid = "white", high = "red")
plot.new()
plot.window(ylim = c(0, 1), xlim = range(levels), xaxs = "i", yaxs = "i")
rect(levels[-length(levels)],0, levels[-1L],0.3, col = col, border = NA)
axis(3,tick = FALSE, at = do.breaks(c(start.point, end.point), end.point/10))
title(plot.title)
par(mar = c(5.1,4.1,2.1,2.1))
suppressWarnings(par(mgp = c(2,1,0)))
plot.new()
if(!is.null(x.range)) {
row.start <- (x.range[1] - 1) * 5
row.end <- (x.range[2] - 1) * 5
c <- c[row.start:row.end, ]
}
if(!is.null(y.range)) {
col.start <- (y.range[1] - 1) * 5
col.end <- (y.range[2] - 1) * 5
c <- c[, col.start:col.end]
}
x.2D <- (1:dim(c)[1] - 1) / (dim(c)[1] - 1)
y.2D <- (1:dim(c)[2] - 1) / (dim(c)[2] - 1)
plot.window(asp = NA, xlim = range(x.2D), ylim = range(y.2D), "", xaxs = "i", yaxs = "i")
.filled.contour(x.2D, y.2D, z = c, levels, col = col)
#grid(dim(c)[1] - 1, dim(c)[2] - 1, col = "gray", lty = "solid", lwd = 0.3)
box()
mtext(drug.col, 1, cex = 1, padj = 3)
mtext(drug.row, 2, cex = 1, las = 3, padj = -3)
if(!is.null(x.range)) {
x.conc <- x.conc[x.range[1]:x.range[2]]
}
if(!is.null(y.range)) {
y.conc <- y.conc[y.range[1]:y.range[2]]
}
axis(side = 1, at = seq(0, 1, by = 1/(length(x.conc) - 1)), labels = signif(x.conc, 1))
axis(side = 2, at = seq(0, 1, by = 1/(length(y.conc) - 1)), labels = signif(y.conc, 1))
print(syn.3d.plot, position = c(0.5,0, 1, 1), newpage = FALSE)
fig <- recordPlot()
}
plots[[i]] <- fig
if(save.file) {
if (type %in% c("2D", "3D")) {
pdf(file.name, width = 10, height = 10)
} else {
pdf(file.name, width = 12, height = 6)
}
print(fig)
dev.off()
}
}
if(!save.file) {
for(i in num.pairs) {
dev.new(noRStudioGD = TRUE)
replayPlot(plots[[i]])
}
}
}
struckma/synergyfinder documentation built on May 17, 2019, 11:18 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#' f2si2
#'
#' Format numbers with SI units
#'
#' @param numbers vector of numbers
#' @param rounding round numbers
#' @param force_unit NA = auto detect suitable units per number, a unit: convert numbers to that unit
#' @param no_suffix omit the unit prefix in the output
#'
#' @return a vector of converted numbers as characters
#'
#' @export
#'
#' @author Stephan Struckmann \email{stephan.struckmann@expressence.de}
#' @author Roland \url{https://stackoverflow.com/users/1412059/roland)}
#' @author Jonas Stein \email{news@jonasstein.de} \url{https://github.com/jonasstein/sitools}
#'
#' @examples
#' f2si2(c(1000000.1, 10, 100, 1000, 0.001, .000001, .00000001), force_unit = "k", no_suffix = FALSE )
#' f2si2(c(1000000.1, 10, 100, 1000, 0.001, .000001, .00000001), force_unit = "k" )
#' f2si2(c(1000000.1, 10, 100, 1000, 0.001, .000001, .00000001) )
#'
#' @seealso https://stackoverflow.com/a/11341038
f2si2<-function (numbers, rounding=FALSE, force_unit = NA, no_suffix = !is.na(force_unit))
{
sistrings = c()
for (i in 1:length(numbers)) {
number = numbers[[i]]
if (number < 0) {
number <- - number
sgn = "-"
} else
sgn = ""
lut <- c(1e-24, 1e-21, 1e-18, 1e-15, 1e-12, 1e-09, 1e-06,
0.001, 1, 1000, 1e+06, 1e+09, 1e+12, 1e+15, 1e+18, 1e+21,
1e+24)
pre <- c("y", "z", "a", "f", "p", "n", "µ", "m", "", "k",
"M", "G", "T", "P", "E", "Z", "Y")
if (!is.na(force_unit)) {
ix <- which(pre == force_unit)
} else {
if (number == 0)
ix = 9
else
ix <- findInterval(number, lut)
}
if (ix == 0)
ix <- 1
if (lut[ix]!=1) {
if (no_suffix)
unit = ""
else
unit = pre[ix]
if (rounding==T) {
sistring <- paste0(sgn, paste(round(number/lut[ix]), unit, sep = ifelse(unit == "", "", " ")))
}
else {
sistring <- paste0(sgn, paste(number/lut[ix], unit, sep = ifelse(unit == "", "", " ")))
}
} else {
sistring <- paste0(sgn, as.character(number))
}
sistrings = c(sistrings, sistring)
}
return(sistrings)
}
#'Drug interaction landscape
#'
#'A function to visualize the synergy scores for drug combinations as 2D or 3D
#'interaction landscape over the dose-response matrix.
#'
#' @param data a list object generated by function \code{\link{CalculateSynergy}}.
#' @param knitr leave out the dev.off / dev.new calls to make output working in knitr files
#' @param type a parameter to specify the type of the interaction landscape, 2D, 3D
#' or both. By default, 2D interaction landscape is returned.
#' @param save.file a logical parameter to specify if the interaction landscape is saved as a pdf file
#' in the current working directory or returned as an R object. By default, it is FALSE.
#' @param pair.index a parameter to specify which drug combination if there are many drug combinations in the data.
#' By default, it is NULL so that the synergy score visualization of all the drug combinations in the data is returned.
#' @param legend.start a parameter to specify the starting point of the legend. By defualt, it is NULL so the legend
#' starting point is fixed by the data automatically.
#' @param legend.end a parameter to specify the ending point of the legend. By defualt, it is NULL so the legend
#' ending point is fixed by the data automatically.
#' @param x.range a parameter to specify the starting and ending concentration of the drug on x-axis. Use e.g., c(1, 3) to
#' specify that only from 1st to 3rd concentrations of the drug on x-axis are used. By default, it is NULl so all the
#' concentrations are used.
#' @param y.range a parameter to specify the starting and ending concentration of the drug on y-axis. Use e.g., c(1, 3) to
#' specify that only from 1st to 3rd concentrations of the drug on y-axis are used. By default, it is NULl so all the
#' concentrations are used.
#' @param x.unit unit for display
#' @param y.unit unit for display
#' @return a pdf file or the interaction landscapes are only displayed depending on
#' the save.file parameter.
#' @author Liye He \email{liye.he@helsinki.fi}
#' @examples
#' data("mathews_screening_data")
#' data <- ReshapeData(mathews_screening_data)
#' scores <- CalculateSynergy(data)
#' PlotSynergy(scores, "2D")
PlotSynergy <- function(data, knitr = FALSE, type = "2D", save.file = FALSE, pair.index = NULL, legend.start = NULL, legend.end = NULL, x.range = NULL, y.range = NULL,
x.unit = "µ", y.unit = "µ") {
if(!is.list(data)) {
stop("Input data is not a list format!")
}
scores <- data$scores
drug.pairs <- data$drug.pairs
num.pairs <- 1:length(scores)
plots <- list()
if(!is.null(pair.index)) {
num.pairs <- pair.index
}
for (i in num.pairs) {
scores.dose <- t(scores[[i]])
drug.col <- drug.pairs$drug.col[i]
drug.row <- drug.pairs$drug.row[i]
scores.tmp <- scores.dose
if(!is.null(x.range)) {
if(x.range[1] == 1) {
scores.tmp <- scores.tmp[(x.range[1] + 1):x.range[2], ]
} else {
scores.tmp <- scores.tmp[x.range[1]:x.range[2], ]
}
} else {
# delete the first row
scores.tmp <- scores.tmp[-1, ]
}
if(!is.null(y.range)) {
if(y.range[1] == 1) {
scores.tmp <- scores.tmp[, (y.range[1] + 1):y.range[2]]
} else {
scores.tmp <- scores.tmp[, y.range[1]:y.range[2]]
}
} else {
# delete the first column
scores.tmp <- scores.tmp[, -1]
}
#summary.score <- round(mean(scores.dose[-1, -1]), 3)
summary.score <- round(mean(scores.tmp, na.rm = TRUE), 3)
# drug.col: the col-wise drug
# drug.row: the row-wise drug
x.conc <- as.numeric(rownames(scores.dose)) ## col-wise drug concentration
y.conc <- as.numeric(colnames(scores.dose)) ## row-wise drug concentration
# drug.col.conc in index
x <- c(0:(length(x.conc) - 1))
# drug.row.conc in index
y <- c(0:(length(y.conc) - 1))
# smoothing by kriging
tmp <- cbind(expand.grid(x, y),c(as.matrix(scores.dose)))
x.vec <- tmp[, 1]
y.vec <- tmp[, 2]
scores.dose.vec = tmp[, 3]
#kriged = kriging(x.vec,y.vec,delta.dose.vec,lags=3,pixels=50)
pixels.num <- 5 * (length(x.conc) - 1) + 2
if (dim(scores.dose)[1] < 8) {
# kriged <- kriging(x.vec, y.vec, scores.dose.vec,lags = 2, pixels = pixels.num, model = "spherical")
kriged <- kriging(y.vec, x.vec, scores.dose.vec,lags = 2, pixels = pixels.num, model = "spherical")
} else {
# kriged <- kriging(x.vec, y.vec, scores.dose.vec,lags = 3, pixels = pixels.num, model = "spherical")
kriged <- kriging(y.vec, x.vec, scores.dose.vec,lags = 3, pixels = pixels.num, model = "spherical")
}
xseq <- round(kriged[["map"]]$x / kriged$pixel)
yseq <- round(kriged[["map"]]$y / kriged$pixel)
a <- min(xseq):max(xseq)
b <- min(yseq):max(yseq)
z.len <- length(kriged[["map"]]$pred) ## n
na <- length(a)
nb <- length(b)
res1 <- as.double(rep(0, na * nb)) ## c
res2 <- as.integer(rep(0,na * nb)) ## d
for(idx1 in 1:na) {
for(idx2 in 1:nb) {
for(idx3 in 1:z.len) {
if(xseq[idx3] == a[idx1] && yseq[idx3] == b[idx2]) {
res1[idx2+(idx1-1)*nb] <- kriged[["map"]]$pred[idx3]
res2[idx2+(idx1-1)*nb] <- 1
break
}
}
}
}
res1[res2 == 0] <- NA
c <- matrix(res1, na, nb, byrow = TRUE)
plot.title <- paste(data$method,"synergy score:", summary.score, sep = " ")
lut <- c(1e-24, 1e-21, 1e-18, 1e-15, 1e-12, 1e-09, 1e-06,
0.001, 1, 1000, 1e+06, 1e+09, 1e+12, 1e+15, 1e+18, 1e+21,
1e+24)
pre <- c("y", "z", "a", "f", "p", "n", "µ", "m", "", "k",
"M", "G", "T", "P", "E", "Z", "Y")
conc.unit <- drug.pairs$concUnit[i] ## concentration unit
conc.unit.prefix <- substr(conc.unit, 1, 1)
unit.text.x <- paste("(", x.unit, "M)", sep = "")
unit.text.y <- paste("(", y.unit, "M)", sep = "")
file.name <- paste(drug.row, drug.col, "synergy", drug.pairs$blockIDs[i], data$method, "pdf", sep = ".")
drug.row <- paste(drug.row, unit.text.x, sep = " ")
drug.col <- paste(drug.col, unit.text.y, sep = " ")
max.dose <- max(abs(max(scores.dose)), abs(min(scores.dose)))
color.range <- round(max.dose + 5, -1)
if(is.null(legend.start)){
start.point <- -color.range
} else {
start.point <- legend.start
}
if(is.null(legend.end)) {
end.point <- color.range
} else {
end.point <- legend.end
}
if (type == "3D") {
plots.3d <- c
if(!is.null(x.range)) {
row.start <- (x.range[1] - 1) * 5
row.end <- (x.range[2] - 1) * 5
plots.3d <- plots.3d[row.start:row.end, ]
}
if(!is.null(y.range)) {
col.start <- (y.range[1] - 1) * 5
col.end <- (y.range[2] - 1) * 5
plots.3d <- plots.3d[, col.start:col.end]
}
x.conc1 <- x.conc
y.conc1 <- y.conc
if(!is.null(x.range)) {
x.conc1 <- x.conc1[x.range[1]:x.range[2]]
}
if(!is.null(y.range)) {
y.conc1 <- y.conc1[y.range[1]:y.range[2]]
}
# conc1 -> column.values -- ncol(plots.3d)
# conc2 -> row.values -- nrow(plots.3d)
c1 = data.frame(point=seq(from = 1, to = ncol(plots.3d), length.out = length(x.conc1)), conc=x.conc1)
c2 = data.frame(point=seq(from = 1, to = nrow(plots.3d), length.out = length(y.conc1)), conc=y.conc1)
column.values = spline(c1$point, c1$conc, xout = (1:ncol(plots.3d)))
row.values = spline(c2$point, c2$conc, xout = (1:nrow(plots.3d)))
x.at = spline(c1$conc, c1$point, xout = x.conc1)$y
y.at = spline(c2$conc, c2$point, xout = y.conc1)$y
x_axtcks = axisTicks(c(min(y.conc1), max(y.conc1)), FALSE)
y_axtcks = axisTicks(c(min(x.conc1), max(x.conc1)), FALSE)
#axisTicks(c(min(b), max(b)), FALSE, nint=dev.size(units = "px")[[2]]/30)
fig <- wireframe(plots.3d,scales = list(arrows = FALSE,distance=c(0.8,0.8,0.8),col=1,cex=0.8,z = list(tick.number=7)
,
x=list(at = x_axtcks, labels = f2si2(x_axtcks*lut[pre == conc.unit.prefix], force_unit = x.unit)),
y=list(at = y_axtcks, labels = f2si2(y_axtcks*lut[pre == conc.unit.prefix], force_unit = y.unit))
),
drape = TRUE, colorkey = list(space="top",width=0.5),
screen = list(z = 30, x = -55),
zlab=list(expression("Synergy score"),rot=90,cex=1,axis.key.padding = 0),xlab=list(as.character(drug.row),cex=1, rot=20),ylab=list(as.character(drug.col),cex=1,rot=-50),
zlim=c(start.point, end.point),
col.regions=colorRampPalette(c("green","white","red"))(100),
main = plot.title,
at=do.breaks(c(start.point, end.point),100),
par.settings = list(axis.line=list(col="transparent")),
zoom = 1, aspect = 1,
row.values = row.values$y,
column.values = column.values$y
#par.settings=list(layout.widths=list(left.padding=0,right.padding=0), layout.heights=list(top.padding=0, bottom.padding=0)) # margin
)
print(fig)
fig <- recordPlot()
} else if (type == "2D") {
dev.new(noRStudioGD = TRUE)
layout(matrix(c(1, 2), nrow = 2L, ncol = 1L), heights = c(0.1, 1))
par(mar = c(0, 6.1, 2.1, 4.1))
suppressWarnings(par(mgp = c(3, -0.4, 0)))
levels <- seq(start.point, end.point, by = 2)
col <- colorpanel(end.point, low = "green", mid = "white", high = "red")
plot.new()
plot.window(ylim = c(0, 1), xlim = range(levels), xaxs = "i", yaxs = "i")
rect(levels[-length(levels)],0, levels[-1L],0.3, col = col, border = NA)
axis(3,tick = FALSE, at = do.breaks(c(start.point, end.point), end.point/10))
title(plot.title)
par(mar = c(5.1,4.1,2.1,2.1))
suppressWarnings(par(mgp = c(2,1,0)))
plot.new()
if(!is.null(x.range)) {
row.start <- (x.range[1] - 1) * 5
row.end <- (x.range[2] - 1) * 5
c <- c[row.start:row.end, ]
}
if(!is.null(y.range)) {
col.start <- (y.range[1] - 1) * 5
col.end <- (y.range[2] - 1) * 5
c <- c[, col.start:col.end]
}
x.2D <- (1:dim(c)[1] - 1) / (dim(c)[1] - 1)
y.2D <- (1:dim(c)[2] - 1) / (dim(c)[2] - 1)
plot.window(asp = NA, xlim = range(x.2D), ylim = range(y.2D), "", xaxs = "i", yaxs = "i")
.filled.contour(x.2D, y.2D, z = c, levels, col = col)
#grid(dim(c)[1] - 1, dim(c)[2] - 1, col = "gray", lty = "solid", lwd = 0.3)
box()
mtext(drug.col, 1, cex = 1, padj = 3)
mtext(drug.row, 2, cex = 1, las = 3, padj = -3)
if(!is.null(x.range)) {
x.conc <- x.conc[x.range[1]:x.range[2]]
}
if(!is.null(y.range)) {
y.conc <- y.conc[y.range[1]:y.range[2]]
}
axis(side = 1, at = seq(0, 1, by = 1/(length(x.conc) - 1)), labels = signif(x.conc, 1))
axis(side = 2, at = seq(0, 1, by = 1/(length(y.conc) - 1)), labels = signif(y.conc, 1))
fig <- recordPlot()
dev.off()
} else {
plots.3d <- c
if(!is.null(x.range)) {
row.start <- (x.range[1] - 1) * 5
row.end <- (x.range[2] - 1) * 5
plots.3d <- plots.3d[row.start:row.end, ]
}
if(!is.null(y.range)) {
col.start <- (y.range[1] - 1) * 5
col.end <- (y.range[2] - 1) * 5
plots.3d <- plots.3d[, col.start:col.end]
}
x.conc1 <- x.conc
y.conc1 <- y.conc
if(!is.null(x.range)) {
x.conc1 <- x.conc1[x.range[1]:x.range[2]]
}
if(!is.null(y.range)) {
y.conc1 <- y.conc1[y.range[1]:y.range[2]]
}
# conc1 -> column.values -- ncol(plots.3d)
# conc2 -> row.values -- nrow(plots.3d)
c1 = data.frame(point=seq(from = 1, to = ncol(plots.3d), length.out = length(x.conc1)), conc=x.conc1)
c2 = data.frame(point=seq(from = 1, to = nrow(plots.3d), length.out = length(y.conc1)), conc=y.conc1)
column.values = spline(c1$point, c1$conc, xout = (1:ncol(plots.3d)))
row.values = spline(c2$point, c2$conc, xout = (1:nrow(plots.3d)))
x.at = spline(c1$conc, c1$point, xout = x.conc1)$y
y.at = spline(c2$conc, c2$point, xout = y.conc1)$y
x_axtcks = axisTicks(c(min(y.conc1), max(y.conc1)), FALSE)
y_axtcks = axisTicks(c(min(x.conc1), max(x.conc1)), FALSE)
syn.3d.plot <- wireframe(plots.3d,
scales = list(arrows = FALSE,distance=c(0.8,0.8,0.8),col=1,cex=0.8,z = list(tick.number=7)
,
x=list(at = x_axtcks, labels = f2si2(x_axtcks*lut[pre == conc.unit.prefix], force_unit = x.unit)),
y=list(at = y_axtcks, labels = f2si2(y_axtcks*lut[pre == conc.unit.prefix], force_unit = y.unit))
),
drape = TRUE, colorkey = list(space="top",width=0.5),
screen = list(z = 30, x = -55),
zlab=list(expression("Synergy score"),rot=90,cex=1,axis.key.padding = 0),xlab=list(as.character(drug.row),cex=1, rot=20),ylab=list(as.character(drug.col),cex=1,rot=-50),
zlim=c(start.point, end.point),
col.regions=colorRampPalette(c("green","white","red"))(100),
main = plot.title,
at=do.breaks(c(start.point, end.point),100),
par.settings = list(axis.line=list(col="transparent")),
zoom = 1, aspect = 1,
row.values = row.values$y,
column.values = column.values$y
#par.settings=list(layout.widths=list(left.padding=0,right.padding=0), layout.heights=list(top.padding=0, bottom.padding=0)) # margin
)
layout(matrix(c(1, 2, 3, 3), nrow = 2L, ncol = 2L), heights = c(0.1, 1))
par(mar = c(0, 6.1, 2.1, 4.1))
suppressWarnings(par(mgp = c(3, -1, 0)))
levels <- seq(start.point, end.point, by = 2)
col <- colorpanel(end.point, low = "green", mid = "white", high = "red")
plot.new()
plot.window(ylim = c(0, 1), xlim = range(levels), xaxs = "i", yaxs = "i")
rect(levels[-length(levels)],0, levels[-1L],0.3, col = col, border = NA)
axis(3,tick = FALSE, at = do.breaks(c(start.point, end.point), end.point/10))
title(plot.title)
par(mar = c(5.1,4.1,2.1,2.1))
suppressWarnings(par(mgp = c(2,1,0)))
plot.new()
if(!is.null(x.range)) {
row.start <- (x.range[1] - 1) * 5
row.end <- (x.range[2] - 1) * 5
c <- c[row.start:row.end, ]
}
if(!is.null(y.range)) {
col.start <- (y.range[1] - 1) * 5
col.end <- (y.range[2] - 1) * 5
c <- c[, col.start:col.end]
}
x.2D <- (1:dim(c)[1] - 1) / (dim(c)[1] - 1)
y.2D <- (1:dim(c)[2] - 1) / (dim(c)[2] - 1)
plot.window(asp = NA, xlim = range(x.2D), ylim = range(y.2D), "", xaxs = "i", yaxs = "i")
.filled.contour(x.2D, y.2D, z = c, levels, col = col)
#grid(dim(c)[1] - 1, dim(c)[2] - 1, col = "gray", lty = "solid", lwd = 0.3)
box()
mtext(drug.col, 1, cex = 1, padj = 3)
mtext(drug.row, 2, cex = 1, las = 3, padj = -3)
if(!is.null(x.range)) {
x.conc <- x.conc[x.range[1]:x.range[2]]
}
if(!is.null(y.range)) {
y.conc <- y.conc[y.range[1]:y.range[2]]
}
axis(side = 1, at = seq(0, 1, by = 1/(length(x.conc) - 1)), labels = signif(x.conc, 1))
axis(side = 2, at = seq(0, 1, by = 1/(length(y.conc) - 1)), labels = signif(y.conc, 1))
print(syn.3d.plot, position = c(0.5,0, 1, 1), newpage = FALSE)
fig <- recordPlot()
}
plots[[i]] <- fig
if(save.file) {
if (type %in% c("2D", "3D")) {
pdf(file.name, width = 10, height = 10)
} else {
pdf(file.name, width = 12, height = 6)
}
print(fig)
dev.off()
}
}
if(!save.file) {
for(i in num.pairs) {
dev.new(noRStudioGD = TRUE)
replayPlot(plots[[i]])
}
}
}
PlotEffect <- function(data, knitr = FALSE, type = "2D", save.file = FALSE,
effect.label = "Effect",
pair.index = NULL, legend.start = NULL, legend.end = NULL, x.range = NULL, y.range = NULL,
x.unit = "µ", y.unit = "µ") {
if(!is.list(data)) {
stop("Input data is not a list format!")
}
scores <- data$dose.response.mats
drug.pairs <- data$drug.pairs
num.pairs <- 1:length(scores)
plots <- list()
if(!is.null(pair.index)) {
num.pairs <- pair.index
}
for (i in num.pairs) {
scores.dose <- t(scores[[i]])
drug.col <- drug.pairs$drug.col[i]
drug.row <- drug.pairs$drug.row[i]
scores.tmp <- scores.dose
if(!is.null(x.range)) {
if(x.range[1] == 1) {
scores.tmp <- scores.tmp[(x.range[1] + 1):x.range[2], ]
} else {
scores.tmp <- scores.tmp[x.range[1]:x.range[2], ]
}
} else {
# delete the first row
scores.tmp <- scores.tmp[-1, ]
}
if(!is.null(y.range)) {
if(y.range[1] == 1) {
scores.tmp <- scores.tmp[, (y.range[1] + 1):y.range[2]]
} else {
scores.tmp <- scores.tmp[, y.range[1]:y.range[2]]
}
} else {
# delete the first column
scores.tmp <- scores.tmp[, -1]
}
#summary.score <- round(mean(scores.dose[-1, -1]), 3)
summary.score <- round(mean(scores.tmp, na.rm = TRUE), 3)
# drug.col: the col-wise drug
# drug.row: the row-wise drug
x.conc <- as.numeric(rownames(scores.dose)) ## col-wise drug concentration
y.conc <- as.numeric(colnames(scores.dose)) ## row-wise drug concentration
# drug.col.conc in index
x <- c(0:(length(x.conc) - 1))
# drug.row.conc in index
y <- c(0:(length(y.conc) - 1))
# smoothing by kriging
tmp <- cbind(expand.grid(x, y),c(as.matrix(scores.dose)))
x.vec <- tmp[, 1]
y.vec <- tmp[, 2]
scores.dose.vec = tmp[, 3]
#kriged = kriging(x.vec,y.vec,delta.dose.vec,lags=3,pixels=50)
pixels.num <- 5 * (length(x.conc) - 1) + 2
if (dim(scores.dose)[1] < 8) {
kriged <- kriging(y.vec, x.vec, scores.dose.vec,lags = 2, pixels = pixels.num, model = "spherical")
# kriged <- kriging(x.vec, y.vec, scores.dose.vec,lags = 2, pixels = pixels.num, model = "spherical")
} else {
kriged <- kriging(y.vec, x.vec, scores.dose.vec,lags = 3, pixels = pixels.num, model = "spherical")
# kriged <- kriging(x.vec, y.vec, scores.dose.vec,lags = 3, pixels = pixels.num, model = "spherical")
}
xseq <- round(kriged[["map"]]$x / kriged$pixel)
yseq <- round(kriged[["map"]]$y / kriged$pixel)
a <- min(xseq):max(xseq)
b <- min(yseq):max(yseq)
z.len <- length(kriged[["map"]]$pred) ## n
na <- length(a)
nb <- length(b)
res1 <- as.double(rep(0, na * nb)) ## c
res2 <- as.integer(rep(0,na * nb)) ## d
for(idx1 in 1:na) {
for(idx2 in 1:nb) {
for(idx3 in 1:z.len) {
if(xseq[idx3] == a[idx1] && yseq[idx3] == b[idx2]) {
res1[idx2+(idx1-1)*nb] <- kriged[["map"]]$pred[idx3]
res2[idx2+(idx1-1)*nb] <- 1
break
}
}
}
}
res1[res2 == 0] <- NA
c <- matrix(res1, na, nb, byrow = TRUE)
# plot.title <- paste(data$method,"synergy score:", summary.score, "(Effect)", sep = " ")
plot.title <- paste("Dose-response plane", sep = " ")
lut <- c(1e-24, 1e-21, 1e-18, 1e-15, 1e-12, 1e-09, 1e-06,
0.001, 1, 1000, 1e+06, 1e+09, 1e+12, 1e+15, 1e+18, 1e+21,
1e+24)
pre <- c("y", "z", "a", "f", "p", "n", "µ", "m", "", "k",
"M", "G", "T", "P", "E", "Z", "Y")
conc.unit <- drug.pairs$concUnit[i] ## concentration unit
conc.unit.prefix <- substr(conc.unit, 1, 1)
unit.text.x <- paste("(", x.unit, "M)", sep = "")
unit.text.y <- paste("(", y.unit, "M)", sep = "")
file.name <- paste(drug.row, drug.col, "effect", drug.pairs$blockIDs[i], data$method, "pdf", sep = ".")
drug.row <- paste(drug.row, unit.text.x, sep = " ")
drug.col <- paste(drug.col, unit.text.y, sep = " ")
max.dose <- max(abs(max(scores.dose)), abs(min(scores.dose)))
color.range <- round(max.dose + 5, -1)
if(is.null(legend.start)){
start.point <- 0
} else {
start.point <- legend.start
}
if(is.null(legend.end)) {
end.point <- color.range
} else {
end.point <- legend.end
}
if (type == "3D") {
plots.3d <- c
if(!is.null(x.range)) {
row.start <- (x.range[1] - 1) * 5
row.end <- (x.range[2] - 1) * 5
plots.3d <- plots.3d[row.start:row.end, ]
}
if(!is.null(y.range)) {
col.start <- (y.range[1] - 1) * 5
col.end <- (y.range[2] - 1) * 5
plots.3d <- plots.3d[, col.start:col.end]
}
x.conc1 <- x.conc
y.conc1 <- y.conc
if(!is.null(x.range)) {
x.conc1 <- x.conc1[x.range[1]:x.range[2]]
}
if(!is.null(y.range)) {
y.conc1 <- y.conc1[y.range[1]:y.range[2]]
}
# conc1 -> column.values -- ncol(plots.3d)
# conc2 -> row.values -- nrow(plots.3d)
c1 = data.frame(point=seq(from = 1, to = ncol(plots.3d), length.out = length(x.conc1)), conc=x.conc1)
c2 = data.frame(point=seq(from = 1, to = nrow(plots.3d), length.out = length(y.conc1)), conc=y.conc1)
column.values = spline(c1$point, c1$conc, xout = (1:ncol(plots.3d)))
row.values = spline(c2$point, c2$conc, xout = (1:nrow(plots.3d)))
x.at = spline(c1$conc, c1$point, xout = x.conc1)$y
y.at = spline(c2$conc, c2$point, xout = y.conc1)$y
x_axtcks = axisTicks(c(min(y.conc1), max(y.conc1)), FALSE)
y_axtcks = axisTicks(c(min(x.conc1), max(x.conc1)), FALSE)
fig <- wireframe(plots.3d,scales = list(arrows = FALSE,distance=c(0.8,0.8,0.8),col=1,cex=0.8,z = list(tick.number=7)
,
x=list(at = x_axtcks, labels = f2si2(x_axtcks*lut[pre == conc.unit.prefix], force_unit = x.unit)),
y=list(at = y_axtcks, labels = f2si2(y_axtcks*lut[pre == conc.unit.prefix], force_unit = y.unit))
),
drape = TRUE, colorkey = list(space="top",width=0.5),
screen = list(z = 30, x = -55),
zlab=list(effect.label,rot=90,cex=1,axis.key.padding = 0),xlab=list(as.character(drug.row),cex=1, rot=20),ylab=list(as.character(drug.col),cex=1,rot=-50),
zlim=c(start.point, end.point),
col.regions=colorRampPalette(c("green","white","red"))(100),
main = plot.title,
at=do.breaks(c(start.point, end.point),100),
par.settings = list(axis.line=list(col="transparent")),
zoom = 1, aspect = 1,
row.values = row.values$y,
column.values = column.values$y
#par.settings=list(layout.widths=list(left.padding=0,right.padding=0), layout.heights=list(top.padding=0, bottom.padding=0)) # margin
)
print(fig)
fig <- recordPlot()
} else if (type == "2D") {
dev.new(noRStudioGD = TRUE)
layout(matrix(c(1, 2), nrow = 2L, ncol = 1L), heights = c(0.1, 1))
par(mar = c(0, 6.1, 2.1, 4.1))
suppressWarnings(par(mgp = c(3, -0.4, 0)))
levels <- seq(start.point, end.point, by = 2)
col <- colorpanel(end.point, low = "green", mid = "white", high = "red")
plot.new()
plot.window(ylim = c(0, 1), xlim = range(levels), xaxs = "i", yaxs = "i")
rect(levels[-length(levels)],0, levels[-1L],0.3, col = col, border = NA)
axis(3,tick = FALSE, at = do.breaks(c(start.point, end.point), end.point/10))
title(plot.title)
par(mar = c(5.1,4.1,2.1,2.1))
suppressWarnings(par(mgp = c(2,1,0)))
plot.new()
if(!is.null(x.range)) {
row.start <- (x.range[1] - 1) * 5
row.end <- (x.range[2] - 1) * 5
c <- c[row.start:row.end, ]
}
if(!is.null(y.range)) {
col.start <- (y.range[1] - 1) * 5
col.end <- (y.range[2] - 1) * 5
c <- c[, col.start:col.end]
}
x.2D <- (1:dim(c)[1] - 1) / (dim(c)[1] - 1)
y.2D <- (1:dim(c)[2] - 1) / (dim(c)[2] - 1)
plot.window(asp = NA, xlim = range(x.2D), ylim = range(y.2D), "", xaxs = "i", yaxs = "i")
.filled.contour(x.2D, y.2D, z = c, levels, col = col)
#grid(dim(c)[1] - 1, dim(c)[2] - 1, col = "gray", lty = "solid", lwd = 0.3)
box()
mtext(drug.col, 1, cex = 1, padj = 3)
mtext(drug.row, 2, cex = 1, las = 3, padj = -3)
if(!is.null(x.range)) {
x.conc <- x.conc[x.range[1]:x.range[2]]
}
if(!is.null(y.range)) {
y.conc <- y.conc[y.range[1]:y.range[2]]
}
axis(side = 1, at = seq(0, 1, by = 1/(length(x.conc) - 1)), labels = signif(x.conc, 1))
axis(side = 2, at = seq(0, 1, by = 1/(length(y.conc) - 1)), labels = signif(y.conc, 1))
fig <- recordPlot()
dev.off()
} else {
plots.3d <- c
if(!is.null(x.range)) {
row.start <- (x.range[1] - 1) * 5
row.end <- (x.range[2] - 1) * 5
plots.3d <- plots.3d[row.start:row.end, ]
}
if(!is.null(y.range)) {
col.start <- (y.range[1] - 1) * 5
col.end <- (y.range[2] - 1) * 5
plots.3d <- plots.3d[, col.start:col.end]
}
x.conc1 <- x.conc
y.conc1 <- y.conc
if(!is.null(x.range)) {
x.conc1 <- x.conc1[x.range[1]:x.range[2]]
}
if(!is.null(y.range)) {
y.conc1 <- y.conc1[y.range[1]:y.range[2]]
}
# conc1 -> column.values -- ncol(plots.3d)
# conc2 -> row.values -- nrow(plots.3d)
c1 = data.frame(point=seq(from = 1, to = ncol(plots.3d), length.out = length(x.conc1)), conc=x.conc1)
c2 = data.frame(point=seq(from = 1, to = nrow(plots.3d), length.out = length(y.conc1)), conc=y.conc1)
column.values = spline(c1$point, c1$conc, xout = (1:ncol(plots.3d)))
row.values = spline(c2$point, c2$conc, xout = (1:nrow(plots.3d)))
x.at = spline(c1$conc, c1$point, xout = x.conc1)$y
y.at = spline(c2$conc, c2$point, xout = y.conc1)$y
x_axtcks = axisTicks(c(min(y.conc1), max(y.conc1)), FALSE)
y_axtcks = axisTicks(c(min(x.conc1), max(x.conc1)), FALSE)
syn.3d.plot <- wireframe(plots.3d,
scales = list(arrows = FALSE,distance=c(0.8,0.8,0.8),col=1,cex=0.8,z = list(tick.number=7)
,
x=list(at = x_axtcks, labels = f2si2(x_axtcks*lut[pre == conc.unit.prefix], force_unit = x.unit)),
y=list(at = y_axtcks, labels = f2si2(y_axtcks*lut[pre == conc.unit.prefix], force_unit = y.unit))
),
drape = TRUE, colorkey = list(space="top",width=0.5),
screen = list(z = 30, x = -55),
zlab=list(expression("Synergy score"),rot=90,cex=1,axis.key.padding = 0),xlab=list(as.character(drug.row),cex=1, rot=20),ylab=list(as.character(drug.col),cex=1,rot=-50),
zlim=c(start.point, end.point),
col.regions=colorRampPalette(c("green","white","red"))(100),
main = plot.title,
at=do.breaks(c(start.point, end.point),100),
par.settings = list(axis.line=list(col="transparent")),
zoom = 1, aspect = 1,
row.values = row.values$y,
column.values = column.values$y
#par.settings=list(layout.widths=list(left.padding=0,right.padding=0), layout.heights=list(top.padding=0, bottom.padding=0)) # margin
)
layout(matrix(c(1, 2, 3, 3), nrow = 2L, ncol = 2L), heights = c(0.1, 1))
par(mar = c(0, 6.1, 2.1, 4.1))
suppressWarnings(par(mgp = c(3, -1, 0)))
levels <- seq(start.point, end.point, by = 2)
col <- colorpanel(end.point, low = "green", mid = "white", high = "red")
plot.new()
plot.window(ylim = c(0, 1), xlim = range(levels), xaxs = "i", yaxs = "i")
rect(levels[-length(levels)],0, levels[-1L],0.3, col = col, border = NA)
axis(3,tick = FALSE, at = do.breaks(c(start.point, end.point), end.point/10))
title(plot.title)
par(mar = c(5.1,4.1,2.1,2.1))
suppressWarnings(par(mgp = c(2,1,0)))
plot.new()
if(!is.null(x.range)) {
row.start <- (x.range[1] - 1) * 5
row.end <- (x.range[2] - 1) * 5
c <- c[row.start:row.end, ]
}
if(!is.null(y.range)) {
col.start <- (y.range[1] - 1) * 5
col.end <- (y.range[2] - 1) * 5
c <- c[, col.start:col.end]
}
x.2D <- (1:dim(c)[1] - 1) / (dim(c)[1] - 1)
y.2D <- (1:dim(c)[2] - 1) / (dim(c)[2] - 1)
plot.window(asp = NA, xlim = range(x.2D), ylim = range(y.2D), "", xaxs = "i", yaxs = "i")
.filled.contour(x.2D, y.2D, z = c, levels, col = col)
#grid(dim(c)[1] - 1, dim(c)[2] - 1, col = "gray", lty = "solid", lwd = 0.3)
box()
mtext(drug.col, 1, cex = 1, padj = 3)
mtext(drug.row, 2, cex = 1, las = 3, padj = -3)
if(!is.null(x.range)) {
x.conc <- x.conc[x.range[1]:x.range[2]]
}
if(!is.null(y.range)) {
y.conc <- y.conc[y.range[1]:y.range[2]]
}
axis(side = 1, at = seq(0, 1, by = 1/(length(x.conc) - 1)), labels = signif(x.conc, 1))
axis(side = 2, at = seq(0, 1, by = 1/(length(y.conc) - 1)), labels = signif(y.conc, 1))
print(syn.3d.plot, position = c(0.5,0, 1, 1), newpage = FALSE)
fig <- recordPlot()
}
plots[[i]] <- fig
if(save.file) {
if (type %in% c("2D", "3D")) {
pdf(file.name, width = 10, height = 10)
} else {
pdf(file.name, width = 12, height = 6)
}
print(fig)
dev.off()
}
}
if(!save.file) {
for(i in num.pairs) {
dev.new(noRStudioGD = TRUE)
replayPlot(plots[[i]])
}
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.