DrugCombi_AddOmicsData.Rmd

---
title: "CLL drug combinations ex-vivo: Add information on omics data"
author: "Britta Velten"
date: "15 April 2020"
output:
  BiocStyle::html_document:
    toc: true
    toc_depth: 3
params: 
    today: 200415
---

#Introduction

This file adds available molecular data on the patients to the combi screen objects.
```{r, echo=F}
library(knitr )
knitr::opts_chunk$set(warning=FALSE, message=FALSE)

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)
options(stringsAsFactors = FALSE)
```

```{r}
datadir <- "data"
outdir <- "out"
today <- params$today
figdir = paste0("figs", today, "/figuresAddOmicsData/")
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
```{r}
load(file.path(outdir, paste0("CLLCombiDataAfterQC_",today,".RData")))
```


# Add molecular data

## Load data

```{r}
# Sample IGHV status
patMeta <- read.csv(file.path(datadir,"patMeta.csv")) %>% as_tibble()
patMeta %<>% dplyr::rename(PatientID = Patient.ID)
IGHV <- select(patMeta, IGHV.status, PatientID)

# Patient ID conversion (to ID used in BloodCancerMultiOmics2017 package)
encPatientID <- read.csv(file.path(datadir, "encPatientID.csv"), row.names = 1)
PID <- encPatientID$PatientID
names(PID) <- encPatientID$PatientID2

# mutations data
library(BloodCancerMultiOmics2017)
data(mutCOM)
mutData <- channel(mutCOM, "binary")
rownames(mutData) <- PID[rownames(mutData)]
idx <- which(colnames(mutData) %in% c("del13q14_bi", "del13q14_mono"))
mutData <- mutData[,-idx]
colnames(mutData)[which(colnames(mutData)=="del13q14_any")] = "del13q14"

# RNAseq data
library(DESeq2)
data("dds") # RNA profiles contained in BloodCancerMultiOmics2017 (only n = 48 samples, lacking 4 samples used for the analysis in the manuscript)
colnames(dds) <- PID[colData(dds)$PatID]
# load("ddsrna_180717.RData") # updated RNA profiles for all n = 52 samples
dds

# WES data
# methylation
data("methData")
methData
colnames(methData) <- PID[colnames(methData)]
```


## Check for which patients molecular data is available
```{r}
patientsID <- unique(df4ana$PatientID)
stopifnot(length(patientsID) == 52)

OverviewMolData <- data.frame(patientID = patientsID)
OverviewMolData %<>% mutate(hasMeth = patientID %in% colnames(methData),
                            hasRNAseq = patientID %in% colnames(dds),
                            hasGenetics = patientID %in% rownames(mutData),
                            hasIGHV = patientID %in% filter(IGHV, !is.na(IGHV.status))$PatientID)
colSums(OverviewMolData[,-1])
```

## Subset data to patients from this screen
```{r}
minOcc <- 3
methData <- methData[, patientsID]
mutData <- mutData[patientsID, ]
mutData <- mutData[, colSums(exprs(mutData), na.rm=TRUE) >= minOcc]
eData <- dds[,patientsID]
patData <- filter(patMeta, PatientID %in% patientsID)
```

## Construct Mutation DF
```{r}
dfMuts <- mutData %>% exprs() %>% as.data.frame() %>%
  rownames_to_column("PatientID")
dfMuts <- left_join(dfMuts, IGHV, by = "PatientID")
dfMuts %<>% dplyr::rename(IGHV = IGHV.status)
dfMuts %<>% mutate(IGHV = ifelse(IGHV == "M",1,0))

# for testing only take  only consider clinically most relevant:
nOcc <- colSums(dfMuts[,colnames(dfMuts) != "PatientID"])
dfMuts4testing <- select(dfMuts, PatientID, TP53, IGHV) 
dfMuts4testing
```

# Patient overview including genetic background and molecular data
Creates a plot called "PatientOverview.pdf" giving an overview over the patients including infromation about somatic mutations, IGHV status, which omics are available.
```{r MutsOverview}
df4plot <- dfMuts %>%
  gather(key="mutation", value="status", -PatientID)
nOccur <- with(df4plot, tapply(status, mutation, sum, na.rm=T))
df4plot %<>% mutate(status = ifelse(status==1, "mutated", "unmutated"))
df4plot %<>% mutate(nOcc = nOccur[mutation])
df4plot %<>% mutate(mutation = factor(mutation, level = names(sort(nOccur))))
gg1 <- ggplot(df4plot, aes(x=PatientID, y=mutation, fill=status)) + geom_tile(col="black") +
  scale_fill_manual(values=c(mutated="black", unmutated="white")) + 
  theme(axis.text.x = element_blank(),
        axis.ticks = element_blank(), 
        panel.background = element_blank(),
        text = element_text(size=16),
        legend.position = "top",
        legend.key.size = unit(0.4, units="cm"),
        legend.text = element_text(size=12),
        legend.title =  element_text(size=12),
        plot.margin =  unit(c(0,0,-0.2,0), units="cm")) +
  xlab("") + ylab("") + guides(fill=guide_legend(title=""))

 # gg1 

# gg2 <- ggplot(df4plot, aes(x=mutation, y=nOcc)) +
#   geom_bar(stat="identity", fill ="cornflowerblue") +
#   theme(plot.margin =  unit(c(0,0,0,-2), units="cm"),
#         panel.background = element_blank(),
#         text = element_blank(),
#         axis.ticks = element_blank()
#         ) + coord_flip()

gg3 <- ggplot(df4plot, aes(x=PatientID, y=1, fill=PatientID))+ geom_tile(col="black") +
  scale_fill_manual(values=patcol) + 
  theme(axis.text.x = element_blank(),
        axis.ticks = element_blank(), 
        panel.background = element_blank(),
        text = element_text(size=18),
         axis.text.y = element_blank(),
        plot.margin =  unit(c(-0.2,0,0,0), units="cm")) +
  guides(fill=FALSE) + 
  xlab("Patient sample") + ylab("")

cowplot::plot_grid(gg1,gg3, nrow=2, align = "hv", axis = "lb",
                   rel_heights = c(5,0.35), rel_widths = c(5,1))
```

# Data export
```{r}
save(methData, mutData,
     eData,
     patData, dfMuts,dfMuts4testing,
     file = file.path(outdir,paste0("OmicsData",today,".RData")))
```

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

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