DrugCombi_10x10.Rmd

---
title: "CLL drug combinations ex-vivo: 10x10 screens"
author: "Britta Velten"
date: "15 April 2020"
output:
  BiocStyle::html_document:
    toc: true
    toc_depth: 3
params: 
    today: 200415
---

#Introduction

Takes raw data files for 10x10 screen as input.

```{r, echo=F}
library(knitr )
knitr::opts_chunk$set(warning=FALSE, message=FALSE)
options(stringAsFactors = 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)
library(xlsx)
source("plot_utils.R")
```

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


```{r}
load(file.path(outdir, paste0("CLLCombiDataAfterQC_",today,".RData")))
datadir10x10 <- file.path(datadir,"10x10")
```


# Get plate set-up
From Prtocol get set-up of the plate (i.e. which substance in which concentration in which well)
```{r}
FileSetUp10x10 <- as.matrix(read.xlsx(file.path(datadir10x10,
                                                "protocol_10x10_validation.xlsx"),2))

# get set-up for base/test drugs
  tableIdxBase <- which(FileSetUp10x10=="A", arr.ind = T)[1,]
  SetUpBase <- FileSetUp10x10[tableIdxBase[1]+0:15,tableIdxBase[2]+1:24]
  rownames(SetUpBase) <- FileSetUp10x10[tableIdxBase[1]+0:15,tableIdxBase[2]]
  colnames(SetUpBase) <- FileSetUp10x10[tableIdxBase[1]-1,tableIdxBase[2]+1:24]
  SetUpBase <- rownames_to_column(as.data.frame(SetUpBase), var = "WellLetter") %>%
    gather(key="WellNo", value = "base", -WellLetter)
  
  # get concentration values per column of plate
  conc_base_df <- tibble(WellNo = FileSetUp10x10[tableIdxBase[1]-1,tableIdxBase[2]+1:24],
                      base_conc = FileSetUp10x10[tableIdxBase[1]-2,tableIdxBase[2]+1:24])
  
  SetUpBase %<>% left_join(conc_base_df, by = "WellNo")
  SetUpBase %<>% mutate(base_conc = ifelse(base=="DM", NA, base_conc))
  SetUpBase <- filter(SetUpBase, !is.na(base))
  
  
# get set-up for ibrutinib
  tableIdxIbr <- which(FileSetUp10x10=="A", arr.ind = T)[2,]
  SetUpIbr <- FileSetUp10x10[tableIdxIbr[1]+0:15,tableIdxIbr[2]+1:24]
  rownames(SetUpIbr) <- FileSetUp10x10[tableIdxIbr[1]+0:15,tableIdxIbr[2]]
  colnames(SetUpIbr) <- FileSetUp10x10[tableIdxIbr[1]-1,tableIdxIbr[2]+1:24]
  SetUpIbr <- rownames_to_column(as.data.frame(SetUpIbr), var = "WellLetter") %>%
    gather(key="WellNo", value = "combi", -WellLetter)
  
  # get concentration values per row of plate
  conc_combi_df <- tibble(WellLetter = FileSetUp10x10[tableIdxIbr[1]+0:15,tableIdxIbr[2]],
                          combi_conc = FileSetUp10x10[tableIdxIbr[1]+(0:15),tableIdxIbr[2]+25])
  conc_combi_df %<>% mutate(combi_conc = ifelse(combi_conc == "Concentration (Synergie-Substanz)",
                                                NA, combi_conc))

  SetUpIbr %<>% left_join(conc_combi_df, by = "WellLetter")
  
  # set conc for DM to NA
  SetUpIbr %<>% mutate(combi_conc = ifelse(!is.na(combi),
                                           ifelse(combi=="DM", NA, combi_conc), combi_conc))

  #fill in NAs with a conc value by Ibrutinib
  SetUpIbr %<>% mutate(combi = ifelse(combi_conc %in% paste0("c", 1:10), "Ibrutinib", combi))
  SetUpIbr <- filter(SetUpIbr, !is.na(combi))
                       
  # fill in row with changing values
  SetUpIbr %<>% mutate(combi_conc = ifelse(grepl("Ibru c",combi), sub("Ibru ", "", combi), combi_conc))
  SetUpIbr %<>% mutate(combi = ifelse(grepl("Ibru c",combi), "Ibrutinib", combi))

  # merge 
  SetUp <- full_join(SetUpIbr, SetUpBase, by=c("WellLetter", "WellNo"))
  SetUp %<>% mutate(WellLetter = tolower(WellLetter))
```


