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

###################################################
### code chunk number 1: iChip.Rnw:90-93
###################################################
library(iChip)
data(oct4)
head(oct4,n=3L)


###################################################
### code chunk number 2: iChip.Rnw:100-101
###################################################
oct4 = oct4[order(oct4[,1],oct4[,2]),]


###################################################
### code chunk number 3: iChip.Rnw:105-106
###################################################
oct4lmt = lmtstat(oct4[,5:6],oct4[,3:4])


###################################################
### code chunk number 4: iChip.Rnw:120-121
###################################################
oct4Y = cbind(oct4[,1],oct4lmt)


###################################################
### code chunk number 5: iChip.Rnw:131-133
###################################################
set.seed(777)
oct4res2 = iChip2(Y=oct4Y,burnin=2000,sampling=10000,winsize=2,sdcut=2,beta=1.25,verbose=FALSE)


###################################################
### code chunk number 6: iChip.Rnw:142-147
###################################################
par(mfrow=c(2,2), mar=c(4.1, 4.1, 2.0, 1.0))
plot(oct4res2$mu0, pch=".", xlab="Iterations", ylab="mu0")
plot(oct4res2$lambda0, pch=".", xlab="Iterations", ylab="lambda0")
plot(oct4res2$mu1, pch=".", xlab="Iterations", ylab="mu1")
plot(oct4res2$lambda1, pch=".", xlab="Iterations", ylab="lambda1")


###################################################
### code chunk number 7: iChip.Rnw:155-156
###################################################
hist(oct4res2$pp)


###################################################
### code chunk number 8: iChip.Rnw:164-166
###################################################
reg1 = enrichreg(pos=oct4[,1:2],enrich=oct4lmt,pp=oct4res2$pp,cutoff=0.9,method="ppcut",maxgap=500)
print(reg1)


###################################################
### code chunk number 9: iChip.Rnw:172-174
###################################################
reg2 = enrichreg(pos=oct4[,1:2],enrich=oct4lmt,pp=oct4res2$pp,cutoff=0.01,method="fdrcut",maxgap=500)
print(reg2)


###################################################
### code chunk number 10: iChip.Rnw:181-183
###################################################
bed1 = data.frame(chr=paste("chr",reg2[,1],sep=""),reg2[,2:3])
print(bed1[1:3,])


###################################################
### code chunk number 11: iChip.Rnw:189-191
###################################################
bed2 = data.frame(chr=paste("chr",reg2[,1],sep=""),gstart=reg2[,6]-100,gend=reg2[,6]+100)
print(bed2[1:3,])


###################################################
### code chunk number 12: iChip.Rnw:196-197
###################################################
oct4res1 =iChip1(enrich=oct4lmt,burnin=2000,sampling=10000,sdcut=2,beta0=3,minbeta=0,maxbeta=10,normsd=0.1,verbose=FALSE)


###################################################
### code chunk number 13: iChip.Rnw:203-208
###################################################
par(mfrow=c(2,2), mar=c(4.1, 4.1, 2.0, 1.0))
plot(oct4res1$beta, pch=".", xlab="Iterations", ylab="beta")
plot(oct4res1$mu0, pch=".", xlab="Iterations", ylab="mu0")
plot(oct4res1$mu1, pch=".", xlab="Iterations", ylab="mu1")
plot(oct4res1$lambda, pch=".", xlab="Iterations", ylab="lambda")


###################################################
### code chunk number 14: iChip.Rnw:214-216
###################################################
enrichreg(pos=oct4[,1:2],enrich=oct4lmt,pp=oct4res1$pp,cutoff=0.9,method="ppcut",maxgap=500)
enrichreg(pos=oct4[,1:2],enrich=oct4lmt,pp=oct4res1$pp,cutoff=0.01,method="fdrcut",maxgap=500)


###################################################
### code chunk number 15: iChip.Rnw:223-229
###################################################
data(p53)
head(p53,n=3L)
# sort the p53 data by chromosome and genomic position
p53 = p53[order(p53[,1],p53[,2]),]
p53lmt = lmtstat(p53[,9:14],p53[,3:8])
p53Y = cbind(p53[,1],p53lmt)


###################################################
### code chunk number 16: iChip.Rnw:236-237
###################################################
p53res2 = iChip2(Y=p53Y,burnin=2000,sampling=10000,winsize=2,sdcut=2,beta=2.5)


###################################################
### code chunk number 17: iChip.Rnw:241-246
###################################################
par(mfrow=c(2,2), mar=c(4.1, 4.1, 2.0, 1.0))
plot(p53res2$mu0, pch=".", xlab="Iterations", ylab="mu0")
plot(p53res2$lambda0, pch=".", xlab="Iterations", ylab="lambda0")
plot(p53res2$mu1, pch=".", xlab="Iterations", ylab="mu1")
plot(p53res2$lambda1, pch=".", xlab="Iterations", ylab="lambda1")


###################################################
### code chunk number 18: iChip.Rnw:254-255
###################################################
hist(p53res2$pp)


###################################################
### code chunk number 19: iChip.Rnw:259-261
###################################################
enrichreg(pos=p53[,1:2],enrich=p53lmt,pp=p53res2$pp,cutoff=0.9,method="ppcut",maxgap=500)
enrichreg(pos=p53[,1:2],enrich=p53lmt,pp=p53res2$pp,cutoff=0.01,method="fdrcut",maxgap=500)


###################################################
### code chunk number 20: iChip.Rnw:266-267
###################################################
p53res1 =iChip1(enrich=p53lmt,burnin=2000,sampling=10000,sdcut=2,beta0=3,minbeta=0,maxbeta=10,normsd=0.1,verbose=FALSE)


###################################################
### code chunk number 21: iChip.Rnw:271-276
###################################################
par(mfrow=c(2,2), mar=c(4.1, 4.1, 2.0, 1.0))
plot(p53res1$beta, pch=".", xlab="Iterations", ylab="beta")
plot(p53res1$mu0, pch=".", xlab="Iterations", ylab="mu0")
plot(p53res1$mu1, pch=".", xlab="Iterations", ylab="mu1")
plot(p53res1$lambda, pch=".", xlab="Iterations", ylab="lambda")


###################################################
### code chunk number 22: iChip.Rnw:280-282
###################################################
enrichreg(pos=p53[,1:2],enrich=p53lmt,pp=p53res1$pp,cutoff=0.9,method="ppcut",maxgap=500)
enrichreg(pos=p53[,1:2],enrich=p53lmt,pp=p53res1$pp,cutoff=0.01,method="fdrcut",maxgap=500)


###################################################
### code chunk number 23: iChip.Rnw:288-291
###################################################
randomY = cbind(p53[,1],rnorm(10000,0,1))
randomres2 = iChip2(Y=randomY,burnin=1000,sampling=5000,winsize=2,sdcut=2,beta=2.5,verbose=FALSE)
table(randomres2$pp)


###################################################
### code chunk number 24: iChip.Rnw:312-314
###################################################
randomres1 =iChip1(enrich=randomY[,2],burnin=1000,sampling=5000,sdcut=2,beta0=3,minbeta=0,maxbeta=10,normsd=0.1,verbose=FALSE)
table(randomres1$pp)


###################################################
### code chunk number 25: iChip.Rnw:353-354
###################################################
sessionInfo()