###################################################
### chunk number 1: 
###################################################
library(ChIPpeakAnno)
data(myPeakList)
data(TSS.human.NCBI36)
annotatedPeak = annotatePeakInBatch(myPeakList[1:6,], AnnotationData=TSS.human.NCBI36)
as.data.frame(annotatedPeak)


###################################################
### chunk number 2: 
###################################################
myPeak1 = RangedData(IRanges(start=c(967654, 2010897, 2496704, 3075869, 3123260 ,3857501,201089,1543200,1557200,1563000,1569800,167889600), end= c(967754, 2010997, 2496804, 3075969, 3123360 ,3857601,201089,1555199,1560599,1565199,1573799,167893599),names=c("Site1", "Site2", "Site3", "Site4", "Site5", "Site6", "Site7","Site8","Site9","Site10","Site11","Site12")), space=c("1", "2", "3", "4", "5", "6","2","6","6","6","6","5"))

TFbindingSites = RangedData(IRanges(start=c(967659, 2010898, 2496700, 3075866, 3123260 ,3857500, 96765, 201089, 249670, 307586, 312326 ,385750,1549800,1554400,1565000,1569400,167888600), end=c(967869, 2011108, 2496920,3076166,3123470,3857780, 96985, 201299, 249890, 307796,312586,385960,1550599,1560799,1565399,1571199,167888999), names=c("t1", "t2", "t3", "t4", "t5", "t6","t7", "t8", "t9", "t10", "t11", "t12","t13","t14","t15","t16","t17")), space=c("1", "2", "3", "4", "5", "6","1", "2", "3", "4", "5", "6","6","6","6","6","5"), strand=
c(1,1,1,1,1,1,-1,-1,-1,-1,-1,-1,1,1,1,1,1))
annotatedPeak2 = annotatePeakInBatch(myPeak1, AnnotationData=TFbindingSites)
pie(table(as.data.frame(annotatedPeak2)$insideFeature))
as.data.frame(annotatedPeak2)


###################################################
### chunk number 3: 
###################################################
peaks1 = RangedData(IRanges(start=c(967654, 2010897, 2496704, 3075869, 3123260 ,3857501,201089,1543200,1557200,1563000,1569800,167889600), end= c(967754, 2010997, 2496804, 3075969, 3123360 ,3857601,201089,1555199,1560599,1565199,1573799,167893599),names=c("Site1", "Site2", "Site3", "Site4", "Site5", "Site6", "Site7","Site8","Site9","Site10","Site11","Site12")), space=c("1", "2", "3", "4", "5", "6","2","6","6","6","6","5"),strand=as.integer(1))

peaks2 = RangedData(IRanges(start=c(967659, 2010898, 2496700, 3075866, 3123260 ,3857500, 96765, 201089, 249670, 307586, 312326 ,385750,1549800,1554400,1565000,1569400,167888600), end=c(967869, 2011108, 2496920,3076166,3123470,3857780, 96985, 201299, 249890, 307796,312586,385960,1550599,1560799,1565399,1571199,167888999), names=c("t1", "t2", "t3", "t4", "t5", "t6","t7", "t8", "t9", "t10", "t11", "t12","t13","t14","t15","t16","t17")), space=c("1", "2", "3", "4", "5", "6","1", "2", "3", "4", "5", "6","6","6","6","6","5"), strand=
c(1,1,1,1,1,1,-1,-1,-1,-1,-1,-1,1,1,1,1,1))

t1 =findOverlappingPeaks(peaks1, peaks2, maxgap=1000, multiple=F, NameOfPeaks1="TF1", NameOfPeaks2="TF2") 



###################################################
### chunk number 4: 
###################################################
overlappingPeaks = t1$OverlappingPeaks
overlappingPeaks
pie(table(overlappingPeaks$overlapFeature))


###################################################
### chunk number 5: 
###################################################
as.data.frame(t1$MergedPeaks)


###################################################
### chunk number 6: 
###################################################
as.data.frame(t1$Peaks1withOverlaps)


###################################################
### chunk number 7: 
###################################################
as.data.frame(t1$Peaks2withOverlaps)


###################################################
### chunk number 8: 
###################################################
peaks3 = RangedData(IRanges(start=c(967859, 2010868, 2496500, 3075966, 3123460 ,3851500, 96865, 201189, 249600, 307386), end=c(967969, 2011908, 2496720,3076166,3123470,3857680, 96985, 201299, 249890, 307796), names=c("p1", "p2", "p3", "p4", "p5", "p6","p7", "p8", "p9", "p10")), space=c("1", "2", "3", "4", "5", "6","1", "2", "3", "4"), strand=
c(1,1,1,1,1,1,-1,-1,-1,-1))

findOverlappingPeaks(findOverlappingPeaks(peaks1, peaks2, maxgap=1000, multiple=F,NameOfPeaks1="TF1", NameOfPeaks2="TF2")$Peaks1withOverlap, peaks3, maxgap=1000, multiple=F, NameOfPeaks1="TF1TF2", NameOfPeaks2="TF3")$Peaks1withOverlap


###################################################
### chunk number 9: 
###################################################
makeVennDiagram(RangedDataList(peaks1, peaks2), NameOfPeaks=c("TF1", "TF2"), maxgap=0, totalTest=100, cex = 1, counts.col = "red")


###################################################
### chunk number 10: 
###################################################
peaks = RangedData(IRanges(start=c(100, 500), end=c(300, 600), names=c("peak1", "peak2")), space=c("NC_008253", "NC_010468"))
library(BSgenome.Ecoli.NCBI.20080805)
peaksWithSequences = getAllPeakSequence(peaks, upstream = 20,
	downstream = 20, genome = Ecoli)


###################################################
### chunk number 11: 
###################################################
write2FASTA(peaksWithSequences)


###################################################
### chunk number 12: 
###################################################
data(enrichedGO)


###################################################
### chunk number 13: 
###################################################
enrichedGO$bp[1:6,]


###################################################
### chunk number 14: 
###################################################
enrichedGO$mf[1:6,]


###################################################
### chunk number 15: 
###################################################
enrichedGO$cc


###################################################
### chunk number 16: 
###################################################
sessionInfo()