---
title: "sSNAPPY: Singel Sample directioNAl Pathway Perturbation analYsis"
author:
- name: Wenjun Nora Liu
  affiliation: Dame Roma Mitchell Cancer Research Laboratories, Adelaide Medical School, University of Adelaide
  email: wenjun.liu@adelaide.edu.au
pacakge: sSNAPPY
output: 
    BiocStyle::html_document:
    toc: yes
vignette: >
  %\VignetteIndexEntry{Single Sample Directional Pathway Perturbation Analysis}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: inline
---

```{r, echo=FALSE}
knitr::opts_chunk$set(message = FALSE, warning = FALSE, crop = NULL)
```

# Introduction

This vignette demonstrates how to use the package `sSNAPPY` to compute directional single sample pathway perturbation scores by incorporating pathway topologies, utilize sample permutation to test the significance of individual scores and compare average pathway activities across treatments. 

The package also provides a function to visualize overlap between pathway genes contained in perturbed biological pathways as network plots. 

# Installation

The package `sSNAPPY` can be installed using the package `BiocManager`, along with other packages required for this vignette (`tidyverse`, `magrittr`, `ggplot2`, `cowplot`, and `DT`). 

```{r install, eval = FALSE}
if (!"BiocManager" %in% rownames(installed.packages()))
  install.packages("BiocManager")
pkg <- c("tidyverse", "magrittr", "ggplot2", "cowplot", "DT")
BiocManager::install(pkg)
BiocManager::install("sSNAPPY")
```

# Load packages

```{r setup,  results="hide", warning=FALSE}
library(sSNAPPY)
library(tidyverse)
library(magrittr)
library(ggplot2)
library(cowplot)
library(DT)
```

# load data

The example dataset used for this tutorial can loaded with `data()` as shown below. It's a subset of data retrieved from [Singhal H et al. 2016](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4928895/), where ER-positive primary breast cancer tumor tissues collected from 12 patients were split into fragments of equal sizes for different treatments. For the purpose of this tutorial, we only took the RNA-seq data from samples that were treated with vehicle, E2 OR E2 + R5020 for 48 hrs. They were from 5 different patients, giving rise to 15 samples in total. More detailed description of the dataset can be assessed through the help page (`?logCPM_example` and `?metadata_example`).


```{r data}
data(logCPM_example)
data(metadata_example)
# check if samples included in the logCPM matrix and metadata dataframe are identical
setequal(colnames(logCPM_example), metadata_example$sample)
# View sample metadata
datatable(metadata_example,  filter = "top")
```

# Compute weighted single-sample logFCs (ssLogFCs)

It is expected that the logCPM matrix will be filtered to remove undetectable genes and normalised to correct for library sizes or other systematic artefacts, such as gene length or GC contents, prior to applying this method. Filtration and normalisation has been performed on the example dataset. 

Before single-sample logFCs (ssLogFCs) can be computed, row names of the logCPM matrix need to be converted to `entrez ID`. This is because all the pathway topology information retrieved will be in  `entrez ID`. The conversion can be achieved through bioconductor packages `AnnotationHub` and `ensembldb`.  

```{r logCPM_example}
head(logCPM_example)
```

To compute the ssLogFCs, samples must be in matching pairs. In the example, treated samples are matched to the corresponding control sample that were derived from the same patients. So the `factor` parameter of the `weight_ss_fc()` functions needs to be set to be "patient". The function also requires the control treatment level to be specified, which was "Vehicle" in this case. 

`weight_ss_fc()` requires both the logCPM matrix and sample metadata as input. The column names of the logCPM matrix should be sample name, matching to a column called `sample` in the metadata. The metadata must also contain a treatment column, and a column corresponding to the `factor` parameter (i.e.. patient in this case).


```{r ssFC}
#compute weighted single sample logFCs
weightedFC <- weight_ss_fc(logCPM_example, metadata = metadata_example,
factor = "patient", control = "Vehicle")
```


The `weight_ss_fc()` function firstly computes raw ssLogFCs for each gene by subtracting logCPM values of control sample from the logCPM values of treated samples for each patient. 

