###################################################
### chunk number 1: Reading library
###################################################
library(BAC)


###################################################
### chunk number 2: Load the ER data
###################################################
data(ER)


###################################################
### chunk number 3: Calculate joint posterior probabilities
###################################################
# Please uncomment the following line to test the code
# BAConER<-BAC(ER[,5:7], ER[,2:4], B=50,verbose=FALSE,w=5)
# This load the resulting data obtained after executing the line above
load("bac.rda")


###################################################
### chunk number 4: PProbabilities
###################################################
plot(ER[,1],BAConER$jointPP,pch="+",xlab="Genomic Position",ylab="Posterior probabilities")


###################################################
### chunk number 5: Call Regions
###################################################
ERregions<-CallRegions(ER[,1],BAConER$jointPP,cutoff=0.5,maxGap=500)


###################################################
### chunk number 6: Create BED file
###################################################
# Create the BED file
nRegions<-max(ERregions)
BED<-matrix(0,nRegions,4)
for(i in 1:nRegions)
{
BED[i,2:3]<-range(ER[ERregions==i,1])
#The score should be between 0 and 1000
BED[i,4]<-max(BAConER$jointPP[ERregions==i])*1000
}
BED<-data.frame(BED)
# The ER data is a subset of chr 21
BED[,1]<-"chr21"
names(BED)<-c("chrom","chromStart","chromEnd","Score")
# print it
print(BED)