```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,
                      eval = FALSE,
                      message=FALSE,
                      error=FALSE,
                      fig.width = 12)
knitr::opts_knit$set(root.dir = '../')
library("webexercises")
```


#### Prepare environment

Set your working directory and load the libraries we need.

```{r eval=TRUE}
# set working directory
setwd("~/Documents/git_projects/ADDA_taller2/")

# install pacman if not already installed
if (!require("pacman")) install.packages("pacman")

# use pacman to load libraries
pacman::p_load(tidyverse,gprofiler2, scico)

```

## Background

From the previous exercises on differential gene expression, we looked at some comparisons that resulted in >300 genes that were over or under expressed. It is not trivial to make sense of the biological significance of that many genes. One way to try to distill this large list of genes into biological functions is to perform a functional enrichment analysis. Essentially, we want to use curated lists of genes that belong to a specific functional category and test whether some of these are over represented in our case. Perhaps the most widely used knowledge base of this kind, is the [Gene Ontology](http://geneontology.org/), but there are others, such as the [Kyoto Encyclopedia of Genes and Genomes](https://www.genome.jp/kegg/) and the [Reactome](https://reactome.org/). This is also called a `Functional Enrichment` or `Over Representation Analysis`

Let's say we are interested in a single biological process: melanin biosynthesis. Conveniently, a GO term for this exists: [GO:0042438](https://www.ebi.ac.uk/QuickGO/term/GO:0042438), as well as a _Xenopus_-specific Reactome pathway [R-XTR-5662702](https://reactome.org/content/detail/R-XTR-5662702). Essentially, what we want to do is see how many of our differentially expressed genes are either part of this term or not. This would result in a contingency table like this:
 
    
|                              | Melanin biosynthesis | Not melanin biosynthesis |
|------------------------------|:--------------------:|:------------------------:|
| Differentially expressed     | a                    | b                        |
| Not differentially expressed | c                    | d                        |
 
    
We can now perform a Fisher's exact test to get a p-value that would tell us whether genes associated with melanin biosynthesis are over represented in our list of differentially expressed genes. You can imagine however that we don't just want to test this single functional group, but many. We are therefore going to repeat this with all of the terms in the Gene Ontology for example. That many tests therefore requires a multiple-test correction, to get an adjusted p-value.
 
One final thing that we have to keep in mind is the _background_ or _universe_. This refers to the total list of genes that we will compare your set of genes to. Imagine we compare our set of genes that we obtained from skin samples to a full list of genes from the whole genome. In this case, it would be highly likely that we will get a lot of enriched terms or functions related to processes occurring in the skin. While this is not wrong, it is also not informative, because what we are interested in, is in _differences in enriched skin-specific functions_ between two conditions. We therefore ideally want to restrict our background set to only those genes that we realistically expect to find in our samples.
  
To perform such as test, we therefore need:

1. A list of genes/transcripts that are of biological interest (e.g. list of experimentally derived DEGs). These need to be annotated! That is, we need to know if these loci are genes with a known function (i.e. annotated).
2. A 'background' list of genes, that includes any gene that *could* have been differentially expressed in the analysis.
3. A curated database of biologically relevant categories (GO terms, gene pathways etc.), with associated genes.

Lets get to it!


## Load data

Let's load the results files we need.

```{r}
# the list of DEG results from the previous exercises
res<-readRDS("./results/deseq2_results.rds")

# the annotations
xtrop<-read_csv("data/PCU23_annotations_xtr105_genes.csv")
```

## Prepare gene sets and background

Ideally, the annotations of your genes should come from experimental evidence from your organism. If you work with mice, _Drosophila_ or humans for example, many functional enrichment analysis tools will be very easy to implement because they will automatically connect your lists of genes to annotations, background lists etc. This is unlikely going to exist if you work with non-model systems. We are therefore dependent on making our own annotations and lists.  

For our annotations, we are using BLAST results from querying our genes against the proteome of _Xenopus tropicalis_. This is a well studied frog species. It is still only distantly related to our focal species, but this is the best we can do. As we saw in the previous exercise, this results in many genes not having any annotations. These will be excluded unfortunately.

### Make gene sets

First we will use the same bit of code from the previous exercise to make a list of DEGs per comparison.

```{r}
extract_degs<-function(x) {
  return(
    x %>%
      as_tibble(rownames = "gene") %>%
      filter(padj<0.05) %>%
      pull(gene)
  )
}

# now extract all
sig_deg<-lapply(res, FUN=extract_degs)
str(sig_deg)
```

We now have to `convert` those to *Xenopus tropicalis* IDs, by extracting the matching annotations. Because we will do this multiple times over list itmes, it is cleaner to write a function first, then apply it to the list. 

```{r}
# make a function that extracts matching X. tropicalis IDs
extract_xtr<-function(x) {
  return(
      xtrop %>%
        filter(gene_id %in% x) %>%
        pull(xtr_pep_id_x) %>%
        unique()
  )
}

# apply function to list of Pelobates IDs
xtr_deg<-lapply(sig_deg, FUN=extract_xtr)
str(xtr_deg)
```


> _Question:_ Are the gene sets in the list of Xenopus genes the same length as the  original DEG list? How and why are they different?

`r hide("Show Answer")`

```{r}
sapply(sig_deg, length)
sapply(xtr_deg, length)
```

The gene sets are smaller. This is because not every gene for our organism could be annotated.
`r unhide()`

### Make background

There is some discussion about what makes a good background. Ideally, it should be the complete list of genes that *could* be differentially expressed. But what is this?


> _Question:_ A good background is:
> 
> a) All genes in the Xenopus tropicalis proteome
> b) All genes that could be annotated in the Pelobates genome
> c) All genes in the tissue-specific transcriptome