It has been demonstrated previously that in RNA-seq data, lowly expressed genes turn to have larger variance, which is also demonstrated by the plots below. To reduce the impact of this artefact, `weight_ss_fc` also weight each ssLogFCs by estimating the relationship between the variance in ssLogFCs and mean logCPM, and defining the gene-wise weight to be the inverse of the predicted variance of that gene's mean logCPM value.



```{r lowess, fig.width=8,fig.height=4}
perSample_FC <- lapply(levels(metadata_example$patient), function(x){
    temp <- logCPM_example[seq_len(1000),str_detect(colnames(logCPM_example), x)] 
    ratio <- temp[, str_detect(colnames(temp), "Vehicle", negate = TRUE)] - temp[, str_detect(colnames(temp), "Vehicle")] 
}) %>%
    do.call(cbind,.)
aveCPM <- logCPM_example[seq_len(1000),] %>%
    rowMeans() %>%
    enframe(name = "gene_id", 
            value = "aveCPM")
p1 <- perSample_FC %>%
    as.data.frame() %>%
    rownames_to_column("gene_id") %>%
    pivot_longer(cols = -"gene_id",
                 names_to = "name",
                 values_to = "logFC") %>%
    left_join(aveCPM) %>%
    ggplot(aes(aveCPM, logFC)) +
    geom_point() +
    labs(y = "sslogFC", 
         x = "Average logCPM") +
    theme(
        panel.background = element_blank()
    )
p2 <- data.frame(
    gene_id = rownames(perSample_FC),
    variance = perSample_FC %>%
        apply(1,var)) %>%
    left_join(aveCPM) %>%
    ggplot(aes(aveCPM, variance)) +
    geom_point() +
    geom_smooth(method = "loess") +
    labs(y = "Variance in ssLogFCs", 
         x = "Average logCPM") +
    theme(
        panel.background = element_blank()
    )
plot_grid(p1, p2)

```

The output of the `weight_ss_fc()` function is a list with two element, where one is the weighted ssLogFCs matrix and the other is a vector of gene-wise weights.

# Retrieve pathway topologies in required format

*sSNAPPY* adopts the pathway perturbation scoring algorithm proposed in `r Biocpkg("SPIA")`, which makes use of gene-set topologies and gene-gene interaction to propagate pathway genes' logFCs down the topologies to compute pathway perturbation scores, where signs of scores reflect the potential directions of changes. 

Therefore, pathway topology information need to be firstly retrieved from your chosen database and converted to weight adjacency matrices, the format required to apply the scoring algorithm. This step is achieved through a chain of functions that are part of the `r Biocpkg("grapghite")` and has been nested into one simple function in this package: `retrieve_topology()`. Databases that are currently supported are:

```{r pathwayDatabases}
library(graphite)
graphite::pathwayDatabases() %>%
  dplyr::filter(species ==  "hsapiens") %>%
  pander::pander()
```

The retrieved topology information will be a list where each element is a pathway. We recommend you to save it as a file so this step only needs to be performed once for each database.

This vignette chose *KEGG* pathways as an example. To reduce the computation time in the following steps, only half of the randomly sampled *KEGG* pathways were kept. 

```{r retrieve_topology}
gsTopology <- retrieve_topology(database = "kegg")
head(names(gsTopology))
```

If only selected pathways are of interest, it's possible to only retrieve the topologies of those pathways by specifying the pathway names. 

```{r gsTopology_sub}
gsTopology_sub <- retrieve_topology(
  database = "kegg",
  pathwayName = c(
    "Glycolysis / Gluconeogenesis", 
    "Citrate cycle (TCA cycle)",
    "Pentose phosphate pathway"
    ))
names(gsTopology_sub)
```

# Score single sample pathway perturbation

Once the expression matrix, sample metadata and pathway topologies are all ready, single sample pathway perturbation scores (PS) can be computed using function `compute_perturbation_score()`, which returns a data.frame containing the test perturbation scores for each sample each pathway. 

```{r ssPertScore}
ssPertScore <- compute_perturbation_score(weightedFC$logFC, gsTopology)
head(ssPertScore)
```

# Generate null distributions of perturbation scores

