MOVICS: Multi-Omics integration and VIsualization in Cancer Subtyping

Documented in compDrugsen

#' @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)
}