DrugCombi_BaseAlone.Rmd

---
title: "CLL drug combinations ex-vivo: Analyse effect of base drugs alone"
author: "Britta Velten"
date: "15 April 2020"
output:
  BiocStyle::html_document:
    toc: true
    toc_depth: 3
params: 
    today: 200415
---

#Introduction

This file analyses the data from the combination screen for combination of library drugs with DMSO. (Effect of base (library) drug alone)

```{r, echo=F}
library(knitr )
knitr::opts_chunk$set(warning=FALSE, message=FALSE)
```

```{r, echo=F}
library(ggplot2)
library(RColorBrewer)
library(Biobase)
library(abind)
library(grid)
library(gtable)
library(reshape2)
library(gridExtra)
require(pracma)
library(reshape2)
library(dplyr)
library(pheatmap)
library(magrittr)
library(lattice)
library(tidyverse)
source("plot_utils.R")
options(stringsAsFactors = FALSE)
```

```{r}
datadir <- "data"
outdir <- "out"
today <- params$today
figdir = paste0("figs", today, "/figuresBaseAlone/")
if(!dir.exists(figdir)) dir.create(figdir)
knitr::opts_chunk$set(dev = c("png", "pdf"), fig.path = figdir)
```


# Data Import:
Load processed Data from `DrugCombi_QC.Rmd` and `DrugCombi_AddOmics.Rmd` (moelcular data on patient samples)
```{r}
load(file.path(outdir, paste0("CLLCombiDataAfterQC_",today,".RData")))
load(file.path(outdir, paste0("OmicsData",today,".RData")))
```

```{r}
# discard outliers for further analyses, only use averaged values for replicates
df4ana <- df4anaAvreplicatesRmOutliers
rm(df4anaAvreplicatesRmOutliers, df4anaAvreplicates)
```

# Construct dataframe for base effect alone
Note that each experiment consists of 3 combination partner plus DMSO. Therefore, all patient samples with more than 3 combinaiton partners will have several effect values of the DMSO plate from each experiment. In the following we will average those values. For the combination studies comparisons are made within the experiments taken the value from the corresponding DMSO plate.

```{r}
df_effectB <- df4ana %>% select(starts_with("BDrug"), PatientID, effectB)

# this contains multiple values for most patient, drug, concentration pairs:
nrow(df_effectB)
52*32*5
df <- df_effectB %>%
  group_by(BDrugName, BDrugID, BDrugConcId,BDrugConc,PatientID) %>%
  summarize(n_effectB = length(effectB)) 
table(df$n_effectB)

#summarize by mean
df_effectB %<>%  group_by(BDrugName, BDrugID, BDrugConcId,BDrugConc,PatientID) %>%
  summarize(effectB = mean(effectB))%>%
  ungroup()
nrow(df_effectB)
```

# Test for assciations of base effect to mutations
```{r}
df_effectB_vs_muts <- left_join(df_effectB,
                                gather(dfMuts4testing, key="mutation", value="status", -PatientID),
                                by="PatientID")
#remove duplicates, as info is contained in each combi drug
df_effectB_vs_muts <- df_effectB_vs_muts[!duplicated(df_effectB_vs_muts),] 
df_pvals_effectB_vs_muts <- group_by(df_effectB_vs_muts, BDrugName, BDrugID, BDrugConc, BDrugConcId, mutation) %>%
  summarize(pval = t.test(effectB ~ status, var.equal= TRUE)$p.value,
            mean.diff = diff(t.test(effectB ~ status, var.equal= TRUE)$estimate))
df_pvals_effectB_vs_muts$pval.adj <- p.adjust(df_pvals_effectB_vs_muts$pval, method="BH")
nrow(df_pvals_effectB_vs_muts)
hist(df_pvals_effectB_vs_muts$pval)
```

