---

title: "Part 4"

output:

  BiocStyle::html_document

---

```{r, message=FALSE, include=!exists(".standalone"), eval=!exists(".standalone")}

library("BloodCancerMultiOmics2017")

library("Biobase")

library("ggplot2")

library("grid")

```

```{r echo=FALSE}

plotDir = ifelse(exists(".standalone"), "", "part04/")

if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)

```

```{r}

options(stringsAsFactors=FALSE)

```

# Single associations of drug response with gene mutation or type of disease (IGHV included)

We univariantly tested different features (explained in detail below) for their associations with the drug response using Student t-test (two-sided, with equal variance). Each concentration was tested separately. The minimal size of the compared groups was set to 3. p-values were adjusted for multiple testing by applying the Benjamini-Hochberg procedure. Adjusted p-values were then used for setting the significance threshold.

Loading the data.

```{r}

data(list=c("drpar", "patmeta", "drugs", "mutCOM", "conctab"))

```

## Function which test associations of interest

Below is a general function with which all the tests for single associations were performed.

```{r}

testFactors = function(msrmnt, factors, test="student", batch=NA) {

  # cut out the data

  tmp = colnames(factors)

  factors = data.frame(factors[rownames(msrmnt),], check.names=FALSE)

  colnames(factors) = tmp

  for(cidx in 1:ncol(factors))

    factors[,cidx] = factor(factors[,cidx], levels=c(0,1))

  # calculate the group size

  groupSizes = do.call(rbind, lapply(factors, function(tf) {

    tmp = table(tf)

    data.frame(n.0=tmp["0"], n.1=tmp["1"])

  }))

  # remove the factors with less then 2 cases per group

  factors = factors[,names(which(apply(groupSizes, 1,

                                       function(i) all(i>2)))), drop=FALSE]

  # calculate the effect

  effect = do.call(rbind, lapply(colnames(factors), function(tf) {

    tmp = aggregate(msrmnt ~ fac, data=data.frame(fac=factors[,tf]), mean)

    rownames(tmp) = paste("mean", tmp$fac, sep=".")

    tmp = t(tmp[2:ncol(tmp)])

    data.frame(TestFac=tf,

               DrugID=rownames(tmp),

               FacDr=paste(tf, rownames(tmp), sep="."),

               n.0=groupSizes[tf,"n.0"], n.1=groupSizes[tf,"n.1"],

               tmp, WM=tmp[,"mean.0"]-tmp[,"mean.1"])

  }))

  # do the test

  T = if(test=="student") {

    do.call(rbind, lapply(colnames(factors), function(tf) {

      tmp = do.call(rbind, lapply(colnames(msrmnt), function(dr) {

        res = t.test(msrmnt[,dr] ~ factors[,tf], var.equal=TRUE)

        data.frame(DrugID=dr, TestFac=tf,

                   pval=res$p.value, t=res$statistic,

                   conf1=res$conf.int[1], conf2=res$conf.int[2])

      }))

      tmp

    }))

  } else if(test=="anova") {

    do.call(rbind, lapply(colnames(factors), function(tf) {

      tmp = do.call(rbind, lapply(colnames(msrmnt), function(dr) {

        # make sure that the order in batch is the same as in msrmnt

        stopifnot(identical(rownames(msrmnt), names(batch)))

        res = anova(lm(msrmnt[,dr] ~ factors[,tf]+batch))

        data.frame(DrugID=dr, TestFac=tf, pval=res$`Pr(>F)`[1],

                   f=res$`F value`[1], meanSq1=res$`Mean Sq`[1],

                   meanSq2=res$`Mean Sq`[2])

      }))

      tmp

    }))

  } else {

    NA

  }

  enhanceObject = function(obj) {

    # give nice drug names

    obj$Drug = giveDrugLabel(obj$DrugID, conctab, drugs)

    # combine the testfac and drug id

    obj$FacDr = paste(obj$TestFac, obj$DrugID, sep=".")

    # select just the drug name

    obj$DrugID2 = substring(obj$DrugID, 1, 5)

    obj

  }

  list(effect=effect, test=enhanceObject(T))

}

```

## Associations of *ex vivo* drug responses with genomic features in CLL

### Prepare objects for testing

```{r}

## 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)

}))

```

### 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. 

```{r}

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.

```{r}

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. 

