###################################################
### chunk number 1: load
###################################################
#line 85 "vignettes/TEQC/inst/doc/TEQC.Rnw"
library(TEQC)

exptPath <- system.file("extdata", package="TEQC")
targets <- get.targets(targetsfile=paste(exptPath, "ExampleSet_Targets.bed", sep="/"), chrcol=1, startcol=2, endcol=3, skip=0)
targets


###################################################
### chunk number 2: fraction.target
###################################################
#line 108 "vignettes/TEQC/inst/doc/TEQC.Rnw"
ft <- fraction.target(targets, genome="hg19")
ft


###################################################
### chunk number 3: reads
###################################################
#line 125 "vignettes/TEQC/inst/doc/TEQC.Rnw"
reads <- get.reads(paste(exptPath, "ExampleSet_Reads.bed", sep="/"), chrcol=1, startcol=2, endcol=3, idcol=4, zerobased=F, skip=0)
reads


###################################################
### chunk number 4: reads2pairs
###################################################
#line 152 "vignettes/TEQC/inst/doc/TEQC.Rnw"
readpairs <- reads2pairs(reads)
readpairs


###################################################
### chunk number 5: insertsizehist
###################################################
#line 174 "vignettes/TEQC/inst/doc/TEQC.Rnw"
insert.size.hist(readpairs, breaks=10)


###################################################
### chunk number 6: chrombarplot
###################################################
#line 209 "vignettes/TEQC/inst/doc/TEQC.Rnw"
chrom.barplot(reads, targets)


###################################################
### chunk number 7: fraction.reads.target
###################################################
#line 225 "vignettes/TEQC/inst/doc/TEQC.Rnw"
fr <- fraction.reads.target(reads, targets)
fr


###################################################
### chunk number 8: withoffset
###################################################
#line 237 "vignettes/TEQC/inst/doc/TEQC.Rnw"
fraction.reads.target(reads, targets, Offset=100)


###################################################
### chunk number 9: enrichment
###################################################
#line 251 "vignettes/TEQC/inst/doc/TEQC.Rnw"
fr / ft


###################################################
### chunk number 10: fracpairs
###################################################
#line 263 "vignettes/TEQC/inst/doc/TEQC.Rnw"
fraction.reads.target(readpairs, targets)


###################################################
### chunk number 11: coverage.target
###################################################
#line 288 "vignettes/TEQC/inst/doc/TEQC.Rnw"
Coverage <- coverage.target(reads, targets, perTarget=T, perBase=T)
Coverage
targets2 <- Coverage$targetCoverages


###################################################
### chunk number 12: readspertarget
###################################################
#line 307 "vignettes/TEQC/inst/doc/TEQC.Rnw"
targets2 <- readsPerTarget(reads, targets2)
targets2


###################################################
### chunk number 13: covered.k
###################################################
#line 326 "vignettes/TEQC/inst/doc/TEQC.Rnw"
covered.k(Coverage$coverageTarget, k=c(1, 5, 10))


###################################################
### chunk number 14: coveragehist
###################################################
#line 342 "vignettes/TEQC/inst/doc/TEQC.Rnw"
coverage.hist(Coverage$coverageTarget, covthreshold=8)


###################################################
### chunk number 15: coverageuniformity
###################################################
#line 366 "vignettes/TEQC/inst/doc/TEQC.Rnw"
coverage.uniformity(Coverage)


###################################################
### chunk number 16: coveragetargetlength
###################################################
#line 396 "vignettes/TEQC/inst/doc/TEQC.Rnw"
par(mfrow=c(1,2))
coverage.targetlength.plot(targets2, plotcolumn="nReads", pch=16, cex=1.5)
coverage.targetlength.plot(targets2, plotcolumn="avgCoverage", pch=16, cex=1.5)


###################################################
### chunk number 17: get.baits
###################################################
#line 419 "vignettes/TEQC/inst/doc/TEQC.Rnw"
baitsfile <- paste(exptPath, "ExampleSet_Baits.txt", sep="/")
baits <- get.baits(baitsfile, chrcol=3, startcol=4, endcol=5, seqcol=2)