`r hide("Show Answer")`

It is not always clear, but in this case, I would argue that the tissue-specific transcriptome is the closest right answer.

We will use the full set of genes that were returned by `DESeq2`. This set should have filtered out genes that have low counts (i.e. unlikely to be expressed across any of our tissues/conditions).

`r unhide()`

We can use the same function from earlier to convert our list of _Pelobates_ IDs to _Xenopus_ peptide IDs.

```{r}
xtr_bg<-extract_xtr(rownames(res$bD_bV))
str(xtr_bg)
```

## Functional Enrichment Analysis

We are now ready to go! There are a number of software and R packages that let you perform functional enrichment analysis. Here, we will use [g:Profiler])(https://biit.cs.ut.ee/gprofiler/), because it plays particularly well with R and with Ensembl gene/peptide annotations.  

The analysis can be performed with a single command, even if our query is a list of multiple gene sets!  

An important thing to remember is that the associated R package `gprofiler2` is just an API, and the actual analysis will be performed on the g:Profiler server. This version is continuously updated, to match the updates Because of this, it is important to tell gprofiler which version of Ensembl you would like to use. Our annotations came from version 105, which is a slightly older version, so we have to make sure we set this correctly. We can do that by looking up the [archives](https://biit.cs.ut.ee/gprofiler/page/archives) on the g:Profiler homepage.

```{r}
# set base url:
set_base_url("https://biit.cs.ut.ee/gprofiler_archive3/e105_eg52_p16/")


# run the analysis
res_ora<-gost(multi_query = FALSE, # returns separate results tables for multiquery
              custom_bg = xtr_bg, # our background
              query=xtr_deg, # our list of gene sets
              organism="xtropicalis", # the organism our annotations belong to
              exclude_iea = FALSE, # include GO terms that were electronically assigned
              correction_method = "gSCS", # the recommended multiple testing correction.
              sources=c("GO:BP","GO:CC","GO:MF", "KEGG","REAC"), # the functional sets we are interested in 
              evcodes=FALSE, ## evcodes TRUE needed for downstream analysis like enrichment maps in Cytoscape, but this takes longer.
              significant= FALSE) # return all terms, not just the significant ones


# the results are stored as a "results" dataframe 
head(res_ora$result)
```


> _Question:_ What do these resuls look like and how many functionally enriched terms are there for each gene set?:

`r hide("Show Answer")`

The results returned tell us a few things:

```{r}
colnames(res_ora$result)
```

most importantly:
  
* the `term_id` and associated `term_name` which are the functional terms from the sources we requested (GO, KEGG, Reactome).
* the `p_value`, which is actually the adjusted p-value (misleading!), that the given term is enriched (over represented) in our query (gene set).
* the `term_size`, `query_size` and `intersect_size` which tell you how many genes make up the given term, how many genes were in your query and how many genes from both are overlapping.

`r unhide()`

> _Question:_ How many terms have been significantly enriched for each of the comparisons?

`r hide("Show Answer")`

We can use a p_value cutoff of 0.05 to see how many terms have been functionally enriched in each term.

```{r}
res_ora$result %>%
  filter(p_value<0.05) %>%
  group_by(query) %>%
  dplyr::count(query, sort=TRUE)
```

We see that the dorsal-ventral comparisons have the most enriched terms (`26` and `24`) and the dorsal-dorsal comparison has the fewest (`8`) This is not surprising given the small number of genes in that last set.

`r unhide()`

gprofiler also has a few visualization tools as well. For example an interactive Manhattan-style plot:

```{r fig.height=12, fig.fullwidth=TRUE}
gostplot(res_ora)
```

You should see that the Dorsal-Ventral comparisons look fairly similar. Let's take some time to explore this plot.

> _Question:_ What REACTOME pathway is enriched for all of the comparisons?

`r hide("Show Answer")`
* Melanin biosynthesis
`r unhide()`

Another useful plot to show enrichment results are dot plots. These are slightly more informative and make comparing across multiple sets more intuitive.  
  
We can easily make our custom dot plot using the gprofiler results tables and ggplot.

```{r fig.height=12, fig.fullwidth=TRUE}
res_ora$result %>%
  select(query,term_name, p_value, intersection_size, query_size,source) %>%
  filter(p_value<0.05) %>%
  mutate(GeneRatio=intersection_size/query_size) %>%
  arrange(GeneRatio) %>%
  mutate(term_name = factor(term_name, levels=unique(term_name))) %>%
  ggplot(aes(x=GeneRatio, y=term_name)) +
  geom_point(aes(color=p_value, size=intersection_size)) +
  ylab("") +
  scale_color_scico(palette = "batlow", direction = 1) +
  facet_grid(source~query,scales = "free_y",space = "free") +
  theme_bw()
```

> _Question:_ After having done the differential gene expression analysis and performed a functional enrichment analysis, what have we learned about the pigmentation plasticity of tadpoles?

`r hide("Show Answer")`

We have learned for example that:
  
* Melanin differential biosynthesis is relevant for all comparisons (even white ventral and black ventral?).
* Dorsal-ventral differentiation involves other, non-pigment related processes such as "system development" and "extracellular matrix".

`r unhide()`


## Concluding remarks

* A functional enrichment analysis is useful when you have many differentially expressed genes and you would like to make biological sense of this.
* These are 100% dependent on a) accurate annotations and b) curated functional databases. Both of these can be very problematic when working with non-model organisms.
* The analysis also tends to be very sensitive to a number of factors, such as what the background set is, whether electronically annotations are included and what significance thresholds to use.

