---
title: "DrDimont: Drug Response Prediction from Differential Multi-Omics Networks"
author: "Katharina Baum, Pauline Hiort, Julian Hugo, Spoorthi Kashyap, Nataniel Mueller, and Justus Zeinert"
date: "`r Sys.Date()`"
output: 
    rmarkdown::html_vignette: 
    toc: true 
vignette: >
  %\VignetteIndexEntry{DrDimont: Drug Response Prediction from Differential Multi-Omics Networks}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)
```
## Introduction

The main purpose of the DrDimont pipeline is to easily and efficiently generate, reduce, and combine molecular networks from two groups or conditions (e.g., of patients) to compute a differential drug interaction score based on drug targets. This allows for improved predictions of the effect of drugs (e.g., for cancer) on two groups with different characteristics. 

![DrDimont Pipeline Overview](220401_overview.png){width=80%}

(Figure taken from Figure 1A by <i>Hiort et al. (2022)</i>)

## Installation

The R package `DrDimont` can be installed via CRAN. The R dependencies of DrDimont will also be installed. The complete source code is available at https://gitlab.com/PHiort/DrDimont.

You can install the package from CRAN with:
```{r Package installation, message=FALSE, warning=FALSE, eval=FALSE}

### please also install these dependencies of WGCNA (used in DrDimont) explicitly if not already installed
if (!require('BiocManager', quietly = TRUE))
    install.packages('BiocManager')
BiocManager::install(c('GO.db', 'preprocessCore', 'impute'))

install.packages('DrDimont')
```


After installation, you can load the package into your session with:
```{r Loading DrDimont, message=FALSE}
library(DrDimont)
```

### Installation of Python and its dependencies
The pipeline uses a Python script for one of the intermediate steps. For differential drug response score computation, Python (>=3.8) must be installed on the system before running DrDimont. Please install either the standalone Python 3.9 (or 3.8) (https://www.python.org/downloads/), or Python via Anaconda (https://www.anaconda.com/download) or via miniconda (https://docs.conda.io/en/latest/miniconda.html) before running DrDimont.

If Python is installed on your system, you can use the `install_python_dependencies()` function of DrDimont to install the required Python dependencies. The Python packages will be installed in a virtual Python or Conda environment called 'r-DrDimont'. Depending on which Python package manager (standalone or Anaconda/miniconda Python) is installed on your system, the dependencies can be installed using pip (default; if standalone Python is installed) or or Conda (if Anaconda/miniconda Python is installed).

For pip run:
```{r Install python with pip, echo=TRUE, warning=FALSE, eval=FALSE}
install_python_dependencies(package_manager="pip")
```

For Conda run:
```{r Install python with conda, echo=TRUE, warning=FALSE, eval=FALSE}
install_python_dependencies(package_manager="conda")
```
<b> ATTENTION </b>: When using pip, Python version >= 3.9 must be installed on the system. When using conda Python 3.9 will be automatically installed in the conda environment. If the installation does not work with DrDimont's internal function, please refer to a manual installation of the libraries as described below.

To manually create and install the python packages in the `r-DrDimont` environment, please run the following on your command line (outside of R):

With conda run:
```{bash, echo=TRUE, warning=FALSE, eval=FALSE}
conda create -n r-DrDimont -c conda-forge --yes python=3.9
conda activate r-DrDimont
pip install numpy tqdm igraph ray
```

With pip run:
```{bash, echo=TRUE, warning=FALSE, eval=FALSE}
#on Windows run the following in your home folder:
mkdir  .\Documents\.virtualenvs\
python -m venv .\Documents\.virtualenvs\r-DrDimont
.\Documents\.virtualenvs\r-DrDimont\Scripts\activate
pip install --upgrade pip numpy tqdm igraph ray

