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

###################################################
### code chunk number 1: loadTA
###################################################
require(IRanges)
library(nucleR)
data(nucleosome_tiling)
head(nucleosome_tiling, n=25)


###################################################
### code chunk number 2: import
###################################################
data(nucleosome_htseq)
class(nucleosome_htseq)
nucleosome_htseq


###################################################
### code chunk number 3: processReads
###################################################
# Process the paired end reads, but discard those with length > 200
reads_pair <- processReads(nucleosome_htseq, type="paired", fragmentLen=200)

# Process the reads, but now trim each read to 40bp around the dyad
reads_trim <- processReads(nucleosome_htseq, type="paired",
                           fragmentLen=200, trim=40)


###################################################
### code chunk number 4: coverage
###################################################
# Calculate the coverage, directly in reads per million (r.p.m)
cover_pair <- coverage.rpm(reads_pair)
cover_trim <- coverage.rpm(reads_trim)


###################################################
### code chunk number 5: figcover
###################################################
# Compare both coverages
t1 <- as.vector(cover_pair[[1]])[1:2000]
t2 <- as.vector(cover_trim[[1]])[1:2000]
t1 <- (t1 - min(t1)) / max(t1 - min(t1)) # Normalization
t2 <- (t2 - min(t2)) / max(t2 - min(t2)) # Normalization
par(mar=c(5, 4, 1, 1), cex=0.7)
plot(t1, type="l", lwd="1", col="blue", ylab="norm. coverage", xlab="position")
lines(t2, lwd="1", col="orange")


###################################################
### code chunk number 6: figcover
###################################################
# Compare both coverages
t1 <- as.vector(cover_pair[[1]])[1:2000]
t2 <- as.vector(cover_trim[[1]])[1:2000]
t1 <- (t1 - min(t1)) / max(t1 - min(t1)) # Normalization
t2 <- (t2 - min(t2)) / max(t2 - min(t2)) # Normalization
par(mar=c(5, 4, 1, 1), cex=0.7)
plot(t1, type="l", lwd="1", col="blue", ylab="norm. coverage", xlab="position")
lines(t2, lwd="1", col="orange")


###################################################
### code chunk number 7: mnase
###################################################
# Toy example
map <- syntheticNucMap(as.ratio=TRUE, wp.num=50, fuz.num=25)

exp <- coverage(map$syn.reads)
ctr <- coverage(map$ctr.reads)

corrected <- controlCorrection(exp, ctr)


###################################################
### code chunk number 8: figmnase
###################################################
# Toy example
par(mar=c(5,4,1,1), cex=0.7)
plot(exp, col="darkgrey", type="l", ylab="coverage", xlab="position")
lines(corrected, col="darkblue", lty=2)
legend("top", c("normal", "corrected"), hor=TRUE, bty="n",
       fill=c("darkgrey", "blue"))


###################################################
### code chunk number 9: fignoise
###################################################
par(mar=c(5,4,1,1), mfcol=c(1,4), cex=0.4)
plot(nucleosome_tiling[1:2000], type="l", main="original", ylim=c(-3,2),
     xlab="position", ylab="intensity")
plot(filter(nucleosome_tiling[1:2000], c(rep(1,20))/20), type="l",
     main="sliding w. 20bp", ylim=c(-3,2), xlab="position", ylab="")
plot(filter(nucleosome_tiling[1:2000], c(rep(1,50))/50), type="l",
     main="sliding w. 50bp", ylim=c(-3,2), xlab="position", ylab="")
plot(filter(nucleosome_tiling[1:2000], c(rep(1,100))/100), type="l",
     main="sliding w. 100bp", ylim=c(-3,2), xlab="position", ylab="")


###################################################
### code chunk number 10: fignoise
###################################################
par(mar=c(5,4,1,1), mfcol=c(1,4), cex=0.4)
plot(nucleosome_tiling[1:2000], type="l", main="original", ylim=c(-3,2),
     xlab="position", ylab="intensity")
plot(filter(nucleosome_tiling[1:2000], c(rep(1,20))/20), type="l",
     main="sliding w. 20bp", ylim=c(-3,2), xlab="position", ylab="")
plot(filter(nucleosome_tiling[1:2000], c(rep(1,50))/50), type="l",
     main="sliding w. 50bp", ylim=c(-3,2), xlab="position", ylab="")
plot(filter(nucleosome_tiling[1:2000], c(rep(1,100))/100), type="l",
     main="sliding w. 100bp", ylim=c(-3,2), xlab="position", ylab="")


###################################################
### code chunk number 11: fft
###################################################
fft_ta = filterFFT(nucleosome_tiling, pcKeepComp=0.01, showPowerSpec=TRUE)


###################################################
### code chunk number 12: figfft
###################################################
par(mar=c(6,4,1,1), cex=0.7)
fft_ta <- filterFFT(nucleosome_tiling, pcKeepComp=0.01, showPowerSpec=TRUE)


###################################################
### code chunk number 13: figfft
###################################################
par(mar=c(6,4,1,1), cex=0.7)
fft_ta <- filterFFT(nucleosome_tiling, pcKeepComp=0.01, showPowerSpec=TRUE)


