5.2.1 Effect of drugs with similar target

Here we compare the effect of drugs designed to target components of the same signalling pathway.

# select the data
mtcll = as.data.frame(t(giveViabMatrix(diag="CLL",
                                       screen="main",
                                       chnnl="viaraw.4_5")))
colnames(mtcll) = drugs[colnames(mtcll),"name"]

# function which plots the scatter plot
scatdr = function(drug1, drug2, coldot, mtNEW, min){
  
  dataNEW = mtNEW[,c(drug1, drug2)]
  colnames(dataNEW) = c("A", "B")

  p = ggplot(data=dataNEW,  aes(A,  B)) + geom_point(size=3, col=coldot, alpha=0.8) +
    labs(x = drug1, y = drug2) + ylim(c(min, 1.35)) +  xlim(c(min, 1.35)) +
    theme(panel.background = element_blank(),
          axis.text = element_text(size = 15),
          axis.title = element_text(size = rel(1.5)),
          axis.line.x = element_line(colour = "black", size = 0.5),
          axis.line.y = element_line(colour = "black", size = 0.5),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank()) +
    geom_smooth(method=lm) +
    geom_text(x=1, y=min+0.1,
              label=paste0("Pearson-R = ",
                           round(cor(dataNEW$A, dataNEW$B ), 2)),
              size = 5)
  
  return(p)
}

8.1.1 Additional functions

someMatch <- function(...) {
  rv <- match(...)
  if (all(is.na(rv)))
    stop(sprintf("`match` failed to match any of the following: %s",
                 paste(x[is.na(rv)], collapse=", ")))
  rv  
}

8.1.2 Preprocess IGHV gene usage

Find the largest categories, which we will represent by colours, and merge the others in a group other.

colnames(pData(lpdAll))

## [1] "Diagnosis"               "Gender"                 
## [3] "IGHV Uppsala gene usage" "IGHV Uppsala % SHM"     
## [5] "IGHV Uppsala U/M"        "ATPday0"                
## [7] "ATP48h"

gu    <- pData(lpdAll)$`IGHV Uppsala gene usage`
tabgu <- sort(table(gu), decreasing = TRUE)
biggestGroups <- names(tabgu)[1:5]
gu[ is.na(match(gu, biggestGroups)) & !is.na(gu) ] <- "other"
pData(lpdAll)$`IGHV gene usage` <- factor(gu, levels = c(biggestGroups, "other"))

8.1.3 Some drugs that we are particularly interested in

stopifnot(is.null(drugs$id))
drugs$id <- rownames(drugs)

targetedDrugNames <- c("ibrutinib", "idelalisib",  "PRT062607 HCl", 
   "duvelisib", "spebrutinib", "selumetinib", "MK-2206",  
   "everolimus", "encorafenib")
id1 <- safeMatch(targetedDrugNames, drugs$name)
targetedDrugs <- paste( rep(drugs[id1, "id"], each = 2), 4:5, sep="_" )

chemoDrugNames <- c("fludarabine", "doxorubicine",  "nutlin-3")
id2 <- safeMatch(chemoDrugNames, drugs$name)
chemoDrugs <- paste( rep(drugs[id2, "id"], each = 5), 3:5, sep="_" )

tzselDrugNames <- c("ibrutinib", "idelalisib", "duvelisib", "selumetinib", 
    "AZD7762", "MK-2206", "everolimus", "venetoclax", "thapsigargin", 
    "AT13387", "YM155", "encorafenib", "tamatinib", "ruxolitinib", 
    "PF 477736", "fludarabine", "nutlin-3")
id3 <- safeMatch(tzselDrugNames, drugs$name)
tzselDrugs <- unlist(lapply(seq(along = tzselDrugNames), function(i)
  paste(drugs[id3[i], "id"], 
      if (tzselDrugNames[i] %in% c("fludarabine", "nutlin-3")) 2:3 else 4:5, 
    sep = "_" )))

8.1.4 Feature weighting and score for dendrogram reordering

The weights are used for weighting the similarity metric used in heatmap clustering. weights$hclust affects the clustering of the patients. weights$score affects the dendrogram reordering of the drugs.

bcrDrugs <- c("ibrutinib", "idelalisib", "PRT062607 HCl", "spebrutinib")
everolID <- drugs$id[ safeMatch("everolimus",  drugs$name)]
bcrID    <- drugs$id[ safeMatch(bcrDrugs, drugs$name)]

is_BCR  <- 
  (fData(lpdAll)$id %in% bcrID)    & (fData(lpdAll)$subtype %in% paste(4:5))
is_mTOR <- 
  (fData(lpdAll)$id %in% everolID) & (fData(lpdAll)$subtype %in% paste(4:5))

myin <- function(x, y) as.numeric( (x %in% y) & !is.na(x) )

weights1 <- data.frame(
  hclust = rep(1, nrow(lpdAll)) + 1.75 * is_mTOR, 
  score  = as.numeric( is_BCR ), 
  row.names = rownames(lpdAll))

weights2 <- data.frame(
  row.names = tzselDrugs,
  hclust = myin(drugs$target_category[id3], "B-cell receptor") * 0.3 + 0.7,
  score  = rep(1, length(tzselDrugs))) 

Remove drugs that failed quality control: NSC 74859, bortezomib.

badDrugs <- c(bortezomib = "D_008", `NSC 74859` = "D_025")
stopifnot(identical(drugs[ badDrugs, "name"], names(badDrugs)))

candDrugs <- rownames(lpdAll)[  
  fData(lpdAll)$type=="viab" & !(fData(lpdAll)$id %in% badDrugs) &
  fData(lpdAll)$subtype %in% paste(2:5)
]

Threshold parameters: drugs are accepted that for at least effectNum samples have a viability effect less than or equal to effectVal. On the other hand, the mean viability effect should not be below viab.

thresh <- list(effectVal = 0.7, effectNum = 4, viab = 0.6, maxval = 1.1)
overallMean  <- rowMeans(Biobase::exprs(lpdAll)[candDrugs, ])
nthStrongest <- apply(Biobase::exprs(lpdAll)[candDrugs, ], 1,
                      function(x) sort(x)[thresh$effectNum])

par(mfrow = c(1, 3))
hist(overallMean,  breaks = 30, col = "pink")
abline(v = thresh$viab,      col="blue")
hist(nthStrongest, breaks = 30, col = "pink")
abline(v = thresh$effectVal, col="blue")
plot(overallMean, nthStrongest)
abline(h = thresh$effectVal, v = thresh$viab, col = "blue")

seldrugs1 and d1: as in the version of Figure 3A we had in the first submission to JCI. d2: two concentrations for each drug in tzselDrugNames.

seldrugs1 <- candDrugs[ overallMean >= thresh$viab &
                        nthStrongest <= thresh$effectVal ] %>%
            union(targetedDrugs) %>% 
            union(chemoDrugs)

d1 <- Biobase::exprs(lpdAll[seldrugs1,, drop = FALSE ]) %>%
  deckel(lower = 0, upper = thresh$maxval)

d2 <- Biobase::exprs(lpdAll[tzselDrugs,, drop = FALSE ]) %>%
  deckel(lower = 0, upper = thresh$maxval)

We are going to scale the data. But was is the best measure of scale (or spread)? Let’s explore different measures of spread. We’ll see that it does not seem to matter too much which one we use. We’ll use median centering and scaling by mad.

spreads <- sapply(list(mad = mad, `Q95-Q05` = function(x)
  diff(quantile(x, probs = c(0.05, 0.95)))), function(s) apply(d1, 1, s))
plot( spreads )

jj <- names(which( spreads[, "mad"] < 0.15 & spreads[, "Q95-Q05"] > 0.7))
jj

## [1] "D_041_2"

drugs[ stripConc(jj), "name" ]

## [1] "BAY 11-7085"

medianCenter_MadScale <- function(x) {
  s <- median(x)
  (x - s) / deckel(mad(x, center = s), lower = 0.05, upper = 0.2)
}

scaleDrugResp  <- function(x) t(apply(x, 1, medianCenter_MadScale)) 

scd1 <- scaleDrugResp(d1)
scd2 <- scaleDrugResp(d2)

8.3.1 Matrix row / column clustering

The helper function matClust clusters a matrix x, whose columns represent samples and whose rows represent drugs. Its arguments control how the columns are clustered:

• weights: a data.frame with a weight for each row of x. The weights are used in the computation of distances between columns and thus for column sorting. The data.frame’s column hclust contains the weights for hclust(dist()). The column score contains the weights for computing the score used for dendrogram reordering. See weights1 and weights2 defined above.
• colgroups: a factor by which to first split the columns before clustering
• reorderrows: logical. FALSE for previous behaviour (old Fig. 3A), TRUE for reordering the row dendrogram, too.

matClust <- function(x, 
                     rowweights,
                     colgroups   = factor(rep("all", ncol(x))),
                     reorderrows = FALSE) {
 
  stopifnot(is.data.frame(rowweights),
            c("hclust", "score") %in% colnames(rowweights),
            !is.null(rownames(rowweights)),
            !is.null(rownames(x)), 
            all(rownames(x) %in% rownames(rowweights)),
            is.factor(colgroups), 
            !any(is.na(colgroups)), 
            length(colgroups) == ncol(x))

  wgt <- rowweights[ rownames(x), ]
  
  columnsClust <- function(xk) {
    score   <- -svd(xk * wgt[, "score"])$v[, 1]
    cmns    <- colSums(xk * wgt[, "score"])
    ## make sure that high score = high sensitivity  
    if (cor(score, cmns) > 0) score <- (-score)
    
    ddraw  <- as.dendrogram(hclust(dist(t(xk * wgt[, "hclust"]), 
                                        method = "euclidean"), 
                                   method = "complete"))
    dd  <- reorder(ddraw, wts = -score, agglo.FUN = mean)
    ord <- order.dendrogram(dd)
    list(dd = dd, ord = ord, score = score)
  }

  sp <- split(seq(along = colgroups), colgroups)
  cc <- lapply(sp, function(k) columnsClust(x[, k, drop=FALSE]))
  cidx <- unlist(lapply(seq(along = cc), function (i)
    sp[[i]][ cc[[i]]$ord ]))
  csc  <- unlist(lapply(seq(along = cc), function (i)
    cc[[i]]$score[ cc[[i]]$ord ]))

  rddraw <- as.dendrogram(hclust(dist(x, method = "euclidean"),
                                   method = "complete"))
  ridx  <- if (reorderrows) {
    ww <- (colgroups == "CLL")
    stopifnot(!any(is.na(ww)), any(ww))
    rowscore <- svd(t(x) * ww)$v[, 1]
    dd <- reorder(rddraw, wts = rowscore, agglo.FUN = mean)
    order.dendrogram(dd)
  } else {
    rev(order.dendrogram(dendsort(rddraw)))
  }

  res <- x[ridx, cidx]
  
  stopifnot(identical(dim(res), dim(x)))
  attr(res, "colgap") <- cumsum(sapply(cc, function(x) length(x$score)))
  res
}

8.3.2 Prepare sample annotations

I.e. the right hand side color bar. IGHV Uppsala U/M is implied by IGHV Uppsala % SHM (see sensi2.Rmd). cut(..., right=FALSE) will use intervals that are closed on the left and open on the right.

translation = list(IGHV=c(U=0, M=1),
                   Methylation_Cluster=c(`LP-CLL`=0, `IP-CLL`=1, `HP-CLL`=2))

make_pd <- function(cn, ...) {
  df <- function(...) data.frame(..., check.names = FALSE)
  
  x <- lpdAll[, cn]
  pd <- df(
    t(Biobase::exprs(x)[    c("del17p13", "TP53", "trisomy12"), , drop = FALSE]) %>% 
      `colnames<-`(c("del 17p13", "TP53", "trisomy 12")))
  
  # pd <- df(pd,
  #    t(Biobase::exprs(x)[  c("SF3B1", "del11q22.3", "del13q14_any"),, drop = FALSE]) %>% 
  #    `colnames<-`(c("SF3B1", "del11q22.3", "del13q14"))) 
  
  pd <- df(pd,
      cbind(as.integer(Biobase::exprs(x)["KRAS",] | Biobase::exprs(x)["NRAS",])) %>% 
      `colnames<-`("KRAS | NRAS")) 
  
  pd <- df(pd,
    # IGHV = Biobase::exprs(x)["IGHV Uppsala U/M", ],
    `IGHV (%)` = cut(x[["IGHV Uppsala % SHM"]],
                     breaks = c(0, seq(92, 100, by = 2), Inf), right = FALSE),
    `Meth. Cluster` = names(translation$Methylation_Cluster)[
          someMatch(paste(Biobase::exprs(x)["Methylation_Cluster", ]),
                    translation$Methylation_Cluster)],
    `Gene usage` = x$`IGHV gene usage`)
  
  if(length(unique(x$Diagnosis)) > 1)  
    pd <- df(pd, Diagnosis = x$Diagnosis)
  
  pd <- df(pd, 
    pretreated  = ifelse(patmeta[colnames(x),"IC50beforeTreatment"],"no","yes"), 
    alive       = ifelse(patmeta[colnames(x),"died"]>0, "no", "yes"),
    sex         = factor(x$Gender))
  
  rownames(pd) <- colnames(Biobase::exprs(x))
  
  for (i in setdiff(colnames(pd), "BCR score")) {
    if (!is.factor(pd[[i]]))
      pd[[i]] <- factor(pd[[i]])
    if (any(is.na(pd[[i]]))) {
      levels(pd[[i]]) <- c(levels(pd[[i]]), "n.d.")
      pd[[i]][ is.na(pd[[i]]) ] <- "n.d."   
    }
  }
  
  pd
}