# Read in values from 10x10 plates
```{r}
import10x10Data <- function(filename, SetUp){
  print(filename)
  data <- read.xlsx2(file.path(datadir10x10,filename),1, header = F, startRow=3)
  
  #info from filename
  NameComp <- strsplit(filename, "_")[[1]]
  patID <- NameComp[grep("P0",NameComp)] 
  
  # in files named PatientID.... first named drug is second in screen
  # in files named 10x10-validation..... first named drug is second in screen.
  if(grepl("^10x10-validation", filename)){
  TestDrugA <- NameComp[3] 
  TestDrugB <- sub(".xlsx","", NameComp[4])
  } else {
      TestDrugB <- NameComp[3] 
      TestDrugA <- sub(".xlsx","", NameComp[4])
  }
    
  
  # table from excel file (ignore meta data about screen)
  tableIdx <- which(data=="<>", arr.ind = T)
  stopifnot(nrow(tableIdx)==1)
  plateValues <- data[tableIdx[1]+1:11,tableIdx[2]+1:24]
  colnames(plateValues) <- 1:ncol(plateValues)
  rownames(plateValues) <- letters[2 + 1:nrow(plateValues)]  #labels start with C
  
  # turn into data.frame
  df <- rownames_to_column(plateValues, var = "WellLetter")
  df %<>% gather(key="WellNo", value="rawValue", -WellLetter)
  df %<>% mutate(rawValue = as.numeric(rawValue))
  
  # add setup info
  df %<>% left_join(SetUp,  by=c("WellLetter", "WellNo"))
  
  # fill in test drug names: test drug 1 is first of file name
  df %<>% mutate(base = ifelse(base =="DM", "DM", ifelse(base==1,TestDrugA, TestDrugB)))
  df %<>% dplyr::rename(BaseDrug = base)
  df %<>% dplyr::rename(CombiDrug = combi)

  # normalized by pure DMSO median
  dfDMSO <- filter(df, BaseDrug=="DM" , CombiDrug=="DM")
  df %<>% mutate(pureDMSOMedian = median(dfDMSO$rawValue))
  df %<>% mutate(pureDMSOMean = mean(dfDMSO$rawValue))
  df %<>% mutate(normalizedValue = rawValue/pureDMSOMedian)

  # get table of concentrations (for Ibrutinib, drug A and drug B)
  # (bug in reaad.xlsx but dilution factor is always 2 and c10 is zero, replace by hand)
  conctableIdx <- which(data=="µM"| data=="nM", arr.ind = T)
  concArray <- sapply(1:3, function(i) {
    concvec <- rep(0,10)
    names(concvec) <- paste("c", 1:10, sep="")
    concvec[1] <- as.numeric(as.character(data[conctableIdx[i,1]+0,conctableIdx[i,2]+1]))
    if(data[conctableIdx[i,1],conctableIdx[i,2]]=="nM")   concvec[1] <- concvec[1]/1000
    for(j in 2:9) concvec[j] <- concvec[j-1]/2
    concvec
  })
  colnames(concArray) <- data[conctableIdx-c(1,1)]
  concArray <- cbind(concArray,DM=rep(NA,10))
  colnames(concArray)[colnames(concArray)=="CAL-101"] <- "CAL101"
  colnames(concArray)[colnames(concArray)=="Abt-199"] <- "Abt199"
  colnames(concArray)[colnames(concArray)=="ABT-263"] <- "ABT263"
  concValueDf <- concArray %>% as.data.frame() %>% 
    rownames_to_column(var = "concId") %>% gather(key="drug", value = "concValue", -concId)
  
  # add true conc values to df
  df %<>% left_join(concValueDf, by = c("CombiDrug" = "drug", "combi_conc" = "concId"))
  df %<>% dplyr::rename(concCvalue = concValue)

  df %<>% left_join(concValueDf, by = c("BaseDrug" = "drug", "base_conc" = "concId"))
  df %<>% dplyr::rename(concBvalue = concValue)
  df %<>% mutate(PatientID = patID)

  return(df)
  }
```

```{r}
allFiles <- setdiff(list.files(file.path(datadir10x10), "*.xlsx"), list.files(file.path(datadir10x10),"protocol"))
listOfTables <- lapply(allFiles, function(file) import10x10Data(file, SetUp))
for(l in 1:length(listOfTables)) listOfTables[[l]]$PlateID <- sprintf("10x10PL%03d", l)
df10x10 <- do.call(rbind,listOfTables)
```

```{r}
# use nice drug names
df10x10 %<>% mutate(BaseDrugName = ifelse(BaseDrug == "CAL101", "Idelalisib", 
                           ifelse(BaseDrug == "Abt199", "Venetoclax",
                                  ifelse(BaseDrug == "ABT263","Navitoclax", BaseDrug))))
```

