```{r setup, include=FALSE, echo=FALSE, eval=TRUE}

options(width = 300)
# options(encoding = 'UTF-8')
knitr::opts_chunk$set(
  fig.width = 7, fig.height = 5, 
  fig.path = 'figures/enrich-',
  fig.align = "center", 
  size = "tiny", 
  echo = TRUE, 
  eval = TRUE, 
  warning = FALSE, 
  message = FALSE, 
  results = TRUE, 
  comment = "")

options(scipen = 12) ## Max number of digits for non-scientific notation
# knitr::asis_output("\\footnotesize")

```


```{r libraries, echo=FALSE, eval=TRUE}

#### Load libraries ####
message("Loading libraries")
requiredLib <- c("knitr","gprofiler2")
for (lib in requiredLib) {
  if (!require(lib, character.only = TRUE)) {
    install.packages(lib, )
  }
  require(lib, character.only = TRUE)
}

library(gprofiler2)

```

## Goal of this practical

This practical aims at showing how to analyze functional enrichment on a list of genes with `gProfiler`.

## Study case : The Cancer Gene Atlas (TCGA) Breast Invasive Cancer (BIC) data

Description : <https://du-bii.github.io/study-cases/Homo_sapiens/TCGA_study-case/>

## References

- R package : <https://cran.r-project.org/web/packages/gprofiler2>

- Reference manual: <https://cran.r-project.org/web/packages/gprofiler2/gprofiler2.pdf>

- Vignette: <https://cran.r-project.org/web/packages/gprofiler2/vignettes/gprofiler2.html>

## Load (clusters of differentially) expressed genes

First, you need to load the gene list, which can be found on the github repository, in the following file. 

- 1000 top DEG genes: <https://raw.githubusercontent.com/DU-Bii/module-3-Stat-R/master/stat-R_2021/data/TCGA_BIC_subset/BIC_edgeR_DEG_top_1000_geneIDs.txt>


```{r data_loading}

#### Load the gene list ####
message("Loading gene list")
## URL of the data folder on github
data_folder_url <- 'https://raw.githubusercontent.com/DU-Bii/module-3-Stat-R/master/stat-R_2021/data/TCGA_BIC_subset/'

gene_list_file <- "BIC_edgeR_DEG_top_1000_geneIDs.txt"

## Trick : the file is read as a one-column table, and we transform it to a vector with unlist()
gene_list <- unlist(read.delim(file = file.path(data_folder_url, gene_list_file), header = FALSE))

## Check the length of the gene list
length(gene_list)

```


## Gene list functional enrichment analysis with gost

1. Use the `gost()` function to get the enrichment results.

- Note: we use FDR as multiple testing correction.

2. Visualize the results as an interactive plot

We can now tune the parameters to get some additional information about the evidence code and the list of genes at the intersection between the query set and the functional classes. Note that this generates a table with very large columns. 

3. Get a web link to share the results with your coworkers

4. Make a static plot 

5. Highlight the `GO:0071840` and `HPA:0300000` terms

6. Produce a custom table with choosen results

If you feel adventurous, you can try to do the same with multiple gene lists (for instance, clusters from the previous practical).

## Save your session info

For the sake of traceability, store the specifications of your R environment in the report, with the command `sessionInfo()`. This will indicate the version of R as well as of all the libraries used in this notebook. 

```{r session_info}

sessionInfo()

```