# on Linux and Mac run the following in your home folder:
mkdir .virtualenvs/
python -m venv .virtualenvs/r-DrDimont
source .virtualenvs/r-DrDimont/bin/activate
pip install --upgrade pip numpy tqdm igraph ray
```

 <b> ATTENTION </b>: The Python dependencies have to be installed into a virtual or conda environment with the name `r-DrDimont`, otherwise, the execution of the Python script will not work. 

A text file (`requirments_pip.txt` or `requirements_conda.txt`) with all required packages can also be downloaded from gitlab in the `inst/` directory or you can find them in your R package directory folder in the `DrDimont/` folder.

## Example Data Set Description

The following exemplary pipeline application showcases the usage of molecular breast cancer data with ER+ (Estrogen receptor-positive) patient samples as group A and ER- (Estrogen receptor-negative) as group B. A reduced exemplary data set is included within the package.

The breast cancer data by <i>Krug et al. (2020)</i> used for this tutorial is already preprocessed and only includes samples with tumor purity > 0.5 and known ER status. Metabolite data was sampled randomly to generate distributions similar to those reported, e.g., in <i>Terunuma et al. (2014)</i>.

The data set contains observations from:

* 78 ER+ samples
* 34 ER- samples

|   |Number of genes, etc.|Preprocessing|Identifier|
|---|---|---|---|
|mRNA|13915|quantified mRNA expression; log2-transformed FPKM values, NAs set to -11, removed mRNAs with > 90% of zero measurements, reduced|gene name|
|Protein|5809 (ER+) and 5845 (ER-)|quantified proteomics data; normalized, standardized, removed proteins with > 20% NAs, reduced|NCBI RefSeq ID, gene name|
|phosphosites|10272 (ER+) and 11318 (ER-)|quantified phosphoproteomics data; normalized, removed phosphosites with > 20% NAs, reduced|phosphosite, gene name, NCBI RefSeq ID|
|Metabolite| 275 from 33 (ER+) and 34 (ER-) samples|randomly sampled metabolomics data; removed metabolites with > 50% NAs|biochemical name, PubChem ID, metabolon ID|

To limit the runtime and space requirements of the example, we reduced the mRNA, protein, and phosphosite data to a random set of 50 genes. The 50 genes were randomly selected from the set of genes with known drug targets from The Drug Gene Interaction Database (https://dgidb.org/). The metabolite data were also randomly reduced to 50 metabolites.

### Load the data

First, you load the preprocessed data. This data is included in the package and does not need to be manually loaded, but can be directly accessed once `library(DrDimont)` is called.

```{r Load data}
data("mrna_data")
data("protein_data")
data("phosphosite_data")
data("metabolite_data")
data("metabolite_protein_interactions")
data("drug_gene_interactions")
```

### Transform the data to the required input format

After loading the data, you can use formatting functions to bring your data into the required input formats:

* make_layer() - creates individual molecular layers from raw data and unique identifiers
* make_connection() - specifies connections between two individual layers
* make_drug_target() - formats drug target interactions

#### Create the individual layers data structure from the molecular data
Before running the pipeline, you can create individual layer objects using `make_layer()`. Please supply raw data stratified over two patient groups and unique identifiers for the molecular entities, e.g, genes. The function `make_layer()` requires the following input parameters: `name`, `data_groupA`, `data_groupB`, `identifiers_groupA`, and `identifiers_groupB`. Please give each layer a unique name with the `name` argument. The `identifiers_groupA` and `identifiers_groupB` parameters are given dataframes which should contain one or more uniquely named columns with identifiers of the molecular entities in the rows, e.g., gene names. You can supply the raw data with the `data_groupA` and `data_groupB` parameters with the <b>molecular entities (e.g, genes) as rows</b> and the <b>samples as columns</b>. Please ensure that the identifiers of the molecular entities are in the same order as the columns in the raw data. If you have only one group to analyse, then you can set the parameters `data_groupB=NULL` and `identifiers_groupB=NULL`.


Run the code below for exemplary raw dataframes:
```{r Data inspection}
# Data inspection
mrna_data$groupA[1:3, 1:5]
protein_data$groupA[1:3, 1:5]
phosphosite_data$groupA[1:3, 1:5]
metabolite_data$groupA[1:3, 1:5]

```

Run the code below to create the individual layers:
```{r Create layers}
# Create individual layers
mrna_layer <- make_layer(name="mrna",
                         data_groupA=mrna_data$groupA[,-1],
                         data_groupB=mrna_data$groupB[,-1],
                         identifiers_groupA=data.frame(gene_name=mrna_data$groupA$gene_name),
                         identifiers_groupB=data.frame(gene_name=mrna_data$groupB$gene_name))

