---

title: "Part 5"

output:

  BiocStyle::html_document

---

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

library("BloodCancerMultiOmics2017")

library("Biobase")

library("dplyr")

library("tidyr")

library("broom")

library("ggplot2")

library("grid")

library("gridExtra")

library("reshape2")

library("foreach")

library("doParallel")

library("scales")

library("knitr")

registerDoParallel()

```

```{r echo=FALSE}

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

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

```

```{r}

options(stringsAsFactors=FALSE)

```

# Associations of drug responses with mutations in CLL (IGHV not included)

In this part, we use both gene mutations and chromosome aberrations to test for gene-drug response associations. In contrast to the analysis done previously, we exclude IGHV status from testing. Additionally, we use information on patient treatment status to account for its effect on drug response screening.

## Additional functions

Accessor functions:

```{r}

# get drug responsee data

get.drugresp <- function(lpd) {

  drugresp = t(Biobase::exprs(lpd[fData(lpd)$type == 'viab'])) %>%

    dplyr::tbl_df() %>% dplyr::select(-ends_with(":5")) %>%

    dplyr::mutate(ID = colnames(lpd)) %>%

    tidyr::gather(drugconc, viab, -ID) %>%

    dplyr::mutate(drug = drugs[substring(drugconc, 1, 5), "name"],

           conc = sub("^D_([0-9]+_)", "", drugconc)) %>%

    dplyr::mutate(conc = as.integer(gsub("D_CHK_", "", conc)))

  drugresp

}

# extract mutations and IGHV status

get.somatic <- function(lpd) {

  somatic = t(Biobase::exprs(lpd[Biobase::fData(lpd)$type == 'gen' | 

                                   Biobase::fData(lpd)$type == 'IGHV']))

  ## rename IGHV Uppsala to 'IGHV' (simply)

  colnames(somatic)[grep("IGHV", colnames(somatic))] = "IGHV"

  ## at least 3 patients should have this mutation

  min.samples = which(Matrix::colSums(somatic, na.rm = T) > 2)

  somatic = dplyr::tbl_df(somatic[, min.samples]) %>%

    dplyr::select(-one_of("del13q14_bi", "del13q14_mono", 

                          "Chromothripsis", "RP11-766F14.2")) %>%

    dplyr::rename(del13q14 = del13q14_any) %>% 

    dplyr::mutate(ID = colnames(lpd)) %>%

    tidyr::gather(mutation, mut.value, -ID)

  somatic

}

```

Define the ggplot themes

```{r ggTheme}

t1<-theme(                              

  plot.background = element_blank(), 

  panel.grid.major = element_line(),

  panel.grid.major.x = element_line(linetype = "dotted", colour = "grey"),

  panel.grid.minor = element_blank(), 

  panel.border = element_blank(), 

  panel.background = element_blank(),

  axis.line = element_line(size=.4),

  axis.line.x = element_line(),

  axis.line.y = element_line(),

  axis.text.x  = element_text(angle=90, size=12, 

                              face="bold", hjust = 1, vjust = 0.4),

  axis.text.y = element_text(size = 14),

  axis.ticks.x = element_line(linetype = "dotted"),

  axis.ticks.length = unit(0.3,"cm"),

  axis.title.x = element_text(face="bold", size=16), 

  axis.title.y = element_text(face="bold", size=20),

  plot.title = element_text(face="bold", size=16, hjust = 0.5)

)

## theme for the legend

t.leg <-  theme(legend.title = element_text(face='bold', 

                                            hjust = 1, size=11),

                legend.position = c(0, 0.76),

                legend.key = element_blank(),

                legend.text = element_text(size=12),

                legend.background = element_rect(color = "black"))

```

Define the main color palette: 

```{r colorPalette}

colors= c("#015872","#3A9C94","#99977D","#ffbf00","#5991C7","#99cc00",

          "#D5A370","#801416","#B2221C","#ff5050","#33bbff","#5c5cd6",

          "#E394BB","#0066ff","#C0C0C0")

```

Get pretreatment status:

```{r}

get.pretreat <- function(patmeta, lpd) {

  patmeta = patmeta[rownames(patmeta) %in% colnames(lpd),]

  data.frame(ID=rownames(patmeta), pretreat=!patmeta$IC50beforeTreatment) %>% 

    mutate(pretreat = as.factor(pretreat))

}

```

Merge drug response, pretreatment information and somatic mutation data sets

```{r}

make.dr <- function(resp, features, patmeta, lpd) {

  treat = get.pretreat(patmeta, lpd)

  dr = full_join(resp, features) %>% 

    inner_join(treat) 

}

```

Summarize viabilities using Tukey's medpolish

```{r}

