Introduction.R
In ChIPseqR: Identifying Protein Binding Sites in High-Throughput Sequencing Data

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

###################################################
### code chunk number 1: distances
###################################################
	set.seed(1)
	dist <- sample(c(170:300, 400), 10100, prob=c(dnbinom(0:130, mu=30, size=5), 0.2), replace=TRUE)
	pos <- cumsum(dist)
	pos <- pos[pos < 2e6 - 200]


###################################################
### code chunk number 2: readPos1
###################################################
	fwdRegion <- unlist(lapply(pos, function(x) (x-88):(x-68)))
	revRegion <- unlist(lapply(pos, function(x) (x+68):(x+88)))


###################################################
### code chunk number 3: readPos2
###################################################
	fwd <- sample(fwdRegion, 5e4, replace=TRUE)
	rev <- sample(revRegion, 5e4, replace=TRUE)


###################################################
### code chunk number 4: readPos3
###################################################
	fwd <- c(fwd, sample(25:(2e6-25), 2e5, replace=TRUE))
	rev <- c(rev, sample(25:(2e6-25), 2e5, replace=TRUE))


###################################################
### code chunk number 5: formatData
###################################################
	reads <- data.frame(chromosome="chr1", position=c(fwd, rev), length=25, 
			strand=rep(c("+", "-"), each=250000))


###################################################
### code chunk number 6: library
###################################################
	library(ChIPseqR)


###################################################
### code chunk number 7: simple
###################################################
	counts <- strandPileup(reads, chrLen=2e6, extend=1, plot=FALSE)
	nucs <- simpleNucCall(counts, bind=147, support=15, plot=FALSE)


###################################################
### code chunk number 8: plotNucPrint (eval = FALSE)
###################################################
## 	plot(nucs, type="density")
## 	plot(nucs, type="qqplot")


###################################################
### code chunk number 9: plotNuc
###################################################
	par(mfrow=c(1,2))
	plot(nucs, type="density")
	plot(nucs, type="qqplot")


###################################################
### code chunk number 10: plotWindowPrint (eval = FALSE)
###################################################
## 	predicted <- peaks(nucs)[[1]][911]
## 	plot(counts, chr="chr1", center=predicted, score=nucs, width=1000, type="window")


###################################################
### code chunk number 11: plotWindowPrint2 (eval = FALSE)
###################################################
## 	abline(v=pos[pos < predicted + 1000 & pos > predicted - 1000], col=3, lty=3)


###################################################
### code chunk number 12: plotWindow
###################################################
	predicted <- peaks(nucs)[[1]][911]
	plot(counts, chr="chr1", center=predicted, score=nucs, width=1000, type="window")
	abline(v=pos[pos < predicted + 1000 & pos > predicted - 1000], col=3, lty=3)


###################################################
### code chunk number 13: noOver
###################################################
	calls <- peaks(nucs)[[1]][c(1,which(diff(peaks(nucs)[[1]]) >= 170)+1)]
	length(calls)


###################################################
### code chunk number 14: spec
###################################################
	table(sapply(calls, function(x) any((x-20):(x+20) %in% pos)))


###################################################
### code chunk number 15: getBindLen
###################################################
	bLen <- getBindLen(counts, bind = c(100,200), support = c(5, 50))


###################################################
### code chunk number 16: getBindLenPrint (eval = FALSE)
###################################################
## 	bLen <- getBindLen(counts, bind = c(100,200), support = c(5, 50))


###################################################
### code chunk number 17: decompress
###################################################
	counts2 <- decompress(counts)


###################################################
### code chunk number 18: sessionInfo
###################################################
	sessionInfo()