```{r}
if(!dir.exists(file.path(outdir, "tables"))) dir.create(file.path(outdir, "tables"))
write.csv(df_pvals_effectB_vs_muts, file= paste0(outdir,"/tables/pvalues_mut_BaseDrugs",today,".csv"))
```


```{r}
# FDR threshold
alpha <- 0.05
```

## Volcanoe plot
```{r volcanoe_base, fig.height=4, fig.width=8}
# take only most significant concentration for plotting (FDR control is on level of single concentrations)
df_pvals_max <- df_pvals_effectB_vs_muts %>%
  group_by(BDrugName, BDrugID, mutation) %>%
  summarize(idx = which.min(pval), pval = pval[idx],
            pval.adj= pval.adj[idx], mean.diff = mean.diff[idx],
            conc = BDrugConcId[idx])
df_pvals_max %<>% mutate(sig = pval.adj<alpha, col = mean.diff>0)

ggplot(df_pvals_max, aes(x=mean.diff, y=-log10(pval.adj), col=col)) +
  geom_point(alpha=0.7) +
  ggrepel::geom_label_repel(aes(label=ifelse(sig, BDrugName, "")),
                            label.size= NA, fill = NA, size=4.5) +
  scale_color_manual(values= c("deeppink", "navy")) +
  facet_wrap(~mutation)+ guides(col=FALSE) +
  theme_bw(base_size = 14) +
  geom_hline(yintercept = -log10(alpha), lty= "dashed", alpha=0.3) +
  xlab ("Viability effect") +
  theme(strip.background = element_blank(),
        strip.text = element_text(size=14))
```

## List of significant pharmacogenomic interactions
```{r}
df_sign <- df_pvals_effectB_vs_muts %>% filter(pval.adj < alpha)
df_nconc <- df_sign %>% group_by(BDrugName, mutation) %>% summarize(nconc =n())

# TP53
df_signTP53 <- df_sign %>% filter(mutation == "TP53")
paste(unique(df_signTP53$BDrugName), collapse = ", ")

# IGHV
df_signIGHV <- df_sign %>% filter(mutation == "IGHV")
paste(unique(df_signIGHV$BDrugName), collapse = ", ")
```

## Boxplots and Curves
Draw boxplots for TP53 and IGHV including all drugs that have significant concentrations at the specified FDR threshold.
```{r boxplot_base_TP53, fig.width=8, fig.height=3}
plotBaseBoxplots(df_effectB_vs_muts, unique(df_signTP53$BDrugName), gene = "TP53")
```
```{r curves_base_TP53, fig.width=7, fig.height=2.8}
plotBaseResponseCurves(df_effectB_vs_muts, unique(df_signTP53$BDrugName), gene = "TP53")
```

```{r, boxplot_base_IGHV, fig.width=8, fig.height=7.5}
selDrugs <- c("Ibrutinib", "R406", "Spebrutinib", "Dasatinib",
              "AZD7762", "PF 477736", "Navitoclax", "Venetoclax")
stopifnot(all(selDrugs %in% df_signIGHV$BDrugName))
plotBaseBoxplots(df_effectB_vs_muts, selDrugs, gene = "IGHV")
```

```{r curves_base_IGHV, fig.width=7, fig.height =6.5}
plotBaseResponseCurves(df_effectB_vs_muts, selDrugs, gene = "IGHV")
```