get.medp <- function(drugresp) {

  tab = drugresp %>% group_by(drug, conc) %>% 

    do(data.frame(v = .$viab, ID = .$ID)) %>% spread(ID, v)

  med.p = foreach(n=unique(tab$drug), .combine = cbind) %dopar% {

    tb = filter(tab, drug == n) %>% ungroup() %>% dplyr::select(-(drug:conc)) %>% 

      as.matrix %>% `rownames<-`(1:5)

    mdp = stats::medpolish(tb)

    df = as.data.frame(mdp$col) + mdp$overall

    colnames(df) <- n

    df

  }

  medp.viab = dplyr::tbl_df(med.p) %>% dplyr::mutate(ID = rownames(med.p)) %>%

    tidyr::gather(drug, viab, -ID) 

  medp.viab

}

```

Process labels for the legend:

```{r}

get.labels <- function(pvals) {

  lev = levels(factor(pvals$mutation))

  lev = gsub("^(gain)([0-9]+)([a-z][0-9]+)$", "\\1(\\2)(\\3)", lev)

  lev =  gsub("^(del)([0-9]+)([a-z].+)$", "\\1(\\2)(\\3)", lev)

  lev = gsub("trisomy12", "trisomy 12", lev)

  lev

}

```

Get order of mutations

```{r}

get.mutation.order <- function(lev) {

  ord = c("trisomy 12", "TP53",

          "del(11)(q22.3)", "del(13)(q14)",

          "del(17)(p13)",

          "gain(8)(q24)",

          "BRAF", "CREBBP", "PRPF8",

          "KLHL6", "NRAS", "ABI3BP", "UMODL1")

  mut.order = c(match(ord, lev),

                grep("Other", lev), grep("Below", lev))

  mut.order

}

```

Group drugs by pathway/target

```{r}

get.drug.order <- function(pvals, drugs) {

  ## determine drug order by column sums of log-p values

  dr.order = pvals %>% 

    mutate(logp = -log10(p.value)) %>% 

    group_by(drug) %>% summarise(logsum = sum(logp)) 

  dr.order = inner_join(dr.order, pvals %>%

                          group_by(drug) %>% 

                          summarise(n = length(unique(mutation)))) %>% 

    arrange(desc(n), desc(logsum))

  dr.order = inner_join(dr.order, drugs %>% rename(drug = name))

  dr.order = left_join(dr.order, dr.order %>% 

                         group_by(`target_category`) ) %>%

    arrange(`target_category`, drug) %>%

    filter(! `target_category` %in% c("ALK", "Angiogenesis", "Other")) %>%

    filter(!is.na(`target_category`))

  dr.order

}

```

Add pathway annotations for selected drug classes

```{r}

make.annot <- function(g, dr.order) {

  # make a color palette for drug pathways

  drug.class = c("#273649", "#647184", "#B1B2C8",

                 "#A7755D", "#5D2E1C", "#38201C")

  pathways = c("BH3","B-cell receptor","DNA damage",

               "MAPK", "PI3K", "Reactive oxygen species")

  names(pathways) = c("BH3", "BCR inhibitors", "DNA damage",

                      "MAPK", "PI3K", "ROS")

  for (i in 1:6) {

    prange = grep(pathways[i], dr.order$`target_category`)

    path.grob <- grobTree(rectGrob(gp=gpar(fill=drug.class[i])),

                          textGrob(names(pathways)[i], 

                                   gp = gpar(cex =0.8, col = "white")))

    g = g + 

      annotation_custom(path.grob, 

                        xmin = min(prange) -0.25 - 0.1 * ifelse(i == 2, 1, 0), 

                        xmax = max(prange) + 0.25 + 0.1 * ifelse(i == 2, 1, 0), 

                        ymin = -0.52, ymax = -0.2)

  }

  g

}

```

Define a function for `glegend`

```{r}

g_legend<-function(a.gplot){

  tmp <- ggplot_gtable(ggplot_build(a.gplot))

  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")

  legend <- tmp$grobs[[leg]]

  legend

} ## end define

```

## Data setup

Load the data.

```{r}

data(list=c("conctab", "drugs", "lpdAll", "patmeta"))

```

Get drug response data.

```{r preprocesslpd}

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

## extract viability data for 5 different concentrations

drugresp = get.drugresp(lpdCLL)

```

Get somatic mutations and structural variants.

```{r preprocessMuts}

## extract somatic variants

somatic = get.somatic(lpdCLL) %>% 

  mutate(mut.value = as.factor(mut.value))

```

## Test for drug-gene associations

Summarize drug response using median polish.

```{r medPolish, warning=FALSE, results='hide'}

## compute median polish patient effects and recalculate p-values

medp.viab = get.medp(drugresp)

dr = make.dr(medp.viab, somatic, patmeta, lpdCLL)

```

Calculate $p$ values and FDR (10%).

```{r pvalsAndFDR}

pvals = dr %>% group_by(drug, mutation) %>%

  do(tidy(t.test(viab ~ mut.value, data = ., var.equal = T))) %>%

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

```

Remove unnecessary mutations and bad drugs.

```{r subsetPvals}

# selected mutations

select.mutations = c("trisomy12", "TP53",

        "del11q22.3", "del13q14",

        "del17p13",

        "gain8q24",

        "BRAF", "CREBBP", "PRPF8",

        "KLHL6", "NRAS", "ABI3BP", "UMODL1")

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

