###################################################
### chunk number 1: 
###################################################
#line 81 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
library(BayesPeak)


###################################################
### chunk number 2: ex1 eval=FALSE
###################################################
## #line 89 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
## raw.output <- bayespeak("H3K4me3-chr16.bed", "Input-chr16.bed",
## 			chr = "chr16", start = 9.2E7, end = 9.5E7, job.size = 6E6)
## output <- summarize.peaks(raw.output, method = "lowerbound")


###################################################
### chunk number 3: ex2 eval=FALSE
###################################################
## #line 107 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
## raw.output <- bayespeak("H3K4me3-chr16.bed", "Input-chr16.bed")
## output <- summarize.peaks(raw.output, method = "lowerbound")


###################################################
### chunk number 4: savetodisk eval=FALSE
###################################################
## #line 117 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
## library(rtracklayer)
## export(output, "H3K4me3output.bed")
## export(output, "H3K4me3output.gff")


###################################################
### chunk number 5: savetodisk2 eval=FALSE
###################################################
## #line 125 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
## write.table(as.data.frame(output), file = "H3K4me3output.txt", quote = FALSE)
## write.csv(as.data.frame(output), file = "H3K4me3output.csv", quote = FALSE)


###################################################
### chunk number 6: Help eval=FALSE
###################################################
## #line 152 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
## raw.output <- bayespeak("H3K4me3-chr16.bed", "Input-chr16.bed")


###################################################
### chunk number 7: Help eval=FALSE
###################################################
## #line 186 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
## library(multicore)
## raw.output <- bayespeak("H3K4me3-chr16.bed", "Input-chr16.bed",
## 	use.multicore = TRUE, mc.cores = 4)
## output <- summarize.peaks(raw.output, method = "lowerbound")


###################################################
### chunk number 8: Help eval=FALSE
###################################################
## #line 201 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
## library(snow)
## cl <- makeCluster(4, type="SOCK") ##adapt this line as appropriate for your cluster
## clusterEvalQ(cl, library(BayesPeak)) ##load BayesPeak on the cluster
## raw.output <- bayespeak("H3K4me3-chr16.bed", "Input-chr16.bed",
## 	snow.cluster = cl)
## output <- summarize.peaks(raw.output, method = "lowerbound")


###################################################
### chunk number 9: PP1
###################################################
#line 241 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
data(raw.output)
raw.output$call


###################################################
### chunk number 10: PP3a
###################################################
#line 276 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
min.job <- min(raw.output$peaks$job)
max.job <- max(raw.output$peaks$job)

par(mfrow = c(2,2), ask = TRUE)
for(i in min.job:max.job) {plot.PP(raw.output, job = i, ylim = c(0,50))}


###################################################
### chunk number 11: job324
###################################################
#line 295 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
i <- 324
plot.PP(raw.output, job = i, ylim = c(0,50))


###################################################
### chunk number 12: job324fig
###################################################
#line 304 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
#line 295 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw#from line#304#"
i <- 324
plot.PP(raw.output, job = i, ylim = c(0,50))
#line 305 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"


###################################################
### chunk number 13: job325
###################################################
#line 322 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
i <- 325
plot.PP(raw.output, job = i, ylim = c(0,50))


###################################################
### chunk number 14: job325fig
###################################################
#line 331 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
#line 322 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw#from line#331#"
i <- 325
plot.PP(raw.output, job = i, ylim = c(0,50))
#line 332 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"


###################################################
### chunk number 15: OF1
###################################################
#line 351 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
data(raw.output)


###################################################
### chunk number 16: OF3
###################################################
#line 386 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"

plot.overfitdiag(raw.output, whatX = "calls", whatY = "lambda1",
	logX = TRUE, logY = TRUE)


###################################################
### chunk number 17: OF3b
###################################################
#line 391 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
plot.overfitdiag(raw.output, whatX = "calls", whatY = "score",
	logX = TRUE, logY = TRUE)


###################################################
### chunk number 18: OF3fig
###################################################
#line 401 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
#line 386 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw#from line#401#"

plot.overfitdiag(raw.output, whatX = "calls", whatY = "lambda1",
	logX = TRUE, logY = TRUE)
#line 402 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"


###################################################
### chunk number 19: OF3bfig
###################################################
#line 411 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
#line 391 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw#from line#411#"
plot.overfitdiag(raw.output, whatX = "calls", whatY = "score",
	logX = TRUE, logY = TRUE)
#line 412 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"


###################################################
### chunk number 20: OF3
###################################################
#line 423 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
plot.overfitdiag(raw.output, whatX = "calls", whatY = "score",
	logX = TRUE, logY = FALSE)


###################################################
### chunk number 21: OF3b
###################################################
#line 430 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
plot.overfitdiag(raw.output, whatX = "lambda1", whatY = "score",
	logX = TRUE, logY = TRUE)


###################################################
### chunk number 22: OF4
###################################################
#line 444 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
unreliable.jobs <- log(raw.output$QC$lambda1) < 1.5
output <- summarize.peaks(raw.output, method = "lowerbound",
	exclude.jobs = unreliable.jobs)


###################################################
### chunk number 23: OF5
###################################################
#line 452 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
unreliable.jobs2 <- log(raw.output$QC$lambda1) < 1.5 | log(raw.output$QC$calls) > 5
output.2 <- summarize.peaks(raw.output, method = "lowerbound",
	exclude.jobs = unreliable.jobs2)


###################################################
### chunk number 24: region.overfitdiag.help
###################################################
#line 460 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
unreliable.jobs3 <- region.overfitdiag(raw.output, whatX = "lambda1", whatY = "score",
	logX = TRUE, logY = TRUE)
##user defines a polygon on the resulting plot
##left-click to place each vertex, right-click to close polygon 
output.3 <- summarize.peaks(raw.output, method = "lowerbound",
	exclude.jobs = unreliable.jobs3)


###################################################
### chunk number 25: QCcommand
###################################################
#line 483 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
log(raw.output$QC[,c("calls", "score", "lambda1")])


###################################################
### chunk number 26: coda eval=FALSE
###################################################
## #line 523 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
## data(raw.output)
## library(coda)
## mcmc.job1 <- mcmc(raw.output$p.samples[[316]], thin = 10)
## geweke.diag(mcmc.job1)


###################################################
### chunk number 27: 
###################################################
#line 570 "vignettes/BayesPeak/inst/doc/BayesPeak.Rnw"
sessionInfo()