###################################################
### chunk number 18: coverageGC
###################################################
#line 438 "vignettes/TEQC/inst/doc/TEQC.Rnw"
coverage.GC(Coverage$coverageAll, baits, pch=16, cex=1.5)


###################################################
### chunk number 19: coverageplot
###################################################
#line 459 "vignettes/TEQC/inst/doc/TEQC.Rnw"
coverage.plot(Coverage$coverageAll, targets, Offset=100, chr="chr1", Start=11157524, End=11158764)


###################################################
### chunk number 20: make.wigfiles eval=FALSE
###################################################
## #line 477 "vignettes/TEQC/inst/doc/TEQC.Rnw"
## make.wigfiles(Coverage$coverageAll)


###################################################
### chunk number 21: duplicatesbarplot
###################################################
#line 505 "vignettes/TEQC/inst/doc/TEQC.Rnw"
duplicates.barplot(reads, targets)


###################################################
### chunk number 22: duplicatesbarplot2
###################################################
#line 534 "vignettes/TEQC/inst/doc/TEQC.Rnw"
duplicates.barplot(readpairs, targets, ylab="Fraction of read pairs")


###################################################
### chunk number 23: collapse1
###################################################
#line 552 "vignettes/TEQC/inst/doc/TEQC.Rnw"
dupfun <- function(x) duplicated(x$ranges)
params <- RDApplyParams(rangedData=reads, applyFun=dupfun)
dups <- unlist(rdapply(params))
reads.collapsed <- reads[!dups,,drop=T]


###################################################
### chunk number 24: collapse2
###################################################
#line 565 "vignettes/TEQC/inst/doc/TEQC.Rnw"
params2 <- RDApplyParams(rangedData=readpairs, applyFun=dupfun)
dups2 <- unlist(rdapply(params2))
ID.nondups <- readpairs$ID[!dups2]
sel <- reads$ID %in% ID.nondups
reads.collapsed.pairs <- reads[sel,,drop=T]


###################################################
### chunk number 25: cov1
###################################################
#line 579 "vignettes/TEQC/inst/doc/TEQC.Rnw"
coverage.target(reads.collapsed, targets, perBase=F, perTarget=F)


###################################################
### chunk number 26: cov2
###################################################
#line 589 "vignettes/TEQC/inst/doc/TEQC.Rnw"
coverage.target(reads.collapsed.pairs, targets, perBase=F, perTarget=F)


###################################################
### chunk number 27: newsample
###################################################
#line 606 "vignettes/TEQC/inst/doc/TEQC.Rnw"
r <- sample(nrow(reads), 0.1 * nrow(reads))
reads2 <- reads[-r,,drop=T]
Coverage2 <- coverage.target(reads2, targets, perBase=T)


###################################################
### chunk number 28: coveragedensity
###################################################
#line 629 "vignettes/TEQC/inst/doc/TEQC.Rnw"
covlist <- list(Coverage, Coverage2)
par(mfrow=c(1,2))
coverage.density(covlist)
coverage.density(covlist, normalized=F)


###################################################
### chunk number 29: coverageuniformity2
###################################################
#line 656 "vignettes/TEQC/inst/doc/TEQC.Rnw"
coverage.uniformity(Coverage, addlines=F)
coverage.uniformity(Coverage2, addlines=F, add=T, col="blue", lty=2)


###################################################
### chunk number 30: coverageplot2
###################################################
#line 671 "vignettes/TEQC/inst/doc/TEQC.Rnw"
coverage.plot(Coverage$coverageAll, targets, Offset=100, chr="chr1", Start=11157524, End=11158764)
coverage.plot(Coverage2$coverageAll, add=T, col.line=2, chr="chr1", Start=11157524, End=11158764)


###################################################
### chunk number 31: covcor
###################################################
#line 698 "vignettes/TEQC/inst/doc/TEQC.Rnw"
coverage.correlation(covlist, plotfrac=0.1, cex.pch=4)


###################################################
### chunk number 32: sessioninfo
###################################################
#line 719 "vignettes/TEQC/inst/doc/TEQC.Rnw"
toLatex(sessionInfo())