---
title: "Differential Expression Analysis"
author: 
- name: Lis Arend
bibliography: references.bib
biblio-style: apalike
link-citation: yes
colorlinks: yes
output: 
  bookdown::html_document2:
    toc: true
    toc_depth: 2
    number_sections: true
    fig_caption: true
pkgdown:
  as_is: true
vignette: >
  %\VignetteIndexEntry{5. Differential Expression Analysis}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---


```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>", message = TRUE, warning = FALSE,
  fig.width=8,
  fig.height =6
)
```

# Load PRONE Package

```{r setup, message = FALSE}
library(PRONE)
```

# Load Data (TMT)

Here, we are directly working with the [SummarizedExperiment](https://bioconductor.org/packages/release/bioc/html/SummarizedExperiment.html) data. For more information on how to create the SummarizedExperiment from a proteomics data set, please refer to the ["Get Started"](PRONE.html) vignette.

The example TMT data set originates from [@biadglegne_mycobacterium_2022].

```{r load_real_tmt}
data("tuberculosis_TMT_se")
se <- tuberculosis_TMT_se
```

# Normalize Data

In order to compare the performance of different normalization methods on their ability to detect differentially expressed proteins, we performed some normalization with batch effect correction technqiues here. For more details about how to normalize data and evaluate the normalization approaches quantitatively and qualitatively using PRONE, please refer to the ["Normalization"](Normalization.html) vignette.

```{r}
se_norm <- normalize_se(se, c("IRS_on_RobNorm", "IRS_on_Median", 
                              "IRS_on_LoessF", "IRS_on_Quantile"), 
                        combination_pattern = "_on_")
```

# Differential Expression Analysis

After having performed normalization and evaluated the different normalization methods via qualitative and quantitative analysis, differential expression analysis can be used to further analyze the differences of the normalization methods and assess the impact of normalization on downstream analyses.

However before, you need to remove the reference samples in case of a TMT experiment. This can be easily done with the function `remove_reference_samples()`.

```{r}
se_norm <- remove_reference_samples(se_norm)
```

## Run DE Analysis

### Specific Comparisons

First, you need to specify the comparisons you want to perform in DE analysis. For this, the function `specify_comparisons()` was developed which helps to build the right comparison strings. 

However, you can also just simply create a vector of comparisons to ensure the correct order and handle this vector over to the DE analysis method.

```{r}
comparisons <- specify_comparisons(se_norm, condition = "Group", 
                                   sep = NULL, control = NULL)

comparisons <- c("PTB-HC", "TBL-HC", "TBL-PTB", "Rx-PTB")
```

The function `specify_comparisons()` becomes handy when having a lot of sample groups and foremost samples of multiple conditions (not the case for this dataset). 

For instance, you have a data set with diabetic and healthy samples (condition) that were measured at day 3 (3d), day 7 (7d), and day 17 (d14) after operation (timepoints. In DE analysis, you want to perform comparisons between the diabetic and healthy samples at each time point but also between the different time points for a fixed condition (e.g. diabetic). Instead of writing all comparisons manually, the function helps to build the right comparison strings by only considering comparisons where at least one of the groups (timepoint or condition) remains static. For this you only need to have a column in your metadata named for instance "Condition" that combines the two groups by "_", with values diabetic_3d, diabetic_7d, diabetic_14d, etc. And in the function you specify the column name "Condition" and the separator (sep = "_").

To ensure clarity, a comparison in the form of Condition-Reference, where the logFC is positive, signifies that the protein intensity is greater in the Condition compared to the Reference.

### Perform DE Analysis

The function `run_DE()` performs the DE analysis on the selected SummarizedExperiment and comparisons. A novel feature of PRONE is to compare the DE results of different normalization methods as it has been shown that normalization has an impact on downstream analyses. 

Therefore, DE analysis can be performed on multiple assays (normalization methods) at once using the already known "ain" parameter.

In the following, a more detailed explanation of the other parameters is given:

* The condition of the SummarizedExperiment object, specified at the beginning, can be used (condition = NULL) or any other column of the meta data can specified. 
* Three methods are available for DE analysis: limma [@ritchie_limma_2015], DEqMS [@zhu_deqms:_2020], and ROTS [@suomi_rots:_2017]. 
* The logFC parameter specifies if the log2 fold change should be calculated.
* The logFC_up and logFC_down parameters specify the log2 fold change thresholds for up- and down-regulated proteins.
* The p_adj parameter specifies if the p-values should be adjusted.
* The alpha parameter specifies the significance level for the p-values.

In addition, the following parameters are specific for the DE analysis methods:

* Limma (DOI: 10.18129/B9.bioc.limma):
  * The covariate parameter can be used to include a covariate in the DE analysis in case limma is performed. This can be any column of the meta data.
  * The trend parameter specifies if the trend test should be performed in case limma is used.
  * The robust parameter specifies if the robust method should be used in case limma is used.
  
* ROTS (DOI: 10.18129/B9.bioc.ROTS): 
  * B is the number of bootstrapping and K is the number of top-ranked features for reproducibility optimization.
* DEqMS (DOI: 10.18129/B9.bioc.DEqMS)
  * The DEqMS_PSMs_column need to be specified. This can be any column in the `rowData(se)` that represents the number of quantified peptides or PSMs. It is used in DEqMS to estimate prior variance for proteins quantified by different number of PSMs. 

```{r}
de_res <- run_DE(se = se_norm, 
                     comparisons = comparisons,
                     ain = NULL, 
                     condition = NULL, 
                     DE_method = "limma", 
                     logFC = TRUE, 
                     logFC_up = 1, 
                     logFC_down = -1, 
                     p_adj = TRUE, 
                     alpha = 0.05,
                     covariate = NULL, 
                     trend = TRUE,
                     robust = TRUE,
                     B = 100,
                     K = 500
                 ) 
```

If you want to apply other logFC or p-value threshold, there is no need to re-run the DE analysis again. With `apply_thresholds()`, you can simply change the threshold values.

For instance, if you want to not apply a logFC threshold and only consider proteins with an adjusted p-value of 0.05 as DE, just set the logFC parameter to FALSE. In this case, the proteins are not classified into up- and down-regulated but only as significant change or no change.

```{r}
new_de_res <- apply_thresholds(de_res = de_res, logFC = FALSE, p_adj = TRUE, 
                               alpha = 0.05)
```

However, if you still want to see if a protein has a positive or negative logFC, you can set the logFC parameter to TRUE with the logFC_up and logFC_down parameters to 0.
 
```{r}
new_de_res <- apply_thresholds(de_res = de_res, logFC = TRUE, 
                               logFC_up = 0, logFC_down = 0, 
                               p_adj = TRUE, alpha = 0.05)
```


## Visualize DE Results

### Barplot
To get an overview of the DE results of the different normalization methods, you can visualize the number of significant DE proteins per normalization method in a barplot using `plot_overview_DE_bar()`. This plot can be generated in different ways by specifying the "plot_type" parameter:

* single: plot a single barplot for each comparison (default)
* facet_comp: facet the barplot by comparison
* stacked: stack the number of DE per comparison
* facet_regulation: stack the number of DE per comparison but facet by up- and down-regulated

```{r, fig.cap = "Barplot showing the number of DE proteins per normalization method colored by comparison and facetted by up- and down-regulation."}
plot_overview_DE_bar(de_res, ain = NULL, comparisons = comparisons, 
                     plot_type = "facet_regulation")
```

You can also just visualize two specific comparisons:

```{r, fig.cap = "Barplot showing the number of DE proteins per normalization method faceted by comparison."}
plot_overview_DE_bar(de_res, ain = NULL, comparisons = comparisons[seq_len(2)], 
                     plot_type = "facet_comp")
```


#### Tile Plot

You can also get an overview of the DE results in form of a heatmap using the `plot_overview_DE_tile()`.

```{r, fig.cap = "Heatmap showing the number of DE proteins per comparison and per normalization method."}
plot_overview_DE_tile(de_res)
```

#### Volcano Plots

Another option is to generate volcano plots for each comparison. The function `plot_volcano_DE()` generates a grid of volcano plots for each normalization techniques (facet_norm = TRUE) or for each comparison (facet_comparison = TRUE). A list of volcano plots is returned.

```{r, fig.cap = "Volcano plots per normalization method for a single comparison."}
plot_volcano_DE(de_res, ain = NULL, comparisons = comparisons[1], 
                facet_norm = TRUE)
```

#### Heatmap of DE Results

Furthermore, you can visualize the DE results in form of a heatmap. The function `plot_heatmap_DE()` generates a heatmap of the DE results for a specific comparison and normalization method.

```{r, fig.cap = "Individual heatmap of the DE results for a specific comparison and a selection of normalization methods. The adjusted p-values are added as row annotation, while the condition of each sample is shown as column annotation."}
plot_heatmap_DE(se_norm, de_res, ain = c("RobNorm", "IRS_on_RobNorm"), 
                comparison = "PTB-HC", condition = NULL, label_by = NULL, 
                pvalue_column = "adj.P.Val")
```

### Intersection Analysis of DE Results

Moreover, you can also intersect the DE results of different normalization methods to see how many DE proteins overlap. You can either plot for each requested comparison an individual upset plot (plot_type = "single") or stack the number of overlapping DE proteins per comparison ("stacked"). Not only the upset plot(s) are returned, but also a table with the intersections is provided by the `plot_upset_DE()`function.

```{r, fig.cap = "Upset plot showing the overlapping DE proteins of different normalization methods colored by comparison.", fig.height = 12}
intersections <- plot_upset_DE(de_res, ain = NULL, 
                               comparisons = comparisons[seq_len(3)], 
                               min_degree = 6, plot_type = "stacked")
if(!is.null(intersections)){
  # put legend on top due to very long comparisons
intersections$upset[[2]] <- intersections$upset[[2]] + 
  ggplot2::theme(legend.position = "top", legend.direction = "vertical")
intersections$upset
}
```

Additionally, the Jaccard similarity index can be calculated to quantify the similarity of the DE results between the different normalization methods. A individual heatmap can be generated for each comparison ("plot_type = single"), a single heatmap facetted by comparison ("facet_comp") or a single heatmap taking all comparisons into account ("all") can be generated.

```{r, fig.cap = "Heatmap showing the Jaccard similarity indices of the DE results between different normalization methods for all comparisons."}
plot_jaccard_heatmap(de_res, ain = NULL, comparisons = comparisons, 
                     plot_type = "all")
```

PRONE offers the functionality to extract a consensus set of DEPs based on a selection of normalization methods and a threshold for the number of methods that need to agree on the DE status of a protein. The function `get_consensus_DE()` returns a list of consensus DE proteins for either each comparison separately or for all comparisons combined.

```{r}
DT::datatable(extract_consensus_DE_candidates(de_res, ain = NULL, 
                                              comparisons = comparisons, 
                                              norm_thr = 0.8, 
                                              per_comparison = TRUE), 
              options = list(scrollX = TRUE))
```

# Functional enrichment analysis

Finally, to biologically interpret the DE results of the different normalization methods, PRONE offers the ability to perform functional enrichment on the different sets of DE proteins originating from different normalization methods. The function `plot_intersection_enrichment()` can be used to visualize the functional enrichment results of the DE proteins of different normalization methods. This function returns two types of plot: a heatmap of enriched terms visualizing whether a term in enriched in the DE set of a specific normalization method and a heatmap showing Jaccard similarity coefficients to quantify the similarity of the enriched terms between the different normalization methods.

For this analysis, the `gprofiler2` package is used to perform the functional enrichment analysis. The following parameters are required:

* se: SummarizedExperiment object
* de_res: DE results
* comparison: Specific comparison
* ain: Vector of normalization methods to perform functional enrichment
* id_column: Column name of the rowData of the SummarizedExperiment object which includes the gene names. If MaxQuant data is provided with multiple gene names per protein group, the entries in the id_column are split by ";".
* organism: Organism name (gprofiler parameter, example: hsapiens, mmusculus)
* source: Data source to use (gprofiler parameter, example: KEGG, GO:BP, GO:MF, GO:CC, REAC, HP, ...)

```{r}
enrich_results <- plot_intersection_enrichment(se, de_res, comparison = "PTB-HC", 
                             ain = c("IRS_on_Median", "IRS_on_RobNorm", "RobNorm"), 
                             id_column = "Gene.Names", organism = "hsapiens", 
                             source = "KEGG", signif_thr = 0.05)
```

```{r, fig.cap = "Heatmap showing enriched terms of the DE sets between different normalization methods for a specific comparison."}
enrich_results$enrichment_term_plot

```

```{r, fig.cap = "Heatmap showing the Jaccard similarity indices of the enriched terms between different normalization methods for a specific comparison."}
enrich_results$jaccard_intersection_plot
```

# Session Info

```{r}
utils::sessionInfo()
```

# References