```

Reshape names of genomic rearrangements.

```{r renameMuts}

## order of mutations

lev = get.labels(pvals)

folge = get.mutation.order(lev)

```

Set order of drugs.

```{r pvalsForGGplot}

drugs = drugs[,c("name", "target_category")]

# get the drug order

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

```

## Plot results

### Main Figure

Function for generating the figure.

```{r}

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

}

```

```{r pvalsMain, fig.path=plotDir, dev=c("png", "pdf"), fig.width=14, fig.height=10}

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

```

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

```{r}

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

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

```

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

```{r pvalsSupp, fig.path=plotDir, dev=c("png", "pdf"), fig.width=14, fig.height=10}

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

```

## Comparison of $P$-Values

```{r}

pvals.main = dr %>% group_by(drug, mutation) %>%

  do(tidy(t.test(viab ~ mut.value, data = ., var.equal = T))) %>%

  dplyr::select(drug, mutation, p.value)

p.main.adj = p.adjust(pvals.main$p.value, method = "BH")

fdr.main = max(pvals.main$p.value[which(p.main.adj <= fd.thresh/100)])

pvals.main = filter(pvals.main, mutation != "IGHV") %>%

             rename(p.main = p.value)

## lm(viab ~ mutation + pretreatment.status)

pvals.sup = dr %>% group_by(drug, mutation) %>%

  do(tidy(lm(viab ~ mut.value + pretreat, data = .))) %>%

  filter(term == 'mut.value1') %>%

  dplyr::select(drug, mutation, p.value)

p.sup.adj = p.adjust(pvals.sup$p.value, method = "BH")

fdr.sup = max(pvals.sup$p.value[which(p.sup.adj <= fd.thresh/100)])

pvals.sup = filter(pvals.sup, mutation != "IGHV") %>%

  rename(p.sup = p.value)

pvals = inner_join(pvals.main, pvals.sup)

pvals = mutate(pvals, signif = ifelse(p.main > fdr.main, 

                                      ifelse(p.sup > fdr.sup, 

                                             "Below 10% FDR in both models", 

                                             "Significant with pretreatment accounted"), 

                                      ifelse(p.sup > fdr.sup, 

                                             "Significant without pretreatment in the model", 

                                             "Significant in both models")))

t2<-theme(                              

  plot.background = element_blank(), 

  panel.grid.major = element_line(),

  panel.grid.major.x = element_line(),

  panel.grid.minor = element_blank(), 

  panel.border = element_blank(), 

  panel.background = element_blank(),

  axis.line = element_line(size=.4),

  axis.line.x = element_line(),

  axis.line.y = element_line(),

  axis.text.x  = element_text(size=12),

  axis.text.y = element_text(size = 12),

  axis.title.x = element_text(face="bold", size=12), 

  axis.title.y = element_text(face="bold", size=12),

  legend.title = element_text(face='bold', 

                                            hjust = 1, size=10),

  legend.position = c(0.78, 0.11),

  legend.key = element_blank(),

  legend.text = element_text(size=10),

  legend.background = element_rect(color = "black")

)

```

```{r pvalComparisonScatterplot, fig.path=plotDir, dev=c("png", "pdf"), fig.width=10, fig.height=7}

#FIG# S19

ggplot(pvals, aes(-log10(p.main), -log10(p.sup), colour = factor(signif))) + 

  geom_point() + t2 + labs(x = expression(paste(-log[10], "p, pretreatment not considered", sep = "")),

                                   y = expression(paste(-log[10], "p, accounting for pretreatment", sep = ""))) +

  coord_fixed() +

  scale_x_continuous(breaks = seq(0,9,by = 3)) + 

  scale_y_continuous(breaks = seq(0,9,by = 3)) + 

  geom_abline(intercept = 0, slope = 1, linetype = "dashed") + 

  scale_color_manual(name = "Statistical Significance",

                     values = c("#F1BB7B","#669999", "#FD6467", "#5B1A18"))

```

What are the drug-mutation pairs that are significant only in one or another model (i.e. only without pretreatment or with pretreatment included)? 

```{r}

signif.in.one = filter(pvals, 

                       signif %in% c("Significant with pretreatment accounted",

                            "Significant without pretreatment in the model")) %>%

        arrange(signif)

kable(signif.in.one, digits = 4,

      align = c("l", "l", "c", "c", "c"),

      col.names = c("Drug", "Mutation", "P-value (Main)",

                    "P-value (Supplement)", "Statistical significance"),

      format.args = list(width = 14))

```

Produce LaTeX output for the Supplement:

```{r, comment=NA, eval=FALSE}

print(kable(signif.in.one, format = "latex", digits = 4,

       align = c("l", "l", "c", "c", "c"),

      col.names = c("Drug", "Mutation", "P-value (Main)",

                    "P-value (Supplement)", "Statistical significance")))

```

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

sessionInfo()

```

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

rm(list=ls())

```