## QC
Check whether the concentration in the xls files are plausible: They are different in each experiment, consistent would be 2µM for navitoclax and venetoclax and 20µM for afatinib, idelalisib and ibrutinib in 2-fold serial dilutions, ending at 0.
```{r}
conc_BCL2 <- sapply(0:9, function(i) 2 * 0.5^i)
conc_BCL2[10] <- 0
names(conc_BCL2) <- paste0("c", 1:10)
conc_other <- sapply(0:9, function(i) 20 * 0.5^i)
conc_other[10] <- 0
names(conc_other) <- paste0("c", 1:10)
df10x10corr <- mutate(df10x10, concBvalue = ifelse(BaseDrugName %in% c("Afatinib", "Idelalisib"), conc_other[base_conc],
                                 ifelse(BaseDrugName %in% c("Navitoclax", "Venetoclax"),  conc_BCL2[base_conc], NA)))
df10x10corr <-mutate(df10x10corr, concCvalue = ifelse(CombiDrug == "Ibrutinib", conc_other[combi_conc], NA))
```

Check that concentration values for B are correct (different per patient) in the spreadsheeds. YES
```{r, fig.width=10}
gg1 <- filter(df10x10, concCvalue == 0, CombiDrug == "Ibrutinib", BaseDrugName!="DM") %>%
  ggplot(aes(x= concBvalue, y=normalizedValue, col= PatientID)) +
  geom_line() + facet_wrap(~BaseDrugName, scales = "free_x")+ggtitle("concentrations as spreadsheets")

gg2 <- filter(df10x10corr, concCvalue == 0, CombiDrug == "Ibrutinib", BaseDrugName!="DM") %>%
  ggplot(aes(x= concBvalue, y=normalizedValue, col= PatientID)) +
  geom_line() + facet_wrap(~BaseDrugName, scales = "free_x") +ggtitle("concentrations as in manuscript")

cowplot::plot_grid(gg1,gg2, ncol=2)

filter(df10x10,  CombiDrug == "Ibrutinib", BaseDrugName!="DM") %>%
  ggplot(aes(x= concBvalue, y=normalizedValue, col= PatientID)) +
  geom_line() + facet_grid(concCvalue ~BaseDrugName, scales = "free_x")
```

## Overview
```{r}
# data overview
df10x10 %>% group_by(BaseDrugName, CombiDrug) %>%
  summarize(nPats = length(unique(PatientID)),
            pats = paste(unique(PatientID), collapse =", "))

df10x10 %>% filter(base_conc== "c1") %>%
  group_by(BaseDrugName, concBvalue) %>%
  summarize(nP = length(unique(PatientID)))
```

## Data export
```{r}
save(df10x10, file = file.path(outdir, paste0("df10x10_",today,".Rdata")))

if(!dir.exists(file.path(outdir, "tables"))) dir.create(file.path(outdir, "tables"))
write.csv(df10x10, file= paste0(outdir,"/tables/10x10_",today,".csv"))
```


# Tile-Plots of normalized viabilties

```{r}
knitr::opts_chunk$set(fig.width=8)
```

## Navitoclax
```{r 10x10_Navitoclax_P0016}
plotTiles10x10(df10x10, "Navitoclax", "P0016")
```
```{r 10x10_Navitoclax_P0347}
plotTiles10x10(df10x10, "Navitoclax", "P0347")
```
```{r 10x10_Navitoclax_P0445}
plotTiles10x10(df10x10, "Navitoclax", "P0445")
```
```{r 10x10_Navitoclax_P0309}
plotTiles10x10(df10x10, "Navitoclax", "P0309")
```
```{r 10x10_Navitoclax_P0369}
plotTiles10x10(df10x10, "Navitoclax", "P0369")
```

## Afatinib
```{r 10x10_Afatinib_P0010}
plotTiles10x10(df10x10, "Afatinib", "P0010")
```
```{r 10x10_Afatinib_P0050}
plotTiles10x10(df10x10, "Afatinib", "P0050")
```
```{r 10x10_Afatinib_P0445}
plotTiles10x10(df10x10, "Afatinib", "P0069")
```
```{r 10x10_Afatinib_P0480}
plotTiles10x10(df10x10, "Afatinib", "P0480")
```
```{r 10x10_Afatinib_P0484}
plotTiles10x10(df10x10, "Afatinib", "P0484")
```

## Venetoclax
```{r 10x10_Venetoclax_P0016}
plotTiles10x10(df10x10, "Venetoclax", "P0016")
```
```{r 10x10_Venetoclax_P0347}
plotTiles10x10(df10x10, "Venetoclax", "P0347")
```
```{r 10x10_Venetoclax_P0445}
plotTiles10x10(df10x10, "Venetoclax", "P0445")
```