8.3.3 Define the annotation colors

gucol <- rev(brewer.pal(nlevels(lpdAll$`IGHV gene usage`), "Set3")) %>% 
             `names<-`(sort(levels(lpdAll$`IGHV gene usage`)))
gucol["IGHV3-21"] <- "#E41A1C"

make_ann_colors <- function(pd) {
  bw <- c(`TRUE` = "darkblue", `FALSE` = "#ffffff")
  res <- list(
  Btk = bw, Syk = bw, PI3K = bw, MEK = bw)
  
  if ("exptbatch" %in% colnames(pd))
    res$exptbatch <- brewer.pal(nlevels(pd$exptbatch), "Set2") %>%
    `names<-`(levels(pd$exptbatch))
  
  if ("IGHV (%)" %in% colnames(pd))
    res$`IGHV (%)` <- 
       c(rev(colorRampPalette(
         brewer.pal(9, "Blues"))(nlevels(pd$`IGHV (%)`)-1)), "white") %>% 
                        `names<-`(levels(pd$`IGHV (%)`))
  
  if ("CD38" %in% colnames(pd))
    res$CD38 <- colorRampPalette(
      c("blue", "yellow"))(nlevels(pd$CD38)) %>% `names<-`(levels(pd$CD38))
  
  if("Gene usage" %in% colnames(pd)) 
    res$`Gene usage` <- gucol 

  if("Meth. Cluster" %in% colnames(pd)) 
    res$`Meth. Cluster` <- brewer.pal(9, "Blues")[c(1, 5, 9)] %>% 
        `names<-`(names(translation$Methylation_Cluster))
  
  res <- c(res, BloodCancerMultiOmics2017:::sampleColors)   # from addons.R  
  
  if("Diagnosis" %in% colnames(pd)) 
    res$Diagnosis <- BloodCancerMultiOmics2017:::colDiagS[
      names(BloodCancerMultiOmics2017:::colDiagS) %in% levels(pd$Diagnosis) ] %>% 
         (function(x) x[order(names(x))])
  
  for(i in names(res)) {
    whnd <- which(names(res[[i]]) == "n.d.")
    if(length(whnd)==1)
      res[[i]][whnd] <- "#e0e0e0" else 
        res[[i]] <- c(res[[i]], `n.d.` = "#e0e0e0")
    stopifnot(all(pd[[i]] %in% names(res[[i]])))
  }
  
  res
}

8.3.4 Heatmap drawing function

theatmap <- function(x, cellwidth = 7, cellheight = 11) {
  
  stopifnot(is.matrix(x))
  patDat <- make_pd(colnames(x))
  
  bpp <- brewer.pal(9, "Set1")
  breaks <- 2.3 * c(seq(-1, 1, length.out = 101)) %>% `names<-`(
      colorRampPalette(c(rev(brewer.pal(7, "YlOrRd")), 
                         "white", "white", "white", 
                         brewer.pal(7, "Blues")))(101))

  if (!is.null(attr(x, "colgap")))
    stopifnot(last(attr(x, "colgap")) == ncol(x))

  pheatmapwh(deckel(x, lower = first(breaks), upper = last(breaks)), 
         cluster_rows = FALSE,  
         cluster_cols = FALSE,
         gaps_col = attr(x, "colgap"),
         gaps_row = attr(x, "rowgap"),
         scale = "none",
         annotation_col = patDat,
         annotation_colors = make_ann_colors(patDat),
         color = names(breaks),
         breaks = breaks, 
         show_rownames = TRUE, 
         show_colnames = !TRUE, 
         cellwidth = cellwidth, cellheight = cellheight, 
         fontsize = 10, fontsize_row = 11, fontsize_col = 8,
         annotation_legend = TRUE, drop_levels = TRUE)
  }

8.4.1 Identify places where gaps between the rows should be

setGapPositions <- function(x, gapAt) {
  rg <- if (missing(gapAt)) c(0) else {
    s <- strsplit(gapAt, split = "--")
    stopifnot(all(listLen(s) == 2L))
    s <- strsplit(unlist(s), split = ":")
    spname <- drugs[safeMatch(sapply(s, `[`, 1), drugs$name), "id"]
    spconc <- as.numeric(sapply(s, `[`, 2))
    spi <- mapply(function(d, cc) {
      i <- which(cc == conctab[d, ])
      stopifnot(length(i) == 1)
      i
     }, spname, spconc)
    spdrug <- paste(spname, spi, sep = "_")
    mt <- safeMatch(spdrug, rownames(x))
    igp <- seq(1, length(mt), by = 2L)
    stopifnot(all( mt[igp] - mt[igp + 1] == 1))
    #stopifnot(all( mt[igp] - mt[igp + 1] == 1))
    c(mt[igp + 1], 0)
  }
  attr(x, "rowgap") <- rg
  x
}

clscd1 %<>% setGapPositions(gapAt = c(
  "PF 477736:10--idelalisib:10",
  "spebrutinib:2.5--PF 477736:2.5",
  "PRT062607 HCl:10--selumetinib:2.5",
  "selumetinib:10--MK-2206:2.5",
  "MK-2206:0.156--tipifarnib:10",
  "AT13387:0.039--encorafenib:10",
  "encorafenib:2.5--SD07:1.111",
  "doxorubicine:0.016--encorafenib:0.625",
  "encorafenib:0.156--rotenone:2.5",
  "SCH 900776:0.625--everolimus:0.625",
  "everolimus:10--afatinib:1.667",
  "arsenic trioxide:1--thapsigargin:5",
  "thapsigargin:0.313--fludarabine:0.156"
))

clscd2 %<>% setGapPositions(gapAt = c(
  "AT13387:0.039--everolimus:0.156",
  "everolimus:0.625--nutlin-3:10", 
  "fludarabine:10--thapsigargin:0.078",
  "thapsigargin:0.313--encorafenib:0.625", 
  "encorafenib:0.156--ruxolitinib:0.156"
))

#FIG# S8
rownames(clscd1) <- with(fData(lpdAll)[ rownames(clscd1),, drop = FALSE], 
  paste0(drugs[id, "name"], " ", conctab[cbind(id, paste0("c", subtype))], "uM")) 
rownames(clscd1)

##   [1] "vorinostat 0.313uM"    "BAY 11-7085 10uM"      "fludarabine 0.156uM"  
##   [4] "thapsigargin 0.313uM"  "thapsigargin 1.25uM"   "thapsigargin 5uM"     
##   [7] "arsenic trioxide 1uM"  "nutlin-3 0.156uM"      "chaetoglobosin A 5uM" 
##  [10] "nutlin-3 0.625uM"      "nutlin-3 2.5uM"        "vorinostat 1.25uM"    
##  [13] "fludarabine 2.5uM"     "fludarabine 0.625uM"   "tofacitinib 10uM"     
##  [16] "rigosertib 10uM"       "KX2-391 0.625uM"       "KX2-391 2.5uM"        
##  [19] "afatinib 1.667uM"      "everolimus 10uM"       "everolimus 0.156uM"   
##  [22] "everolimus 0.625uM"    "SCH 900776 0.625uM"    "SCH 900776 2.5uM"     
##  [25] "silmitasertib 10uM"    "chaetocin 0.031uM"     "navitoclax 0.004uM"   
##  [28] "orlistat 2.5uM"        "orlistat 10uM"         "venetoclax 0.004uM"   
##  [31] "navitoclax 0.016uM"    "thapsigargin 0.078uM"  "rotenone 0.156uM"     
##  [34] "rotenone 0.625uM"      "rotenone 2.5uM"        "encorafenib 0.156uM"  
##  [37] "encorafenib 0.625uM"   "doxorubicine 0.016uM"  "doxorubicine 0.25uM"  
##  [40] "doxorubicine 0.063uM"  "YM155 0.008uM"         "SD51 0.37uM"          
##  [43] "SD51 1.111uM"          "MIS-43 0.37uM"         "MIS-43 1.111uM"       
##  [46] "SD07 1.111uM"          "encorafenib 2.5uM"     "encorafenib 10uM"     
##  [49] "AT13387 0.039uM"       "tipifarnib 2.5uM"      "tipifarnib 10uM"      
##  [52] "MK-2206 0.156uM"       "MK-2206 0.625uM"       "MK-2206 2.5uM"        
##  [55] "selumetinib 10uM"      "selumetinib 0.156uM"   "selumetinib 0.625uM"  
##  [58] "selumetinib 2.5uM"     "PRT062607 HCl 10uM"    "BIX02188 10uM"        
##  [61] "rabusertib 10uM"       "enzastaurin 10uM"      "VE-821 10uM"          
##  [64] "sotrastaurin 2.5uM"    "sotrastaurin 10uM"     "saracatinib 2.5uM"    
##  [67] "saracatinib 10uM"      "MK-2206 10uM"          "TAE684 2.5uM"         
##  [70] "sunitinib 10uM"        "CCT241533 2.5uM"       "SGI-1776 2.5uM"       
##  [73] "SGX-523 10uM"          "KU-60019 10uM"         "NU7441 10uM"          
##  [76] "ibrutinib 10uM"        "KU-60019 2.5uM"        "NU7441 2.5uM"         
##  [79] "PF 477736 0.625uM"     "PF 477736 2.5uM"       "spebrutinib 2.5uM"    
##  [82] "spebrutinib 10uM"      "PRT062607 HCl 0.156uM" "spebrutinib 0.156uM"  
##  [85] "spebrutinib 0.625uM"   "MK-1775 2.5uM"         "idelalisib 0.156uM"   
##  [88] "idelalisib 0.625uM"    "idelalisib 2.5uM"      "ibrutinib 0.156uM"    
##  [91] "ibrutinib 0.625uM"     "ibrutinib 2.5uM"       "tamatinib 2.5uM"      
##  [94] "tamatinib 10uM"        "MK-1775 10uM"          "AZD7762 0.156uM"      
##  [97] "AZD7762 0.625uM"       "BX912 2.5uM"           "dasatinib 0.156uM"    
## [100] "PRT062607 HCl 0.625uM" "PRT062607 HCl 2.5uM"   "duvelisib 0.156uM"    
## [103] "duvelisib 0.625uM"     "duvelisib 2.5uM"       "duvelisib 10uM"       
## [106] "idelalisib 10uM"       "PF 477736 10uM"        "AZD7762 2.5uM"        
## [109] "AZD7762 10uM"          "BX912 10uM"            "dasatinib 10uM"       
## [112] "dasatinib 0.625uM"     "dasatinib 2.5uM"       "AT13387 0.156uM"      
## [115] "AT13387 0.625uM"       "AT13387 2.5uM"

theatmap(clscd1)

#FIG# 3A
rownames(clscd2) <- with(fData(lpdAll)[ rownames(clscd2),, drop = FALSE], 
  paste0(drugs[id, "name"], " ", conctab[cbind(id, paste0("c", subtype))], "uM")) 

rownames(clscd2)

##  [1] "venetoclax 0.016uM"   "venetoclax 0.004uM"   "YM155 0.031uM"       
##  [4] "YM155 0.008uM"        "ruxolitinib 0.625uM"  "ruxolitinib 0.156uM" 
##  [7] "encorafenib 0.156uM"  "encorafenib 0.625uM"  "thapsigargin 0.313uM"
## [10] "thapsigargin 0.078uM" "fludarabine 10uM"     "fludarabine 2.5uM"   
## [13] "nutlin-3 2.5uM"       "nutlin-3 10uM"        "everolimus 0.625uM"  
## [16] "everolimus 0.156uM"   "AT13387 0.039uM"      "AT13387 0.156uM"     
## [19] "PF 477736 0.156uM"    "PF 477736 0.625uM"    "MK-2206 0.156uM"     
## [22] "MK-2206 0.625uM"      "selumetinib 0.156uM"  "selumetinib 0.625uM" 
## [25] "AZD7762 0.625uM"      "AZD7762 0.156uM"      "ibrutinib 0.156uM"   
## [28] "ibrutinib 0.625uM"    "duvelisib 0.156uM"    "duvelisib 0.625uM"   
## [31] "idelalisib 0.625uM"   "idelalisib 0.156uM"   "tamatinib 0.156uM"   
## [34] "tamatinib 0.625uM"