protein_layer <- make_layer(name="protein",
                            data_groupA=protein_data$groupA[, c(-1,-2)],
                            data_groupB=protein_data$groupB[, c(-1,-2)],
                            identifiers_groupA=data.frame(gene_name=protein_data$groupA$gene_name, 
                                ref_seq=protein_data$groupA$ref_seq),
                            identifiers_groupB=data.frame(gene_name=protein_data$groupB$gene_name, 
                                ref_seq=protein_data$groupB$ref_seq))

phosphosite_layer <- make_layer(name="phosphosite",
                                data_groupA=phosphosite_data$groupA[, c(-1,-2, -3)],
                                data_groupB=phosphosite_data$groupB[, c(-1,-2, -3)],
                                identifiers_groupA=data.frame(phosphosite_data$groupA[, 1:3]),
                                identifiers_groupB=data.frame(phosphosite_data$groupB[, 1:3]))

metabolite_layer <- make_layer(name="metabolite",
                               data_groupA=metabolite_data$groupA[, c(-1,-2, -3)],
                               data_groupB=metabolite_data$groupB[, c(-1,-2, -3)],
                               identifiers_groupA=data.frame(metabolite_data$groupA[, 1:3]),
                               identifiers_groupB=data.frame(metabolite_data$groupB[, 1:3]))

```

Run the code below to create a list of all individual layers for the pipeline input:
```{r Make layers list}
all_layers <- list(mrna_layer, protein_layer, phosphosite_layer, metabolite_layer)
```

#### Create inter-layer connections data structure
You can supply inter-layer connections with `make_connection()`. The parameters `from` and `to` have to match a name given in the previously created layers by `make_layer()`. The established connection will result in an undirected combined graph. The parameter `group` indicates whether the connection will be applied to `both` groups (default) or only group `A` or `B`. There are two options to connect layers: (i) based on identical identifiers of entities, or (ii) based on a given interaction table. 

For (i), two layers should contain one matching column name in their `identifiers_groupA`/`identifiers_groupB` dataframes that is passed as the parameter `connect_on`. Two entities in the different layers with the same ID therein are connected with an edge of fixed weight (indicated by the `weight` parameter, default 1).

For example:
```{r Connections, eval=FALSE}
# (i) make inter-layer connection 
make_connection(from='mrna', to='protein', connect_on='gene_name', weight=1, group="both")
```

For (ii), an interaction table containing three columns is required. Two columns should contain entity IDs that are also given in the respective identifiers `identifiers_groupA`/`identifiers_groupB` of the two layers to be connected. One column of those should have the same name as a column name given in the `identifiers_groupA`/`identifiers_groupB` dataframes of one layer, and the second column should have the same name for the second layer. The third column should contain the weights with which the respective entities of the two layers are to be connected. 

See `data(metabolite_protein_interactions)` for an exemplary interaction table. The table contains the columns "pubchem_id" also given for the metabolite layer, "gene_name" also given for the protein layer, and "combined_score" containing the weights for the respective interactions:
```{r Data inspection interactions}
# Data inspection
metabolite_protein_interactions[1:3, ]
```

The interaction table is passed to the `connect_on` parameter of `make_connection()` and the column name of the column containing the weights to the `weight` parameter. 
For example:
```{r Inter-layer connection, eval=FALSE}
# (ii) make inter-layer connection 
make_connection(from='protein', to='metabolite', 
                connect_on=metabolite_protein_interactions, 
                weight='combined_score', group="both")