# Visulaization of response curves 
## Single-Dose Response Curves
```{r}
ggplot(filter(df10x10, BaseDrugName == "Navitoclax", CombiDrug == "Ibrutinib"),
       aes(x = (concBvalue), y= normalizedValue, col = PatientID,
           shape = PatientID %in% c("P0369", "P0309"))) +
  geom_point() + facet_wrap(~concCvalue) +scale_x_log10() 

ggplot(filter(df10x10, BaseDrugName == "Navitoclax", CombiDrug == "Ibrutinib"),
       aes(x = (concCvalue), y= normalizedValue, col = PatientID,
           shape = PatientID %in% c("P0369", "P0309"))) +
  geom_point() + facet_wrap(~(concBvalue*1000)) +scale_x_log10() 
```

## Dose Response curves of combination drug
```{r}
# B drug curves
filter(df10x10, BaseDrugName == "Navitoclax", CombiDrug == "Ibrutinib", concCvalue==0) %>%
ggplot(aes(x = (concBvalue), y= normalizedValue,
           col = PatientID, shape = PatientID %in% c("P0369", "P0309"))) +
  geom_point() + facet_wrap(~(concCvalue*1000)) + scale_x_log10() + geom_line()

filter(df10x10, BaseDrugName == "Venetoclax", CombiDrug == "Ibrutinib", concCvalue==0) %>%
ggplot(aes(x = (concBvalue), y= normalizedValue,
           col = PatientID, shape = PatientID %in% c("P0369", "P0309"))) +
  geom_point() + facet_wrap(~(concCvalue*1000)) + scale_x_log10() + geom_line()
ggplot(filter(df10x10, BaseDrugName == "Afatinib", CombiDrug == "Ibrutinib", concCvalue==0),
       aes(x = (concBvalue), y= normalizedValue, col = PatientID,
           shape = PatientID %in% c("P0369", "P0309"))) +
  geom_point() + facet_wrap(~(concCvalue*1000)) +scale_x_log10() +geom_line()
```

## Dose Response curves of Ibrutinib
```{r}
# C Drug curves
ggplot(filter(df10x10, BaseDrugName == "Navitoclax", CombiDrug == "Ibrutinib", concBvalue==0),
       aes(x = (concCvalue), y= normalizedValue, col = PatientID,
           shape = PatientID %in% c("P0369", "P0309"))) +
  geom_point() + facet_wrap(~(concBvalue*1000)) +scale_x_log10() +geom_line()

ggplot(filter(df10x10, BaseDrugName == "Venetoclax", CombiDrug == "Ibrutinib", concBvalue==0),
       aes(x = (concCvalue), y= normalizedValue, col = PatientID,
           shape = PatientID %in% c("P0369", "P0309"))) +
  geom_point() + facet_wrap(~(concBvalue*1000)) +scale_x_log10() +geom_line()

ggplot(filter(df10x10, BaseDrugName == "Afatinib", CombiDrug == "Ibrutinib", concBvalue==0),
       aes(x = (concCvalue), y= normalizedValue, col = PatientID,
           shape = PatientID %in% c("P0369", "P0309"))) +
  geom_point() + facet_wrap(~(concBvalue*1000)) +scale_x_log10() +geom_line()
```

# Assesing synergy using SynergyFinder
To assess synergy we use the ZIP score proposed (here)[https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4759128/] and implemented in the R package synergyfinder. The results are overall similar to using the Loewe score. Note that the implementation of the Loewe score provided in the package uses the maximal effect of single compounds at the additive concentration value if the Loewe equation cannot be solved resulting in artifacts with strong antagonism, when the dose-response relationship differs between two drugs. These cases are should be replaced by the NA or based on the value that is closest to solving the Loewe equation instead if the Loewe score was used.

```{r}
library("synergyfinder")
```

```{r ZIP_surface}
dfZIP <- plotCITiles(df10x10, "ZIP")
```

# Summarize scores
## Navitoclax

```{r summary_Navitoclax_ZIP, fig.width = 10 }
ggNavi <- plotSummaryCI10x10(dfZIP, "Navitoclax", type = "row")
ggNavi
```

## Afatinib

```{r summary_Afatinib_ZIP, fig.width = 10 }
ggAfa <- plotSummaryCI10x10(dfZIP, "Afatinib", type ="row")
ggAfa
```

## Venetoclax

### ZIP score
```{r summary_Venetoclax_ZIP, fig.width = 10 }
ggVen <- plotSummaryCI10x10(dfZIP, "Venetoclax", type ="row")
ggVen
```

## Idelalisib

```{r summary_Idelalisib_ZIP, fig.height = 8 }
ggIdel <- plotSummaryCI10x10(dfZIP, "Idelalisib")
ggIdel
```


```{r}
sessionInfo()
```