library(ChIPseeker)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
library(clusterProfiler)
library(org.Hs.eg.db)


peak <- readPeakFile("Mec1_ATM_shRNA_c_deseq2_sig_up_only.bed")

peakAnno <- annotatePeak(peak, tssRegion=c(-2000, 2000),
                         TxDb=txdb, annoDb="org.Hs.eg.db")

promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
tagMatrix <- getTagMatrix(peak, windows=promoter)
tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")

plotAvgProf(tagMatrix, xlim=c(-3000, 3000),
            xlab="Genomic Region (5'->3')", ylab = "Read Count Frequency")

plotAnnoPie(peakAnno)

plotAnnoBar(peakAnno)

vennpie(peakAnno)

upsetplot(peakAnno)

upsetplot(peakAnno, vennpie=TRUE)

plotDistToTSS(peakAnno,
              title="Distribution of transcription factor-binding loci\nrelative to TSS")



genes = (as.data.frame(peakAnno)$geneId)
#genes = na.omit(genes)
#names(genes) = sub("_", "\n", names(genes))


kk <- enrichKEGG(genes, organism = "hsa", pvalueCutoff=0.05, pAdjustMethod="none", 
                 qvalueCutoff=0.2)

head(summary(kk))
dotplot(kk, showCategory=25)


ego <- clusterProfiler::enrichGO(gene          = genes,
                                 OrgDb         = org.Hs.eg.db,
                                 ont           = "BP",
                                 pAdjustMethod = "BH",
                                 pvalueCutoff  = 0.05,
                                 qvalueCutoff  = 0.2, 
                                 readable      = TRUE)
head(summary(ego)[,-8])
dotplot(ego, showCategory=20)
clusterProfiler::dotplot(ego, showCategory=15)
cnetplot(ego, categorySize="geneNum", foldChange=geneList)