theatmap(clscd2)

options(stringsAsFactors=TRUE)

9.3.1 BCR drugs

#FIG# 3B
makeTernary(lpdCLL, c("PI3K", "BTK", "SYK"), ighv=NA)

## 
## FALSE  TRUE 
##    52   131

#FIG# 3B
makeTernary(lpdCLL, c("PI3K", "BTK", "SYK"), ighv="M")

## 
## FALSE  TRUE 
##    14    83

#FIG# 3B
makeTernary(lpdCLL, c("PI3K", "BTK", "SYK"), ighv="U")

## 
## FALSE  TRUE 
##    33    41

9.3.2 BTK & MEK & MTOR

#FIG# 3BC
makeTernary(lpdCLL, c("MTOR", "BTK", "MEK"), ighv=NA)

## 
## FALSE  TRUE 
##    52   131

#FIG# 3BC
makeTernary(lpdCLL, c("MTOR", "BTK", "MEK"), ighv="M")

## 
## FALSE  TRUE 
##    14    83

#FIG# 3BC
makeTernary(lpdCLL, c("MTOR", "BTK", "MEK"), ighv="U")

## 
## FALSE  TRUE 
##    33    41

9.3.3 PI3K & MEK & MTOR

All CLL samples included.

#FIG# S9 left
makeTernary(lpdCLL, c("MTOR", "PI3K", "MEK"), ighv=NA)

## 
## FALSE  TRUE 
##    52   131

9.3.4 SYK & MEK & MTOR

All CLL samples included.

#FIG# S9 right
makeTernary(lpdCLL, c("MTOR", "SYK", "MEK"), ighv=NA)

## 
## FALSE  TRUE 
##    52   131

16.3.1 Functions for enrichment analysis and plots

Define cut-off.

pCut = 0.05

Function to run GSEA or PAGE in R.

runGSEA <- function(inputTab, gmtFile, GSAmethod="gsea", nPerm=1000){
    inGMT <- loadGSC(gmtFile,type="gmt")
     #re-rank by score
    rankTab <- inputTab[order(inputTab[,1],decreasing = TRUE),,drop=FALSE] 
    if (GSAmethod == "gsea"){
        #readin geneset database
        #GSEA analysis
        res <- runGSA(geneLevelStats = rankTab,
                      geneSetStat = GSAmethod,
                      adjMethod = "fdr", gsc=inGMT,
                      signifMethod = 'geneSampling', nPerm = nPerm)
        GSAsummaryTable(res)
    } else if (GSAmethod == "page"){
        res <- runGSA(geneLevelStats = rankTab,
                      geneSetStat = GSAmethod,
                      adjMethod = "fdr", gsc=inGMT,
                      signifMethod = 'nullDist')
        GSAsummaryTable(res)
    }
}

Function which run the GSE for each response group.

runGSE = function(gmt) {
  
  Res <- list()
  for (i in names(DEres)) {
    dataTab <- data.frame(DEres[[i]])
    dataTab$ID <- rownames(dataTab)
    #filter using pvalues
    dataTab <- filter(dataTab, pvalue <= pCut) %>%
                  arrange(pvalue) %>% 
                  mutate(Symbol = rowData(ddsCLL[ID,])$symbol)
    dataTab <- dataTab[!duplicated(dataTab$Symbol),]
    statTab <- data.frame(row.names = dataTab$Symbol, stat = dataTab$stat)
    resTab <- runGSEA(inputTab=statTab, gmtFile=gmt, GSAmethod="page")
    Res[[i]] <- arrange(resTab,desc(`Stat (dist.dir)`))
  }
  Res
}

Function to get the list of genes enriched in a set.

getGenes <- function(inputTab, gmtFile){
  geneList <- loadGSC(gmtFile,type="gmt")$gsc
  enrichedUp <- lapply(geneList, function(x) 
    intersect(rownames(inputTab[inputTab[,1] >0,,drop=FALSE]),x))
  enrichedDown <- lapply(geneList, function(x)
    intersect(rownames(inputTab[inputTab[,1] <0,,drop=FALSE]),x))
  return(list(up=enrichedUp, down=enrichedDown))
}

A function to plot the heat map of intersection of genes in different gene sets.

plotSetHeatmap <- 
  function(geneTab, enrichTab, topN, gmtFile, tittle="",
           asterixList = NULL, anno=FALSE) {

    if (nrow(enrichTab) < topN) topN <- nrow(enrichTab)
    enrichTab <- enrichTab[seq(1,topN),]

    geneList <- getGenes(geneTab,gmtFile)
    
    geneList$up <- geneList$up[enrichTab[,1]]
    geneList$down <- geneList$down[enrichTab[,1]]
    
    #form a table 
    allGenes <- unique(c(unlist(geneList$up),unlist(geneList$down)))
    allSets <- unique(c(names(geneList$up),names(geneList$down)))
    plotTable <- matrix(data=NA,ncol = length(allSets),
                        nrow = length(allGenes),
                        dimnames = list(allGenes,allSets))
    for (setName in names(geneList$up)) {
      plotTable[geneList$up[[setName]],setName] <- 1
    }
    for (setName in names(geneList$down)) {
      plotTable[geneList$down[[setName]],setName] <- -1
    }

    if(is.null(asterixList)) {
      #if no correlation table specified, order by the number of
      # significant gene
      geneOrder <- rev(
        rownames(plotTable[order(rowSums(plotTable, na.rm = TRUE),
                                 decreasing =FALSE),]))
    } else {
      #otherwise, order by the p value of correlation
      asterixList <- arrange(asterixList, p)
      geneOrder <- filter(
        asterixList, symbol %in% rownames(plotTable))$symbol
      geneOrder <- c(
        geneOrder, rownames(plotTable)[! rownames(plotTable) %in% geneOrder])
    }
    
    plotTable <- melt(plotTable)
    colnames(plotTable) <- c("gene","set","value")
    plotTable$gene <- as.character(plotTable$gene)
    
    if(!is.null(asterixList)) {
      #add + if gene is positivily correlated with sensitivity, else add "-"
      plotTable$ifSig <- asterixList[
        match(plotTable$gene, asterixList$symbol),]$coef
      plotTable <- mutate(plotTable, ifSig =
                            ifelse(is.na(ifSig) | is.na(value), "",
                                   ifelse(ifSig > 0, "-", "+")))
    }
    plotTable$value <- replace(plotTable$value,
                               plotTable$value %in% c(1), "Up")
    plotTable$value <- replace(plotTable$value,
                               plotTable$value %in%  c(-1), "Down")
    
    allSymbols <- plotTable$gene
    
    geneSymbol <- geneOrder
    
    if (anno) { #if add functional annotations in addition to gene names
      annoTab <- tibble(symbol = rowData(ddsCLL)$symbol, 
                        anno = sapply(rowData(ddsCLL)$description,
                                      function(x) unlist(strsplit(x,"[[]"))[1]))
      annoTab <- annoTab[!duplicated(annoTab$symbol),]
      annoTab$combine <- sprintf("%s (%s)",annoTab$symbol, annoTab$anno)
      plotTable$gene <- annoTab[match(plotTable$gene,annoTab$symbol),]$combine
      geneOrder <- annoTab[match(geneOrder,annoTab$symbol),]$combine
      geneOrder <- rev(geneOrder)
    }
    
    plotTable$gene <- factor(plotTable$gene, levels =geneOrder)
    plotTable$set <- factor(plotTable$set, levels = enrichTab[,1])

    g <- ggplot(plotTable, aes(x=set, y = gene)) +
      geom_tile(aes(fill=value), color = "black") +
      scale_fill_manual(values = c("Up"="red","Down"="blue")) +
      xlab("") + ylab("") + theme_classic() +
      theme(axis.text.x=element_text(size=7, angle = 60, hjust = 0),
            axis.text.y=element_text(size=7),
            axis.ticks = element_line(color="white"),
            axis.line = element_line(color="white"),
            legend.position = "none") +
      scale_x_discrete(position = "top") +
      scale_y_discrete(position = "right")
    
    if(!is.null(asterixList)) {
      g <- g + geom_text(aes(label = ifSig), vjust =0.40)
    }
    
    # construct the gtable
    wdths = c(0.05, 0.25*length(levels(plotTable$set)), 5)
    hghts = c(2.8, 0.1*length(levels(plotTable$gene)), 0.05)
    gt = gtable(widths=unit(wdths, "in"), heights=unit(hghts, "in"))
    ## make grobs
    ggr = ggplotGrob(g)
    ## fill in the gtable
    gt = gtable_add_grob(gt, gtable_filter(ggr, "panel"), 2, 2)
    gt = gtable_add_grob(gt, ggr$grobs[[5]], 1, 2) # top axis
    gt = gtable_add_grob(gt, ggr$grobs[[9]], 2, 3) # right axis
    
    return(list(list(plot=gt,
                     width=sum(wdths),
                     height=sum(hghts),
                     genes=geneSymbol)))
}    

Prepare stats per gene for plotting.

statTab = setNames(lapply(c("mTORnone","BTKnone","MEKnone"), function(gr) {
  dataTab <- data.frame(DEres[[gr]])
  dataTab$ID <- rownames(dataTab)
  #filter using pvalues
  dataTab <- filter(dataTab, pvalue <= pCut) %>%
    arrange(pvalue) %>%
    mutate(Symbol = rowData(ddsCLL[ID,])$symbol) %>%
    filter(log2FoldChange > 0)
  dataTab <- dataTab[!duplicated(dataTab$Symbol),]
  data.frame(row.names = dataTab$Symbol, stat = dataTab$stat)
}), nm=c("mTORnone","BTKnone","MEKnone"))

16.3.2 Geneset enrichment based on Hallmark set (H)

Perform enrichment analysis using PAGE method.

hallmarkRes = runGSE(gmt=gmts[["H"]])

## Checking arguments...done!
## Calculating gene set statistics...done!
## Calculating gene set significance...done!
## Adjusting for multiple testing...done!
## Checking arguments...done!
## Calculating gene set statistics...done!
## Calculating gene set significance...done!
## Adjusting for multiple testing...done!
## Checking arguments...done!
## Calculating gene set statistics...done!
## Calculating gene set significance...done!
## Adjusting for multiple testing...done!

16.3.3 Geneset enrichment based on oncogenic signature set (C6)

Perform enrichment analysis using PAGE method.

c6Res = runGSE(gmt=gmts[["C6"]])

## Checking arguments...done!
## Calculating gene set statistics...done!
## Calculating gene set significance...done!
## Adjusting for multiple testing...done!
## Checking arguments...done!
## Calculating gene set statistics...done!
## Calculating gene set significance...done!
## Adjusting for multiple testing...done!
## Checking arguments...done!
## Calculating gene set statistics...done!
## Calculating gene set significance...done!
## Adjusting for multiple testing...done!

16.4.1 Correlation test

Prepare drug sensitivity vector and gene expression matrix

ddsCLL.mTOR <- ddsCLL.norm[,ddsCLL.norm$group %in% "mTOR"]
viabMTOR <- Biobase::exprs(lpdCLL["D_063_4:5", ddsCLL.mTOR$PatID])[1,]
stopifnot(all(ddsCLL.mTOR$PatID == colnames(viabMTOR)))  

Filtering and applying variance stabilizing transformation on RNAseq data

#only keep genes that have counts higher than 10 in any sample
keep <- apply(assay(ddsCLL.mTOR), 1, function(x) any(x >= 10)) 
ddsCLL.mTOR <- ddsCLL.mTOR[keep,]
dim(ddsCLL.mTOR)

## [1] 9775   15

Association test using Pearson correlation

corResult <- tibble(ID = rownames(ddsCLL.mTOR), 
                    symbol = rowData(ddsCLL.mTOR)$symbol,
                    coef = NA,
                    p = NA)

for (i in seq(nrow(corResult))) {
  res <- cor.test(viabMTOR, assay(ddsCLL.mTOR[i,])[1,], method = "pearson")
  corResult[i,c(3,4)] <- c(res$estimate[[1]], res$p.value)
}

corResult <- arrange(corResult, p) %>% mutate(p.adj = p.adjust(p, method="BH"))

16.4.2 Enrichment heatmaps for mTOR group with overlapped genes indicated

The genes that are positively correlated with everolimus sensitivity are labeled as “+” in the heatmap and the negatively correlated genes are labeled as “-”.

Plot for C6 gene sets

pCut = 0.05
corResult.sig <- filter(corResult, p <= pCut)
c6Plot <- plotSetHeatmap(geneTab=statTab[["mTORnone"]],
                enrichTab=c6Res[["mTORnone"]],
                topN=5, gmtFile=gmts[["C6"]],
                #add asterix in front of the overlapped genes
                asterixList = corResult.sig, 
                anno=TRUE, i)