To derive the empirical p-values for each single sample PS or normalize the raw scores for comparing overall treatment effects, null distributions of scores for each pathway is generated through a sample-label permutation approach. 

For each round of permutation, sample labels are randomly shuffled to derive the permuted ssLogFCs, which are then used to score pathway perturbation. We recommend to perform a minimum of 1000 rounds of permutation, which means at least 8 samples are required. The `generate_permuted_scores()` function does not require sample metadata but the number of treatments in the study design, including the control treatment, need to be specified. In this example data, the number of treatment was 3.  

Output of the `generate_permuted_scores()` function is a list where each element is a vector of permuted perturbation scores for a pathway. 

The permutation step relies on the parallel computing feature provided by `r Biocpkg("BiocParallel")`. User can choose to customize the parallel back-end or stick with the default one returned by `BiocParallel::bpparam()`. Depending on the size of the data, this step can take some time to complete. If the sample size is large, we recommend users to consider performing this step on a HPC. 

```{r permutedScore, eval=FALSE}
permutedScore <- generate_permuted_scores(
  expreMatrix  = logCPM_example, 
  numOfTreat = 3, NB = 1000, 
  gsTopology = gsTopology, 
  weight = weightedFC$weight
)
```

# Test Significances

## Of individal score

After the empirical null distributions are generated, the median and mad will be calculated for each pathway to convert the test single-sample perturbation scores derived from the `compute_perturbation_score()` to robust z-scores: $$ (Score - Median)/MAD$$ Two-sided p-values associated with each robust z-scores are also computed and will be corrected for multiple-testing using a user-define approach. The default is `fdr`.

The `normalise_by_permu` function requires the test perturbation score and permuted perturbation scores as input. 
```{r normalisedScores, eval=FALSE}
normalisedScores <- normalise_by_permu(permutedScore, ssPertScore)
```

Since the permutation step takes a long time to run and the output is too large to be included as part of the package, results of the `normalise_by_permu` step has been pre-computed and can be loaded with:

```{r}
load(system.file("extdata", "normalisedScores.rda", package = "sSNAPPY"))
```

The pathways that were significant perturbed within individual samples are:
```{r DT_indi}
normalisedScores %>%
    dplyr::filter(adjPvalue < 0.05) %>%
    left_join(metadata_example) %>%
    mutate_at(vars(c("sample", "gs_name")), as.factor) %>%
    mutate_if(is.numeric, sprintf, fmt = '%#.4f') %>%
    mutate(Direction = ifelse(robustZ < 0, "Inhibited", "Activation")) %>%
    dplyr::select(
        sample, patient, Treatment = treatment,
        `Perturbation Score` = robustZ, Direction,
        `Gene-set name` = gs_name, 
        `P-value` = pvalue, 
        FDR = adjPvalue
    ) %>%
    datatable(
        filter = "top", 
        options = list(
            columnDefs = list(list(targets = "Direction", visible = FALSE))
        )
    ) %>% 
    formatStyle(
        'Perturbation Score', 'Direction',
        backgroundColor = styleEqual(c("Inhibited", "Activation"), c('lightblue', 'indianred'))
    )
```

### Visualisation

We can use the `plot_gs_network` function to visualise the significantly perturbed biological pathways as networks, where edges between gene-sets reflect how much overlap those two gene-sets share. The function can take `normalise_by_permu`'s output, or a subset of it as its direct input. 

Nodes in the network plots could be coloured by the predicted direction of perturbation (i.e.. sign of robust z-score):

```{r sigGS_nt_zscore, fig.width= 10, fig.height=4}
pl <- normalisedScores %>%
  dplyr::filter(adjPvalue < 0.05) %>%
  split(f = .$sample) %>%
  lapply(
    plot_gs_network, 
    # layout = "dh",
    gsTopology = gsTopology, 
    colorBy = "robustZ"
  ) %>%
    lapply(function(x){
        x + theme(
        panel.grid = element_blank(), 
        panel.background = element_blank()
        ) })
plot_grid(
    plotlist = pl, 
    nrow = 1)
```

Or p-values:

