---
title: "Using the depmap data"
author:
- name: Theo Killian
  affiliation: VIB Center for Cancer Biology, KU Leuven
- name: Laurent Gatto
  affiliation: Computational Biology, UCLouvain
date: "`r Sys.Date()`"
output:
  BiocStyle::html_document:
    toc_float: true
vignette: >
  %\VignetteIndexEntry{depmap use cases}
  %\VignetteEngine{knitr::rmarkdown}
  %\usepackage[utf8]{inputenc}
---

```{r style, echo = FALSE, results = 'asis'}
BiocStyle::markdown()
```

This vignette illustrates use cases and visualizations of the data
found in the [depmap](https://bioconductor.org/packages/devel/data/experiment/html/depmap.html)
package. See the *depmap* vignette for details about the datasets.

# Introduction

The [depmap](https://bioconductor.org/packages/devel/data/experiment/html/depmap.html)
package aims to provide a reproducible research framework to cancer dependency
data described by [Tsherniak, Aviad, et al. "Defining a cancer dependency map."
Cell 170.3 (2017): 564-576.](https://www.ncbi.nlm.nih.gov/pubmed/28753430). The
data found in the [depmap](https://bioconductor.org/packages/devel/data/experiment/html/depmap.html)
package has been formatted to facilitate the use of common R packages such as
`dplyr` and `ggplot2`. We hope that this package will allow researchers to more
easily mine, explore and visually illustrate dependency data taken from the
Depmap cancer genomic dependency study.

# Use cases

Perhaps the most interesting datasets found within the
[depmap](https://bioconductor.org/packages/devel/data/experiment/html/depmap.html)
package are those that relate to the cancer gene dependency score, such as
`rnai` and `crispr`. These datasets contain a score expressing how vital a
particular gene is in terms of how lethal the knockout/knockdown of that gene is
on a target cell line. For example, a highly negative dependency score implies
that a cell line is highly dependent on that gene.

Load necessary libaries.

```{r load_libraries, message=FALSE, warning=FALSE, echo=TRUE}
library("dplyr")
library("ggplot2")
library("viridis")
library("tibble")
library("gridExtra")
library("stringr")
library("depmap")
library("ExperimentHub")
```

Load the `rnai`, `crispr` and `copyNumber` datasets for visualization. Note: the
datasets listed below are from the 19Q3 release. Newer datasets, such as those
from the 20Q1 release are available.

```{r load_data, message=FALSE, warning=FALSE, echo=TRUE}
## create ExperimentHub query object
eh <- ExperimentHub()
query(eh, "depmap")
rnai <- eh[["EH3080"]]
mutationCalls <- eh[["EH3085"]]
metadata <- eh[["EH3086"]]
TPM <- eh[["EH3084"]]
copyNumber <- eh[["EH3082"]]
# crispr <- eh[["EH3081"]]
# drug_sensitivity <- eh[["EH3087"]]
```

By importing the `depmap` data into the R environment, the data can be mined
more effectively. For example, if one interested researching soft tissue
sarcomas and wanted to search all such cancer cell lines for the gene with the
greatest dependency, it is possible to accomplish this task by using data
manipulation and visualization tools `dplyr` and `ggplot2`. Below, the `rnai`
dataset is selected for cell lines with *"SOFT_TISSUE"* in the CCLE name, and
displaying a list of the highest dependency scores.

```{r, soft_tissue_cell_lines, echo=TRUE}
## list of dependency scores
rnai %>% dplyr::select(cell_line, gene_name, dependency) %>%
         dplyr::filter(stringr::str_detect(cell_line, "SOFT_TISSUE")) %>%
         dplyr::arrange(dependency) %>% 
         head(10)
```

As the gene `RPL14` appears several times in the top dependencies scores, it may
make an interesting candidate target. Below, a plot of the `rnai` data is
displayed as a histogram showing the distribution of dependency scores for gene
`RPL14`. 

```{r message=FALSE, warning=FALSE}
## Basic histogram
rnai %>% dplyr::select(gene, gene_name, dependency) %>% 
         dplyr::filter(gene_name == "RPL14") %>% 
         ggplot(aes(x = dependency)) +
         geom_histogram() +
         geom_vline(xintercept = mean(rnai$dependency, na.rm = TRUE),
                    linetype = "dotted", color = "red") +
         ggtitle("Histogram of dependency scores for gene RPL14")
```

A more complex plot of the `rnai` data, as shown below involves plotting the
distribution of dependency scores for gene `RPL14` for each major type of
cancer, while highlighting the nature of mutations of this gene in such cancer
cell lines (e.g. if such are COSMIC hotspots, damaging, etc.). Notice that the
plot above reflects the same overall distribution in two dimensions.

```{r message=FALSE, warning=FALSE}
meta_rnai <- metadata %>%
             dplyr::select(depmap_id, lineage) %>%
             dplyr::full_join(rnai, by = "depmap_id") %>%
             dplyr::filter(gene_name == "RPL14") %>% 
             dplyr::full_join((mutationCalls %>%
                              dplyr::select(depmap_id, entrez_id,
                                            is_cosmic_hotspot, var_annotation)),
                                            by = c("depmap_id", "entrez_id"))
p1 <- meta_rnai %>%
      ggplot(aes(x = dependency, y = lineage)) +
      geom_point(alpha = 0.4, size = 0.5) +
      geom_point(data = subset(
         meta_rnai, var_annotation == "damaging"), color = "red") +
      geom_point(data = subset(
         meta_rnai, var_annotation == "other non-conserving"), color = "blue") +
      geom_point(data = subset(
         meta_rnai, var_annotation == "other conserving"), color = "cyan") +
      geom_point(data = subset(
         meta_rnai, is_cosmic_hotspot == TRUE), color = "orange") +
      geom_vline(xintercept=mean(meta_rnai$dependency, na.rm = TRUE),
                 linetype = "dotted", color = "red") +
      ggtitle("Scatterplot of dependency scores for gene RPL14 by lineage")
p1
```

Below is a boxplot displaying expression values for gene `RPL14` by lineage:

```{r message=FALSE, warning=FALSE}
metadata %>%
      dplyr::select(depmap_id, lineage) %>%
      dplyr::full_join(TPM, by = "depmap_id") %>%
      dplyr::filter(gene_name == "RPL14") %>% 
      ggplot(aes(x = lineage, y = expression, fill = lineage)) +
      geom_boxplot(outlier.alpha = 0.1) +
      ggtitle("Boxplot of expression values for gene RPL14 by lineage") +
      theme(axis.text.x = element_text(angle = 45, hjust=1)) +
      theme(legend.position = "none")
```

High dependency, high expression genes are more likely to interesting research
targets. Below is a plot of expression vs rnai gene dependency for
Rhabdomyosarcoma Sarcoma:

```{r  message=FALSE, warning=FALSE}
## expression vs rnai gene dependency for Rhabdomyosarcoma Sarcoma
sarcoma <- metadata %>%
           dplyr::select(depmap_id, cell_line,
                         primary_disease, subtype_disease) %>%
           dplyr::filter(primary_disease == "Sarcoma",
                         subtype_disease == "Rhabdomyosarcoma")

rnai_sub <- rnai %>% dplyr::select(depmap_id, gene, gene_name, dependency)
tpm_sub <- TPM %>% dplyr::select(depmap_id, gene, gene_name, expression)

sarcoma_dep <- sarcoma %>%
               dplyr::left_join(rnai_sub, by = "depmap_id") %>%
               dplyr::select(-cell_line, -primary_disease,
                             -subtype_disease, -gene_name)

sarcoma_exp <- sarcoma %>% dplyr::left_join(tpm_sub, by = "depmap_id")

sarcoma_dat_exp <- dplyr::full_join(sarcoma_dep, sarcoma_exp,
                             by = c("depmap_id", "gene")) %>%
                             dplyr::filter(!is.na(expression))

p2 <- ggplot(data = sarcoma_dat_exp, aes(x = dependency, y = expression)) +
      geom_point(alpha = 0.4, size = 0.5) +
      geom_vline(xintercept=mean(sarcoma_dat_exp$dependency, na.rm = TRUE),
                 linetype = "dotted", color = "red") +
      geom_hline(yintercept=mean(sarcoma_dat_exp$expression, na.rm = TRUE),
                 linetype = "dotted", color = "red") +
      ggtitle("Scatterplot of rnai dependency vs expression values for gene")
p2 + theme(axis.text.x = element_text(angle = 45))
```

A selection of the genes shown above with the lowest depenency scores, also
displaying gene expression in TPM in the last column.

```{r}
sarcoma_dat_exp %>%
    dplyr::select(cell_line, gene_name, dependency, expression) %>%
    dplyr::arrange(dependency) %>% 
    head(10)
```

Below is a boxplot displaying log genomic copy number for gene `RPL14` by
lineage:

```{r message=FALSE, warning=FALSE}
metadata %>%
    dplyr::select(depmap_id, lineage) %>%
    dplyr::full_join(copyNumber, by = "depmap_id") %>%
    dplyr::filter(gene_name == "RPL14") %>%
    ggplot(aes(x = lineage, y = log_copy_number, fill = lineage)) +
    geom_boxplot(outlier.alpha = 0.1) +
    ggtitle("Boxplot of log copy number for gene RPL14 by lineage") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
    theme(legend.position = "none")
```

# Session information

```{r echo = FALSE}
sessionInfo()
```