```

If you have only one layer, you can skip the next step and set the parameter `inter_layer_connections=NULL` later on.

Run the code below to create a list of all inter-layer connections for pipeline input:
```{r Inter-layer connections}
all_inter_layer_connections = list(
    make_connection(from='mrna', to='protein', connect_on='gene_name', weight=1, group="both"),
    make_connection(from='protein', to='phosphosite', connect_on='gene_name', weight=1, group="both"),
    make_connection(from='protein', to='metabolite', 
        connect_on=metabolite_protein_interactions, weight='combined_score', group="both")
)
```

#### Create drug-target interaction data structure
Drug-target interactions are required to run the entire pipeline. For that you need an interaction table mapping drugs to their targets, e.g, proteins. The table should contain two columns: one column containing the drug IDs with the name `drug_name` and another column containing the drug targets with a name matching a column name in the `identifiers_groupA`/`identifiers_groupB` dataframes of the target layer. The example data includes a table from The Drug Gene Interaction Database providing interactions of drugs with genes. The exemplary dataframe has three columns (gene_name, drug_name, drug_chembl_id), one containing the gene names also given for the target protein layer, the second containing the drug names, which are used to identify the drugs, and a third column containing the ChEMBL IDs of drugs, which will be ignored in the pipeline. The dataframe of the drug-target interactions should have a column named `drug_name` containing drug identifiers.

Example:
```{r Data inspection drug-target}
# Data inspection
drug_gene_interactions[1:3, ]
```

The function `make_drug_target()` generates the required format for the drug-target interactions in the pipeline. The parameter `target_molecules` should match one of the layer names, e.g., `protein`. The dataframe supplied with the parameter `interaction_table` should map drugs to their target as described above. The column in the interaction table containing the targets should be given with the `match_on` parameter, e.g, `match_on=gene_name` for `protein` as targets.

Run the code below to create a list containing the drug-target input for the pipeline:
```{r Make drug-target interaction}
all_drug_target_interactions <- make_drug_target(
                                    target_molecules='protein',
                                    interaction_table=drug_gene_interactions,
                                    match_on='gene_name')
```

#### Check input data structures
When the input data structures of the individual layers, the inter-layer connections, and the drug target interactions are created, they are checked automatically for validity. Additionally, the function below checks for various possible input formatting and connection errors and reports registered data set sizes (samples, entities) for the user to compare with the intended input.

```{r Check for errors}
return_errors(check_input(layers=all_layers, 
                          inter_layer_connections=all_inter_layer_connections, 
                          drug_target_interactions=all_drug_target_interactions))
```


## Run the complete pipeline

The pipeline can be run entirely or in individual steps. To set global pipeline , you can create a settings list using the `drdimont_settings()` function. This function contains default parameters that can be modified as shown below. For a detailed explanation of all possible settings and parameters, please refer to the function documentation by calling `?drdimont_settings()`.

Please be aware of the Python script used in one of the pipeline steps (see <em>Requirements</em> above). If you have installed Python and the required packages via pip, then you should set the `drdimont_settings()` parameter `conda=FALSE`. If you have installed Python and the required packages via conda, then set the `drdimont_settings()` parameters `conda=TRUE`. `drdimont_settings()` will automatically check if Python can be found and prints a warning if not.

The intermediate pipeline and drug response scores output (parameter `save_data`) is deactivated (default), but especially for large data files consider turning it on. You can specify the output location of files with the `saving_path` parameter. If not specified, all files will be written to a temporary file created by R. For this example, the data will be saved in a temporary directory. If you want to save the data elsewhere, change the parameter `saving_path` below. The intermediate output data includes RData files of the correlation matrices, the individual graphs, the combined graphs, the drug target edges, the interaction score graphs, and the differential score graph. The drug response scores are saved in a TSV file in the specified output directory if `save_data=TRUE`.

See <em>Running the individual pipeline steps</em> and call `?drdimont_settings()` for further explanations of settings parameters.

Run the following code to create a settings list for the example:
```{r Settings}
example_settings <- drdimont_settings(
                        ### saving
                        saving_path = tempdir(),
                        save_data = FALSE,
                        ### network generation
                        correlation_method = "spearman",
                        handling_missing_data = list(
                            default = "pairwise.complete.obs",
                            mrna = "all.obs"),
                        ### network reduction
                        reduction_method = "pickHardThreshold",
                        ### pickHardThreshold
                        r_squared=list(default=0.65, metabolite=0.1),
                        cut_vector=list(default=seq(0.2, 0.65, 0.01)),
                        mean_number_edges = NULL,
                        edge_density = NULL,
                        ### p-value (not used in this example)
                        p_value_adjustment_method = "BH",
                        reduction_alpha = 0.05,
                        ### interaction_score
                        conda = FALSE,
                        max_path_length = 3,
                        num_cpus = 1,
                        int_score_mode = "auto",
                        ### drug response score
                        median_drug_response=FALSE,
                        absolute_difference=FALSE
                        )
