---
title: "CLL example: predicting drug responses using 3 omics (drug/mRNA/meth) on original scales"
author: "Britta Velten"
date: "11/05/2019"
output:
  BiocStyle::html_document:
    toc: true
---

```{r, message=FALSE, warning=FALSE}
library(ggplot2)
library(dplyr)
library(magrittr)
library(plyr)
library(cowplot)
library(reshape2)
options(stringsAsFactors = FALSE)
library(graper)
library(tidyverse)
library("wesanderson")
library(ggbeeswarm)
source("../util_defs.R") # contains color schemes and defines which methods to include
```

# Perparations
Set input/output paths.
```{r}
datadir <- "data"
outdir <- "2018-11-06"
knitr::opts_chunk$set(fig.path = "figs/", dev = c('png',"pdf"))
```

Load data used for fitting.
```{r}
load(file.path(datadir, "dataCLL.RData"))
dim(data$X)
```

Load summarised results.
```{r}
load(file.path(outdir,"result_CLL.RData"))
```

Colors for omics
```{r}
cols4groups <- c(wes_palette("GrandBudapest1"),
                 wes_palette("GrandBudapest2")[1])
names(cols4groups) <- c("Drugs", "Methylation","mRNA")
```

# Make plots
## RMSE
Compare prediciton performance between the methods in terms of root mean squared error.
```{r RMSE, fig.widht=8, fig.height=6}
df_RMSE <- melt(sapply(resultList, function(l) l$RMSE),
                varnames = c("method", "run"), value.name="RMSE")
df_RMSE %<>% mutate(method = make_nicenames(method))
df_RMSE %<>% mutate(method_type=ifelse(method%in% methods2compare_sparse,
                                       "sparse", "dense"))  
df_RMSE %<>% filter(method %in% methods2compare_sparse |
                      method %in% methods2compare_dense) 

# plot 
ggRMSE <- ggplot(df_RMSE, aes(x=method, y=RMSE, fill=method)) +
  geom_boxplot(alpha=0.5, outlier.shape = NA) + 
  # ggtitle ("Prediction of ibrutinib response")  +
  geom_beeswarm(cex=2.5) + facet_wrap(~method_type, scale="free_x") +
  scale_fill_manual(values = cols4methods) +
  theme(axis.text.x = element_text(angle=60, vjust=1, hjust=1),
        axis.title.x = element_blank()) + guides(fill= FALSE) +
    coord_cartesian(ylim = c(0.03,0.16))
ggRMSE
```

## Penalty Factors (in sparse graper)
```{r hyperparameters, fig.widht=8, fig.height=3}
df_pf <- melt(lapply(resultList, function(l) l$pf_mat[, "graper_SS", drop=F]),
              varnames = c("omic", "method"))
df_pf$omic <-  c("Drugs", "Methylation","mRNA")[df_pf$omic]
df_pf$value <- as.numeric(df_pf$value)

gg1 <- ggplot(df_pf,aes(x=omic, y=value)) + geom_boxplot(alpha=0.5, outlier.shape = NA) +
  ylab(expression(hat(gamma))) + geom_beeswarm(aes(col=omic), cex=2.5) +
  guides(col=FALSE) + xlab("omic type") + 
  theme_bw(base_size = 15) + scale_color_manual(values = cols4groups) +
  guides(col=FALSE)

df_pf <- melt(lapply(resultList, function(l) l$sparsity_mat[, "graper_SS", drop=F]),
              varnames = c("omic", "method"))
df_pf$omic <-  c("Drugs", "Methylation","mRNA")[df_pf$omic]

gg2 <- ggplot(df_pf,aes(x=omic, y=value)) + geom_boxplot(alpha=0.5, outlier.shape = NA) +
  ylab(expression(hat(pi))) + geom_beeswarm(aes(col=omic), cex = 2.5)+
    xlab("omic type") + 
  theme_bw(base_size = 15) +
  scale_color_manual(values = cols4groups) +
  guides(col=FALSE)

ggPF <- plot_grid(gg1,gg2, rel_widths = c(1,1), nrow=1, align = "hv", axis = "lb")
ggPF
```

# Joint plot
```{r joint, fig.widht=8, fig.height=9}
plot_grid(ggRMSE, ggPF, ncol=1, labels = letters[1:2],
          label_size = 18, rel_heights = c(2,1))
```


```{r selected, eval = FALSE, echo =FALSE}
## Selected features
getSelected <- function(resultList, methodnm){
  nfolds <- length(resultList)
  chosen_per_fold <- lapply(1:nfolds, function(i) {
    names(which(resultList[[i]]$beta_mat[,methodnm]!=0))
    })
  df <- data.frame(times_selected = as.numeric(table(Reduce(c,chosen_per_fold))),
                 feature = names( table(Reduce(c,chosen_per_fold))))
  df$omic = sapply(df$feature, function(s) {
    resultList[[1]]$annot[which(rownames(resultList[[1]]$beta_mat)==s)]
    })
  df$omic <-  c("Drugs", "Methylation","mRNA")[df$omic]
  
  ggplot(df, aes(x = feature, y = times_selected, fill = omic)) +
  geom_bar(stat="identity") + 
  theme(axis.text.x = element_text(angle=60, hjust=1, vjust=1),
        axis.title.x = element_blank()) +
  ylab("# selected") + geom_hline(yintercept=nfolds, lty="dashed") +
  ggtitle(methodnm)
}

## for graper_SS
gg_graper <- getSelected(resultList, "graper_SScutoff")

## for Lasso
gg_Lasso <- getSelected(resultList, "Lasso") +
  theme(axis.text.x = element_blank())

## for EN
gg_EN <- getSelected(resultList, "ElasticNet") +
  theme(axis.text.x = element_blank())

## for IPF-Lasso
gg_IPF <- getSelected(resultList, "IPFLasso") +
  theme(axis.text.x = element_blank())

# joint plot
plot_grid(gg_graper, gg_Lasso, gg_IPF, gg_EN,
          align="hv", axis = "lb")
```

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