```{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)
```

We load selected age-associated peaks from BED file, and convert them to RangedData
```{r}
suppressMessages(library(ChIPpeakAnno))
peaks <- read.table("data/mmc2_Model_All_data.bed", sep="\t", header=F, stringsAsFactors=F)
# peaks <- read.table("data/Alisch_et_al_Sup_Table2_all.bed", sep="\t", header=F, stringsAsFactors=F)
# peaks <- read.table("data/Supp_Table_3.bed", sep="\t", header=F, stringsAsFactors=F)
peaks.bed <- BED2RangedData(peaks) # Convert to RangedData
```

We pull transcript annotations from BiomaRt, and annotate our peaks in relation to them. We also peek inside how the annotation looks like.
```{r}
library(biomaRt)
mart<-useMart(biomart="ensembl",dataset="hsapiens_gene_ensembl")
# myFullAnnotation = getAnnotation(mart, featureType=c("TSS","miRNA", "Exon", "5utr", "3utr", "ExonPlusUtr", "transcript")) # An example in relation to what we can get annotations
myFullAnnotation = getAnnotation(mart, featureType=c("TSS")) # Getting annotations is time consuming
annotatedPeak = annotatePeakInBatch(peaks.bed, AnnotationData=myFullAnnotation)
pander(head(as.data.frame(annotatedPeak))) # An overview of how annotation results look like
colnames(as.data.frame(annotatedPeak)) # Check column names in the annotation results
```

So, we have "distancetoFeature" column, and can subset the annotated peaks by it. For convenience, we use designated function for it, finding peaks within 10K around transcripts.  Let us also check where they are located.

```{r}
# t1 <- findOverlappingPeaks(group1.bed, annotatedPeak, maxgap=2000, select="first", NameOfPeaks1 = "peaks", NameOfPeaks2 = "TSS", annotate=1) # Very time consuming for lagre list of peaks
t1 <- readRDS("group1.t1.RDS")
names(t1) # Names of the columns

overlappingPeaks <- t1$OverlappingPeaks # overlapping peaks with maximum gap 2kb
pie(table(overlappingPeaks$overlapFeature)) # and a pie graph describing the distribution of relative position of peaks1 to peaks2 for overlapping peaks

head(as.data.frame(t1$MergedPeaks)) # merged overlapping peaks
head(as.data.frame(t1$Peaks1withOverlaps)) # the peaks in peaks1 that overlaps with peaks in peaks2
head(as.data.frame(t1$Peaks2withOverlaps)) # the peaks in peaks2 that overlap with peaks in peaks1

group1.venn <- makeVennDiagram(RangedDataList(group1.bed, myFullAnnotation), NameOfPeaks = c("peaks", "TSS"), maxgap=0, minoverlap=1, totalTest = 100, cex=1, counts.col= "red", useFeature=FALSE)
```

The most interesting is to understand what, if anything, those genes are doing. We check their biological processes and molecular functions.

```{r}
library(org.Hs.eg.db)
enrichedGO = getEnrichedGO(annotatedPeak, orgAnn = 'org.Hs.eg.db', maxP = 0.01, multiAdj = FALSE, minGOterm = 10, multiAdjMethod = 'none' ) # GO enrichment analysis is time consuming
pander(unique(enrichedGO$bp[, c(-1, -4, -8:-10)][1:10, ]))
pander(unique(enrichedGO$mf[, c(-1, -4, -8:-10)][1:10, ]))
```

And, KEGG canonical pathways.

```{r}
library(KEGG.db)
enrichedPATH <- getEnrichedPATH(annotatedPeak, orgAnn= 'org.Hs.eg.db', pathAnn='KEGG.db', feature_id_type="ensembl_gene_id", maxP=0.01, minPATHterm=10, multiAdjMethod="none")
pander(unique(enrichedPATH[, c(9,3:5)]))
```

Let us also check what genes are actually there:

```{r}
annotatedPeakGene <- addGeneIDs(annotatedPeak,"org.Hs.eg.db",c("symbol"))
pander(head(as.data.frame(annotatedPeakGene)))
```