Plot for Hallmark gene sets

hallmarkPlot <- plotSetHeatmap(geneTab=statTab[["mTORnone"]],
                enrichTab=hallmarkRes[["mTORnone"]],
                topN=5, gmtFile=gmts[["H"]],
                asterixList = corResult.sig,
                anno=TRUE, i)

16.5.1 Correlation test

Prepare drug sensitivity vector and gene expression matrix

ddsCLL.BTK <- ddsCLL.norm[,ddsCLL.norm$group %in% "BTK"]
viabBTK <- Biobase::exprs(lpdCLL["D_002_4:5", ddsCLL.BTK$PatID])[1,]
stopifnot(all(ddsCLL.BTK$PatID == colnames(viabBTK)))  

Filtering and applying variance stabilizing transformation on RNAseq data

#only keep genes that have counts higher than 10 in any sample
keep <- apply(assay(ddsCLL.BTK), 1, function(x) any(x >= 10)) 
ddsCLL.BTK <- ddsCLL.BTK[keep,]
dim(ddsCLL.BTK)

## [1] 10416    30

Association test using Pearson correlation

corResult <- tibble(ID = rownames(ddsCLL.BTK), 
                    symbol = rowData(ddsCLL.BTK)$symbol,
                    coef = NA,
                    p = NA)

for (i in seq(nrow(corResult))) {
  res <- cor.test(viabBTK, assay(ddsCLL.BTK[i,])[1,], method = "pearson")
  corResult[i,c(3,4)] <- c(res$estimate[[1]], res$p.value)
}

corResult <- arrange(corResult, p) %>% mutate(p.adj = p.adjust(p, method="BH"))

16.5.2 Enrichment heatmaps for BTK group with overlapped genes indicated

Plot for C6 gene sets

pCut = 0.05
corResult.sig <- filter(corResult, p <= pCut)
c6Plot <- plotSetHeatmap(geneTab=statTab[["BTKnone"]],
                enrichTab=c6Res[["BTKnone"]],
                topN=5, gmtFile=gmts[["C6"]],
                #add asterix in front of the overlapped genes
                asterixList = corResult.sig, 
                anno=TRUE, i)

Plot for Hallmark gene sets

hallmarkPlot <- plotSetHeatmap(geneTab=statTab[["BTKnone"]],
                enrichTab=hallmarkRes[["BTKnone"]],
                topN=5, gmtFile=gmts[["H"]],
                asterixList = corResult.sig,
                anno=TRUE, i)

16.5.3 Selumetinib response VS gene expression (within MEK group)

Correlation test was performed to identify genes whose expressions are correlated with the sensitivity to the MEK inhibitor (selumetinib) sensitivity within the MEK group.

16.5.4 Correlation test

Prepare drug sensitivity vector and gene expression matrix

ddsCLL.MEK <- ddsCLL.norm[,ddsCLL.norm$group %in% "MEK"]
viabMEK <- Biobase::exprs(lpdCLL["D_012_4:5", ddsCLL.MEK$PatID])[1,]
stopifnot(all(ddsCLL.MEK$PatID == colnames(viabMEK)))  

Filtering and applying variance stabilizing transformation on RNAseq data

#only keep genes that have counts higher than 10 in any sample
keep <- apply(assay(ddsCLL.MEK), 1, function(x) any(x >= 10)) 
ddsCLL.MEK <- ddsCLL.MEK[keep,]
dim(ddsCLL.MEK)

## [1] 10174    18

Association test using Pearson correlation

corResult <- tibble(ID = rownames(ddsCLL.MEK), 
                    symbol = rowData(ddsCLL.MEK)$symbol,
                    coef = NA,
                    p = NA)

for (i in seq(nrow(corResult))) {
  res <- cor.test(viabMEK, assay(ddsCLL.MEK[i,])[1,], method = "pearson")
  corResult[i,c(3,4)] <- c(res$estimate[[1]], res$p.value)
}

corResult <- arrange(corResult, p) %>% mutate(p.adj = p.adjust(p, method="BH"))

Within MEK group, no gene expression was correlated with Selumetinib response

16.5.5 Enrichment heatmaps for MEK group

Plot for C6 gene sets

pCut = 0.05
corResult.sig <- filter(corResult, p <= pCut)
c6Plot <- plotSetHeatmap(geneTab=statTab[["MEKnone"]],
                enrichTab=c6Res[["MEKnone"]],
                topN=5, gmtFile=gmts[["C6"]],
                anno=TRUE, i)

Plot for Hallmark gene sets

hallmarkPlot <- plotSetHeatmap(geneTab=statTab[["MEKnone"]],
                enrichTab=hallmarkRes[["MEKnone"]],
                topN=5, gmtFile=gmts[["H"]],
                asterixList = corResult.sig,
                anno=TRUE, i)

options(stringsAsFactors=FALSE)

17.2.1 Prepare objects for testing

## VIABILITIES
## list of matrices; one matrix per screen/day
## each matrix contains all CLL patients
measurements=list()

### Main Screen
patM = colnames(drpar)[which(patmeta[colnames(drpar),"Diagnosis"]=="CLL")]
measurements[["main"]] =
  do.call(cbind,
          lapply(list("viaraw.1","viaraw.2","viaraw.3","viaraw.4","viaraw.5"),
                 function(viac) {
  tmp = t(assayData(drpar)[[viac]][,patM])
  colnames(tmp) = paste(colnames(tmp), conctab[colnames(tmp),
            paste0("c",substring(viac,8,8))], sep="-")
  tmp
}))

pats = sort(unique(patM))

## TESTING FACTORS
testingFactors = list()
# ighv
ighv = setNames(patmeta[pats, "IGHV"], nm=pats)
# mutations
tmp = cbind(IGHV=ifelse(ighv=="U",1,0), assayData(mutCOM)$binary[pats,])
testingFactors[["mutation"]] = tmp[,-grep("Chromothripsis", colnames(tmp))]

# BATCHES
j = which(pData(drpar)[patM, "ExpDate"] < as.Date("2014-01-01"))
k = which(pData(drpar)[patM, "ExpDate"] < as.Date("2014-08-01") &
            pData(drpar)[patM, "ExpDate"] > as.Date("2014-01-01"))
l = which(pData(drpar)[patM, "ExpDate"] > as.Date("2014-08-01"))

measurements[["main"]] = measurements[["main"]][c(patM[j], patM[k], patM[l]),]
batchvec = factor(
  setNames(c(rep(1, length(j)), rep(2, length(k)), rep(3, length(l))),
           nm=c(patM[j], patM[k], patM[l])))

# LABELS FOR GROUPING
beelabs = t(sapply(colnames(testingFactors[["mutation"]]), function(fac) {
  if(fac=="IGHV")
    c(`0`="IGHV mut", `1`="IGHV unmut")
  else if(grepl("[[:upper:]]",fac)) # if all letters are uppercase
    c(`0`=paste(fac, "wt"),`1`=paste(fac, "mt"))
  else
    c(`0`="wt",`1`=fac)
}))

17.2.2 Assesment of importance of batch effect

We first used the approach explained in the introduction section to test for associations between drug viability assay results and genomic features, which comprised: somatic mutations (aggregated at the gene level), copy number aberrations and IGHV status.

allresultsT = testFactors(msrmnt=measurements[["main"]],
                          factors=testingFactors[["mutation"]],
                          test="student", batch=NA)

resultsT = allresultsT$test
resultsT$adj.pval = p.adjust(resultsT$pval, method="BH")

However, we ware aware that the main screen was performed in three groups of batches over a time period of 1.5 years; these comprise, respectively, the samples screened in 2013, in 2014 before August and in 2014 in August and September. Therefore, to control for confounding by the different batch groups we repeated the drug-feature association tests using batch group as a blocking factor and a two-way ANOVA test.

allresultsA = testFactors(msrmnt=measurements[["main"]],
                          factors=testingFactors[["mutation"]],
                          test="anova", batch=batchvec)

resultsA = allresultsA$test
resultsA$adj.pval = p.adjust(resultsA$pval, method="BH")

We then compared the p-values from both tests.

Only one drug, bortezomib, showed discrepant p-values, and exploration of its data suggested that it lost its activity during storage. The data for this drug and NSC 74859 were discarded from further analysis.

badrugs = c("D_008", "D_025")

measurements = lapply(measurements, function(drres) {
  drres[,-grep(paste(badrugs, collapse="|"), colnames(drres))]
})

For all remaining associations, testing with and without batch as a blocking factor yielded equivalent results. Therefore, all reported p-values for associations come from the t-tests without using blocking for batch effects.

17.2.3 Associations of drug response with mutations in CLL

We tested for associations between drug viability assay results and genomic features (43 features for the pilot screen and 63 for the main screen). p-values were adjusted for multiple testing by applying the Benjamini-Hochberg procedure, separately for the main screen and for each of the two incubation times of the pilot screen.

allresults1 = lapply(measurements, function(measurement) {
  testFactors(msrmnt=measurement, factors=testingFactors[["mutation"]],
              test="student", batch=NA)
})

effects1 = lapply(allresults1, function(res) res[["effect"]])
results1 = lapply(allresults1, function(res) res[["test"]])

results1 = lapply(results1, function(res) {
  res$adj.pval = p.adjust(res$pval, method="BH")
  res
})

measurements1 = measurements
testingFactors1 = testingFactors
beelabs1 = beelabs

17.2.4 Volcano plots: summary of the results

In this section we summarize all significant associations for a given mutation in a form of volcano plots. The pink color spectrum indicates a resistant phenotype and the blue color spectrum a sensitive phenotype in the presence of the tested mutation. FDR of 10 % was used.

IGHV.

Trisomy 12.

19.4.1 Main Figure

Function for generating the figure.

plot.pvalues <- function(pvals, dr.order, folge, colors, shapes) {
  g = ggplot(data = filter(pvals, drug %in% dr.order$drug)) +
  geom_point(aes(x = factor(drug, levels = dr.order$drug), y = -log10(p.value), 
                 colour = factor(mutation, levels(factor(mutation))[folge]),
                 shape  = factor(mutation, levels(factor(mutation))[folge])),
             size=5, show.legend = T)  + 
  scale_color_manual(name = "Mutations",
                     values = colors,
                     labels = lev[folge]) + 
  scale_shape_manual(name = "Mutations",
                     values = shapes,
                     labels = lev[folge]) + t1 + 
  labs(x = "", y = expression(paste(-log[10], "p")), title = "") +
    scale_y_continuous(expression(italic(p)*"-value"),
                       breaks=seq(0,10,5),
                       labels=math_format(expr=10^.x)(-seq(0,10,5))) 
  g
}

#FIG# 4A
## plot the p-values 
g = plot.pvalues(pvals, dr.order, folge, 
                 colors, shapes = c(rep(16,13), c(1,1)))

## add FDR threshold
g = g + geom_hline(yintercept = -log10(fdr),
                   linetype="dashed", size=0.3)

g = g + 
  annotation_custom(grob = textGrob(label = paste0("FDR", fd.thresh, "%"), 
                                    hjust = 1, vjust = 1, 
                                    gp = gpar(cex = 0.5,
                                              fontface = "bold",
                                              fontsize = 25)),
                    ymin = -log10(fdr) - 0.2, 
                    ymax = -log10(fdr) + 0.5, 
                    xmin = -1.3, xmax = 1.5) + 
  theme(legend.position = "none")

# generate pathway/target annotations for certain drug classes
#g = make.annot(g, dr.order)

# legend guide
leg.guides <- guides(colour = guide_legend(ncol = 1, 
                                           byrow = TRUE,
                                           override.aes = list(size = 3),
                                           title = "Mutations",
                                           label.hjust = 0,
                                           keywidth = 0.4,
                                           keyheight = 0.8), 
                      shape = guide_legend(ncol = 1, 
                                           byrow = TRUE,
                                           title = "Mutations",
                                           label.hjust = 0,
                                           keywidth = 0.4,
                                           keyheight = 0.8))

# create a legend grob
legend = g_legend(g + t.leg +  leg.guides)

## arranget the main plot and the legend
# using grid graphics
gt <- ggplot_gtable(ggplot_build(g + theme(legend.position = 'none')))
gt$layout$clip[gt$layout$name == "panel"] <- "off"

grid.arrange(gt, legend,
             ncol=2, nrow=1, widths=c(0.92,0.08))

19.4.2 Supplementary Figure (incl. pretreatment)

In the supplementary figure we use pretreatment status as a blocking factor, i.e. we model drug sensitivity - gene variant associations as: lm(viability ~ mutation + pretreatment)

