---
title: "Using fgseaMultilevel function"
date: "2018-11-14"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Using fgseaMultilevel function}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---
`fgseaMultilevel` is a new function in `fgsea` package. The basis of this function is the adaptive multilevel splitting Monte Carlo approach. This approach allows us to exceed the results of simple sampling and calculate arbitrarily small P-value.

## Loading necessary libraries

```{r message=FALSE}
library(fgsea)
library(ggplot2)
library(BiocParallel)
register(SerialParam())
```

## Quick run

Loading example pathways and gene-level statistics:

```{r}
data(examplePathways)
data(exampleRanks)
```

Running fgseaMultilevel:
```{r}
fgseaMultilevelRes <- fgseaMultilevel(pathways = examplePathways, 
                                      stats = exampleRanks,
                                      minSize=15,
                                      maxSize=500)
```

The resulting table contains enrichment scores and p-values:
```{r}
head(fgseaMultilevelRes[order(pval), ])
```

The output table is similar in structure to the table obtained using the `fgsea` function, except for the `log2err` column. This column corresponds to the standard deviation of the P-value logarithm.

## Function argument

 - `pathways` - List of gene sets to check.
 - `stats` - Named vector of gene-level stats. Names should be the same as in `pathways`
 - `sampleSize` - The size of randomly generated gene sets, where each set has pathway size.
 - `minSize` - Minimal size of a gene set to test. All pathways below the threshold are excluded.
 - `maxSize` - Maximal size of a gene set to test. All pathways above the threshold are excluded.
 - `absEps` - This parameter sets the boundary for calculating the p value.
 - `nproc` - If not equal to zero sets BPPARAM to use nproc workers (default = 0).
 - `BPPARAM` - Parallelization parameter used in bplapply.
 
## Briefly about the approach
Let $q$ be the initial set of genes of size $k$ for which we calculate the P-value.
Then the general scheme of the algorithm is as follows:

 - Random sets of $k$-size genes are generated in an amount equal to `sampleSize`. For each set of genes, an enrichment score (ES) is calculated. The median value of the ES for a given set of genes is then determined. 
 - The next step generates sets of genes with ES exceeding the median value obtained in the previous step. For newly generated samples, the median ES value is determined, and so on.
 
Let $m_{k,j}$ be the median value of ES for set of $k$-size genes obtained at the $j$ step. As a result, at some step, we obtain that the following relation will be satisfied for the initial set of genes:
$$m_{k,j} \leq ES(q) < m_{k,j+1}$$
Then the P-value is in the range from $2^{-j-1}$ to $2^{-j}$. In addition, we can achieve even better accuracy by using not only the medians themselves, but also the intermediate data.

## Compare `fgsea` vs `fgseaMultilevel`
Compare `fgsea` with `fgseaMultilevel` on example data.

```{r}
set.seed(42)

fgseaRes <- fgsea(pathways=examplePathways, stats=exampleRanks,
                  minSize=15, maxSize=500, nperm=10000)

fgseaMultilevelRes <- fgseaMultilevel(pathways=examplePathways, 
                                      stats=exampleRanks,
                                      minSize=15, maxSize=500, 
                                      sampleSize=100)

logPvalsDf <- data.frame(fgseaData=-log10(fgseaRes$pval),
                         fgseaMultilevelData=-log10(fgseaMultilevelRes$pval))
```

The following figure compares the P-value calculated using two approaches.

```{r fig.width=7, fig.height=5}
ggplot(logPvalsDf, aes(x=fgseaData, y=fgseaMultilevelData)) +
  geom_point(color='royalblue4', size=3) +
  xlab("Fgsea: -log10(P-value)") +
  ylab("FgseaMultilevel: -log10(P-value)") +
  geom_abline(size=1, color='orangered2') +
  labs(title="Fgsea Vs FgseaMultilevel")
```

It can be seen that up to a certain value, which is characterized by the parameter $\dfrac{1}{nperm}$, both algorithms give similar results. In this case, the accuracy of the P-value determination by the gsea algorithm is limited by the parameter  $\dfrac{1}{nperm}$, while the fgseaMultilevel does not have these restrictions.