# to disable multi-threading for example run: (not recommended for actual data processing)
WGCNA::disableWGCNAThreads()

```


To run the entire pipeline from beginning-to-end the `run_pipeline()` function can be used:
```{r Run pipeline, eval=FALSE}
run_pipeline(layers=all_layers, 
             inter_layer_connections=all_inter_layer_connections, 
             drug_target_interactions=all_drug_target_interactions, 
             settings=example_settings)
```

## Run the individual pipeline steps

You can also use the pipeline in a modular fashion. The modules then refer to the different steps:

<ol>
<li>Compute correlation matrices</li>
<li>Generate individual graphs</li>
<li>Combine graphs</li>
<li>Identify drug targets and their edges</li>
<li>Calculate integrated interaction score</li>
<li>Generate differential graph</li>
<li>Calculate differential drug response score</li>
</ol>

### Step 1: Compute correlation matrices

In step one, correlation matrices are computed for the layers created above. The parameter `correlation_method` in `drdimont_settings()` can be set to "spearman" (default), "pearson", or "kendall" as the correlation methods. The list of layers and the settings list are passed to `compute_correlation_matrices()`. 

To reduce runtime, the following example will only analyze the first 10 genes and patients of the mRNA layer (to compute all layers, set `layers=all_layers` in `compute_correlation_matrices()`):
```{r Correlation matrices, message=FALSE, results='hide'}
reduced_mrna_layer <- make_layer(name="mrna",
                          data_groupA=t(mrna_data$groupA[1:10,2:11]),
                          data_groupB=t(mrna_data$groupB[1:10,2:11]),
                          identifiers_groupA=data.frame(gene_name=mrna_data$groupA$gene_name[1:10]),
                          identifiers_groupB=data.frame(gene_name=mrna_data$groupB$gene_name[1:10]))

example_correlation_matrices <- compute_correlation_matrices(
                                    layers=list(reduced_mrna_layer), 
                                    settings=example_settings)
```

The resulting data structure `example_correlation_matrices` is a nested named list with three levels containing the correlation matrices and annotation dataframes. The first level is `correlation_matrices` and `annotations`. The correlation matrices are separated on the second level by group (`groupA` and `groupB`) and on the third level by layer name (e.g., `mrna`, `protein`, etc.). The annotations element contains annotations for each group and a combination of both at the second level (`groupA`, `groupB`, and `both`). The third level consists of named dataframes, for each layer one dataframe, which contains the pipeline-internal mapping of the `identifiers_groupA`/`identifiers_groupB` dataframes of the individual layers to layer-specific node IDs.

Example:
```{r}
# Data inspection
data("correlation_matrices_example")
correlation_matrices_example$annotations$groupA$protein[1:3, ]
```

### Step 2: Generate individual graphs

Next, the individual graphs are generated. In this step, edge weights are established based on the correlation computation, and the edges are reduced by the specified reduction method. Reduction can be done based on maximizing scale-freeness employing `WGCNA::pickHardThreshold` (`reduction_method="pickHardThreshold"` in `drdimont_settings()`; default) or based on the significance of the correlation (`reduction_method="p_value"`). Please call `?drdimont_settings()` for more information on additional "pickHardThreshold" and "p_value" settings. With "pickHardThreshold" the networks can also be reduced to at most a given number of mean edges or a given density with the parameters `mean_number_edges` and `edge_density` respectively (default of both: `NULL`).

Run the following code to generate the individual graphs:
```{r Individual graphs, message=FALSE, results='hide'}
data("correlation_matrices_example")
example_individual_graphs <- generate_individual_graphs(
                                 correlation_matrices=correlation_matrices_example, 
                                 layers=all_layers, 
                                 settings=example_settings)