## lm(viab ~ mutation + pretreatment.status)
pvals = dr %>% group_by(drug, mutation) %>%
  do(tidy(lm(viab ~ mut.value + pretreat, data = .))) %>%
  filter(term == 'mut.value1') %>%
  dplyr::select(drug, mutation, p.value)

# compute the FDR threshold
fd.thresh = 10
padj = p.adjust(pvals$p.value, method = "BH")
fdr = max(pvals$p.value[which(padj <= fd.thresh/100)])


pvals = filter(pvals, mutation != 'IGHV')
pvals = pvals %>% ungroup() %>%
  mutate(mutation = ifelse(p.value > fdr,
                           paste0("Below ", fd.thresh,"% FDR"),
                           mutation)) %>%
  mutate(mutation = ifelse(!(mutation %in% select.mutations) &
                             !(mutation == paste0("Below ",
                                                  fd.thresh,"% FDR")), 
                           "Other", mutation)) %>%
  filter(drug != "bortezomib" & drug != "NSC 74859")


lev = get.labels(pvals)
folge = get.mutation.order(lev)

# get the drug order
dr.order = get.drug.order(pvals, drugs)

## Joining, by = "drug"
## Joining, by = "drug"

## Warning: Factor `target_category` contains implicit NA, consider using
## `forcats::fct_explicit_na`

## Joining, by = c("drug", "logsum", "n", "target_category")

mut.order = folge[!is.na(folge)]

After recomputing the \(p\)-values (using a linear model that accounts for pretreatment status), plot the figure:

#FIG# S19
## plot the p-values 
g = plot.pvalues(pvals, dr.order, mut.order, 
                 colors[which(!is.na(folge))], shapes = c(rep(16,9), c(1,1)))

## add FDR threshold
g = g + geom_hline(yintercept = -log10(fdr),
                   linetype="dashed", size=0.3)

g = g + 
  annotation_custom(grob = textGrob(label = paste0("FDR", fd.thresh, "%"), 
                                    hjust = 1, vjust = 1, 
                                    gp = gpar(cex = 0.5,
                                              fontface = "bold",
                                              fontsize = 25)),
                    ymin = -log10(fdr) - 0.2, 
                    ymax = -log10(fdr) + 0.5, 
                    xmin = -1.3, xmax = 1.5) + 
  theme(legend.position = "none")

# generate pathway/target annotations for certain drug classes
#g = make.annot(g, dr.order)

# legend guide
leg.guides <- guides(colour = guide_legend(ncol = 1, 
                                           byrow = TRUE,
                                           override.aes = list(size = 3),
                                           title = "Mutations",
                                           label.hjust = 0,
                                           keywidth = 0.4,
                                           keyheight = 0.8), 
                      shape = guide_legend(ncol = 1, 
                                           byrow = TRUE,
                                           title = "Mutations",
                                           label.hjust = 0,
                                           keywidth = 0.4,
                                           keyheight = 0.8))

# create a legend grob
legend = g_legend(g + t.leg +  leg.guides)

## arranget the main plot and the legend
# using grid graphics
gt <- ggplot_gtable(ggplot_build(g + theme(legend.position = 'none')))
gt$layout$clip[gt$layout$name == "panel"] <- "off"

grid.arrange(gt, legend,
             ncol=2, nrow=1, widths=c(0.92,0.08))

22.1.1 Heatmap plot of differentially expressed genes

A heat map plot was used to show the normalized expression value (Z-score) of the differentially expressed genes in samples with and without trisomy 12.

Prepare matrix for heat map plot.

#filter genes
fdrCut <- 0.1
fcCut <- 1.5

allDE <- data.frame(DEres.shr) %>%
  rownames_to_column(var = "ID") %>% 
  mutate(Symbol = rowData(ddsCLL[ID,])$symbol,
         Chr = rowData(ddsCLL[ID,])$chromosome) %>% 
  filter(padj <= fdrCut & abs(log2FoldChange) > fcCut) %>% 
  arrange(pvalue) %>% filter(!duplicated(Symbol)) %>%
  mutate(Chr12 = ifelse(Chr == 12, "yes", "no"))

#get the expression matrix
plotMat <- assay(ddsCLL.norm[allDE$ID,])
colnames(plotMat) <- ddsCLL.norm$PatID
rownames(plotMat) <- allDE$Symbol

#sort columns of plot matrix based on trisomy 12 status
plotMat <- plotMat[,order(ddsCLL.norm$trisomy12)]

#calculate z-score and scale
plotMat <- t(scale(t(plotMat)))
plotMat[plotMat >= 4] <- 4
plotMat[plotMat <= -4] <- -4

Plot the heat map.

#FIG# S23 B
#prepare colums and row annotations
annoCol <- data.frame(row.names=ddsCLL.norm$PatID, Tris12=ddsCLL.norm$trisomy12)
levels(annoCol$Tris12) <- list(wt = 0, mut =1)
annoRow <- data.frame(row.names = allDE$Symbol, Chr12 = allDE$Chr12)
annoColor <- list(Tris12 = c(wt = "grey80", mut = "black"),
                  Chr12 = c(yes="red", no = "grey80"))


pheatmap(plotMat,
         color=colorRampPalette(rev(brewer.pal(n=7, name="RdBu")))(100),
         cluster_cols = FALSE,
         annotation_row = annoRow, annotation_col = annoCol,
         show_colnames = FALSE, fontsize_row = 3,
         breaks = seq(-4,4, length.out = 101),
         annotation_colors = annoColor, border_color = NA)

According to the gene expression heat map, the samples with trisomy 12 show distinct expression pattern. 84 genes are significantly up-regulated in trisomy 12 samples and 37 genes are down-regulated in trisomy 12 samples (FDR =0.1 and log2FoldChange > 1.5). Among the 84 up-regulated genes, only 12 genes are from chromosome 12, suggested the distinct expression pattern of trisomy 12 samples can not be merely explained by gene dosage effect.

22.2.1 Perform enrichment analysis

Function to run PAGE in R.

runGSEA <- function(inputTab, gmtFile, GSAmethod="gsea", nPerm=1000){
    inGMT <- loadGSC(gmtFile,type="gmt")
    #re-rank by score
    rankTab <- inputTab[order(inputTab[,1],decreasing = TRUE),,drop=FALSE] 
    if (GSAmethod == "gsea"){
        #readin geneset database
        #GSEA analysis
        res <- runGSA(geneLevelStats = rankTab,
                      geneSetStat = GSAmethod,
                      adjMethod = "fdr", gsc=inGMT,
                      signifMethod = 'geneSampling', nPerm = nPerm)
        GSAsummaryTable(res)
    } else if (GSAmethod == "page"){
        res <- runGSA(geneLevelStats = rankTab,
                      geneSetStat = GSAmethod,
                      adjMethod = "fdr", gsc=inGMT,
                      signifMethod = 'nullDist')
        GSAsummaryTable(res)
    }
}

Function for plotting enrichment bar.

plotEnrichmentBar <- function(resTab, pCut=0.05, ifFDR=FALSE,
                              setName="Signatures") {
  pList <- list()
  rowNum <- c()
  for (i in names(resTab)) {
    plotTab <- resTab[[i]]
    if (ifFDR) {
      plotTab <- dplyr::filter(
        plotTab, `p adj (dist.dir.up)` <= pCut | `p adj (dist.dir.dn)` <= pCut)
    } else {
      plotTab <- dplyr::filter(
        plotTab, `p (dist.dir.up)` <= pCut | `p (dist.dir.dn)` <= pCut)
    }
    if (nrow(plotTab) == 0) {
      print("No sets passed the criteria")
      next
    } else {
      #firstly, process the result table
      plotTab <- apply(plotTab, 1, function(x) {
        statSign <- as.numeric(x[3])
        data.frame(Name = x[1],
                   p = as.numeric(ifelse(statSign >= 0, x[4], x[6])),
                   geneNum = ifelse(statSign >= 0, x[8], x[9]),
                   Direction = ifelse(statSign > 0, "Up", "Down"),
                   stringsAsFactors = FALSE)
      }) %>% do.call(rbind,.)

      plotTab$Name <- sprintf("%s (%s)",plotTab$Name,plotTab$geneNum)
      plotTab <- plotTab[with(plotTab,order(Direction, p, decreasing=TRUE)),]
      plotTab$Direction <- factor(plotTab$Direction, levels = c("Down","Up"))
      plotTab$Name <- factor(plotTab$Name, levels = plotTab$Name)
      #plot the barplot
      pList[[i]] <- ggplot(data=plotTab, aes(x=Name, y= -log10(p),
                                             fill=Direction)) +
        geom_bar(position="dodge",stat="identity", width = 0.5) +
        scale_fill_manual(values=c(Up = "blue", Down = "red")) +
        coord_fixed(ratio = 0.5) + coord_flip() + xlab(setName) +
        ggtitle(i) + theme_bw() + theme(
          plot.title = element_text(face = "bold", hjust =0.5),
          axis.title = element_text(size=15))
      rowNum <-c(rowNum,nrow(plotTab))
    }
  }

  if (length(pList) == 0) {
    print("Nothing to plot")
  } else {
    rowNum <- rowNum
    grobList <- lapply(pList, ggplotGrob)
    grobList <- do.call(rbind,c(grobList,size="max"))
    panels <- grobList$layout$t[grep("panel", grobList$layout$name)]
    grobList$heights[panels] <- unit(rowNum, "null")
  }
  return(grobList)
}

Prepare input table for gene set enrichment analysis. A cut-off of raw p value < 0.05 was used to select genes for the analysis.

pCut <- 0.05

dataTab <- data.frame(DEres)
dataTab$ID <- rownames(dataTab)

#filter using raw pvalues
dataTab <- filter(dataTab, pvalue <= pCut) %>%
  arrange(pvalue) %>%
  mutate(Symbol = rowData(ddsCLL[ID,])$symbol)
dataTab <- dataTab[!duplicated(dataTab$Symbol),]
statTab <- data.frame(row.names = dataTab$Symbol, stat = dataTab$stat)

Gene set enrichment analysis using Hallmarks gene set from MsigDB.

hallmarkRes <- list()

#run PAGE
resTab <- runGSEA(statTab, gmts$H ,GSAmethod = "page")

## Checking arguments...done!
## Calculating gene set statistics...done!
## Calculating gene set significance...done!
## Adjusting for multiple testing...done!

#remove the HALLMARK_
resTab$Name <- gsub("HALLMARK_","",resTab$Name)

hallmarkRes[["Gene set enrichment analysis"]] <- 
  arrange(resTab,desc(`Stat (dist.dir)`))

hallBar <- plotEnrichmentBar(hallmarkRes, pCut = 0.01, ifFDR = TRUE,
                             setName = "Hallmark gene sets")

Gene set enrichment analysis using kegg gene set from MsigDB.

keggRes <- list()

resTab <- runGSEA(statTab,gmts$KEGG,GSAmethod = "page")

## Checking arguments...done!
## Calculating gene set statistics...done!
## Calculating gene set significance...done!
## Adjusting for multiple testing...done!

#remove the KEGG_
resTab$Name <- gsub("KEGG_","",resTab$Name)

keggRes[["Gene set enrichment analysis"]] <- resTab

keggBar <- plotEnrichmentBar(keggRes, pCut = 0.01, ifFDR = TRUE,
                             setName = "KEGG gene sets")

22.2.2 Heatmap for selected gene sets

Heatmap plots were used to show the expression values of differentially expressed genes from KEGG_CHEMOKINE_SIGNALING_PATHWAY gene set

Prepare the matrix for heatmap plot.

#select differentially expressed genes
fdrCut <- 0.05
cytoDE <- data.frame(DEres) %>% rownames_to_column(var = "ID") %>% 
  mutate(Symbol = rowData(ddsCLL[ID,])$symbol,
         Chr=rowData(ddsCLL[ID,])$chromosome) %>% 
  filter(padj <= fdrCut, log2FoldChange > 0) %>% 
  arrange(pvalue) %>% filter(!duplicated(Symbol)) %>%
  mutate(Chr12 = ifelse(Chr == 12, "yes", "no"))

#get the expression matrix
plotMat <- assay(ddsCLL.norm[cytoDE$ID,])
colnames(plotMat) <- ddsCLL.norm$PatID
rownames(plotMat) <- cytoDE$Symbol

#sort columns of plot matrix based on trisomy 12 status
plotMat <- plotMat[,order(ddsCLL.norm$trisomy12)]

#calculate z-score and sensor
plotMat <- t(scale(t(plotMat)))
plotMat[plotMat >= 4] <- 4
plotMat[plotMat <= -4] <- -4

annoCol <- data.frame(row.names = ddsCLL.norm$PatID,
                      Tris12 = ddsCLL.norm$trisomy12)
levels(annoCol$Tris12) <- list(wt = 0, mut =1)
annoRow <- data.frame(row.names = cytoDE$Symbol, Chr12 = cytoDE$Chr12)

