---
title: "pathfindR - An R Package for Pathway Enrichment Analysis Utilizing Active Subnetworks"
author: "Ege Ulgen"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{pathfindR - An R Package for Pathway Enrichment Analysis Utilizing Active Subnetworks}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

`pathfindR` is an R package for pathway enrichment analysis of gene-level omics data utilizing active subnetworks. The package also enables hierarchical clustering of the enriched pathways. The method is described in detail in Ulgen E, Ozisik O, Sezerman OU. 2018. pathfindR: An R Package for Pathway Enrichment Analysis Utilizing Active Subnetworks. bioRxiv. [https://doi.org/10.1101/272450](https://doi.org/10.1101/272450)

Our motivation to develop this package was that direct pathway enrichment analysis of differential RNA/protein expression or DNA methylation results may not provide the researcher with the full picture. That is to say; pathway enrichment of only the list of significant genes may not be informative enough to explain the underlying disease mechanisms. 

An active subnetwork is defined as a group of interconnected genes in a protein-protein interaction network (PIN) that contains most of the significant genes. Therefore, these active subnetworks define distinct disease-associated sets of genes, whether discovered through differential expression analysis or discovered because of being in interaction with a significant gene.

Therefore, we propose to leverage information from a PIN to identify distinct active subnetworks and then perform pathway enrichment analyses on these subnetworks. Briefly, this workflow first maps the significant genes onto a PIN and finds active subnetworks. Next, pathway enrichment analyses are performed using each gene set of the identified active subnetworks. Finally, these enrichment results are summarized and returned as a data frame. This workflow is implemented as the function `run_pathfindR()` and further described in the "Enrichment Workflow" section of this vignette.

This process usually yields a great number of enriched pathways with related biological functions. We therefore implemented a pairwise distance metric based on kappa statistics (as proposed by Huang et al. [^1]) between terms and based on this distance metric, also implemented hierarchical clustering and fuzzy clustering [^1] of the pathways. Details of clustering and partitioning of pathways are presented in the "Pathway Clustering" section of this vignette.

# Enrichment Workflow

The overview of the enrichment workflow is presented in the figure below:

<img src="./pathfindr.png" width="650" height="139.81" />

For this workflow, the wrapper function `run_pathfindR()` can be used. This function takes in a data frame consisting of `Gene Symbol`, `log-fold-change` (optional) and `adjusted-p values`. The first 6 rows of an example input dataset (of rheumatoid arthritis differential-expression) can be found below:

```{r}
suppressPackageStartupMessages(library(pathfindR))
data("RA_input")
knitr::kable(head(RA_input))
```

Executing the workflow is straightforward (but takes several minutes):
```{r eval=FALSE}
RA_output <- run_pathfindR(RA_input)
```

The user may want to change certain arguments of the function:
```{r eval=FALSE}
# to change the output directory
RA_output <- run_pathfindR(RA_input, output = "new_directory")

# to change the PIN (default = Biogrid)
RA_output <- run_pathfindR(RA_input, pin_name = "IntAct")
# to use an external PIN of user's choice
RA_output <- run_pathfindR(RA_input, pin_name = "/path/to/myPIN.sif")

# available gene sets are KEGG, Reactome, BioCarta, GO-BP, GO-CC and GO-MF
# default is KEGG
# to change the gene sets used for enrichment analysis
RA_output <- run_pathfindR(RA_input, gene_sets = "BioCarta")

# to change the active subnetwork search algorithm (default = "GR", i.e. greedy algorithm)
# for simulated annealing:
RA_output <- run_pathfindR(RA_input, search_method = "SA")

# to change the number of iterations (default = 10)
RA_output <- run_pathfindR(RA_input, iterations = 5) 

# to manually specify the number processes used during parallel loop by foreach
# defaults to the number of detected cores 
RA_output <- run_pathfindR(RA_input, n_processes = 2)

# to report the non-DEG active subnetwork genes
RA_output <- run_pathfindR(RA_input, list_active_snw_genes = TRUE)
```
For a full list of arguments, see `?run_pathfindR`.

The workflow consists of the following steps :

After input testing, the program attempts to convert any gene symbol that is not in the PIN to an alias symbol that is in the PIN. Next, active subnetwork search is performed via the selected algorithm. The available algorithms for active subnetwork search are: 

- Greedy Algorithm (based on Ideker et al. [^2]),
- Simulated Annealing Algorithm (based on Ideker et al. [^2]) and 
- Genetic Algorithm (based on Ozisik et al. [^3]).

Next, pathway enrichment analyses are performed using the genes in each of the active subnetworks. For this, up-to-date information on human gene sets from KEGG, Reactome, BioCarta and Gene Ontology were retrieved and is available for use within the package. The user may specify custom gene sets, including gene sets for non-human organisms, as described in the section "pathfindR Analysis with Custom Gene Sets".

During enrichment analyses, pathways with adjusted-p values larger than the `enrichment_threshold` (an argument of `run_pathfindR()`, defaults to 0.05) are discarded. The results of enrichment analyses over all active subnetworks are combined by keeping only the lowest adjusted-p value for each pathway. 

This process of active subnetwork search and enrichment analyses is repeated for a selected number of iterations (indicated by the `iterations` argument of `run_pathfindR()`), which is performed in parallel via the R package `foreach`. 

The wrapper function returns a data frame that contains the lowest and the highest adjusted-p values for each enriched pathway, as well as the numbers of times each pathway is encountered over all iterations. The first two rows of the example output of the pathfindR-enrichment workflow (performed on the rheumatoid arthritis data `RA_output`) is shown below:

```{r}
data("RA_output")
knitr::kable(head(RA_output, 2))
```

The function also creates an HTML report `results.html` that is saved in a directory, by default named `pathfindr_Results` but can be changed by changing the argument `output_dir`, under the current working directory. This report contains links to two other HTML files:

## 1. all_pathways.html
This document contains a table of the active subnetwork-oriented pathway enrichment results. Each enriched pathway name is linked to the visualization of that pathway, with the gene nodes colored according to their log-fold-change values. This table contains the same information as the returned data frame. Columns are:

- ID: KEGG ID of enriched pathway
- Pathway: Description the pathway
- Fold_Enrichment: Fold enrichment value for the pathway.
- occurrence: The number of times the pathway was found to be enriched over all iterations
- lowest_p: the lowest adjusted-p value of the pathway over all iterations
- highest_p: the highest adjusted-p value of the pathway over all iterations
- Up_regulated: the up-regulated genes involved in the pathway
- Down_regulated: the down-regulated genes involved in the pathway

## 2. genes_table.html
This document contains a table of converted gene symbols. Columns are:

- Old Symbol: the original gene symbol
- Converted Symbol: the alias symbol that was found in the PIN
- Change: the provided change value
- p-value: the provided adjusted p value

The document contains a second table of genes for which no interactions were identified (after checking for alias symbols).

# Pathway Clustering

<img src="./pw_clustering.png" width="650" height="221.464" />

For this workflow, the wrapper function `cluster_pathways()` is used. This function first calculates the pairwise kappa statistics between the terms in the input data frame. By default, the function performs hierarchical clustering of the terms using this kappa matrix, automatically determines the optimal number of clusters by maximizing the average silhouette width and returns a data frame with cluster assignments:

```{r hierarchical, fig.height=4, fig.width=8}
data("RA_output")
RA_clustered <- cluster_pathways(RA_output)
## First 2 rows of clustered terms data frame
knitr::kable(head(RA_clustered, 2))
## The 8 representative terms
knitr::kable(RA_clustered[RA_clustered$Status == "Representative", ])

# to display the heatmap of kappa statistics
RA_clustered <- cluster_pathways(RA_output, plot_hmap = TRUE, plot_clusters_graph = FALSE)

# to display the dendrogram and optimal clusters
RA_clustered <- cluster_pathways(RA_output, plot_dend = TRUE, plot_clusters_graph = FALSE)

# to change agglomeration method (default = "average")
RA_clustered <- cluster_pathways(RA_output, hclu_method = "centroid")
```

Alternatively, the `fuzzy` clustering method (as described by Huang et al.[^1]) can be used:

```{r fuzzy, fig.height=4, fig.width=8}
RA_clustered <- cluster_pathways(RA_output, method = "fuzzy")
```

# Pathway Scores per Sample

The function `calculate_pw_scores` can be used to calculate the pathway scores per sample. This allows the user to individually examine the scores and infer whether a pathway is activated or repressed in a given sample.

For a set of pathways $P = \{P_1, P_2, ... , P_k\}$, where each $P_i$ contains a set of genes, i.e. $P_i = \{g_1, g_2, ...\}$, the pathway score matrix $PS$ is defined as:

$PS_{p,s} = \frac{1}{k} \sum_{g \in P_p} GS_{g,s}$ for each pathway $p$ and for each sample $s$.

$GS$ is the gene score per sample matrix and is defined as:
$GS_{g,s} = (EM_{g,s} - \bar{x}_g) / sd_g$ where $EM$ is the expression matrix (columns are samples, rows are genes), $\bar{x}_g$ is the mean expression value of the gene and $sd_g$ is the standard deviation of the expression values for the gene.

An example application is provided below:

```{r, fig.height=4, fig.width=8}
## Pathway data frame
pws_table <- pathfindR::RA_clustered
# selecting "Representative" pathways for clear visualization
pws_table <- pws_table[pws_table$Status == "Representative", ]

## Expression matrix
exp_mat <- pathfindR::RA_exp_mat

## Vector of "Case" IDs
cases <- c("GSM389703", "GSM389704", "GSM389706", "GSM389708", 
           "GSM389711", "GSM389714", "GSM389716", "GSM389717", 
           "GSM389719", "GSM389721", "GSM389722", "GSM389724", 
           "GSM389726", "GSM389727", "GSM389730", "GSM389731", 
           "GSM389733", "GSM389735")

## Calculate pathway scores and plot heatmap
score_matrix <- calculate_pw_scores(pws_table, exp_mat, cases)
```

# pathfindR Analysis with Custom Gene Sets

It is possible to use `run_pathfindR()` with custom gene sets. Here we provide an example application with a gene set consisting of targets of two transcription factors.

We first load and prepare the gene sets:
```{r}
## CREB target genes
CREB_targets <- normalizePath(system.file("extdata/CREB.txt",
                                      package = "pathfindR"))
CREB_targets <- read.delim(CREB_targets, skip = 2, header = FALSE, stringsAsFactors = FALSE)
CREB_targets <- CREB_targets$V1

## MYC target genes
MYC_targets <- normalizePath(system.file("extdata/MYC.txt",
                                      package = "pathfindR"))
MYC_targets <- read.delim(MYC_targets, skip = 2, header = FALSE, stringsAsFactors = FALSE)
MYC_targets <- MYC_targets$V1

## Prep for use
custom_genes <- list(CREB_targets, MYC_targets)
names(custom_genes) <- c("TF1", "TF2")

custom_pathways <- c("CREB target genes", "MYC target genes")
names(custom_pathways) <- c("TF1", "TF2")
```

We next prepare an example input, because of the way we choose genes we expect significant enrichment for MYC targets (50 MYC target genes + 5 CREB target genes). Because this is only an example, we set the logFC values to 1.5 and the adj.P.Vals to 0.001:

``` {r}
set.seed(123)
selected_genes <- sample(MYC_targets, 50)
selected_genes <- c(selected_genes,
                    sample(CREB_targets, 5))
input <- data.frame(Gene.symbol = selected_genes,
                    logFC = 1.5,
                    adj.P.Val = 0.001)
head(input)
```

Finally, we run the wrapper script `run_pathfindR` using the custom genes as the gene sets:

``` {r eval=FALSE}
custom_result <- run_pathfindR(input, 
                               gene_sets = "Custom", 
                               custom_genes = custom_genes, 
                               custom_pathways = custom_pathways, 
                               iterations = 1)
```

``` {r include = FALSE}
data(custom_result)
```

``` {r}
knitr::kable(custom_result)
```

The most significantly enriched gene set is, unsurprisingly, "MYC target genes". Interestingly, through interaction information, pathfindR was also able to identify significant enrichment for "CREB target genes" as well.

**It is also possible to perform pathfindR analysis using non-human organism annotation. For this, custom gene sets and a custom PIN belonging to the species-of-interest must be specified.**


[^1]: Huang DW, Sherman BT, Tan Q, et al. The DAVID Gene Functional Classification Tool: a novel biological module-centric algorithm to functionally analyze large gene lists. Genome Biol. 2007;8(9):R183.

[^2]: Ideker T, Ozier O, Schwikowski B, Siegel AF. Discovering regulatory and signalling circuits in molecular interaction networks. Bioinformatics. 2002;18 Suppl 1:S233-40.

[^3]: Ozisik O, Bakir-Gungor B, Diri B, Sezerman OU. Active Subnetwork GA: A Two Stage Genetic Algorithm Approach to Active Subnetwork Search. Current Bioinformatics. 2017; 12(4):320-8. \doi 10.2174/1574893611666160527100444
