## ----style, echo=FALSE, results='asis', message=FALSE--------------------
BiocStyle::markdown()
knitr::opts_chunk$set(tidy         = FALSE,
                      warning      = FALSE,
                      message      = FALSE)

## ----echo=FALSE, results='hide', message=FALSE---------------------------
library(GenomicFeatures)
library(GenomicRanges)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(org.Hs.eg.db)
library(clusterProfiler)
library(ReactomePA)
library(ChIPseeker)

## ------------------------------------------------------------------------
## loading packages
library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
library(clusterProfiler)

## ------------------------------------------------------------------------
files <- getSampleFiles()
print(files)
peak <- readPeakFile(files[[4]])
peak

## ----fig.height=8, fig.width=10------------------------------------------
covplot(peak, weightCol="V5")

## ----fig.height=4, fig.width=10------------------------------------------
covplot(peak, weightCol="V5", chrs=c("chr17", "chr18"), xlim=c(4.5e7, 5e7))

## ------------------------------------------------------------------------
## promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
## tagMatrix <- getTagMatrix(peak, windows=promoter)
##
## to speed up the compilation of this vignettes, we use a precalculated tagMatrix
data("tagMatrixList")
tagMatrix <- tagMatrixList[[4]]

## ----fig.cap="Heatmap of ChIP peaks binding to TSS regions", fig.align="center", fig.height=12, fig.width=4----
tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")

## ----eval=FALSE----------------------------------------------------------
#  peakHeatmap(files[[4]], TxDb=txdb, upstream=3000, downstream=3000, color="red")

## ----eval=TRUE, fig.cap="Average Profile of ChIP peaks binding to TSS region", fig.align="center", fig.height=4, fig.width=7----
plotAvgProf(tagMatrix, xlim=c(-3000, 3000), 
            xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency")

## ----eval=FALSE----------------------------------------------------------
#  plotAvgProf2(files[[4]], TxDb=txdb, upstream=3000, downstream=3000,
#               xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency")

## ----fig.cap="Average Profile of ChIP peaks binding to TSS region", fig.align="center", fig.height=4, fig.width=7, eval=F----
#  plotAvgProf(tagMatrix, xlim=c(-3000, 3000), conf = 0.95, resample = 1000)

## ------------------------------------------------------------------------
peakAnno <- annotatePeak(files[[4]], tssRegion=c(-3000, 3000), 
                         TxDb=txdb, annoDb="org.Hs.eg.db")

## ----fig.cap="Genomic Annotation by pieplot", fig.align="center", fig.height=6, fig.width=8----
plotAnnoPie(peakAnno)

## ----fig.cap="Genomic Annotation by barplot", fig.align="center", fig.height=4, fig.width=10----
plotAnnoBar(peakAnno)

## ----fig.cap="Genomic Annotation by vennpie", fig.align="center", fig.height=8, fig.width=11----
vennpie(peakAnno)

## ----eval=F, fig.cap="Genomic Annotation by upsetplot", fig.align="center", fig.height=8, fig.width=12----
#  upsetplot(peakAnno)

## ----eval=F, fig.cap="Genomic Annotation by upsetplot", fig.align="center", fig.height=8, fig.width=12----
#  upsetplot(peakAnno, vennpie=TRUE)

## ----fig.cap="Distribution of Binding Sites", fig.align="center", fig.height=2, fig.width=6----
plotDistToTSS(peakAnno, 
              title="Distribution of transcription factor-binding loci\nrelative to TSS")

## ----fig.width=8, fig.height=5-------------------------------------------
library(ReactomePA)

pathway1 <- enrichPathway(as.data.frame(peakAnno)$geneId)
head(pathway1, 2)

gene <- seq2gene(peak, tssRegion = c(-1000, 1000), flankDistance = 3000, TxDb=txdb)
pathway2 <- enrichPathway(gene)
head(pathway2, 2)
dotplot(pathway2)

## ----eval=TRUE, fig.cap="Average Profiles of ChIP peaks among different experiments", fig.align="center", fig.height=4, fig.width=6----
## promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
## tagMatrixList <- lapply(files, getTagMatrix, windows=promoter)
## 
## to speed up the compilation of this vigenette, we load a precaculated tagMatrixList
data("tagMatrixList")
plotAvgProf(tagMatrixList, xlim=c(-3000, 3000))

## ----eval=FALSE, fig.cap="Average Profiles of ChIP peaks among different experiments", fig.align="center", fig.height=7, fig.width=6----
#  plotAvgProf(tagMatrixList, xlim=c(-3000, 3000), conf=0.95,resample=500, facet="row")

## ----eval=TRUE, fig.cap="Heatmap of ChIP peaks among different experiments", fig.align="center", fig.height=8, fig.width=7----
tagHeatmap(tagMatrixList, xlim=c(-3000, 3000), color=NULL)

## ------------------------------------------------------------------------
peakAnnoList <- lapply(files, annotatePeak, TxDb=txdb, 
                       tssRegion=c(-3000, 3000), verbose=FALSE)

## ----fig.cap="Genomic Annotation among different ChIPseq data", fig.align="center", fig.height=4, fig.width=6----
plotAnnoBar(peakAnnoList)

## ----fig.cap="Distribution of Binding Sites among different ChIPseq data", fig.align="center", fig.height=5, fig.width=8----
plotDistToTSS(peakAnnoList)

## ----fig.width=8.5, fig.height=8.5---------------------------------------
genes = lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)
names(genes) = sub("_", "\n", names(genes))
compKEGG <- compareCluster(geneCluster   = genes, 
                         fun           = "enrichKEGG",
                         pvalueCutoff  = 0.05, 
                         pAdjustMethod = "BH")
plot(compKEGG, showCategory = 15, title = "KEGG Pathway Enrichment Analysis")

## ----fig.cap="Overlap of annotated genes", fig.align="center", fig.height=7, fig.width=7----
genes= lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)
vennplot(genes)

## ------------------------------------------------------------------------
p <- GRanges(seqnames=c("chr1", "chr3"), 
             ranges=IRanges(start=c(1, 100), end=c(50, 130)))
shuffle(p, TxDb=txdb)

## ------------------------------------------------------------------------
enrichPeakOverlap(queryPeak     = files[[5]], 
                  targetPeak    = unlist(files[1:4]), 
                  TxDb          = txdb, 
                  pAdjustMethod = "BH", 
                  nShuffle      = 50, 
                  chainFile     = NULL,
                  verbose       = FALSE)

## ------------------------------------------------------------------------
getGEOspecies()

## ------------------------------------------------------------------------
getGEOgenomeVersion()

## ------------------------------------------------------------------------
hg19 <- getGEOInfo(genome="hg19", simplify=TRUE)
head(hg19)

## ----eval=FALSE----------------------------------------------------------
#  downloadGEObedFiles(genome="hg19", destDir="hg19")

## ----eval=FALSE----------------------------------------------------------
#  gsm <- hg19$gsm[sample(nrow(hg19), 10)]
#  downloadGSMbedFiles(gsm, destDir="hg19")

## ----echo=FALSE----------------------------------------------------------
sessionInfo()