```

The resulting data structure `example_individual_graphs` is a nested named list with three levels containing the graphs and annotation dataframes. The element `graphs` contains, similar to the correlation matrices, the two groups on the second level and the graphs as iGraphs objects for each molecular layer on the third level. The `annotations` element is a copy from the `annotations` element of the `example_correlation_matrices` data (see <em>Step 1</em>).


### Step 3: Combine graphs

In this step, the individual graphs are combined into a single combined graph per group based on the above created inter-layer connections. The function creates the disjoint union of the individual graphs and adds inter-layer edges with the specified weight.

Run the following code to combine the individual layers:
```{r Combine graphs, message=FALSE, results='hide'}
example_combined_graphs <- generate_combined_graphs(
                               graphs=example_individual_graphs[["graphs"]], 
                               annotations=example_individual_graphs[["annotations"]], 
                               inter_layer_connections=all_inter_layer_connections, 
                               settings=example_settings)
```

The resulting data structure `example_combined_graphs` is a nested named list with two levels containing the combined graphs and a combined annotation dataframe. The element `graphs` includes the two groups on the second level with the combined graphs as iGraphs objects. The `annotations` element consists of a dataframe on the second level named `both`, which contains the mapping of the identifier dataframes of the individual layers to the layer-specific node IDs for all layers together.

Example:
```{r}
# Data inspection
example_combined_graphs$annotations$both[1:3, ]
```


### Step 4: Identify drug targets and their edges

Next, the drug targets are identified in the combined graph for each group to extract the list of relevant drugs. Here, the node IDs of the specified drug targets are found in the combined graph, and the drugs are mapped to their target nodes. Additionally, edge lists are returned containing the incident edges of drug target nodes for which integrated interaction scores are computed in the next step.

Run the following code to extract the drug targets and their edges:
```{r Drug targets and their edges, message=FALSE, results='hide'}
example_drug_target_edges <- determine_drug_targets(
                                 graphs=example_combined_graphs[["graphs"]], 
                                 annotations=example_combined_graphs[["annotations"]], 
                                 drug_target_interactions=all_drug_target_interactions, 
                                 settings=example_settings)
```

The resulting data structure `example_drug_target_edges` is a nested named list with two levels. The first level comprises the elements `targets` and `edgelists`. The element `targets` contains the dataframe `target_nodes` and the dictionary-like list `drugs_to_target_nodes`. The dataframe `target_nodes` contains the node IDs of the nodes that are drug targets and TRUE/FALSE values if they are present in the graph of each group. The list `drugs_to_target_nodes` maps the drugs to the node IDs of their targets. The element `edgelists` consists of a dataframe for each of the groups (`groupA` and `groupB`), which, respectively, contain the incident edges of the drug targets and their weights (columns `from`, `to`, and `weight`).


### Step 5: Calculate integrated interaction score

In this step, the combined graphs for each group and the edge list from `drug_target_edges` are used to calculate the integrated interaction scores for all edges incident to drug targets. The pipeline uses a Python script to compute the integrated interaction scores in the function `generate_interaction_score_graphs()`. The input data for the Python script (combined graphs for both groups in GML format and relevant edges lists for both groups in TSV format) are written to disk, and the script is called to calculate the scores. Output files written by the Python script are two graphs in GML format containing the interaction score as an additional edge attribute called `interactionweight`. These are then loaded into R and returned in a named list containing the graphs for `groupA` and `groupB`, respectively.

<b> ATTENTION </b>: Data exchange via files is necessary and can take a long time for large datasets. Additionally, the interaction score computation can be slow. Therefore, do not set `max_path_length` in `drdimont_settings()` to a large value (default: 3). If `max_path_length=1`, then the integrated interaction scores will be the same as the correlation-based edge weights. 

The Python script for integrated interaction score computation is parallelized using ray (https://www.ray.io/). Refer to the Ray documentation if you encounter problems with running the Python script in parallel. Use the setting `int_score_mode="sequential"` in `drdimont_settings()` for forced sequential computation or `int_score_mode="ray"` and `num_cpus=X` (use an appropriate number of cpus for X) for parallel computation otherwise one of the two will be automatically chosen based on the size of the data.

<b> Running Ray on a compute cluster </b>: Please specify the number of cpus you want the process to use with the parameter `num_cpus` in `drdimont_settings()`. (Deprecated: If you want to run DrDimont on a cluster, run `ray start --head --num-cpus 12` (change number of CPUs as fit) on the command line, before starting the R session.)

Run the following code to calculate the integrated interaction scores:
```{r Calculate interaction score, eval=FALSE, message=FALSE, results='hide'}
example_interaction_score_graphs <- generate_interaction_score_graphs(
                                        graphs=example_combined_graphs[["graphs"]], 
                                        drug_target_edgelists=example_drug_target_edges[["edgelists"]], 
                                        settings=example_settings)
