###################################################
### chunk number 1: 
###################################################
#line 34 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
options(width=60)
options(continue=" ")


###################################################
### chunk number 2: preliminaries
###################################################
#line 47 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
library(R453Plus1Toolbox)


###################################################
### chunk number 3: createRocheAVASet
###################################################
#line 71 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
avaDir = system.file("extdata", "AVASet", package = "R453Plus1Toolbox") 
avaSet = AVASet(avaDir)


###################################################
### chunk number 4: showAVASet
###################################################
#line 82 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
avaSet


###################################################
### chunk number 5: assayDataAVA
###################################################
#line 104 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
assayData(avaSet)$totalForwCount[1:3, ]


###################################################
### chunk number 6: fDataAVA
###################################################
#line 116 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
fData(avaSet)[1:3, ]


###################################################
### chunk number 7: pDataAVA
###################################################
#line 122 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
pData(avaSet)


###################################################
### chunk number 8: assayDataAmpAVA
###################################################
#line 136 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
assayDataAmp(avaSet)$forwCount


###################################################
### chunk number 9: fDataAmpAVA
###################################################
#line 147 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
fDataAmp(avaSet)


###################################################
### chunk number 10: referenceSequences
###################################################
#line 157 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
library(ShortRead)
referenceSequences(avaSet)
sread(referenceSequences(avaSet))


###################################################
### chunk number 11: AVASubSet1
###################################################
#line 206 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
avaSubSet = avaSet[1:10, "Sample_1"]


###################################################
### chunk number 12: AVASubSet2
###################################################
#line 214 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
avaSubSet = subset(avaSet, subset=1:10, dimension="variants")


###################################################
### chunk number 13: AVASubSet3
###################################################
#line 218 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
avaSubSet = subset(subset(avaSet, subset=1:10, dimension="variants"), subset="Sample_1", dimension="samples")


###################################################
### chunk number 14: AVASubSet4
###################################################
#line 223 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
avaSubSet = subset(avaSet, subset=c("TET2_E11.04", "TET2_E06"), dimension="amplicons")


###################################################
### chunk number 15: filterAVA1
###################################################
#line 240 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
avaSetFiltered1 = setVariantFilter(avaSet, filter=0.05)


###################################################
### chunk number 16: filterAVA2
###################################################
#line 247 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
avaSetFiltered2 = setVariantFilter(avaSet, filter=c(0.1, 0.05))


###################################################
### chunk number 17: filterAVA4
###################################################
#line 254 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
avaSet = setVariantFilter(avaSetFiltered1, filter=0)


###################################################
### chunk number 18: filterAVA5
###################################################
#line 258 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
avaSet = setVariantFilter(avaSetFiltered2)


###################################################
### chunk number 19: varPercentages1
###################################################
#line 267 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
getVariantPercentages(avaSet, direction="both")[20:25, 1:4]


###################################################
### chunk number 20: varPercentages2
###################################################
#line 271 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
(assayData(avaSet)[[1]] + assayData(avaSet)[[3]]) / (assayData(avaSet)[[2]] + assayData(avaSet)[[4]])


###################################################
### chunk number 21: alignShortReads eval=FALSE
###################################################
## #line 286 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
## library(BSgenome.Hsapiens.UCSC.hg19)
## seqNames = names(Hsapiens)[1:24]
## avaSet = alignShortReads(avaSet, bsGenome=Hsapiens, seqNames=seqNames, ensemblNotation=TRUE)


###################################################
### chunk number 22: annotateVariants eval=FALSE
###################################################
## #line 298 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
## avaSet = setVariantFilter(avaSet, filter=0.05)
## avaAnnot = annotateVariants(avaSet)


###################################################
### chunk number 23: htmlReport eval=FALSE
###################################################
## #line 324 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
## blocks = as.character(sapply(annotatedVariants(avaAnnot), function(x) x$genes$external_gene_id))
## htmlReport(avaSet, annot=avaAnnot, blocks=blocks, dir="htmlReportExampleAVA", title="htmlReport Example", minMut=3)


###################################################
### chunk number 24: plotAmpCov1
###################################################
#line 337 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
plotAmpliconCoverage(avaSet[, 2], type="amplicon")


###################################################
### chunk number 25: plotAmpCov2
###################################################
#line 346 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
plotAmpliconCoverage(avaSet, bothDirections=TRUE, type="amplicon")


###################################################
### chunk number 26: loadXLS
###################################################
#line 365 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
file = system.file("extdata", "AVAVarFreqExport", "AVAVarFreqExport.xls", package="R453Plus1Toolbox")


###################################################
### chunk number 27: plotVarFreq
###################################################
#line 368 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
plotVariationFrequency(file, plotRange=c(50, 150))


###################################################
### chunk number 28: loadXLS
###################################################
#line 391 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
data(avaSetFiltered_annot)
data(transcriptdf)


###################################################
### chunk number 29: plotVariants
###################################################
#line 395 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
plotVariants(avaSetFiltered_annot, transcript=transcriptdf, legend=TRUE, regions=c(700,1600))