# Test for associaiton to methylation data
```{r methylation_cluster_base, fig.height=8, fig.width=10}
df_base_vs_meth <- left_join(df_effectB, patData, by = "PatientID")
df_pvals_effB_vs_meth <- df_base_vs_meth %>% group_by(BDrugName, BDrugConcId) %>%
  dplyr::summarise(pval = anova(lm(effectB ~ Methylation_Cluster))$'Pr(>F)'[1])

# adjust for multiple testing and filter to FDR < alpha
df_pvals_effB_vs_meth$pval.adj <- p.adjust(df_pvals_effB_vs_meth$pval, method = "BH")
hist(df_pvals_effB_vs_meth$pval)
dfsig_effB_vs_meth <- df_pvals_effB_vs_meth %>% filter(pval.adj < alpha)
unique(dfsig_effB_vs_meth$BDrugName)
length(unique(dfsig_effB_vs_meth$BDrugName))

# plot signifcant base drugs
sel_drugs <- dfsig_effB_vs_meth %>% select(BDrugName)

# plot curves
df4plot <- df_base_vs_meth %>%
  filter(BDrugName %in% sel_drugs$BDrugName, !is.na(Methylation_Cluster)) %>% 
    mutate(BDrugConcId = factor(BDrugConcId, levels = paste0("c", 5:1))) %>%
    mutate(BDrugConc = factor(round(1000*BDrugConc, 0)))
  
  gg <- ggplot(df4plot, aes(x=BDrugConc, y=effectB, group=Methylation_Cluster)) +
    stat_summary(fun.data = "mean_se", aes(col=factor(Methylation_Cluster)), geom="line") +
    facet_wrap(~ BDrugName, ncol=3, scales = "free_x") +
    stat_summary(fun.data = "mean_se", aes(col=factor(Methylation_Cluster)),
                 geom="errorbar", width=0.2, fun.args = list(mult =2)) +
    ggpubr::stat_compare_means(method = "t.test", aes(group = Methylation_Cluster, x = BDrugConc,
                                                      label =  ..p.signif..),
                               hide.ns = TRUE, label.y = 1.3) +
    guides(col = guide_legend(title="Methylation Cluster")) +
    xlab("Concentration (nM)") + ylab("viability")+theme_bw(base_size = 14) +
    theme(strip.background = element_blank(),
          axis.text.x = element_text(angle=90, vjust=1, hjust =1)) +
    scale_color_manual(values= c("LP" = "navy", "IP" = "orange", "HP"="darkred")) +
    ylim(c(0,1.4))
  
  gg
```



# Heatmap on base viability values
```{r}
df_effectB %<>% mutate(BDrugNameConc = paste(BDrugName, BDrugConcId, sep="_"))
effectB_mat <- df_effectB %>% 
  select(BDrugNameConc, PatientID, effectB)  %>%
  spread(key=BDrugNameConc, value=effectB) %>% 
  column_to_rownames("PatientID") %>%
  as.matrix()
```

## Base-DrugxDrug
```{r heatmap_base_drug_drug}
pheatmap(cor(effectB_mat, use = "complete.obs"),
         show_rownames = FALSE,treeheight_row = 10, treeheight_col = 10)
```


Average over concentrations
```{r heatmap_base_drug_drug_av, fig.height=7, fig.width=6.1}
effectB_mat_av <- df_effectB %>% 
  select(BDrugName, PatientID, effectB)  %>%
  group_by(BDrugName,PatientID) %>%
  summarise(effectB = mean(effectB, na.rm=TRUE)) %>%
  ungroup() %>%
  spread(key=BDrugName, value=effectB) %>% 
  column_to_rownames("PatientID") %>%
  as.matrix()
pheatmap(cor(effectB_mat_av, use = "complete.obs"), show_rownames = FALSE,treeheight_row = 10,
         treeheight_col = 10, fontsize = 14, breaks = seq(-1,1,0.01),
         color = colorRampPalette(rev(brewer.pal(n=7, name = "RdBu")))(200))
```

