inst/doc/iSeq.R

In iSeq: Bayesian Hierarchical Modeling of ChIP-seq Data Through Hidden Ising Models

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

###################################################
### code chunk number 1: iSeq.Rnw:88-91
###################################################
library(iSeq)
data(nrsf)
names(nrsf)


###################################################
### code chunk number 2: iSeq.Rnw:96-100
###################################################
chip = rbind(nrsf$chipFC1592,nrsf$chipFC1862,nrsf$chipFC2002)
mock = rbind(nrsf$mockFC1592,nrsf$mockFC1862,nrsf$mockFC2002)
print(chip[1:2,])
print(mock[1:2,])


###################################################
### code chunk number 3: iSeq.Rnw:107-111
###################################################
chip[,2] = chip[,2]+12
mock[,2] = mock[,2]+12
print(chip[1:2,])
print(mock[1:2,])


###################################################
### code chunk number 4: iSeq.Rnw:128-130
###################################################
tagct = mergetag(chip=chip,control=mock,maxlen=80,minlen=10,ntagcut=10)
print(tagct[1:3,])


###################################################
### code chunk number 5: iSeq.Rnw:144-149
###################################################
tagct22 = tagct[tagct[,1]=="chr22",]  #chr22 data 
reg0 = peakreg(chrpos=tagct22[,1:3],count=(tagct22[,5:6]-tagct22[,7:8]),
               pp=(tagct22[,4]>10),cutoff=0,method="ppcut",maxgap=200)
print(dim(reg0)) 
print(reg0[1:3,]) # each row contain the information for an enriched region


###################################################
### code chunk number 6: iSeq.Rnw:158-160
###################################################
reg0 = reg0[order(reg0$sym,reg0$ct12,decreasing = TRUE),]
print(reg0[1:5,]) # Top five regions


###################################################
### code chunk number 7: iSeq.Rnw:169-172
###################################################
set.seed(777)
res1 = iSeq1(Y=tagct22[,1:4],gap=200,burnin=200,sampling=1000,ctcut=0.95,
     a0=1,b0=1,a1=5,b1=1,k0=3,mink=0,maxk=10,normsd=0.1,verbose=FALSE)


###################################################
### code chunk number 8: iSeq.Rnw:182-187
###################################################
par(mfrow=c(2,2), mar=c(4.1, 4.1, 2.0, 1.0))
hist(res1$pp)
plot(res1$kappa, pch=".", xlab="Iterations", ylab="kappa")
plot(res1$lambda0, pch=".", xlab="Iterations", ylab="lambda0")
plot(res1$lambda1, pch=".", xlab="Iterations", ylab="lambda1")


###################################################
### code chunk number 9: iSeq.Rnw:201-205
###################################################
reg1 = peakreg(chrpos=tagct22[,1:3],count=(tagct22[,5:6]-tagct22[,7:8]),pp=res1$pp,
     cutoff=0.5,method="ppcut",maxgap=200)
print(dim(reg1))
print(reg1[1:3,])


###################################################
### code chunk number 10: iSeq.Rnw:211-216
###################################################
par(mfrow=c(2,2), mar=c(4.1, 4.1, 2.0, 1.0))
for(i in 1:4){
  ID = (reg1[i,4]):(reg1[i,5])
  plotreg(tagct22[ID,2:3],tagct22[ID,5:6],tagct22[ID,7:8],peak=reg1[i,6])
}


###################################################
### code chunk number 11: iSeq.Rnw:223-228
###################################################
reg2 = peakreg(tagct22[,1:3],tagct22[,5:6]-tagct22[,7:8],res1$pp,0.05,
     method="fdrcut",maxgap=200)
print(dim(reg2))
print(reg2[1:3,])



###################################################
### code chunk number 12: iSeq.Rnw:235-236
###################################################
bed = data.frame(chr=reg2[,1],gstart=reg2[,6]-100,gend=reg2[,6]+100)


###################################################
### code chunk number 13: iSeq.Rnw:252-254
###################################################
res2 = iSeq2(Y=tagct22[,1:4],gap=200, burnin=100,sampling=500,winsize=2,
     ctcut=0.95,a0=1,b0=1,a1=5,b1=1,k=1.0,verbose=FALSE)


###################################################
### code chunk number 14: iSeq.Rnw:261-265
###################################################
par(mfrow=c(2,2), mar=c(4.1, 4.1, 2.0, 1.0))
hist(res2$pp)
plot(res2$lambda0, pch=".", xlab="Iterations", ylab="lambda0")
plot(res2$lambda1, pch=".", xlab="Iterations", ylab="lambda1")


###################################################
### code chunk number 15: iSeq.Rnw:270-273
###################################################
reg2 = peakreg(tagct22[,1:3],tagct22[,5:6],res2$pp,0.5,method="ppcut",maxgap=200)
print(dim(reg2))
print(reg2[1:3,])


###################################################
### code chunk number 16: iSeq.Rnw:279-284
###################################################
par(mfrow=c(2,2), mar=c(4.1, 4.1, 2.0, 1.0))
for(i in 1:4){
  ID = (reg2[i,4]):(reg2[i,5])
  plotreg(tagct22[ID,2:3],tagct22[ID,5:6],tagct22[ID,7:8],peak=reg2[i,6])
}


###################################################
### code chunk number 17: iSeq.Rnw:292-298
###################################################
tagctY = tagct[tagct[,1]=="chrY",]
print(table(tagctY[,4]))
res1 = iSeq1(Y=tagctY[,1:4],gap=200,burnin=1000,sampling=5000,ctcut=0.95,
     a0=1,b0=1,a1=5,b1=1,k0=3,mink=0,maxk=10,normsd=0.5,verbose=FALSE)
res2 = iSeq2(Y=tagctY[,1:4],gap=200, burnin=1000,sampling=5000,winsize=2,ctcut=0.95,
  a0=1,b0=1,a1=5,b1=1,k=3.0,verbose=FALSE)


###################################################
### code chunk number 18: iSeq.Rnw:302-307
###################################################
par(mfcol=c(2,2), mar=c(4.1, 4.1, 2.0, 1.0))
hist(res1$pp)
plot(res1$lambda1, pch=".", xlab="Iterations", ylab="lambda1")
hist(res2$pp)
plot(res2$lambda1, pch=".", xlab="Iterations", ylab="lambda1")


###################################################
### code chunk number 19: iSeq.Rnw:360-361
###################################################
sessionInfo()