Heatmap for genes from KEGG_CHEMOKINE_SIGNALING_PATHWAY geneset.

gsc <- loadGSC(gmts$KEGG)
geneList <- gsc$gsc$KEGG_CHEMOKINE_SIGNALING_PATHWAY
plotMat.chemo <- plotMat[rownames(plotMat) %in% geneList,]
keggHeatmap <- pheatmap(plotMat.chemo,
                        color = colorRampPalette(
                          rev(brewer.pal(n=7, name="RdBu")))(100),
         cluster_cols = FALSE, clustering_method = "ward.D2",
         annotation_row = annoRow, annotation_col = annoCol,
         show_colnames = FALSE, fontsize_row = 8,
         breaks = seq(-4,4, length.out = 101), treeheight_row = 0,
         annotation_colors = annoColor, border_color = NA,
         main = "CHEMOKINE_SIGNALING_PATHWAY",silent = TRUE)$gtable

Combine enrichment plot and heatmap plot.

#FIG# S24 ABC
ggdraw() + 
  draw_plot(hallBar, 0, 0.7, 0.5, 0.3) + 
  draw_plot(keggBar, 0.5, 0.7, 0.5, 0.3) +
  draw_plot(keggHeatmap, 0.1, 0, 0.8, 0.65) +
  draw_plot_label(c("A","B","C"), c(0, 0.5, 0.1), c(1, 1, 0.7), 
                  fontface = "plain", size=20)

Based on the gene set enrichment analysis results, genes from PI3K_ATK_MTOR pathway are significantly up-regulated in the samples with trisomy 12, which partially explained the increased sensitivity of trisomy 12 samples to PI3K and MTOR inhibitors. In addition, genes that are up-regulated in trisomy 12 are enrichment in chemokine signaling pathway.

options(stringsAsFactors=FALSE)

23.1.1 Data pre-processing

Filtering steps and transformations.

e<-dds
colnames(e)<-colData(e)$PatID


#only consider CLL patients
CLLPatients<-rownames(patmeta)[which(patmeta$Diagnosis=="CLL")]

#Methylation Data
methData = t(assay(methData)) 


#RNA Data
eCLL<-e[,colnames(e) %in% CLLPatients]
###
#filter out genes without gene namce
AnnotationDF<-data.frame(EnsembleId=rownames(eCLL),stringsAsFactors=FALSE)
AnnotationDF$symbol <- mapIds(org.Hs.eg.db,
                     keys=rownames(eCLL),
                     column="SYMBOL",
                     keytype="ENSEMBL",
                     multiVals="first")

## 'select()' returned 1:many mapping between keys and columns

eCLL<-eCLL[AnnotationDF$EnsembleId[!is.na(AnnotationDF$symbol)],]

#filter out low count genes
###
minrs <- 100
rs  <- rowSums(assay(eCLL))
eCLL<-eCLL[ rs >= minrs, ]
#variance stabilize the data
#(includes normalizing for library size and dispsersion estimation) 
vstCounts<-varianceStabilizingTransformation(eCLL)
vstCounts<-assay(vstCounts)
#no NAs in data
any(is.na(vstCounts))

## [1] FALSE

#filter out low variable genes
ntop<-5000
vstCountsFiltered<-vstCounts[order(apply(vstCounts, 1, var, na.rm=T),
                                   decreasing = T)[1:ntop],]
eData<-t(vstCountsFiltered)
#no NAs
any(is.na(eData))

## [1] FALSE

#genetics
#remove features with less than 5 occurences
mutCOMbinary<-channel(mutCOM, "binary")
mutCOMbinary<-mutCOMbinary[featureNames(mutCOMbinary) %in% CLLPatients,]
genData<-Biobase::exprs(mutCOMbinary)
idx <- which(colnames(genData) %in% c("del13q14_bi", "del13q14_mono"))
genData <- genData[,-idx]
colnames(genData)[which(colnames(genData)=="del13q14_any")] = "del13q14"
minObs <- 5
#remove feutes with less than 5 occurecnes
genData<-genData[,colSums(genData, na.rm=T)>=minObs]

#IGHV
translation <- c(`U` = 0, `M` = 1)
stopifnot(all(patmeta$IGHV %in% c("U","M", NA)))
IGHVData <- matrix(translation[patmeta$IGHV], 
                   dimnames = list(rownames(patmeta), "IGHV"), ncol = 1)
IGHVData<-IGHVData[rownames(IGHVData) %in% CLLPatients,,drop=F]
#remove patiente with NA IGHV status
IGHVData<-IGHVData[!is.na(IGHVData), ,drop=F]
any(is.na(IGHVData))

## [1] FALSE

#demographics (age and sex)
patmeta<-subset(patmeta, Diagnosis=="CLL")
gender <- ifelse(patmeta[,"Gender"]=="m",0,1)


# impute missing values in age by mean
ImputeByMean <- function(x) {x[is.na(x)] <- mean(x, na.rm=TRUE); return(x)}
age<-ImputeByMean(patmeta[,"Age4Main"])


demogrData <- cbind(age=age,gender=gender)
rownames(demogrData) <- rownames(patmeta)

#Pretreatment
pretreated<-patmeta[,"IC50beforeTreatment", drop=FALSE]

##### drug viabilites
summaries <- c(paste("viaraw", 1:5, sep=".") %>% `names<-`(paste(1:5)), 
               `4:5` = "viaraw.4_5", `1:5` = "viaraw.1_5")
a <- do.call( abind, c( along=3, lapply( summaries, 
                                         function(x) assayData(drpar)[[x]]))) 
dimnames(a)[[3]] <- names(summaries)
names(dimnames(a)) <- c( "drug", "patient", "summary" )
viabData <- acast( melt(a), patient ~ drug + summary )
rownames(viabData)<-c(substr(rownames(viabData),1,4)[1:3],
                      substr(rownames(viabData),1,5)[4:nrow(viabData)])

Check overlap of data and take care of missing values present in methylation and genetic data.

# common patients 
Xlist<-list(RNA=eData, meth=methData, gen=genData, IGHV=IGHVData,
            demographics=demogrData, drugs=viabData, pretreated=pretreated)
PatientsPerOmic<-lapply(Xlist, rownames)
sapply(PatientsPerOmic, length)

##          RNA         meth          gen         IGHV demographics 
##          136          196          200          188          200 
##        drugs   pretreated 
##          249          200

allPatients<-Reduce(union, PatientsPerOmic)
PatientOverview<-sapply(Xlist, function(M) allPatients %in% rownames(M))
Patients <- (1:nrow(PatientOverview))
Omics <- (1:ncol(PatientOverview))
image(Patients,Omics, PatientOverview*1, axes=F, col=c("white", "black"),
      main="Sample overview across omics")
axis(2, at = 1:ncol(PatientOverview), labels=colnames(PatientOverview), tick=F)

commonPatients<-Reduce(intersect, PatientsPerOmic)
length(commonPatients)

## [1] 112

XlistCommon<-lapply(Xlist, function(data) data[commonPatients,, drop=F])

#Take care of missing values (present in  genetic data)
ImputeByMean <- function(x) {x[is.na(x)] <- mean(x, na.rm=TRUE); return(x)}

#NAs in genetic
#remove feauters with less 90% completeness
RarlyMeasuredFeautres<-
  which(colSums(is.na(XlistCommon$gen))>0.1*nrow(XlistCommon$gen))
XlistCommon$gen<-XlistCommon$gen[,-RarlyMeasuredFeautres]
#remove patients with less than 90% of genetic feautres measured
IncompletePatients<-
  rownames(XlistCommon$gen)[
    (rowSums(is.na(XlistCommon$gen))>0.1*ncol(XlistCommon$gen))]
commonPatients<-commonPatients[!commonPatients %in% IncompletePatients]
XlistCommon<-lapply(XlistCommon, function(data) data[commonPatients,, drop=F])
#replace remaining NA by mean and round to 0 or 1
XlistCommon$gen<-round(apply(XlistCommon$gen, 2, ImputeByMean))

#NAs in methylation
#remove feauters with less 90% completeness
XlistCommon$meth<-
  XlistCommon$meth[,colSums(is.na(XlistCommon$meth))<0.1*nrow(methData)]
#impute remainin missing values by mean for each feautre across patients
XlistCommon$meth<-(apply(XlistCommon$meth, 2, ImputeByMean))

#final dimensions of the data
sapply(XlistCommon, dim)

##       RNA meth gen IGHV demographics drugs pretreated
## [1,]  102  102 102  102          102   102        102
## [2,] 5000 5000  11    1            2   448          1

Use top 20 PCs of methylation and expression as predictors.

pcaMeth<-prcomp(XlistCommon$meth, center=T, scale. = F)
XlistCommon$MethPCs<-pcaMeth$x[,1:20]
colnames(XlistCommon$MethPCs)<-
  paste("meth",colnames(XlistCommon$MethPCs), sep="")

pcaExpr<-prcomp(XlistCommon$RNA, center=T, scale. = F)
XlistCommon$RNAPCs<-pcaExpr$x[,1:20]
colnames(XlistCommon$RNAPCs)<-paste("RNA",colnames(XlistCommon$RNAPCs), sep="")

Choose drug viabilites of interest as response variables.

DOI <- c("D_006_1:5", "D_010_1:5", "D_159_1:5","D_002_4:5", "D_003_4:5",
         "D_012_4:5", "D_063_4:5", "D_166_4:5")
drugviab<-XlistCommon$drugs
drugviab<-drugviab[,DOI, drop=F]
colnames(drugviab) <- drugs[substr(colnames(drugviab),1,5),"name"]

Construct list of designs used and scale all predictors to mean zero and unit variance.

ZPCs<-list(expression=XlistCommon$RNAPCs,
        genetic=XlistCommon$gen, 
        methylation= XlistCommon$MethPCs,
        demographics=XlistCommon$demographics, 
        IGHV=XlistCommon$IGHV,
        pretreated = XlistCommon$pretreated)
ZPCs$all<-do.call(cbind, ZPCs)
ZPCsunscaled<-ZPCs
ZPCsscaled<-lapply(ZPCs, scale)
lapply(ZPCsscaled, colnames)

## $expression
##  [1] "RNAPC1"  "RNAPC2"  "RNAPC3"  "RNAPC4"  "RNAPC5"  "RNAPC6"  "RNAPC7" 
##  [8] "RNAPC8"  "RNAPC9"  "RNAPC10" "RNAPC11" "RNAPC12" "RNAPC13" "RNAPC14"
## [15] "RNAPC15" "RNAPC16" "RNAPC17" "RNAPC18" "RNAPC19" "RNAPC20"
## 
## $genetic
##  [1] "del6q21"    "gain8q24"   "del11q22.3" "trisomy12"  "del13q14"  
##  [6] "del17p13"   "BRAF"       "MYD88"      "NOTCH1"     "SF3B1"     
## [11] "TP53"      
## 
## $methylation
##  [1] "methPC1"  "methPC2"  "methPC3"  "methPC4"  "methPC5"  "methPC6" 
##  [7] "methPC7"  "methPC8"  "methPC9"  "methPC10" "methPC11" "methPC12"
## [13] "methPC13" "methPC14" "methPC15" "methPC16" "methPC17" "methPC18"
## [19] "methPC19" "methPC20"
## 
## $demographics
## [1] "age"    "gender"
## 
## $IGHV
## [1] "IGHV"
## 
## $pretreated
## [1] "IC50beforeTreatment"
## 
## $all
##  [1] "expression.RNAPC1"    "expression.RNAPC2"    "expression.RNAPC3"   
##  [4] "expression.RNAPC4"    "expression.RNAPC5"    "expression.RNAPC6"   
##  [7] "expression.RNAPC7"    "expression.RNAPC8"    "expression.RNAPC9"   
## [10] "expression.RNAPC10"   "expression.RNAPC11"   "expression.RNAPC12"  
## [13] "expression.RNAPC13"   "expression.RNAPC14"   "expression.RNAPC15"  
## [16] "expression.RNAPC16"   "expression.RNAPC17"   "expression.RNAPC18"  
## [19] "expression.RNAPC19"   "expression.RNAPC20"   "genetic.del6q21"     
## [22] "genetic.gain8q24"     "genetic.del11q22.3"   "genetic.trisomy12"   
## [25] "genetic.del13q14"     "genetic.del17p13"     "genetic.BRAF"        
## [28] "genetic.MYD88"        "genetic.NOTCH1"       "genetic.SF3B1"       
## [31] "genetic.TP53"         "methylation.methPC1"  "methylation.methPC2" 
## [34] "methylation.methPC3"  "methylation.methPC4"  "methylation.methPC5" 
## [37] "methylation.methPC6"  "methylation.methPC7"  "methylation.methPC8" 
## [40] "methylation.methPC9"  "methylation.methPC10" "methylation.methPC11"
## [43] "methylation.methPC12" "methylation.methPC13" "methylation.methPC14"
## [46] "methylation.methPC15" "methylation.methPC16" "methylation.methPC17"
## [49] "methylation.methPC18" "methylation.methPC19" "methylation.methPC20"
## [52] "demographics.age"     "demographics.gender"  "IGHV"                
## [55] "IC50beforeTreatment"