## Patient x Drugs
```{r heatmap_base_patient_drug}
# outlying values cut off at 1.4
effectB_mat[is.na(effectB_mat)] <- 1.4

# add genetic background
dfanno <- dfMuts4testing %>% column_to_rownames("PatientID") %>% select(which(colSums(.) > 10))
dfanno <- as.data.frame(ifelse(dfanno==1, "mut", "wt"))

# set colors for mutation status
cols_mut        <- c("white", "black")
names(cols_mut) <- c("wt", "mut")
anno_colors <- lapply(colnames(dfanno), function(x) cols_mut)
names(anno_colors) <- colnames(dfanno)

callbackIGHV = function(hc, mat){
      ighv <- dfMuts$IGHV
      names(ighv) <- dfMuts$PatientID
      dend = reorder(as.dendrogram(hc), wts = ighv[rownames(mat)])
      as.hclust(dend)
}

 pheatmap(effectB_mat, na_col="gray", clustering_distance_rows="euclidean",
         clustering_distance_cols = "correlation",
         color = colorRampPalette(rev(brewer.pal(n = 7, name ="RdYlBu")))(140),
         breaks=seq(0,filter_th,0.01), show_rownames = FALSE, show_colnames = FALSE,
         treeheight_row = 10, treeheight_col = 10,
         annotation_row = dfanno, annotation_colors =anno_colors, annotation_legend = FALSE,
         clustering_callback = callbackIGHV)
```

Average over concentrations
```{r heatmap_base_patient_drug_av}
effectB_mat_av <- df_effectB %>% 
  select(BDrugName, PatientID, effectB)  %>%
  group_by(BDrugName,PatientID) %>%
  summarise(effectB = mean(effectB, na.rm=TRUE)) %>%
  ungroup() %>%
  spread(key=BDrugName, value=effectB) %>% 
  column_to_rownames("PatientID") %>%
  as.matrix()

# add genetic background
dfanno <- dfMuts4testing %>% column_to_rownames("PatientID") %>% select(which(colSums(.) > 10))
dfanno <- as.data.frame(ifelse(dfanno==1, "mut", "wt"))

# set colors for mutation status
cols_mut        <- c("white", "black")
names(cols_mut) <- c("wt", "mut")
anno_colors <- lapply(colnames(dfanno), function(x) cols_mut)
names(anno_colors) <- colnames(dfanno)

 # pheatmap(effectB_mat_av, na_col="gray", clustering_distance_rows="euclidean",
 #         clustering_distance_cols = "correlation",
 #         color = colorRampPalette(rev(brewer.pal(n = 7, name ="RdYlBu")))(140),
 #         breaks=seq(0,filter_th,0.01), show_rownames = FALSE, show_colnames = TRUE,
 #         treeheight_row = 10, treeheight_col = 10, annotation_row = dfanno, annotation_colors =anno_colors,
 #         annotation_legend = FALSE, fontsize = 12,  clustering_callback = callbackIGHV)
 pheatmap(effectB_mat_av, na_col="gray", clustering_distance_rows="euclidean",
         clustering_distance_cols = "correlation",
    show_rownames = FALSE, show_colnames = TRUE,
         treeheight_row = 10, treeheight_col = 10, annotation_row = dfanno, annotation_colors =anno_colors,
         annotation_legend = FALSE, fontsize = 12,  clustering_callback = callbackIGHV,
         breaks = seq(0, filter_th, 0.01),
         color = colorRampPalette((brewer.pal(n=7, name = "PRGn")))(200)[1:140])
```

columns-wise Z-score
```{r}
mat <- apply(effectB_mat_av,2,scale)
rownames(mat)  <- rownames(effectB_mat_av)
 pheatmap(mat, na_col="gray", show_rownames = FALSE, show_colnames = TRUE,
         treeheight_row = 10, treeheight_col = 10, annotation_row = dfanno, annotation_colors =anno_colors,
         annotation_legend = FALSE, fontsize = 12,  clustering_callback = callbackIGHV)
```

# Single correlation values
```{r}
# navitoclax - venetoclax
cor(as.numeric(effectB_mat[,grep("Navitoclax", colnames(effectB_mat))]), as.numeric(effectB_mat[,grep("Venetoclax", colnames(effectB_mat))]))

cor(as.numeric(effectB_mat[,grep("YM155", colnames(effectB_mat))]),
    as.numeric(effectB_mat[,grep("SNS-032", colnames(effectB_mat))]))
```

#SessionInfo
```{r}
sessionInfo()
```