### R code from vignette source 'vignettes/htSeqTools/inst/doc/htSeqTools.Rnw'

###################################################
### code chunk number 1: import1
###################################################
options(width=70)
library(htSeqTools)
data(htSample)
htSample
sapply(htSample,nrow)


###################################################
### code chunk number 2: getmds
###################################################
cmds1 <- cmds(htSample,k=2)
cmds1


###################################################
### code chunk number 3: figmds
###################################################
plot(cmds1)


###################################################
### code chunk number 4: figmds
###################################################
plot(cmds1)


###################################################
### code chunk number 5: qcontrol1
###################################################
ssdCoverage(htSample)
giniCoverage(htSample,mc.cores=1)


###################################################
### code chunk number 6: figgini1
###################################################
giniCoverage(htSample[['ctrlBatch2']],mc.cores=1,mk.plot=TRUE)


###################################################
### code chunk number 7: figgini2
###################################################
giniCoverage(htSample[['ipBatch2']],mc.cores=1,mk.plot=TRUE)


###################################################
### code chunk number 8: figgini1
###################################################
giniCoverage(htSample[['ctrlBatch2']],mc.cores=1,mk.plot=TRUE)


###################################################
### code chunk number 9: figgini2
###################################################
giniCoverage(htSample[['ipBatch2']],mc.cores=1,mk.plot=TRUE)


###################################################
### code chunk number 10: qcontrol2
###################################################
contamSample <- RangedData(IRanges(rep(1,200),rep(36,200)),space=rep('chr2L',200),strand='+')
contamSample <- rbind(htSample[['ctrlBatch1']],contamSample)
nrepeats <- tabDuplReads(contamSample)
nrepeats


###################################################
### code chunk number 11: qcontrol3
###################################################
q <- which(cumsum(nrepeats/sum(nrepeats))>.999)[1]
q
fdrest <- fdrEnrichedCounts(nrepeats, use=1:q, components=1)
numRepeArtif <- rownames(fdrest[fdrest$fdrEnriched<0.01,])[1]
numRepeArtif


###################################################
### code chunk number 12: figrepeats1
###################################################
plot(fdrest$fdrEnriched,type='l',ylab='',xlab='Number of repeats (r)')
lines(fdrest$pdfOverall,col=4)
lines(fdrest$pdfH0,col=2,lty=2)
legend('topright',c('FDR','P(r | no artifact)','P(r)'),lty=c(1,2,1),col=c(1,2,4))


###################################################
### code chunk number 13: figrepeats1
###################################################
plot(fdrest$fdrEnriched,type='l',ylab='',xlab='Number of repeats (r)')
lines(fdrest$pdfOverall,col=4)
lines(fdrest$pdfH0,col=2,lty=2)
legend('topright',c('FDR','P(r | no artifact)','P(r)'),lty=c(1,2,1),col=c(1,2,4))


###################################################
### code chunk number 14: figpreprocess1
###################################################
covbefore <- coverage(htSample[['ipBatch2']])
covbefore <- seqselect(covbefore[['chr2L']],295108,297413)
plot(as.integer(covbefore),type='l',ylab='Coverage')


###################################################
### code chunk number 15: preprocess1
###################################################
ip2Aligned <- alignPeaks(htSample[['ipBatch2']],strand='strand', npeaks=100)
covafter <- coverage(ip2Aligned)
covafter <- seqselect(covafter[['chr2L']],295108,297413)
covplus <- coverage(htSample[['ipBatch2']][htSample[['ipBatch2']][['strand']]=='+',])
covplus <- seqselect(covplus[['chr2L']],295108,297413)
covminus <- coverage(htSample[['ipBatch2']][htSample[['ipBatch2']][['strand']]=='-',])
covminus <- seqselect(covminus[['chr2L']],295108,297413)


###################################################
### code chunk number 16: figpreprocess2
###################################################
plot(as.integer(covafter),type='l',ylab='Coverage')
lines(as.integer(covplus),col='blue',lty=2)
lines(as.integer(covminus),col='red',lty=2)


###################################################
### code chunk number 17: figpreprocess1
###################################################
covbefore <- coverage(htSample[['ipBatch2']])
covbefore <- seqselect(covbefore[['chr2L']],295108,297413)
plot(as.integer(covbefore),type='l',ylab='Coverage')


