### R code from vignette source 'vignettes/PICS/inst/doc/PICS.Rnw'

###################################################
### code chunk number 1: Loading PICS
###################################################
library(PICS)
citation(package="PICS")


###################################################
### code chunk number 2: Loading PICS
###################################################
library(PICS)


###################################################
### code chunk number 3: 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")


###################################################
### code chunk number 4: 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]


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


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


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


###################################################
### code chunk number 8: 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)))


###################################################
### code chunk number 9: 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")


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


###################################################
### code chunk number 11: 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))))


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


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


###################################################
### code chunk number 14: PICS.Rnw:179-180
###################################################
sfStop()