```{r batchEffect, results='asis', echo=FALSE, fig.path=plotDir, dev=c('png','pdf'), fig.height=20, fig.width=14}

#FIG# S30

xylim = 1e-8

tmp = merge(resultsT[,c("FacDr","Drug","DrugID","DrugID2","FacDr","pval")],

            resultsA[,c("FacDr","pval")], by.x="FacDr", by.y="FacDr")

tmp$DrugName = toCaps(drugs[tmp$DrugID2, "name"])

tmp$Shape = ifelse(tmp$pval.x < xylim | tmp$pval.y < xylim, "tri", "dot")

tmp$pval.cens.x = ifelse(tmp$pval.x < xylim, xylim, tmp$pval.x)

tmp$pval.cens.y = ifelse(tmp$pval.y < xylim, xylim, tmp$pval.y)

ggplot(tmp) + geom_abline(intercept=0, slope=1, colour="hotpink") +

  geom_point(aes(x=pval.cens.x, y=pval.cens.y, shape=Shape), alpha=0.6) +

  facet_wrap(~DrugName, ncol=7) +

  scale_x_log10(labels=scientific_10, limits=c(xylim,1)) +

  scale_y_log10(labels=scientific_10, limits=c(xylim,1)) +

  theme_bw() + coord_fixed() + xlab("P-value from Student t-test") +

  ylab("P-value from ANOVA (including batch group factor)") +

  theme(axis.text.x=element_text(angle=0, hjust=1, vjust=0.5, size=rel(1.5)),

        axis.title=element_text(size=18)) + guides(shape=FALSE)

```

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.

```{r}

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.

### 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.

```{r}

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

})

```

```{r}

measurements1 = measurements

testingFactors1 = testingFactors

beelabs1 = beelabs

```

### 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.

```{r echo=FALSE}

## CREATE THE PLOTS

plotmp01 = BloodCancerMultiOmics2017:::run.ggvolcGr2(results=results1, effects=effects1,

                                screen="main", mts="IGHV", fdr=0.1, maxX=0.75,

                                maxY=NA, expY=0.05, hghBox=0.15, axisMarkY=4,

                                breaksX=c(-0.75,-0.5,-0.25,0,0.25,0.5,0.75),

                                arrowLength=0.5, Xhang=0.3, minConc=3)

plotmp02 = BloodCancerMultiOmics2017:::run.ggvolcGr2(results=results1, effects=effects1,

                                screen="main", mts="trisomy12", fdr=0.1,

                                maxX=0.75, maxY=NA, expY=0.05, hghBox=0.15,

                                axisMarkY=4,

                                breaksX=c(-0.75,-0.5,-0.25,0,0.25,0.5,0.75),

                                arrowLength=0.5, Xhang=0.3, minConc=1)

```

IGHV.

```{r volc_IGHV, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=plotmp01$IGHV$figure$width, fig.height=plotmp01$IGHV$figure$height}

#FIG# 4B

grid.draw(plotmp01$IGHV$figure$plot)

```

```{r volc_IGHV_legend, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=plotmp01$IGHV$legend$width, fig.height=plotmp01$IGHV$legend$height, out.height=300, out.width=150}

#FIG# 4B legend

grid.draw(plotmp01$IGHV$legend$plot)

```

Trisomy 12.

```{r volc_trisomy12, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=plotmp02$trisomy12$figure$width, fig.height=plotmp02$trisomy12$figure$height}

#FIG# 4C

grid.draw(plotmp02$trisomy12$figure$plot)

```

```{r volc_trisomy12_legend, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=plotmp02$trisomy12$legend$width, fig.height=plotmp02$trisomy12$legend$height, out.height=300, out.width=150}

#FIG# 4C legend

grid.draw(plotmp02$trisomy12$legend$plot)

```

## Associations of drug responses with genomic features in CLL independently of IGHV status

To assess associations between drug effects and genomic features independently of IGHV status, we performed the analyses separately within U-CLL and M-CLL samples. These analyses were only performed if 3 or more samples carried the analyzed feature within both M-CLL and U-CLL subgroups.

Find out which factors we will be testing (with threshold >2 patients in each of the four groups).

```{r}

fac2test = lapply(measurements, function(mea) {

  tf = testingFactors[["mutation"]][rownames(mea),]

  names(which(apply(tf,2,function(i) {

    if(length(table(i,tf[,1]))!=4)

      FALSE

    else

      all(table(i,tf[,1])>2)

  })))

})

```

Construct the table with drug responses.

```{r}

measurements2 = setNames(lapply(names(measurements), function(mea) {

  ig = testingFactors[["mutation"]][rownames(measurements[[mea]]),"IGHV"]

  patU = names(which(ig==1))

  patM = names(which(ig==0))

  list(U=measurements[[mea]][patU,], M=measurements[[mea]][patM,])

}), nm=names(measurements))

