#' @name compDrugsen
#' @title Comparison of drug sensitivity
#' @description This function estimates the IC50 of specific drugs for each Subtype by developing a ridge regression predictive model based on all/specific cell lines derived from Genomics of Drug Sensitivity in Cancer (GDSC, \url{https://www.cancerrxgene.org/}).
#' @param moic.res An object returned by `getMOIC()` with one specified algorithm or `get\%algorithm_name\%` or `getConsensusMOIC()` with a list of multiple algorithms.
#' @param norm.expr A matrix of normalized expression data with rows for genes and columns for samples; FPKM or TPM without log2 transformation is recommended.
#' @param drugs A string vector to indicate the names of the drugs for which you would like to predict sensitivity, one of Erlotinib, Rapamycin, Sunitinib, PHA-665752, MG-132, Paclitaxel, Cyclopamine, AZ628, Sorafenib, VX-680, Imatinib, TAE684, Crizotinib, Saracatinib, S-Trityl-L-cysteine, Z-LLNle-CHO, Dasatinib, GNF-2, CGP-60474, CGP-082996, A-770041, WH-4-023, WZ-1-84, BI-2536, BMS-536924, BMS-509744, CMK, Pyrimethamine, JW-7-52-1, A-443654, GW843682X, MS-275, Parthenolide, KIN001-135, TGX221, Bortezomib, XMD8-85, Roscovitine, Salubrinal, Lapatinib, GSK269962A, Doxorubicin, Etoposide, Gemcitabine, Mitomycin C, Vinorelbine, NSC-87877, Bicalutamide, QS11, CP466722, Midostaurin, CHIR-99021, AP-24534, AZD6482, JNK-9L, PF-562271, HG-6-64-1, JQ1, JQ12, DMOG, FTI-277, OSU-03012, Shikonin, AKT inhibitor VIII, Embelin, FH535, PAC-1, IPA-3, GSK-650394, BAY 61-3606, 5-Fluorouracil, Thapsigargin, Obatoclax Mesylate, BMS-754807, Lisitinib, Bexarotene, Bleomycin, LFM-A13, GW-2580, AUY922, Phenformin, Bryostatin 1, Pazopanib, LAQ824, Epothilone B, GSK1904529A, BMS345541, Tipifarnib, BMS-708163, Ruxolitinib, AS601245, Ispinesib Mesylate, TL-2-105, AT-7519, TAK-715, BX-912, ZSTK474, AS605240, Genentech Cpd 10, GSK1070916, KIN001-102, LY317615, GSK429286A, FMK, QL-XII-47, CAL-101, UNC0638, XL-184, WZ3105, XMD14-99, AC220, CP724714, JW-7-24-1, NPK76-II-72-1, STF-62247, NG-25, TL-1-85, VX-11e, FR-180204, Tubastatin A, Zibotentan, YM155, NSC-207895, VNLG/124, AR-42, CUDC-101, Belinostat, I-BET-762, CAY10603, Linifanib , BIX02189, CH5424802, EKB-569, GSK2126458, KIN001-236, KIN001-244, KIN001-055, KIN001-260, KIN001-266, Masitinib, MP470, MPS-1-IN-1, BHG712, OSI-930, OSI-027, CX-5461, PHA-793887, PI-103, PIK-93, SB52334, TPCA-1, TG101348, Foretinib, Y-39983, YM201636, Tivozanib, GSK690693, SNX-2112, QL-XI-92, XMD13-2, QL-X-138, XMD15-27; two common chemodrugs (i.e., Cisplatin and Paclitaxel) will be analyzed by default if no indication.
#' @param tissueType A string value to specify if you would like to train the models on only a subset of the CGP cell lines (based on the tissue type from which the cell lines originated); Allowed values contain c("all", "aero_digestive_tract", "blood", "bone", "breast", "digestive_system", "lung", "nervous_system", "skin", "urogenital_system") and "all" by default.
#' @param clust.col A string vector storing colors for annotating each Subtype.
#' @param prefix A string value to indicate the prefix of the output plot.
#' @param fig.path A string value to indicate the output path for storing the boxviolin plot.
#' @param seed A integer value to indicate the seed for reproducing ridge regression.
#' @param width A numeric value to indicate the width of boxviolin plot.
#' @param height A numeric value to indicate the height of boxviolin plot.
#' @param test.method A string value to indicate the method for statistical testing. Allowed values contain c('nonparametric', 'parametric'); nonparametric means two-sample wilcoxon rank sum test for two subtypes and Kruskal-Wallis rank sum test for multiple subtypes; parametric means two-sample t-test when only two subtypes are identified, and anova for multiple subtypes comparison; "nonparametric" by default.
#' @return Data.frame(s) storing the estimated IC50 of specified drugs per sample within each Subtype.
#' @export
#' @import ggplot2
#' @importFrom ggpubr stat_compare_means
#' @references Geeleher P, Cox N, Huang R S. (2014). pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One, 9(9):e107468.
#'
#' Geeleher P, Cox N J, Huang R S. (2014). Clinical drug response can be predicted using baseline gene expression levels and in vitro drug sensitivity in cell lines. Genome Biol, 15(3):1-12.
#' @examples # There is no example and please refer to vignette.
compDrugsen <- function(moic.res = NULL,
norm.expr = NULL,
drugs = c("Cisplatin", "Paclitaxel"),
tissueType = "all",
test.method = "nonparametric",
clust.col = c("#2EC4B6","#E71D36","#FF9F1C","#BDD5EA","#FFA5AB","#011627","#023E8A","#9D4EDD"),
prefix = NULL,
seed = 123456,
fig.path = getwd(),
width = 5,
height = 5) {
if(!is.element(test.method, c("nonparametric","parametric"))) {
stop("test.method can be one of nonparametric or parametric.")
}
# data processing
comsam <- intersect(moic.res$clust.res$samID, colnames(norm.expr))
# check data
if(length(comsam) == nrow(moic.res$clust.res)) {
message("--all samples matched.")
} else {
message(paste0("--",(nrow(moic.res$clust.res)-length(comsam))," samples mismatched from current subtypes."))
}
moic.res$clust.res <- moic.res$clust.res[comsam,,drop = FALSE]
norm.expr <- norm.expr[,comsam]
n.moic <- length(unique(moic.res$clust.res$clust))
sam.order <- moic.res$clust.res[order(moic.res$clust.res$clust, decreasing = FALSE), "samID"]
colvec <- clust.col[1:length(unique(moic.res$clust.res$clust))]
names(colvec) <- paste0("CS",unique(moic.res$clust.res$clust))
annCol <- data.frame("Subtype" = paste0("CS",moic.res$clust.res[sam.order,"clust"]),
samID = sam.order,
row.names = sam.order,
stringsAsFactors = FALSE)
if(max(norm.expr) < 25 | (max(norm.expr) >= 25 & min(norm.expr) < 0)) {
message("--expression profile seems to have veen standardised (z-score or log transformation), no more action will be performed.")
gset <- norm.expr
}
if(max(norm.expr) >= 25 & min(norm.expr) >= 0){
message("--log2 transformation done for expression data.")
gset <- log2(norm.expr + 1)
}
# drug sensitivity prediction
predictedPtype <- predictedBoxdat <- list()
for (drug in drugs) {
set.seed(seed)
predictedPtype[[drug]] <- quiet(pRRopheticPredict(testMatrix = as.matrix(gset[,rownames(annCol)]),
drug = drug,
tissueType = tissueType,
dataset = "cgp2016",
minNumSamples = 5,
selection = 1)) # 1 indicate if multiple genes existed, mean value will be considered
if(!all(names(predictedPtype[[drug]]) == rownames(annCol))) {stop("name mismatched!\n")}
predictedBoxdat[[drug]] <- data.frame("Est.IC50" = predictedPtype[[drug]],
"Subtype" = as.character(annCol$Subtype),
row.names = names(predictedPtype[[drug]]),
stringsAsFactors = FALSE)
message(drug," done...")
# generate boxviolin plot with statistical testing
if(n.moic == 2 & test.method == "nonparametric") {
statistic = "wilcox.test"
ic50.test <- wilcox.test(predictedBoxdat[[drug]]$Est.IC50 ~ predictedBoxdat[[drug]]$Subtype)$p.value
cat(paste0("Wilcoxon rank sum test p value = ", formatC(ic50.test, format = "e", digits = 2), " for ", drug))
}
if(n.moic == 2 & test.method == "parametric") {
statistic = "t.test"
ic50.test <- t.test(predictedBoxdat[[drug]]$Est.IC50 ~ predictedBoxdat[[drug]]$Subtype)$p.value
cat(paste0("Student's t test p value = ", formatC(ic50.test, format = "e", digits = 2), " for ", drug))
}
if(n.moic > 2 & test.method == "nonparametric") {
statistic = "kruskal.test"
ic50.test <- kruskal.test(predictedBoxdat[[drug]]$Est.IC50 ~ predictedBoxdat[[drug]]$Subtype)$p.value
pairwise.ic50.test <- pairwise.wilcox.test(predictedBoxdat[[drug]]$Est.IC50,predictedBoxdat[[drug]]$Subtype,p.adjust.method = "BH")
cat(paste0(drug,": Kruskal-Wallis rank sum test p value = ", formatC(ic50.test, format = "e", digits = 2),"\npost-hoc pairwise wilcoxon rank sum test with Benjamini-Hochberg adjustment presents below:\n"))
print(formatC(pairwise.ic50.test$p.value, format = "e", digits = 2))
}
if(n.moic > 2 & test.method == "parametric") {
statistic = "anova"
ic50.test <- summary(aov(predictedBoxdat[[drug]]$Est.IC50 ~ predictedBoxdat[[drug]]$Subtype))[[1]][["Pr(>F)"]][1]
pairwise.ic50.test <- pairwise.t.test(predictedBoxdat[[drug]]$Est.IC50,predictedBoxdat[[drug]]$Subtype,p.adjust.method = "BH")
cat(paste0(drug,": One-way anova test p value = ", formatC(ic50.test, format = "e", digits = 2),"\npost-hoc pairwise Student's t test with Benjamini-Hochberg adjustment presents below:\n"))
print(formatC(pairwise.ic50.test$p.value, format = "e", digits = 2))
}
p <- ggplot(data = predictedBoxdat[[drug]],
aes(x = Subtype, y = Est.IC50, fill = Subtype)) +
scale_fill_manual(values = colvec) +
geom_violin(alpha = 0.4, position = position_dodge(width = .75),
size = 0.8, color = "black") +
geom_boxplot(notch = TRUE, outlier.size = -1,
color = "black", lwd = 0.8, alpha = 0.7) +
geom_point(shape = 21, size = 2,
position = position_jitterdodge(),
color = "black", alpha = 1) +
theme_classic() +
ylab(bquote("Estimated IC"[50]~"of"~.(drug))) + xlab("") +
theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 12),
axis.ticks = element_line(size = 0.2, color = "black"),
axis.ticks.length = unit(0.2, "cm"),
legend.position = "none",
axis.title = element_text(size = 15),
axis.text = element_text(size = 10)) +
# add statistical inference
stat_compare_means(method = statistic,
hjust = ifelse(n.moic %% 2 == 0, 0.5, 0),
label.x = ifelse(n.moic %% 2 == 0, n.moic / 2 + 0.5, n.moic / 2),
label.y = min(predictedBoxdat[[drug]]$Est.IC50))
# save to pdf
if(is.null(prefix)) {
outFig <- paste0("boxviolin of estimated ic50 for ",drug,".pdf")
} else {
outFig <- paste0(prefix, " for ", drug,".pdf")
}
ggsave(file.path(fig.path, outFig), width = width, height = height)
# print to screen
print(p)
}
return(predictedBoxdat)
}
Datasets
Models
