iChip.R
In iChip: Bayesian Modeling of ChIP-chip Data Through Hidden Ising Models

### R code from vignette source 'iChip.Rnw'

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


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


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


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


###################################################
### code chunk number 5: iChip.Rnw:128-131
###################################################
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:140-145
###################################################
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:153-154
###################################################
hist(oct4res2$pp)


###################################################
### code chunk number 8: iChip.Rnw:162-165
###################################################
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:171-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:2,])


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


###################################################
### code chunk number 12: iChip.Rnw:197-199
###################################################
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:205-210
###################################################
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:216-220
###################################################
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:227-233
###################################################
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:240-241
###################################################
p53res2 = iChip2(Y=p53Y,burnin=2000,sampling=10000,winsize=2,sdcut=2,beta=2.5)


###################################################
### code chunk number 17: iChip.Rnw:245-250
###################################################
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:258-259
###################################################
hist(p53res2$pp)


###################################################
### code chunk number 19: iChip.Rnw:263-267
###################################################
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:272-274
###################################################
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:278-283
###################################################
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:287-291
###################################################
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:297-301
###################################################
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:322-325
###################################################
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:373-374
###################################################
sessionInfo()