```

Testing.

```{r}

allresults2 = setNames(lapply(names(measurements2), function(mea) {

  list(U = testFactors(msrmnt=measurements2[[mea]]$U,

                       factors=testingFactors[["mutation"]][

                         rownames(measurements2[[mea]]$U),fac2test[[mea]]]),

  M = testFactors(msrmnt=measurements2[[mea]]$M,

                  factors=testingFactors[["mutation"]][

                    rownames(measurements2[[mea]]$M),fac2test[[mea]]]))

}), nm=names(measurements2))

```

Divide results to list of effects and list of results.

```{r}

results2 = lapply(allresults2, function(allres) {

  list(U=allres[["U"]][["test"]], M=allres[["M"]][["test"]])

})

effects2 = lapply(allresults2, function(allres) {

  list(U=allres[["U"]][["effect"]], M=allres[["M"]][["effect"]])

})

```

p-values were adjusted for multiple testing by applying the Benjamini-Hochberg procedure to joined results for M-CLL and U-CLL for each screen separately.

```{r}

results2 = lapply(results2, function(res) {

  tmp = p.adjust(c(res$U$pval,res$M$pval), method="BH")

  l = length(tmp)

  res$U$adj.pval = tmp[1:(l/2)]

  res$M$adj.pval = tmp[(l/2+1):l]

  res

})

```

```{r}

testingFactors2 = testingFactors

beelabs2 = beelabs

```

As an example we show the summary of the results for trisomy 12.

```{r echo=FALSE}

## CREATE THE PLOTS

plotmp03 = BloodCancerMultiOmics2017:::run.ggvolcGr2(results=results2$main, effects=effects2$main,

                                screen="U", mts="trisomy12", fdr=0.1,

                                maxX=0.75, maxY=7, expY=0.05, hghBox=NA,

                                axisMarkY=4,

                                breaksX=c(-0.75,-0.5,-0.25,0,0.25,0.5,0.75),

                                arrowLength=0.5, Xhang=0.3, minConc=1,

                                fixedHght=6)

plotmp04 = BloodCancerMultiOmics2017:::run.ggvolcGr2(results=results2$main, effects=effects2$main,

                                screen="M", mts="trisomy12", fdr=0.1,

                                maxX=0.75, maxY=7, expY=0.05, hghBox=NA,

                                axisMarkY=4,

                                breaksX=c(-0.75,-0.5,-0.25,0,0.25,0.5,0.75),

                                arrowLength=0.5, Xhang=0.3, minConc=1,

                                fixedHght=6)

```

Trisomy 12 - IGHV unmutated.

```{r volc_trisomy12_U, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=plotmp03$trisomy12$figure$width, fig.height=plotmp03$trisomy12$figure$height}

#FIG# S21 right

grid.draw(plotmp03$trisomy12$figure$plot)

```

```{r volc_trisomy12_U_legend, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=plotmp03$trisomy12$legend$width, fig.height=plotmp03$trisomy12$legend$height, out.height=300, out.width=150}

#FIG# S21 right legend

grid.draw(plotmp03$trisomy12$legend$plot)

```

Trisomy 12 - IGHV mutated.

```{r volc_trisomy12_M, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=plotmp04$trisomy12$figure$width, fig.height=plotmp04$trisomy12$figure$height}

#FIG# S21 left

grid.draw(plotmp04$trisomy12$figure$plot)

```

```{r volc_trisomy12_M_legend, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=plotmp04$trisomy12$legend$width, fig.height=plotmp04$trisomy12$legend$height, out.height=300, out.width=150}

#FIG# S21 left legend

grid.draw(plotmp04$trisomy12$legend$plot)

```

## Drug response dependance on cell origin of disease

We tested for drug sensitivity differences between different disease entities. The largest group, the CLL samples, was used as the baseline for these comparisons. Here, we compared drug sensitivities across studied diseases entities against all CLL samples Only groups with 3 or more data points were considered (T-PLL, AML, MZL, MCL, B-PLL, HCL, LPL and healthy donor cells hMNC). p-values were adjusted for multiple testing by applying the Benjamini-Hochberg procedure to results for each disease entity separately.

Here we prepare the data for testing the drug response dependence on cell origin of disease.

```{r}

## VIABILITIES

### main

pats = colnames(drpar)

# make the big matrix with viabilities

