---
title: "Annotation_and_Enrichment_REGIONS"
author: "Mikhail Dozmorov"
date: "Spring 2021"
output: html_document
---
```{r setup, echo=FALSE}
# Set up the environment
library(knitr)
opts_chunk$set(cache.path = "cache/", fig.path = "img/", cache = F, tidy = T, fig.keep = "high", echo = T, dpi = 300, out.width = 700)
options(replace.assign = TRUE, width = 220)
suppressMessages(library(pander))
panderOptions("table.split.table", Inf)
set.seed(1)
# Change to your working folder
setwd("/Users/mdozmorov/Documents/Work/Teaching/BIOS691_Cancer_Bioinformatics/static/slides/10_Epigenomics")
```

We load coordinates of age-associated CpG sites from BED file as GenomicRanges. See Steve Horvath, "[DNA methylation age of human tissues and cell types](https://doi.org/10.1186/gb-2013-14-10-r115)" 10 December 2013 for more details about the data.

```{r}
library(ChIPpeakAnno)
library(rtracklayer)

peaks <- import("data/Horvath_cpg_all.bed", format = "BED") # Convert to RangedData
```

Gene annotation data

```{r}
library(EnsDb.Hsapiens.v75)
annoData <- toGRanges(EnsDb.Hsapiens.v75)
annoData[1:2]
```

Annotate CpG sites by nearest/overlapping gene

```{r}
## keep the seqnames in the same style
seqlevelsStyle(peaks) <- seqlevelsStyle(annoData)
## do annotation by nearest TSS
anno <- annotatePeakInBatch(peaks, AnnotationData = annoData)
anno[1:2]
# A pie chart can be used to demonstrate the overlap features of the peaks.
pie1(table(anno$insideFeature))
```

Add gene symbol

```{r}
library(org.Hs.eg.db)
anno <- addGeneIDs(anno,
  orgAnn = "org.Hs.eg.db",
  feature_id_type = "ensembl_gene_id",
  IDs2Add = c("symbol")
)
head(anno)
```

From here, you can already perform gene set enrichment analysis on the identified genes, using, for example, [clusterProfiler](https://bioconductor.org/packages/clusterProfiler/). Try to do, and compare the results with what [GREAT](https://bioconductor.org/packages/rGREAT/) will return

# GREAT

Another approach for peak-to-gene enrichment analysis, directly on peaks.

```{r}
library(rGREAT)
job <- submitGreatJob(peaks)
tb <- getEnrichmentTables(job)
names(tb)
tb[["GO Molecular Function"]]
res <- plotRegionGeneAssociationGraphs(job)
```


