--- a
+++ b/R/compDrugsen.R
@@ -0,0 +1,151 @@
+#' @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)
+}