measureTMP = do.call(cbind,

                     lapply(list("viaraw.1","viaraw.2","viaraw.3",

                                 "viaraw.4","viaraw.5"), function(viac) {

  tmp = t(assayData(drpar)[[viac]][,pats])

  colnames(tmp) = paste(colnames(tmp),

                        conctab[colnames(tmp),

                                paste0("c",substring(viac,8,8))], sep="-")

  tmp

}))

# select diagnosis to work on

pats4diag = tapply(pats, patmeta[pats,"Diagnosis"], function(i) i)

diags = names(which(table(patmeta[pats,"Diagnosis"])>2))

diags = diags[-which(diags=="CLL")]

# there will be two lists: one with CLL and the second with other diagnosis

# (first one is passed as argument to the createObjects function)

pats4diag2 = pats4diag[diags]

# function that creates testingFactors, measurements and beelabs

createObjects = function(pats4diag1, beefix="") {

  measurements=list()

  testingFactors=list()

  # make the list for testing

  for(m in names(pats4diag1)) {

    for(n in names(pats4diag2)) {

      p1 = pats4diag1[[m]]

      p2 = pats4diag2[[n]]

      measurements[[paste(m,n,sep=".")]] = measureTMP[c(p1, p2),]

      testingFactors[[paste(m,n,sep=".")]] = setNames(c(rep(0,length(p1)),

                                                        rep(1,length(p2))),

                                                      nm=c(p1,p2))

    }

  }

  # reformat testingFactors to the df

  pats=sort(unique(c(unlist(pats4diag1),unlist(pats4diag2))))

  testingFactors = as.data.frame(

    do.call(cbind, lapply(testingFactors, function(tf) {

    setNames(tf[pats], nm=pats)

  })))

  # Labels for beeswarms

  beelabs = t(sapply(colnames(testingFactors), function(fac) {

    tmp = unlist(strsplit(fac, ".", fixed=TRUE))

    c(`0`=paste0(tmp[1], beefix),`1`=tmp[2])

  }))

  return(list(msrmts=measurements, tf=testingFactors, bl=beelabs))

}

# all CLL together

res = createObjects(pats4diag1=pats4diag["CLL"])

measurements3 = res$msrmts

testingFactors3 = res$tf

beelabs3 = res$bl

```

Testing.

```{r, results='hide'}

allresults3 = setNames(lapply(names(measurements3), function(mea) {

  tmp = data.frame(testingFactors3[,mea])

  colnames(tmp) = mea

  rownames(tmp) = rownames(testingFactors3)

  testFactors(msrmnt=measurements3[[mea]], factors=tmp)

}), nm=names(measurements3))

```

Divide results to list of effects and list of t-test results.

```{r}

results3 = lapply(allresults3, function(res) res[["test"]])

effects3 = lapply(allresults3, function(res) res[["effect"]])

```

Adjust p-values.

```{r}

results3 = lapply(results3, function(res) {

  res$adj.pval = p.adjust(res$pval, method="BH")

  res

})

```

We summarize the result as a heat map.

```{r echo=FALSE}

tmpheat = BloodCancerMultiOmics2017:::ggheat(results3, effects3)

```

```{r cll.diag, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=tmpheat$figure$width, fig.height=tmpheat$figure$height}

#FIG# S7 plot

grid.draw(tmpheat$figure$plot)

```

```{r cll.diag.legend, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=tmpheat$legend$width, fig.height=tmpheat$legend$height}

#FIG# S7 legend

grid.draw(tmpheat$legend$plot)

```

# Effect of mutation on drug response - examples

```{r}

data(drugs, lpdAll, mutCOM, conctab)

```

```{r}

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

```

Here we highlight the selection of mutation-drug response associations within the different disease subtypes.

```{r beesMutMain, fig.path=plotDir, dev=c("png", "pdf"), fig.width=8, fig.height=10, out.width=280, out.height=350}

#FIG# 4D

BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_010_2", "TP53", cs=T, y1=0, y2=1.2, custc=T)

BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_006_3", "TP53", cs=T,y1=0, y2=1.2, custc=T)

BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_063_5", "CREBBP", cs=T, y1=0, y2=1.2, custc=T)

BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_056_5", "PRPF8", cs=T, y1=0, y2=1.2, custc=T)

BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_012_5", "trisomy12", cs=F,y1=0.6, y2=1.2, custc=T)

```

```{r beesMutSupp, fig.path=plotDir, fig.width = 18, fig.height = 20, dev = c("png", "pdf")}

#FIG# S17

par(mfrow = c(3,4), mar=c(5,4.5,5,2))

