Going through ChIPseeker vignette, will modify as necessary for the actual analysis

ChipSeeker_vignette.R

library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(clusterProfiler)

# Gene data from UCSC hg19
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

# Peak file in bed format (these are MACS files)
files <- getSampleFiles()
files

# Read peak files (produces GRangesList)
peak <- readPeakFile(files[[4]])
peak

# Build peak coverage plot 
covplot(peak, weightCol="V5")

# Plot peak coverage on select chromosomes
covplot(peak, weightCol="V5", chrs=c("chr17", "chr18"), xlim=c(4.5e7, 5e7))

# Profile ChIP peaks binding to TSS regions
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
tagMatrix <- getTagMatrix(peak, windows=promoter)
tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")

# Average peak profile around the TSS
plotAvgProf(tagMatrix, xlim=c(-3000, 3000),
            xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency")

# Estimate confidence interval using Bootstrap method (very slow and can freeze the computer for a bit)
plotAvgProf(tagMatrix, xlim=c(-3000, 3000), conf = 0.95, resample = 1000)

# Peak annotation (by default promoter region is -3kb to +3kb.)
peakAnno <- annotatePeak(files[[4]], tssRegion=c(-3000, 3000),
                          TxDb=txdb, annoDb="org.Hs.eg.db")
peakAnno
# Default region priority for the annotation (this can be changed genomicAnnotationPriority argument)
# Promoter
# 5’ UTR
# 3’ UTR
# Exon
# Intron
# Downstream
# Intergenic

# Pie chart to visualize genomic annotations (do not use piecharts)
plotAnnoPie(peakAnno)

# better visualization of the same data with a bar graph
plotAnnoBar(peakAnno)

# Compare multiple datasets with one command
peakAnnoList <- lapply(files, function(x) annotatePeak(x, TxDb=txdb, tssRegion=c(-3000, 3000),
                                                       annoDb="org.Hs.eg.db"))
peakAnnoList
# Plot annotation bars simultaneusly for multiple datasets
plotAnnoBar(peakAnnoList)

# Compare distribution of distances to TSS between datasets
plotDistToTSS(peakAnnoList)