```{r sigGS_nt_pvalue, fig.width= 10, fig.height=4}
pl <- normalisedScores %>%
    dplyr::filter(adjPvalue < 0.05) %>%
    split(f = .$sample) %>%
    lapply(
        plot_gs_network, 
        gsTopology = gsTopology, 
        colorBy = "pvalue", 
        color_lg_title = "P-value"
    ) %>%
    lapply(function(x){
        x + theme(
        panel.grid = element_blank(), 
        panel.background = element_blank()
        ) })
plot_grid(plotlist = pl, nrow = 1)
```

The function allows you to customize the layout, colour, edge transparency and other aesthetics of the graph. More information can be found in the help page (`?plot_gs_network`). Output of the graph is a `ggplot` object and the theme of it can be changed just as other `ggplot` figures. 

## Of overall treatment effect

Normalised perturbation scores can also be used to model mean treatment effects. An advantage of this method is that it has great flexibility that allows you to incorporate other cofactors or covariate in your modelling. 

For example, in the example dataset, samples were collected from patients with different progesterone receptor (PR) status and we can include it as a cofactor to offset its confounding effect. 

```{r fit}
fit <- normalisedScores %>%
    left_join(metadata_example) %>%
    split(f = .$gs_name) %>%
    #.["Estrogen signaling pathway"] %>%
    lapply(function(x)lm(robustZ ~ 0 + treatment + PR, data = x)) %>%
    lapply(summary)
treat_sig <- lapply(
  names(fit), 
  function(x){
    fit[[x]]$coefficients %>%
      as.data.frame() %>%
      .[seq_len(2),] %>%
      dplyr::select(Estimate, pvalue = `Pr(>|t|)` ) %>%
      rownames_to_column("Treatment") %>%
      mutate(
        gs_name = x, 
        FDR = p.adjust(pvalue, "fdr"), 
        Treatment = str_remove_all(Treatment, "treatment")
      ) 
  }) %>%
  bind_rows() 
```

The pathways that were on average perturbed due to each treatment were:
```{r treat_sig_DT}
treat_sig %>% 
    dplyr::filter(FDR < 0.05) %>%
    mutate_at(vars(c("Treatment", "gs_name")), as.factor) %>%
    mutate_if(is.numeric, sprintf, fmt = '%#.4f') %>%
    mutate(Direction = ifelse(Estimate < 0, "Inhibited", "Activation")) %>%
    dplyr::select(
        Treatment, `Perturbation Score` = Estimate, Direction,
        `Gene-set name` = gs_name, 
        `P-value` = pvalue, 
        FDR
    ) %>%
    datatable(
        filter = "top", 
        options = list(
            columnDefs = list(list(targets = "Direction", visible = FALSE))
        )
    ) %>% 
    formatStyle(
        'Perturbation Score', 'Direction',
        backgroundColor = styleEqual(c("Inhibited", "Activation"), c('lightblue', 'indianred'))
    )
```
Results from this analysis indicate that the estrogen signalling pathway was significantly activated among both E2 and E2+R5020 treated samples, which makes sense biologically.


### Visualisation

Results of `lm` can also be visualised using the `plot_gs_network()` function. We just need to change the name of the `Estimate` column to `robustZ` to  colour the networks by the predicted directionality.

```{r fig.width=12, fig.height=8}
treat_sig %>% 
    dplyr::filter(FDR < 0.05, Treatment == "E2+R5020") %>%
    dplyr::rename(robustZ = Estimate) %>%
    plot_gs_network(
       
        gsTopology = gsTopology, 
        colorBy = "robustZ"
    ) +
    theme(
        panel.grid = element_blank(), 
        panel.background = element_blank()
    )
```

By default, `plot_gs_network()` function does not include nodes that are not connected to any other nodes, which could be turnt of by setting the `plotIsolated` to TURE. 

# References
- Sales G, Calura E, Cavalieri D, Romualdi C (2012). “graphite - a Bioconductor package to convert pathway topology to gene network.” BMC Bioinformatics. https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-13-20.
- Tarca, Adi Laurentiu et al. (2009). “A novel signaling pathway impact analysis.” Bioinformatics vol. 25,1 : 75-82. doi:10.1093/bioinformatics/btn577

# Session Info

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