TEQC.R
In TEQC: Quality control for target capture experiments

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

###################################################
### code chunk number 1: load
###################################################
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


###################################################
### code chunk number 2: fraction.target
###################################################
ft <- fraction.target(targets, genome = "hg38")
ft


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


###################################################
### code chunk number 4: reads2pairs
###################################################
readpairs <- reads2pairs(reads)
readpairs


###################################################
### code chunk number 5: insertsizehist
###################################################
insert.size.hist(readpairs, breaks=10)


###################################################
### code chunk number 6: chrombarplot
###################################################
chrom.barplot(reads, targets)


###################################################
### code chunk number 7: fraction.reads.target
###################################################
fr <- fraction.reads.target(reads, targets)
fr


###################################################
### code chunk number 8: withoffset
###################################################
fraction.reads.target(reads, targets, Offset=100)


###################################################
### code chunk number 9: enrichment
###################################################
fr / ft


###################################################
### code chunk number 10: fracpairs
###################################################
fraction.reads.target(readpairs, targets)


###################################################
### code chunk number 11: coverage.target
###################################################
Coverage <- coverage.target(reads, targets, perTarget=T, perBase=T)
Coverage
targets2 <- Coverage$targetCoverages


###################################################
### code chunk number 12: readspertarget
###################################################
targets2 <- readsPerTarget(reads, targets2)
targets2


###################################################
### code chunk number 13: covered.k
###################################################
covered.k(Coverage$coverageTarget, k=c(1, 5, 10))


###################################################
### code chunk number 14: coveragehist
###################################################
coverage.hist(Coverage$coverageTarget, covthreshold=8)


###################################################
### code chunk number 15: coverageuniformity
###################################################
coverage.uniformity(Coverage)


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


###################################################
### code chunk number 17: get.baits
###################################################
baitsfile <- paste(exptPath, "ExampleSet_Baits.txt", sep="/")
baits <- get.baits(baitsfile, chrcol=3, startcol=4, endcol=5, seqcol=2)


###################################################
### code chunk number 18: coverageGC
###################################################
coverage.GC(Coverage$coverageAll, baits, pch=16, cex=1.5)


###################################################
### code chunk number 19: coverageplot
###################################################
coverage.plot(Coverage$coverageAll, targets, Offset=100, chr="chr1", Start=11157524, End=11158764)


###################################################
### code chunk number 20: make.wigfiles (eval = FALSE)
###################################################
## make.wigfiles(Coverage$coverageAll)


###################################################
### code chunk number 21: duplicatesbarplot
###################################################
duplicates.barplot(reads, targets)


###################################################
### code chunk number 22: duplicatesbarplot2
###################################################
duplicates.barplot(readpairs, targets, ylab="Fraction of read pairs")


###################################################
### code chunk number 23: collapse1
###################################################
reads.collapsed <- unique(reads)


###################################################
### code chunk number 24: collapse2
###################################################
ID.nondups <- names(unique(readpairs))
reads.collapsed.pairs <- reads[names(reads) %in% ID.nondups,,drop=T]


###################################################
### code chunk number 25: cov1
###################################################
coverage.target(reads.collapsed, targets, perBase=F, perTarget=F)


###################################################
### code chunk number 26: cov2
###################################################
coverage.target(reads.collapsed.pairs, targets, perBase=F, perTarget=F)


###################################################
### code chunk number 27: newsample
###################################################
r <- sample(length(reads), 0.1 * length(reads))
reads2 <- reads[-r,,drop=T]
Coverage2 <- coverage.target(reads2, targets, perBase=T)


###################################################
### code chunk number 28: coveragedensity
###################################################
covlist <- list(Coverage, Coverage2)
par(mfrow=c(1,2))
coverage.density(covlist)
coverage.density(covlist, normalized=F)


###################################################
### code chunk number 29: coverageuniformity2
###################################################
coverage.uniformity(Coverage, addlines=F)
coverage.uniformity(Coverage2, addlines=F, add=T, col="blue", lty=2)


###################################################
### code chunk number 30: coverageplot2
###################################################
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)


###################################################
### code chunk number 31: covcor
###################################################
coverage.correlation(covlist, plotfrac=0.1, cex.pch=4)


###################################################
### code chunk number 32: sessioninfo
###################################################
toLatex(sessionInfo())