###################################################
### chunk number 1: 
###################################################
library(BayesPeak)


###################################################
### chunk number 2: ex1 eval=FALSE
###################################################
## raw.output <- bayespeak("H3K4me3-chr16.bed", "Input-chr16.bed",
## 			chr = "chr16", start = 9.2E7, end = 9.5E7, job.size = 6E6)
## output <- summarise.peaks(raw.output, method = "lowerbound")


###################################################
### chunk number 3: ex2 eval=FALSE
###################################################
## raw.output <- bayespeak("H3K4me3-chr16.bed", "Input-chr16.bed")
## output <- summarise.peaks(raw.output, method = "lowerbound")


###################################################
### chunk number 4: Help eval=FALSE
###################################################
## raw.output <- bayespeak("H3K4me3-chr16.bed", "Input-chr16.bed")


###################################################
### chunk number 5: Help eval=FALSE
###################################################
## library(multicore)
## raw.output <- bayespeak("H3K4me3-chr16.bed", "Input-chr16.bed", use.multicore = TRUE, mc.cores = 4)
## output <- combine.peaks(raw.output, method = "lowerbound")


###################################################
### chunk number 6: PP1
###################################################
data(raw.output)
raw.output$call


###################################################
### chunk number 7: PP3
###################################################
PP <- split(raw.output$peaks$PP, raw.output$peaks$job)


###################################################
### chunk number 8: PP3a
###################################################
bin.width <- raw.output$peaks$end[1] - raw.output$peaks$start[1]
job.bins <- (raw.output$QC$end - raw.output$QC$start)/bin.width
job.bins <- job.bins[as.integer(names(PP))]

for(i in 1:length(PP))
{
	PP[[i]] <- c(PP[[i]], rep(0, job.bins[i] - length(PP[[i]])))
}


###################################################
### chunk number 9: PP3a
###################################################
par(mfrow = c(2,2), ask = TRUE)
for(i in 1:length(PP)) hist(PP[[i]], breaks =150, main = names(PP)[i])


###################################################
### chunk number 10: job324
###################################################
i <- 9
hist(PP[[i]], breaks =150, main = names(PP)[i], ylim= c(0,50))


###################################################
### chunk number 11: job324fig
###################################################
i <- 9
hist(PP[[i]], breaks =150, main = names(PP)[i], ylim= c(0,50))


###################################################
### chunk number 12: job325
###################################################
i <- 10
hist(PP[[i]], breaks =150, main = names(PP)[i], ylim= c(0,50))


###################################################
### chunk number 13: job325fig
###################################################
i <- 10
hist(PP[[i]], breaks =150, main = names(PP)[i], ylim= c(0,50))


###################################################
### chunk number 14: OF1
###################################################
data(raw.output)


###################################################
### chunk number 15: OF3
###################################################
plot(log(raw.output$QC$calls), log(raw.output$QC$lambda1),
	main = "Job parameters - enriched bin counts against calls")


###################################################
### chunk number 16: OF3b
###################################################
plot(log(raw.output$QC$calls), log(raw.output$QC$score),
	main = "Job parameters - score against calls")


###################################################
### chunk number 17: OF3fig
###################################################
plot(log(raw.output$QC$calls), log(raw.output$QC$lambda1),
	main = "Job parameters - enriched bin counts against calls")


###################################################
### chunk number 18: OF3bfig
###################################################
plot(log(raw.output$QC$calls), log(raw.output$QC$score),
	main = "Job parameters - score against calls")


###################################################
### chunk number 19: OF4
###################################################
unreliable.jobs <- log(raw.output$QC$lambda1) < 1.5
output.sj <- summarise.peaks(raw.output, method = "lowerbound", exclude.jobs = unreliable.jobs)


###################################################
### chunk number 20: OF5
###################################################
unreliable.jobs2 <- log(raw.output$QC$lambda1) < 1.5 | log(raw.output$QC$calls) > 5
output.sj2 <- summarise.peaks(raw.output, method = "lowerbound", exclude.jobs = unreliable.jobs2)


###################################################
### chunk number 21: 
###################################################
sessionInfo()