###################################################
### chunk number 1: Loading PICS
###################################################
library(PICS)


###################################################
### chunk number 2: Reading the data
###################################################
path<-system.file("extdata",package="PICS")
#Read the experiment : 
dataIP<-read.table(file.path(path,"Treatment_tags_chr21_sort.bed"),header=TRUE,colClasses=c("factor","integer","integer","factor"))
dataIP<-as(dataIP,"RangedData")
dataIP<-as(dataIP,"GenomeData")

dataCont<-read.table(file.path(path,"Input_tags_chr21_sort.bed"),header=TRUE,colClasses=c("factor","integer","integer","factor"))
dataCont<-as(dataCont,"RangedData")
dataCont<-as(dataCont,"GenomeData")


###################################################
### chunk number 3: Reading the mappability profile
###################################################
map<-read.table(file.path(path,"mapProfileShort"),header=TRUE,colClasses=c("factor","integer","integer","NULL"))
map<-as(map,"RangedData")
## Remove the chrM
map<-map[-23]


###################################################
### chunk number 4: Cluster initialization
###################################################
library(snowfall)
sfInit(parallel=TRUE,cpus=2)
sfLibrary(PICS)


###################################################
### chunk number 5: Genome segmentation
###################################################
seg<-segmentReads(dataIP, dataC=dataCont, map=map, minReads=1)


###################################################
### chunk number 6: PICS analysis
###################################################
pics<-PICS(seg,dataType="TF")


###################################################
### chunk number 7: FDR estimation
###################################################
segC<-segmentReads(dataCont,dataC=dataIP,map=map,minReads=1)
picsC<-PICS(segC,dataType="TF")
fdr<-picsFDR(pics,picsC,filter=list(delta=c(50,Inf),se=c(0,50),sigmaSqF=c(0,22500),sigmaSqR=c(0,22500)))


###################################################
### chunk number 8: plot-FDR1
###################################################
plot(pics,picsC,xlim=c(2,8),ylim=c(0,.2),filter=list(delta=c(50,300),se=c(0,50),sigmaSqF=c(0,22500),sigmaSqR=c(0,22500)),type="l")


###################################################
### chunk number 9: plot-FDR2
###################################################
plot(fdr[,c(3,1)])


###################################################
### chunk number 10: Create RangedData data object of enriched regions
###################################################
myFilter=list(delta=c(50,300),se=c(0,50),sigmaSqF=c(0,22500),sigmaSqR=c(0,22500))
rdBed<-makeRangedDataOutput(pics,type="bed",filter=c(myFilter,list(score=c(1,Inf))))


###################################################
### chunk number 11: Create the bed file eval=FALSE
###################################################
## library(rtracklayer)
## export(rdBed,"myfile.bed")


###################################################
### chunk number 12: Create the wig file eval=FALSE
###################################################
## rdBed<-makeRangedDataOutput(pics,type="wig",filter=c(myFilter,list(score=c(1,Inf))))
## export(rdBed,"myfile.wig")