Define colors.

set1 <- brewer.pal(9,"Set1")
colMod<-c(paste(set1[c(4,1,5,3,2,7)],"88",sep=""), "grey")
names(colMod) <-
  c("demographics", "genetic", "IGHV","expression", "methylation", "pretreated",
    "all")

23.1.2 Lasso using multi-omic data

Fit a linear model using Lasso explaining drug response by each one of the omic set separately as well as all together. As measure of explained variance use R2 from linear models for unpenalized models (IGHV) and fraction of variance explained, i.e. 1- cross-validated mean squared error/total sum of squares for others.

To ensure fair treatment of all features they are standardized to mean 0 and unit variance. To study robustness the cross-validation is repeated 100-times to obtain the mean and standard deviation shown in the figure.

Fit model and show resulting omic-prediction profiles.

Optionally, the model can also be fit for all drugs in the study.

nfolds<-10
nrep<-100

DOI <-
  grepl("1:5",colnames(XlistCommon$drugs)) |
  grepl("4:5",colnames(XlistCommon$drugs))
drugviabAll<-XlistCommon$drugs
drugviabAll<-drugviabAll[,DOI]
colnames(drugviabAll) <- 
  paste0(drugs[substr(colnames(drugviabAll),1,5),"name"],
         substr(colnames(drugviabAll),6,9))

R2ForPenReg(Zscaled, drugviabAll, nfold=nfolds, alpha=1, nrep=nrep,
            Parmfrow=c(4,4), ylimMax=0.6)

23.2.1 General definitions

We use the following abbreviations for the different data types.

Color definitions for the groups.

23.2.2 Data pre-processing

Subselect CLL patients.

lpdCLL = lpdAll[ , lpdAll$Diagnosis=="CLL"]

Prepare the data and limit the number of features by subselecting only those which include at least 5 recorded incidences. List the predictors.

lpdCLL = lpdAll[ , lpdAll$Diagnosis=="CLL"]
lpdCLL = BloodCancerMultiOmics2017:::prepareLPD(lpd=lpdCLL, minNumSamplesPerGroup=5)
(predictorList = BloodCancerMultiOmics2017:::makeListOfPredictors(lpdCLL))

## $predictorsM
## [1] "Methylation_Cluster"
## 
## $predictorsG
##  [1] "del6q21"    "gain8q24"   "del11q22.3" "trisomy12"  "del13q14"  
##  [6] "del17p13"   "BRAF"       "NOTCH1"     "SF3B1"      "TP53"      
## 
## $predictorsI
## [1] "IGHV"
## 
## $predictorsP
## [1] "Pretreatment"

23.2.3 Drug response prediction

The prediction will be made for the following drugs and concentrations.

drs = list("1:5"=c("D_006", "D_010", "D_159"),
           "4:5"=c("D_002", "D_003", "D_012", "D_063", "D_166"))

predvar = unlist(BloodCancerMultiOmics2017:::makePredictions(drs=drs,
                                 lpd=lpdCLL,
                                 predictorList=predictorList,
                                 lim=0.15, std=FALSE, adaLasso=TRUE,
                                 colors=coldef),
                 recursive=FALSE)

## [1] "Prediction for: D_006_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_010_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_159_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_002_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_003_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_012_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_063_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_166_4:5; #samples: 168; #features: 13"

Plot the legends.

legends = BloodCancerMultiOmics2017:::makeLegends(legendFor=c("G","I","M", "P"),
                                                  coldef)

Additionaly we plot the prediction for rotenone.

drs_rot = list("4:5"=c("D_067"))
predvar_rot = unlist(BloodCancerMultiOmics2017:::makePredictions(drs=drs_rot,
                                 lpd=lpdCLL,
                                 predictorList=predictorList,
                                 lim=0.23, std=FALSE, adaLasso=TRUE,
                                 colors=coldef),
                 recursive=FALSE)

## [1] "Prediction for: D_067_4:5; #samples: 168; #features: 13"

In a same way the prediction for all the drugs can be made.

alldrs = unique(fData(lpdCLL)[fData(lpdCLL)$type=="viab","id"])
drs = list("1:5"=alldrs, "4:5"=alldrs)

predvar2 = BloodCancerMultiOmics2017:::makePredictions(drs=drs,
                                        lpd=lpdCLL,
                                        predictorList=predictorList,
                                        lim=0.23,
                                        colors=coldef)

## [1] "Prediction for: D_001_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_002_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_003_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_004_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_006_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_007_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_008_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_010_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_011_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_012_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_013_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_015_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_017_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_020_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_021_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_023_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_024_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_025_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_029_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_030_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_032_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_033_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_034_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_035_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_036_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_039_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_040_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_041_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_043_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_045_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_048_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_049_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_050_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_053_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_054_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_056_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_060_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_063_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_066_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_067_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_071_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_074_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_075_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_077_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_078_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_079_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_081_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_082_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_083_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_084_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_127_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_141_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_149_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_159_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_162_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_163_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_164_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_165_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_166_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_167_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_168_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_169_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_172_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_CHK_1:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_001_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_002_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_003_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_004_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_006_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_007_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_008_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_010_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_011_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_012_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_013_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_015_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_017_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_020_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_021_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_023_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_024_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_025_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_029_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_030_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_032_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_033_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_034_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_035_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_036_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_039_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_040_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_041_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_043_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_045_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_048_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_049_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_050_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_053_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_054_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_056_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_060_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_063_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_066_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_067_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_071_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_074_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_075_4:5; #samples: 168; #features: 13"
## [1] "No (0) predictors for given parameters!"
## [1] "Prediction for: D_077_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_078_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_079_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_081_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_082_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_083_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_084_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_127_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_141_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_149_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_159_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_162_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_163_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_164_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_165_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_166_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_167_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_168_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_169_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_172_4:5; #samples: 168; #features: 13"
## [1] "Prediction for: D_CHK_4:5; #samples: 168; #features: 13"

23.2.4 Effect of pre-treatment

In order to find out the effect of pre-treatment on the predictions of drug response we provide the following overview. The summary is made for all drugs, separating however, on ranges of drug concentrations (mean drug effect of: all five (1-5) and two lowest (4-5) concentrations of the drugs).

givePreatreatSum = function(predNum) {
  
  idx = sapply(predvar2[[predNum]], function(x) length(x)==1)
  predvar2[[predNum]] = predvar2[[predNum]][!idx]
  # get model coefficients and reshape
  coeffs <- do.call(cbind,lapply(predvar2[[predNum]], "[[", 'coeffs'))
  coeffs <- coeffs[-1,]
  coeffs <- as.matrix(coeffs)
  # colnames(coeffs) <- unlist(drs["1:5"])
  colnames(coeffs) = names(predvar2[[predNum]])
  colnames(coeffs) <- drugs[colnames(coeffs),"name"]
  coeffDF <- melt(as.matrix(coeffs))
  colnames(coeffDF) <- c("predictor", "drug", "beta")
  coeffDF$selected <- coeffDF$beta !=0
  
  #sort by times being selected
  coeffDF$predictor <- factor(coeffDF$predictor, level=)
  
  # number of drugs a predictor is chosen for
  gg1 <- coeffDF %>% group_by(predictor) %>% 
    dplyr::summarize(selectedSum = sum(selected)) %>%
    mutate(predictor = factor(predictor,
                              levels=predictor[order(selectedSum)])) %>%
    ggplot(aes(x=predictor, y=selectedSum)) + 
    geom_bar(stat="identity")+ylab("# drugs selected for") +
    coord_flip()
  
  # boxplots of non-zero coeffients
  orderMedian <- filter(coeffDF, selected) %>% group_by(predictor) %>% 
    dplyr::summarize(medianBeta = median(abs(beta)))
  coeffDF$predictor <- factor(
    coeffDF$predictor,
    levels=orderMedian$predictor[order(orderMedian$medianBeta)] )
  gg2 <- ggplot(filter(coeffDF, selected), aes(x=predictor, y=abs(beta))) +
    geom_boxplot() +
    coord_flip() + ggtitle("Distribution of non-zero coefficients")
  gridExtra::grid.arrange(gg1,gg2, ncol=1)
  
  # coefficeints per drug
  ggplot(filter(coeffDF, selected), 
         aes(x= drug, y=abs(beta), col= predictor=="Pretreatment")) +
    geom_point() +
    coord_flip()
  
  #drugs pretreatment is selected for
  as.character(filter(coeffDF, predictor=="Pretreatment" & beta!=0)$drug)
  PselDrugs <- as.character(
    filter(coeffDF, predictor=="Pretreatment" & beta!=0)$drug)
  length(PselDrugs)
  # length(drs[[1]])
}

Conc. 1-5

## [1] 26

Conc. 4-5

## [1] 10

options(stringsAsFactors=FALSE)

24.1.1 Define forest functions

forest <- function(Time, endpoint, title, sdrugs, split, sub) {  
  stopifnot(is.character(Time), is.character(title), is.character(split),
            is.character(endpoint), 
            all(c(Time, split, endpoint) %in% colnames(survT)),
            is.logical(survT[[endpoint]]),
            is.character(sdrugs), !is.null(names(sdrugs)))
  
  clrs <- fpColors(box="royalblue",line="darkblue", summary="royalblue")

  res <- lapply(sdrugs, function(g) { 
    drug <- survT[, g] * 10

    suse <- if (identical(sub, "none")) 
      rep(TRUE, nrow(survT)) 
    else 
      (survT[[split]] == sub)
    stopifnot(sum(suse, na.rm = TRUE) > 1)
    
    surv <- coxph(Surv(survT[,Time], survT[,endpoint]) ~ drug, subset=suse)  
    sumsu <- summary(surv) 
    c(p      = sumsu[["coefficients"]][, "Pr(>|z|)"], 
      coef   = sumsu[["coefficients"]][, "exp(coef)"], 
      lower  = sumsu[["conf.int"]][, "lower .95"], 
      higher = sumsu[["conf.int"]][, "upper .95"])
  })
  
  s <- do.call(rbind, res)
  rownames(s) <- names(sdrugs)

  tabletext <- list(c(NA, rownames(s)), append(list("p-value"),
                                               sprintf("%.4f", s[,"p"])))

  forestplot(tabletext, 
           rbind(
             rep(NA, 3), 
             s[, 2:4]), 
           page = new,
           clip = c(0.8,20), 
           xlog = TRUE, xticks = c(0.5,1, 1.5), title = title,
           col = clrs, 
           txt_gp = fpTxtGp(ticks = gpar(cex=1) ),
           new_page = TRUE)
}

Combine OS and TTT in one plot.

com <- function( Time, endpoint, scaleX, sub, d, split, drug_names) {  
  
  res <- lapply(d, function(g)  { 
  
  drug <- survT[,g] * scaleX
  ## all=99, M-CLL=1, U-CLL=0
  if(sub==99) { surv <- coxph(Surv(survT[,paste0(Time)],
                                   survT[,paste0(endpoint)] == TRUE) ~ drug)} 
  if(sub<99)  { surv <- coxph(Surv(survT[,paste0(Time)],
                                   survT[,paste0(endpoint)] == TRUE) ~ drug,
                              subset=survT[,paste0(split)]==sub)}    
  
  c(summary(surv)[[7]][,5], summary(surv)[[7]][,2], 
                summary(surv)[[8]][,3], 
                summary(surv)[[8]][,4])
 })
 s <- do.call(rbind, res)
 colnames(s) <- c("p", "HR", "lower", "higher")
 rownames(s) <- drug_names
 s
}