###################################################
### code chunk number 18: figpreprocess2
###################################################
plot(as.integer(covafter),type='l',ylab='Coverage')
lines(as.integer(covplus),col='blue',lty=2)
lines(as.integer(covminus),col='red',lty=2)


###################################################
### code chunk number 19: islandCounts
###################################################
ip <- rbind(htSample[[2]],htSample[[4]])
ctrl <- rbind(htSample[[1]],htSample[[3]])
pool <- RangedDataList(ip=ip, ctrl=ctrl)
counts <- islandCounts(pool,minReads=10)
head(counts)


###################################################
### code chunk number 20: enrichedRegions
###################################################
mappedreads <- c(ip=nrow(ip),ctrl=nrow(ctrl))
mappedreads
regions <- enrichedRegions(sample1=ip,sample2=ctrl,minReads=10,mappedreads=mappedreads,p.adjust.method='BY',pvalFilter=.05,twoTailed=FALSE)
head(regions)
nrow(regions)


###################################################
### code chunk number 21: enrichedPeaks
###################################################
peaks <- enrichedPeaks(regions, sample1=ip, sample2=ctrl, minHeight=100)
peaks


###################################################
### code chunk number 22: mergeRegions
###################################################
mergeRegions(peaks, maxDist=300, score='height', aggregateFUN='median')


###################################################
### code chunk number 23: analysis2
###################################################
chrLength <- 500000
names(chrLength) <- c('chr2L')
chrregions <- enrichedChrRegions(peaks, chrLength=chrLength, windowSize=10^4-1, fdr=0.05, nSims=1)
chrregions[['chr2L']]


###################################################
### code chunk number 24: figanalysis2
###################################################
chrregions <- RangedData(ranges=IRanges(start=c(100000,200000),end=c(100100,210000)),space=c('chr2L','chr2R'))
plotChrRegions(regions=chrregions, chrLength=c(chr2L=500000,chr2R=500000))


###################################################
### code chunk number 25: figanalysis2
###################################################
chrregions <- RangedData(ranges=IRanges(start=c(100000,200000),end=c(100100,210000)),space=c('chr2L','chr2R'))
plotChrRegions(regions=chrregions, chrLength=c(chr2L=500000,chr2R=500000))


###################################################
### code chunk number 26: figplots1
###################################################
set.seed(1)
peaksanno <- peaks
peaksanno[['start_position']] <- start(peaksanno) + runif(nrow(peaksanno),-500,1000)
peaksanno[['end_position']] <- peaksanno[['start_position']] + 500
peaksanno[['distance']] <- peaksanno[['start_position']] - start(peaksanno)
peaksanno[['strand']] <- sample(c('+','-'),nrow(peaksanno),replace=TRUE)
PeakLocation(peaksanno,peakDistance=1000)


###################################################
### code chunk number 27: figplots1
###################################################
set.seed(1)
peaksanno <- peaks
peaksanno[['start_position']] <- start(peaksanno) + runif(nrow(peaksanno),-500,1000)
peaksanno[['end_position']] <- peaksanno[['start_position']] + 500
peaksanno[['distance']] <- peaksanno[['start_position']] - start(peaksanno)
peaksanno[['strand']] <- sample(c('+','-'),nrow(peaksanno),replace=TRUE)
PeakLocation(peaksanno,peakDistance=1000)


###################################################
### code chunk number 28: regionscov
###################################################
cover <- coverage(ip)
rcov <- regionsCoverage(space(regions),start(regions),end(regions),cover=cover)
names(rcov)
rcov[['views']]
gridcov <- gridCoverage(rcov)
dim(getCover(gridcov))
getViewsInfo(gridcov)


###################################################
### code chunk number 29: figregCoverage
###################################################
ylim <- c(0,max(getViewsInfo(gridcov)[['maxCov']]))
plot(gridcov, ylim=ylim,lwd=2)
for (i in 1:nrow(regions)) lines(getCover(gridcov)[i,], col='gray', lty=2)


###################################################
### code chunk number 30: figregCoverage
###################################################
ylim <- c(0,max(getViewsInfo(gridcov)[['maxCov']]))
plot(gridcov, ylim=ylim,lwd=2)
for (i in 1:nrow(regions)) lines(getCover(gridcov)[i,], col='gray', lty=2)