```


### Step 6: Generate differential graph

To generate a differential graph, the difference of the interaction scores between the two groups is computed by subtracting the values of the edge attributes of `groupB` from `groupA` (i.e., `groupA` - `groupB`). A single differential graph with `differential_score` and `differential_interaction_score` as edge attributes is returned. The edge attribute `differential_score` is the difference of the correlation-based edge weights, and `differential_interaction_score` is the difference in integrated interaction scores. Missing edges in one of the two groups are set to zero before computing the difference. The differential integrated interaction score is set to `NA` if the integrated interaction score was not computed in both groups, i.e., for all edges not incident to a drug target.

```{r Calculate differential score, message=FALSE, results='hide'}
data("interaction_score_graphs_example")
example_differential_graph <- generate_differential_score_graph(
                                  interaction_score_graphs=interaction_score_graphs_example, 
                                  settings=example_settings)

# if interaction score graphs have been computed use the following:
#example_differential_score_graph <- generate_differential_score_graph(
#                                        interaction_score_graphs=example_interaction_score_graphs, 
#                                        settings=example_settings)
```


### Step 7: Calculate differential drug response score

The differential drug response score is calculated based on the differential graph in the last step. The score of a drug is the mean (default) or the median of all differential integrated interaction scores of the edges incident to the drug targets. Drugs that have only targets without any edges have a `NA` as a differential drug response. Only drugs with at least one target present in the network are analysed. The differential drug response score can be calculated in different ways: by computing the mean (default) or the median of the differential integrated interaction scores of the edges incident to a drug's targets (set `median_drug_response=TRUE` in `drdimont_settings()` for median computation). The drug response score can also be computed from the differential (default) or the absolute differential integrated interaction scores (set `absolute_difference=TRUE` in `drdimont_settings()` to use absolute differential interaction scores). The drug response score is reported in absolute values.

```{r Drug response, message=FALSE, results='hide'}
example_drug_response_scores <- compute_drug_response_scores(
                                    differential_graph=example_differential_graph,
                                    drug_targets=example_drug_target_edges[["targets"]],
                                    settings=example_settings)
```

The first few lines of the resulting dataframe are shown below. The dataframe is saved as `drug_response_score.tsv` in the specified output folder (see <em>Running the complete pipeline</em> above). The drug response score is an indirect measure of how the strength of connectivity differs between the groups for the drug targets of the particular drug.

```{r Result Output}
head(dplyr::filter(example_drug_response_scores, !is.na(drug_response_score)))
```

## References
Full citation:

<i>Krug et al. (2020)</i>: 

* Krug, K. et al. (2020) Proteogenomic Landscape of Breast Cancer Tumorigenesis and Targeted Therapy. <i>Cell</i>, 183:1436-1456.e31. https://www.doi.org/10.1016/j.cell.2020.10.036

<i>Terunuma et al. (2014)</i>:

* Terunuma, A. et al. (2014) MYC-driven accumulation of 2-hydroxyglutarate is associated with breast cancer prognosis. <i>J Clin Invest.</i>, 124:398-412. https://www.doi.org/10.1172/JCI71180

<i>Hiort et al. (2022)</i>:

* Hiort, P. et al. (2022) DrDimont: Explainable drug response prediction from differential analysis of multi-omics networks. <i>Bioinformatics</i>, 38:ii113-ii119. https://doi.org/10.1093/bioinformatics/btac477

The package `DrDimont` is an updated version of the previously published `molnet` package (https://github.com/molnet-org/molnet)