BloodCancerMultiOmics2017:::beeF(diag="CLL", drug="D_159_3", mut="TP53", cs=T, y1=0, y2=1.2, custc=T)

BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_006_2", "del17p13", cs=T, y1=0, y2=1.2, custc=T)

BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_159_3", "del17p13", cs=T, y1=0, y2=1.2, custc=T)

BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_010_2", "del17p13", cs=T, y1=0, y2=1.2, custc=T)

BloodCancerMultiOmics2017:::beeF(diag="MCL", "D_006_2", "TP53", cs=T,  y1=0, y2=1.2, custc=T)

BloodCancerMultiOmics2017:::beeF(diag="MCL", "D_010_2", "TP53", cs=T,  y1=0, y2=1.2, custc=T)

BloodCancerMultiOmics2017:::beeF(diag=c("HCL", "HCL-V"), "D_012_3", "BRAF", cs=T, y1=0, y2=1.2,

            custc=F)

BloodCancerMultiOmics2017:::beeF(diag=c("HCL", "HCL-V"), "D_083_4", "BRAF", cs=T, y1=0, y2=1.2,

            custc=F)

BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_012_5", "KRAS", cs=T, y1=0.6, y2=1.2, custc=T)

BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_083_5", "KRAS", cs=T, y1=0.6, y2=1.45, custc=T)

BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_081_4", "UMODL1", cs=T,  y1=0, y2=1.2, custc=T)

BloodCancerMultiOmics2017:::beeF(diag="CLL", "D_001_4", "UMODL1", cs=T,  y1=0, y2=1.2, custc=T)

```

```{r colorbar, echo=FALSE, fig.path=plotDir, dev=c("png", "pdf"), fig.width=2, fig.height=4, out.height=200, out.width=100}

# the below function comes from the CRAN package named monogeneaGM

colorBar <-function(colpalette=NULL, min, max=-min, nticks=11, ticks=seq(min, max, len=nticks), 

 tit="") {

    if(is.null(colpalette)) {

      scale <- (length(colpalette)-1)/(max-min)

    rwb <- c("#99000D","#FB6A4A","white","#6BAED6","#084594")

    colpalette<- colorRampPalette(rwb, space="Lab")(101)

    }

    scale <- (length(colpalette)-1)/(max-min)

    plot(c(0,10), c(min,max), type="n", bty="n", xaxt="n", xlab="", yaxt="n", ylab="", main=tit)

    axis(2, ticks, las=1)

    for (i in 1:(length(colpalette)-1)) {

    y = (i-1)/scale + min

    rect(0,y,10,y+1/scale, col=colpalette[i], border=NA)

    }

}

colorBar(colorRampPalette(c('coral1','blue4'))(100), min=0, max = 1,

         ticks=c(0,0.5,1))

```

Bee swarms for pretreatment.

```{r bee-pretreatment, fig.path=plotDir, fig.width=15, fig.height=10, dev = c("png", "pdf")}

#FIG# S18

par(mfrow = c(2,3), mar=c(5,4.5,2,2)) 

BloodCancerMultiOmics2017:::beePretreatment(lpdCLL, "D_006_1:5", y1=0.2, y2=1.3, fac="TP53",

                       val=c(0,1), name="Fludarabine")

BloodCancerMultiOmics2017:::beePretreatment(lpdCLL, "D_006_1:5", y1=0.2, y2=1.3, fac="TP53",

                       val=c(0),   name="p53 wt:  Fludarabine")

BloodCancerMultiOmics2017:::beePretreatment(lpdCLL, "D_006_1:5", y1=0.2, y2=1.3, fac="TP53",

                       val=c(1),   name="p53 mut: Fludarabine")

BloodCancerMultiOmics2017:::beePretreatment(lpdCLL, "D_002_4:5", y1=0.4, y2=1.3,

                       fac="IGHV Uppsala U/M", val=c(0,1), name="Ibrutinib")

BloodCancerMultiOmics2017:::beePretreatment(lpdCLL, "D_002_4:5", y1=0.4, y2=1.3,

                       fac="IGHV Uppsala U/M", val=c(0), name="U-CLL: Ibrutinib")

BloodCancerMultiOmics2017:::beePretreatment(lpdCLL, "D_002_4:5", y1=0.4, y2=1.3,

                       fac="IGHV Uppsala U/M", val=c(1), name="M-CLL: Ibrutinib")

```

```{r, include=!exists(".standalone"), eval=!exists(".standalone")}

sessionInfo()

```

```{r, message=FALSE, warning=FALSE, include=FALSE}

rm(list=ls())

```