###################################################
### chunk number 1: 
###################################################
library(iSeq)
data(nrsf)
names(nrsf)


###################################################
### chunk number 2: 
###################################################
chip = rbind(nrsf$chipFC1592,nrsf$chipFC1862,nrsf$chipFC2002)
mock = rbind(nrsf$mockFC1592,nrsf$mockFC1862,nrsf$mockFC2002)
print(chip[1:3,])
print(mock[1:3,])


###################################################
### chunk number 3: 
###################################################
tagct = mergetag(chip=chip,control=mock,winsize=50)
print(tagct[1:3,])


###################################################
### chunk number 4: 
###################################################
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=300)
print(dim(reg0)) 
print(reg0[1:3,]) # each row contain the information for an enriched region


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


###################################################
### chunk number 6: 
###################################################
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")


###################################################
### chunk number 7: 
###################################################
reg1 = peakreg(chrpos=tagct22[,1:3],count=(tagct22[,5:6]-tagct22[,7:8]),pp=res1$pp,cutoff=0.5,method="ppcut",maxgap=300)
print(dim(reg1))
print(reg1[1:3,])


###################################################
### chunk number 8: 
###################################################
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])
}


###################################################
### chunk number 9: 
###################################################
reg2 = peakreg(tagct22[,1:3],tagct22[,5:6]-tagct22[,7:8],res1$pp,0.05,method="fdrcut",maxgap=300)
print(dim(reg2))
print(reg2[1:3,])



###################################################
### chunk number 10: 
###################################################
bed = data.frame(chr=reg2[,1],gstart=reg2[,6]-100,gend=reg2[,6]+100)


###################################################
### chunk number 11: 
###################################################
res2 = iSeq2(Y=tagct22[,1:4],gap=300, burnin=100,sampling=500,winsize=2,ctcut=5,
  a0=1,b0=1,a1=5,b1=1,k=1.0,verbose=FALSE)


###################################################
### chunk number 12: 
###################################################
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")


###################################################
### chunk number 13: 
###################################################
reg2 = peakreg(tagct22[,1:3],tagct22[,5:6],res2$pp,0.5,method="ppcut",maxgap=300)
print(dim(reg2))
print(reg2[1:3,])


###################################################
### chunk number 14: 
###################################################
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])
}


###################################################
### chunk number 15: 
###################################################
tagctY = tagct[tagct[,1]=="chrY",]
print(table(tagctY[,4]))
res1 = iSeq1(Y=tagctY[,1:4],gap=300,burnin=1000,sampling=5000,ctcut=3,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=300, burnin=1000,sampling=5000,winsize=2,ctcut=3,
  a0=1,b0=1,a1=5,b1=1,k=3.0,verbose=FALSE)


###################################################
### chunk number 16: 
###################################################
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")


###################################################
### chunk number 17: 
###################################################
sessionInfo()