fp <- function( sub, title, d, split, drug_names, a, b, scaleX) {  
   ttt <- com(Time="T5", endpoint="treatedAfter", sub=sub, d=d,
              split=split, drug_names=drug_names, scaleX=scaleX)
   rownames(ttt) <- paste0(rownames(ttt), "_TTT")
   
   os <-  com(Time="T6", endpoint="died", sub=sub, d=d, split=split,
              drug_names=drug_names, scaleX=scaleX)
   rownames(os) <- paste0(rownames(os), "_OS")
  
   n <- c( p=NA, HR=NA, lower=NA, higher=NA )
   nn <- t( data.frame( n ) )
   for (i in 1:(nrow(ttt)-1) ) { nn <-rbind(nn, n )  }
   rownames(nn) <- drug_names
   
   od <- order( c(seq(nrow(nn)), seq(nrow(ttt)), seq(nrow(os)) ))
      
   s <- data.frame( rbind(nn, ttt, os)[od,  ] )
   s$Name <- rownames(s)
   s$x <- 1:nrow(s)
   s$col <- rep(c("white", "black", "darkgreen"), nrow(ttt) )
   s$Endpoint <- factor( c(rep("nn", nrow(nn) ), rep("TTT", nrow(ttt) ),
                           rep("OS", nrow(os) ) )[od] )
   s$features <- "";  s[ which(s$Endpoint=="OS"),"features"] <- drug_names
   s[which(s$Endpoint=="nn"), "Endpoint"] <- "" 
   s <- rbind(s, rep(NA, 8))

   p <- ggplot(data=s ,aes(x=x, y=HR, ymin=lower, ymax=higher,
                           colour=Endpoint)) +  geom_pointrange() + 
     theme(legend.position="top", legend.text = element_text(size = 20) ) +
     scale_x_discrete(limits=s$x, labels=s$features ) +
     expand_limits(y=c(a,b)) +
     scale_y_log10(breaks=c(0.01,0.1,0.5,1,2,5,10),
                   labels=c(0.01,0.1,0.5,1,2,5,10)) +
     theme(
           panel.grid.minor = element_blank(),
           axis.title.x = element_text(size=16),
           axis.text.x = element_text(size=16, colour="black"),
           axis.title.y  = element_blank(),
           axis.text.y = element_text(size=12, colour="black"),
           legend.key = element_rect(fill = "white"),
           legend.background = element_rect(fill = "white"),
           legend.title = element_blank(),
           panel.background = element_rect(fill = "white", color="black"),
           panel.grid.major = element_blank(),
           axis.ticks.y = element_blank() 
          ) +
     coord_flip() + 
     scale_color_manual(values=c("OS"="darkgreen", "TTT"="black"),
                        labels=c("OS", "TTT", "")) +
     geom_hline(aes(yintercept=1), colour="black", size=1.5,
                linetype="dashed", alpha=0.3)  +
     annotate("text", x = 1:nrow(s)+0.5, y = s$HR+0.003,
              label = ifelse( s$p<0.001, paste0("p<","0.001"), 
                paste0("p=", round(s$p,3) ) ), colour=s$col)
   plot(p)
}

24.1.2 Forest plot for genetic factors

#FIG# S27
d <- c("SF3B1", "NOTCH1", "BRAF", "TP53", "del17p13", "del11q22.3",
       "trisomy12", "IGHV")
drug_names <- c("SF3B1", "NOTCH1", "BRAF", "TP53", "del17p13", "del11q22.3",
                "Trisomy12" ,"IGHV")

fp(sub=99, d=d, drug_names=drug_names, split="IGHV", title="", a=0, b=10,
   scaleX=1)

24.1.3 Forest plot for drug responses

#FIG# 6A
d <- c("flu15", "nut15", "dox15", "ibr45", "ide45", "prt45", "selu45",
       "ever45")
drug_names <- c("Fludarabine",  "Nutlin-3", "Doxorubicine", "Ibrutinib",
                "Idelalisib", "PRT062607 HCl", "Selumetinib" ,"Everolimus")

fp(sub=99, d=d, drug_names=drug_names, split="TP53", title="", a=0, b=5,
   scaleX=10)

24.2.1 Genetics factors

#FIG# S27 left (top+bottom)
par(mfcol=c(1,2))

for (fac in paste(c("IGHV", "TP53"))) {
survplot( Surv(survT$T5, survT$treatedAfter == TRUE)  ~ as.factor(survT[,fac]), 
   snames=c("wt", "mut"),
   lwd=1.5, cex.axis = 1, cex.lab=1, col= c("darkmagenta", "dodgerblue4"),
   show.nrisk = FALSE,
   legend.pos = FALSE, stitle = "", hr.pos= "topright",
   main = paste(fac), 
   xlab = 'Time (Years)', ylab = 'Time to treatment')
}

#FIG# 6B left
#FIG# S27 right (top+bottom)
par(mfcol=c(1,2))

for (fac in paste(c("IGHV", "TP53"))) {
survplot( Surv(survT$T6, survT$died == TRUE)  ~ as.factor(survT[,fac]), 
   snames=c("wt", "mut"),
   lwd=1.5, cex.axis = 1.0, cex.lab=1.0, col= c("darkmagenta", "dodgerblue4"),
   show.nrisk = FALSE,
   legend.pos = FALSE, stitle = "", hr.pos= "bottomleft",
   main = paste(fac), 
   xlab = 'Time (Years)', ylab = 'Overall survival')
}  

24.2.2 Drug responses

Drug responses were dichotomized using maximally selected rank statistics. The analysis is also perforemd within subgroups: TP53 wt/ mut. and IGHV wt/ mut.

Time to next treatment (maxstat).

par(mfrow=c(1,3), mar=c(5,5,2,0.9))
km(drug = "D_006_1:5", split = "TP53", t="TTT",
   title=c("(TP53", "wt)", "mut)"),  hr="tr", c="maxstat")

## Fludarabine cutpoint for TTT:  0.65

km(drug = "D_159_1:5", split = "TP53", t="TTT",
   title=c("(TP53", "wt)", "mut)"), hr="tr",  c="maxstat")

## Doxorubicine cutpoint for TTT:  0.56

km(drug = "D_010_1:5", split = "TP53", t="TTT",
   title=c("(TP53", "wt)", "mut)"), hr="tr",  c="maxstat")  

## Nutlin-3 cutpoint for TTT:  0.85

km(drug = "D_002_4:5", split = "IGHV", t="TTT",
   title=c("(IGHV",  "wt)" , "mut)"), hr="tr", c="maxstat" )

## Ibrutinib cutpoint for TTT:     1

km(drug = "D_003_4:5", split = "IGHV", t="TTT",
   title=c("(IGHV",  "wt)" , "mut)"), hr="tr", c="maxstat" )

## Idelalisib cutpoint for TTT:  0.99

km(drug = "D_166_4:5", split = "IGHV", t="TTT",
   title=c("(IGHV",  "wt)" , "mut)"), hr="tr", c="maxstat" ) 

## PRT062607 HCl cutpoint for TTT:  0.98

Overall survival (maxstat).

par(mfrow=c(1,3), mar=c(5,5,2,0.9))
km(drug = "D_006_1:5", split = "TP53", t="OS",
   title=c("(TP53", "wt)", "mut)"), hr="bl", c="maxstat")

## Fludarabine cutpoint for OS:   0.6

#FIG# 6B right
#FIG# 6C
km(drug = "D_159_1:5", split = "TP53", t="OS", # doxorubicine
   title=c("(TP53", "wt)", "mut)"), hr="bl", c="maxstat" )

## Doxorubicine cutpoint for OS:  0.53

#FIG# 6B middle
km(drug = "D_010_1:5", split = "TP53", t="OS", # nutlin-3
   title=c("(TP53", "wt)", "mut)"), hr="bl", c="maxstat" )

## Nutlin-3 cutpoint for OS:  0.89

km(drug = "D_002_4:5", split = "IGHV", t="OS",
   title=c("(IGHV",  "wt)" , "mut)"), hr="bl", c="maxstat" )

## Ibrutinib cutpoint for OS:  0.95

km(drug = "D_003_4:5", split = "IGHV", t="OS",
   title=c("(IGHV",  "wt)" , "mut)"), hr="bl", c="maxstat" )

## Idelalisib cutpoint for OS:   0.9

## Warning in fitter(X, Y, strats, offset, init, control, weights = weights, :
## Loglik converged before variable 1 ; coefficient may be infinite.

km(drug = "D_166_4:5", split = "IGHV", t="OS",
   title=c("(IGHV",  "wt)" , "mut)"), hr="bl", c="maxstat" ) 

## PRT062607 HCl cutpoint for OS:  0.87

24.3.1 Chemotherapies

surv1 <- coxph(
  Surv(T6, died) ~  
    age +
    as.factor(IC50beforeTreatment) +
    as.factor(trisomy12) +
    as.factor(del11q22.3) +
    as.factor(del17p13) +
    as.factor(TP53) +
    IGHVwt +
    flu15,       # continuous
    #dox15 +     # continuous
    #flu15:TP53,
    #TP53:dox15,
  data = survT )
extractSome(surv1)

##    p.value     HR lower.95..CI upper.95..CI
## 1    0.403    1.2         0.81          1.7
## 2 0.000285   0.11        0.033         0.36
## 3   0.0196    5.3          1.3           22
## 4    0.784    1.2         0.38          3.6
## 5    0.962      1          0.3          3.6
## 6     0.47   0.63         0.18          2.2
## 7   0.0516    2.8         0.99          7.8
## 8    0.085   0.77         0.57            1

cat(sprintf("%s patients considerd in the model; number of events %1g\n", 
            summary(surv1)$n, summary(surv1)[6] ) )

## 156 patients considerd in the model; number of events 24

surv2 <- coxph(
  Surv(T6, died) ~   #as.factor(survT$TP53) , data=survT )
    age +
    as.factor(IC50beforeTreatment) +
    as.factor(trisomy12) +
    as.factor(del11q22.3) +
    as.factor(del17p13) +
    as.factor(TP53) +
    IGHVwt +
    #flu15 +    # continuous
    dox15 ,     # continuous
    #flu15:TP53 ,
    #TP53:dox15,
  data = survT )
extractSome(surv2)

##    p.value     HR lower.95..CI upper.95..CI
## 1    0.121    1.3         0.92            2
## 2 0.000186   0.11        0.036         0.36
## 3   0.0123    6.1          1.5           25
## 4    0.992      1         0.34            3
## 5     0.92    1.1          0.3          3.8
## 6    0.716   0.81         0.26          2.5
## 7   0.0393    2.9          1.1          8.1
## 8   0.0347   0.52         0.28         0.95

cat(sprintf("%s patients considerd in the model; number of events %1g\n", 
            summary(surv2)$n, summary(surv2)[6] ) )

## 156 patients considerd in the model; number of events 24

24.3.2 Targeted therapies

surv4 <- coxph(
  Surv(T5, treatedAfter) ~ 
    age +
    as.factor(IC50beforeTreatment) +
    as.factor(trisomy12) +
    as.factor(del11q22.3) +
    as.factor(del17p13) +
    IGHVwt +
    ibr45 +
    IGHVwt:ibr45,
  data = survT )

extractSome(surv4)

##    p.value     HR lower.95..CI upper.95..CI
## 1    0.284   0.91         0.76          1.1
## 2 9.61e-08   0.22         0.13         0.39
## 3    0.246    1.5         0.76            3
## 4    0.778   0.91         0.48          1.7
## 5     0.31    1.4         0.74          2.6
## 6   0.0111    2.2          1.2          3.9
## 7   0.0289    1.6          1.1          2.5
## 8   0.0417    0.6         0.37         0.98

cat(sprintf("%s patients considerd in the model; number of events %1g\n", 
            summary(surv4)$n, summary(surv4)[6] ) )

## 152 patients considerd in the model; number of events 83

surv6 <- coxph(
  Surv(T5, treatedAfter) ~ 
    age +
    as.factor(IC50beforeTreatment) +
    as.factor(trisomy12) +
    as.factor(del11q22.3) +
    as.factor(del17p13) +
    IGHVwt +
    ide45 +
    IGHVwt:ide45,
  data = survT )

extractSome(surv6)

##    p.value     HR lower.95..CI upper.95..CI
## 1    0.302   0.91         0.76          1.1
## 2 3.44e-08   0.21         0.12         0.37
## 3    0.214    1.5         0.78            3
## 4    0.667   0.87         0.46          1.7
## 5    0.308    1.4         0.74          2.6
## 6  0.00789    2.4          1.3          4.7
## 7   0.0421    1.6            1          2.6
## 8    0.067    0.6         0.35            1

cat(sprintf("%s patients considerd in the model; number of events %1g\n",
            summary(surv6)$n, summary(surv6)[6] ) )

## 152 patients considerd in the model; number of events 83

surv8 <- coxph(
  Surv(T5, treatedAfter) ~ 
    age +
    as.factor(IC50beforeTreatment) +
    as.factor(trisomy12) +
    as.factor(del11q22.3) +
    as.factor(del17p13) +
    IGHVwt +
    prt45 +
    IGHVwt:prt45,
  data = survT )

extractSome(surv8)

##    p.value     HR lower.95..CI upper.95..CI
## 1    0.402   0.93         0.78          1.1
## 2 2.19e-08    0.2         0.12         0.36
## 3    0.397    1.4         0.67          2.7
## 4    0.584   0.83         0.43          1.6
## 5    0.415    1.3         0.69          2.4
## 6  0.00446    2.7          1.4          5.4
## 7    0.012    1.6          1.1          2.4
## 8    0.018   0.58         0.37         0.91

cat(sprintf("%s patients considerd in the model; number of events %1g\n", 
            summary(surv8)$n, summary(surv8)[6] ) )

## 152 patients considerd in the model; number of events 83