

```{r, echo = FALSE}
knitr::opts_chunk$set(
  eval=FALSE
)
```


```{r}
load.libs <- c(
  "ggplot2",
  "EnhancedVolcano",
  "DOSE",
  "GO.db",
  "GSEABase",
  "org.Hs.eg.db",
  "clusterProfiler",
  "dplyr",
  "tidyr",
  "ggplot2",
  "stringr",
  "RColorBrewer",
  "rWikiPathways",
  "RCy3")
options(install.packages.check.source = "no")
options(install.packages.compile.from.source = "never")
if (!require("pacman")) install.packages("pacman"); library(pacman)
p_load(load.libs, update = TRUE, character.only = TRUE)
instatus <- sapply(load.libs,require,character.only = TRUE)
if(all(status)){
    print("SUCCESS: You have successfully installed and loaded all required libraries.")
} else{
    cat("ERROR: One or more libraries failed to install correctly. Check the following list for FALSE cases and try again...\n\n")
    status
}
```

#load your gene list
```{r}
muscle.expr <- read.csv("/Volumes/SSD/predator/DE/SCcav1_cutoff_2.csv", header = T,stringsAsFactors = F)
nrow(muscle.expr)
head(muscle.expr)
```

Now let's prepare up- and down-regulated gene lists using some conventional criteria.
```{r}
up.genes <- muscle.expr[muscle.expr$avg_logFC > cutFC & muscle.expr$p_val < cutP, 1] 
dn.genes <- muscle.expr[muscle.expr$avg_logFC < -cutFC & muscle.expr$p_val < cutP, 1]
#bkgd.genes <- lung.expr[,1]
```

# Enrichment
```{r}
up.genes.entrez <- clusterProfiler::bitr(up.genes,fromType = "SYMBOL",toType = "ENTREZID",OrgDb = org.Hs.eg.db)#specify species change if mouse (Hs)
cat("\n\nWhich column contains my new Entrez IDs?\n")
head(up.genes.entrez)
```

```{r}
keytypes(org.Hs.eg.db)
```

Let's convert our other lists to Entrez IDs:
```{r}
dn.genes.entrez <- bitr(dn.genes,fromType = "SYMBOL",toType = "ENTREZID",OrgDb = org.Hs.eg.db)
```

## Gene Ontology
```{r}
ego.up <- clusterProfiler::enrichGO(
        gene     = up.genes.entrez[[2]], #provide significant genes you want to analyse
        #universe = bkgd.genes.entrez[[2]],
        OrgDb    = org.Hs.eg.db,
        ont      = "BP", #BP biological process we are using
        pAdjustMethod = "fdr",
        pvalueCutoff = 0.05, #p.adjust cutoff
        readable = TRUE)

head(ego.up,10)
```

```{r}
write.csv(ego.up,file = "CAV1up.csv")
```

Tables are fine, but we're in R to see plots. Conveniently, **clusterProfilier** provides a variety of plots with default settings:
```{r}
barplot(ego.up, showCategory = 20)
dotplot(ego.up, showCategory = 20)
emapplot(ego.up, showCategory = 20)
```

#downregulated genes

```{r}
ego.dn <- clusterProfiler::enrichGO(
        gene     = dn.genes.entrez[[2]], #provide significant genes you want to analyse
        #universe = bkgd.genes.entrez[[2]],
        OrgDb    = org.Hs.eg.db,
        ont      = "BP", #BP biological process we are using
        pAdjustMethod = "fdr",
        pvalueCutoff = 0.05, #p.adjust cutoff (https://github.com/GuangchuangYu/clusterProfiler/issues/104)
        readable = TRUE)

head(ego.dn,10)
```

## WikiPathways

```{r}
wp.hs.gmt <- rWikiPathways::downloadPathwayArchive(organism="Homo sapiens", format = "gmt")

# supported organisms (and case-dependent spelling)...
listOrganisms()
```

Now that we've got the latest GMT file for human pathways, we can feed it into **clusterProfiler** and generate the two dataframes we need for *enricher*...
```{r}
wp2gene <- clusterProfiler::read.gmt(wp.hs.gmt)
wp2gene <- wp2gene %>% tidyr::separate(ont, c("name","version","wpid","org"), "%")
wpid2gene <- wp2gene %>% dplyr::select(wpid,gene) #TERM2GENE
wpid2name <- wp2gene %>% dplyr::select(wpid,name) #TERM2NAME
wpid2gene
wpid2name
```

Perform pathway analysis...
```{r}
ewp.up <- clusterProfiler::enricher(
        up.genes.entrez[[2]],
        #universe = bkgd.genes.entrez[[2]],
        pAdjustMethod = "fdr",
        pvalueCutoff = 0.05, #p.adjust cutoff
        TERM2GENE = wpid2gene,
        TERM2NAME = wpid2name)

head(ewp.up)
```

For some reason, *enricher* doesn't automatically add gene symbols to the result object, but there is a handy function in **DOSE** that does...
```{r}
ewp.up <- DOSE::setReadable(ewp.up, org.Hs.eg.db, keyType = "ENTREZID")
head(ewp.up)
```

```{r}
barplot(ewp.up, showCategory = 10)
```

do the same analysis for down-regulated genes...
```{r}
ewp.dn <- enricher(
  dn.genes.entrez[[2]],
  #universe = bkgd.genes[[2]],  #hint: comment out to get any results for demo
  pAdjustMethod = "fdr",
  pvalueCutoff = 0.1, #p.adjust cutoff
  TERM2GENE = wpid2gene,
  TERM2NAME = wpid2name)
ewp.dn <- DOSE::setReadable(ewp.dn, org.Hs.eg.db, keyType = "ENTREZID")
barplot(ewp.dn, showCategory = 10)
```