###################################################
### chunk number 30: gsmDir1
###################################################
#line 426 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
dir_sample01 = system.file("extdata", "MapperSet", "N01", package = "R453Plus1Toolbox")


###################################################
### chunk number 31: gsmDir2
###################################################
#line 429 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
dir_sample03 = system.file("extdata", "MapperSet", "N03", package = "R453Plus1Toolbox")


###################################################
### chunk number 32: gsmDir3
###################################################
#line 432 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
dir_sample04 = system.file("extdata", "MapperSet", "N04", package = "R453Plus1Toolbox")


###################################################
### chunk number 33: gsmDir4
###################################################
#line 434 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
dirs = c(dir_sample01, dir_sample03, dir_sample04)


###################################################
### chunk number 34: createRocheGSMSet
###################################################
#line 437 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
mapperSet = MapperSet(dirs=dirs, samplenames=c("N01", "N03", "N04"))


###################################################
### chunk number 35: showMapperSet
###################################################
#line 472 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
mapperSet


###################################################
### chunk number 36: assayDataMapper1
###################################################
#line 478 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
assayData(mapperSet)$variantForwCount[1:4, ]


###################################################
### chunk number 37: assayDataMapper2
###################################################
#line 481 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
assayData(mapperSet)$totalForwCount[1:4, ]


###################################################
### chunk number 38: fDataMapper
###################################################
#line 484 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
fData(mapperSet)[1:4, ]


###################################################
### chunk number 39: pDataMapper
###################################################
#line 487 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
pData(mapperSet)


###################################################
### chunk number 40: annotateVarMapper eval=FALSE
###################################################
## #line 511 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
## mapperAnnot = annotateVariants(mapperSet)


###################################################
### chunk number 41: htmlReportMapper eval=FALSE
###################################################
## #line 531 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
## htmlReport(mapperSet, annot=mapperAnnot, dir="htmlReportExampleMapper", title="htmlReport Example", minMut=3)


###################################################
### chunk number 42: demultiplex
###################################################
#line 562 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
fnaFile = system.file("extdata", "StructuralVariantDetection", "R_2009_07_30_14_26_52_FLX12080469_Administrator_714003", "D_2009_07_31_08_45_30_flxcluster_fullProcessing_714003", "1.TCA.454Reads.fna", package="R453Plus1Toolbox")
seqs = read.DNAStringSet(fnaFile, format="fasta")
MIDSeqs = genomeSequencerMIDs(c("MID1", "MID2", "MID3"))
dmReads = demultiplexReads(seqs, MIDSeqs, numMismatches=2, trim=TRUE)
length(seqs)
sum(sapply(dmReads, length))


###################################################
### chunk number 43: removeLinker
###################################################
#line 581 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
minReadLength = 15
gSel3 = sequenceCaptureLinkers("gSel3")[[1]]
trimReads = lapply(dmReads, function (reads) {
reads = reads[width(reads) >= minReadLength]
reads = removeLinker(reads, gSel3)
reads = reads[width(reads) >= minReadLength]
readsRev = reverseComplement(reads)
readsRev = removeLinker(readsRev, gSel3)
reads = reverseComplement(readsRev)
reads = reads[width(reads) >= minReadLength]
return(reads)
})


###################################################
### chunk number 44: writeFASTA eval=FALSE
###################################################
## #line 602 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
## write.XStringSet(trimReads[["MID1"]], file="/tmp/N01.fasta", format="fasta")


###################################################
### chunk number 45: readBam
###################################################
#line 630 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
library("Rsamtools")
bamFile = system.file("extdata", "StructuralVariantDetection", "bam", "N01.bam", package="R453Plus1Toolbox")
parameters = ScanBamParam()
bam = scanBam(bamFile, param=parameters)


###################################################
### chunk number 46: filterReads
###################################################
#line 641 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
library("rtracklayer")
bedFile = system.file("extdata", "StructuralVariantDetection", "chip", "CaptureArray_hg19.bed", package="R453Plus1Toolbox")
chip = import.ucsc(bedFile, subformat="bed")
chip = ranges(chip[[1]])
names(chip) = gsub("chr", "", names(chip))

linker = sequenceCaptureLinkers("gSel3")[[1]]

filterReads = filterChimericReads(bam, targetRegion=chip, linkerSeq=linker)
filterReads$log


###################################################
### chunk number 47: detectBreakpoints
###################################################
#line 658 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
bp = detectBreakpoints(filterReads, minClusterSize=1)
bp
table(bp)

mbp = mergeBreakpoints(bp)
summary(mbp)


###################################################
### chunk number 48: plotCR1
###################################################
#line 683 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
plotChimericReads(mbp[1], legend=TRUE)


###################################################
### chunk number 49: plotCR2
###################################################
#line 697 "vignettes/R453Plus1Toolbox/inst/doc/vignette.Rnw"
plotChimericReads(mbp[1], plotBasePairs=TRUE, maxBasePairs=30)