###################################################
### code chunk number 14: figfft3
###################################################
par(mar=c(5, 4, 1, 1), mfcol=c(2, 1), cex=0.7)
tiling_raw <- nucleosome_tiling[1:3000]
tiling_fft <- filterFFT(tiling_raw, pcKeepComp=0.01)
htseq_raw <- as.vector(cover_trim[[1]])[1:3000]
htseq_fft <- filterFFT(htseq_raw, pcKeepComp=0.02)
plot(tiling_raw, type="l", col="darkblue", lwd=3, ylab="intensity", xlab="")
lines(tiling_fft, type="l", col="cyan")
plot(htseq_raw, type="l", col="darkred", lwd=3, ylab="coverage",
     xlab="position")
lines(htseq_fft, type="l", col="pink")


###################################################
### code chunk number 15: figfft3
###################################################
par(mar=c(5, 4, 1, 1), mfcol=c(2, 1), cex=0.7)
tiling_raw <- nucleosome_tiling[1:3000]
tiling_fft <- filterFFT(tiling_raw, pcKeepComp=0.01)
htseq_raw <- as.vector(cover_trim[[1]])[1:3000]
htseq_fft <- filterFFT(htseq_raw, pcKeepComp=0.02)
plot(tiling_raw, type="l", col="darkblue", lwd=3, ylab="intensity", xlab="")
lines(tiling_fft, type="l", col="cyan")
plot(htseq_raw, type="l", col="darkred", lwd=3, ylab="coverage",
     xlab="position")
lines(htseq_fft, type="l", col="pink")


###################################################
### code chunk number 16: corfft
###################################################
tiling_raw <- nucleosome_tiling
tiling_fft <- filterFFT(tiling_raw, pcKeepComp=0.01)
htseq_raw <- as.vector(cover_trim[[1]])
htseq_fft <- filterFFT(htseq_raw, pcKeepComp=0.02)

cor(tiling_raw, tiling_fft, use="complete.obs")
cor(htseq_raw, htseq_fft, use="complete.obs")


###################################################
### code chunk number 17: peaks
###################################################
peaks <- peakDetection(htseq_fft, threshold="25%", score=FALSE)
peaks


###################################################
### code chunk number 18: plotpeak
###################################################
plotPeaks(peaks, htseq_fft, threshold="25%")


###################################################
### code chunk number 19: figpeak
###################################################
par(cex=0.7, mar=c(5,4,1,1))
plotPeaks(peaks, htseq_fft, threshold="25%", yaxt="n", ylab="coverage")


###################################################
### code chunk number 20: peaks2
###################################################
peaks <- peakDetection(htseq_fft, threshold="25%", score=TRUE)
head(peaks)


###################################################
### code chunk number 21: figpeak2
###################################################
par(cex=0.7, mar=c(5,4,1,1))
plotPeaks(peaks, htseq_fft, threshold="25%", yaxt="n", xlim=c(1,4000),
          ylab="coverage")


###################################################
### code chunk number 22: peaks3
###################################################
peaks <- peakDetection(htseq_fft, threshold="25%", score=TRUE, width=140)
peaks


###################################################
### code chunk number 23: figpeak3
###################################################
par(cex=0.7, mar=c(5,4,1,1))
plotPeaks(peaks, htseq_fft, threshold="25%", yaxt="n", xlim=c(1,4000),
          rect.lwd=0.5, rect.thick=2.5, ylab="coverage")


###################################################
### code chunk number 24: ranges
###################################################
nuc_calls <- ranges(peaks[peaks$score > 0.1,])[[1]]
red_calls <- reduce(nuc_calls)
red_class <- RangedData(red_calls, isFuzzy=width(red_calls) > 140)
red_class


###################################################
### code chunk number 25: figranges
###################################################
par(cex=0.7, mar=c(5,4,1,1))
plotPeaks(red_calls, htseq_fft, threshold="25%", yaxt="n", xlim=c(1,4000), rect.lwd=0.5, rect.thick=3, ylab="coverage")


###################################################
### code chunk number 26: syn
###################################################
syntheticNucMap(wp.num=100, wp.del=10, wp.var=30, fuz.num=20, fuz.var=50,
                max.cover=20, nuc.len=147, lin.len=20, rnd.seed=1,
                as.ratio=TRUE, show.plot=TRUE)


###################################################
### code chunk number 27: figsyn
###################################################
par(mar=c(5,4,1,1), cex=0.7)
syn <- syntheticNucMap(wp.num=100, wp.del=10, wp.var=30, fuz.num=20,
                       fuz.var=50, max.cover=20, nuc.len=147, lin.len=20,
                       rnd.seed=1, as.ratio=TRUE, show.plot=TRUE, ylab="",
                       xlab="position")


###################################################
### code chunk number 28: figsyn
###################################################
par(mar=c(5,4,1,1), cex=0.7)
syn <- syntheticNucMap(wp.num=100, wp.del=10, wp.var=30, fuz.num=20,
                       fuz.var=50, max.cover=20, nuc.len=147, lin.len=20,
                       rnd.seed=1, as.ratio=TRUE, show.plot=TRUE, ylab="",
                       xlab="position")