```{r setup, include=FALSE}
library(org.Mm.eg.db)
library(msigdbr)
library(clusterProfiler)
library(pathview)
library(tidyverse)
```

## Exercise 1 - pathview

Load the required packages and data for Day 11 if you have not already done so. 

```{r prepareORAData,}
library(msigdbr)
library(clusterProfiler)
library(pathview)
library(tidyverse)

shrink.d11 <- readRDS("RObjects/Shrunk_Results.d11.rds")
```

> 1. Use `pathview` to export a figure for "mmu04659"or "mmu04658", but this 
> time only use genes that are statistically significant at FDR < 0.01

```{r solution1}
logFC <- shrink.d11 %>% 
  drop_na(FDR, Entrez) %>% 
  filter(FDR < 0.01) %>% 
  pull(logFC, Entrez) 

pathview(gene.data = logFC, 
         pathway.id = "mmu04659", 
         species = "mmu", 
         limit = list(gene=5, cpd=1))
```

mmu04659.pathview.png:

![mmu04659 - Th17 cell differentiation](images/mmu04659_pathview.png)

## Exercise 2 - GO term enrichment analysis

> `clusterProfiler` can also perform over-representation analysis on GO terms.
> using the commmand `enrichGO`. Look at the help page for the command 
> `enrichGO` (`?enrichGO`) and have a look at the instructions in the 
> [clusterProfiler book](http://yulab-smu.top/clusterProfiler-book/chapter5.html#go-over-representation-test).
>
> 1. Run the over-representation analysis for GO terms 
>     - Use genes that have an adjusted p-value (FDR) of less than 0.01 and
> an absolute fold change greater than 2.    
>     - For this analysis you can use Ensembl IDs rather then Entrez
>     - You'll need to provide the background (`universe`) genes, this should be
> all the genes in our analysis.
>     - The mouse database package is called `org.Mm.eg.db`. You'll need to load
>  it using `library` before running the analysis.    
>     - As we are using Ensembl IDs, you'll need to set the `keyType` 
> parameter in the `enrichGO` command to indicate this.
>     - Only test terms in the "Biological Processes" ontology 

```{r solution2}
library(org.Mm.eg.db)

sigGenes <-  shrink.d11 %>% 
    drop_na(FDR) %>% 
    filter(FDR < 0.01 & abs(logFC) > 1) %>% 
    pull(GeneID)

universe <- shrink.d11$GeneID

ego <- enrichGO(gene          = sigGenes, 
                universe      = universe,
                OrgDb         = org.Mm.eg.db,
                keyType       = "ENSEMBL",
                ont           = "BP",
                pvalueCutoff  = 0.01,
                readable      = TRUE)
```

> 2. Use the `dotplot` function to visualise the results.

```{r solution2b, fig.width=12, fig.height=8}
dotplot(ego, font.size = 14)
```



## Exercise 3 - GSEA

> Another common way to rank the genes is to order by pvalue, but also, sorting
> so that upregulated genes are at the start and downregulated at the end - 
> you can do this combining the sign of the fold change and the pvalue. 

First load the pathway details if you have not already done so.

```{r solution3_GSEA_a}
library(msigdbr)
term2gene <- msigdbr(species = "Mus musculus", category = "H") %>% 
  select(gs_name, entrez_gene)
term2name <- msigdbr(species = "Mus musculus", category = "H") %>% 
  select(gs_name, gs_description) %>% 
  distinct()
```

> 1. Rank the genes by statistical significance - you will need to create
> a new ranking value using `-log10({p value}) * sign({Fold Change})`.

```{r solution3_GSEA_1}
# rank genes
rankedGenes.e11 <- shrink.d11 %>%
  drop_na(Entrez, pvalue, logFC) %>%
  mutate(rank = -log10(pvalue) * sign(logFC)) %>%
  arrange(desc(rank)) %>%
  pull(rank, Entrez)
```


> 2. Run `fgsea` using the new ranked genes and the Hallmark pathways.

```{r solution3_GSEA_2}
# conduct analysis:
gseaRes.e11 <- GSEA(rankedGenes.e11,
                TERM2GENE = term2gene,
                TERM2NAME = term2name,
                pvalueCutoff = 1.00, 
                minGSSize = 15,
                maxGSSize = 500)
```

View the results:

```{r top10GseaPrint, eval=FALSE}
as_tibble(gseaRes.e11) %>% 
  arrange(desc(abs(NES))) %>% 
  top_n(10, wt=-p.adjust) %>% 
  select(-core_enrichment) %>%
  mutate(across(c("enrichmentScore", "NES"), round, digits=3)) %>% 
  mutate(across(c("pvalue", "p.adjust", "qvalues"), scales::scientific))
```

```{r top10GseaPrintactual, echo=FALSE}
as_tibble(gseaRes.e11) %>% 
  arrange(desc(abs(NES))) %>% 
  top_n(10, wt=-p.adjust) %>% 
  select(-core_enrichment) %>%
  mutate(across(c("enrichmentScore", "NES"), round, digits=3)) %>% 
  mutate(across(c("pvalue", "p.adjust", "qvalues"), scales::scientific)) %>% 
  DT::datatable(option=list(dom='t'))
```

> 3. Conduct the same analysis for the day 33 Infected vs Uninfected contrast.

```{r solution3_GSEA_3}

# read d33 data in:
shrink.d33 <- readRDS("RObjects/Shrunk_Results.d33.rds")

# rank genes
rankedGenes.e33 <- shrink.d33 %>%
  drop_na(Entrez, pvalue, logFC) %>%
  mutate(rank = -log10(pvalue) * sign(logFC)) %>%
  arrange(desc(rank)) %>%
  pull(rank,Entrez)

# perform analysis
gseaRes.e33 <- GSEA(rankedGenes.e33,
                TERM2GENE = term2gene,
                TERM2NAME = term2name,
                pvalueCutoff = 1.00, 
                minGSSize = 15,
                maxGSSize = 500)
```

View the results:

```{r top10GseaPrint33, eval=FALSE}
as_tibble(gseaRes.e33) %>% 
  arrange(desc(abs(NES))) %>% 
  top_n(10, wt=-p.adjust) %>% 
  select(-core_enrichment) %>%
  mutate(across(c("enrichmentScore", "NES"), round, digits=3)) %>% 
  mutate(across(c("pvalue", "p.adjust", "qvalues"), scales::scientific))
```

```{r top10GseaPrintactual33, echo=FALSE}
as_tibble(gseaRes.e33) %>% 
  arrange(desc(abs(NES))) %>% 
  top_n(10, wt=-p.adjust) %>% 
  select(-core_enrichment) %>%
  mutate(across(c("enrichmentScore", "NES"), round, digits=3)) %>% 
  mutate(across(c("pvalue", "p.adjust", "qvalues"), scales::scientific)) %>% 
  DT::datatable(option=list(dom='t'))
```
