### R code from vignette source 'vignette.Rnw' ################################################### ### code chunk number 1: vignette.Rnw:34-36 ################################################### options(width=60) options(continue=" ") ################################################### ### code chunk number 2: preliminaries ################################################### library(R453Plus1Toolbox) ################################################### ### code chunk number 3: createRocheAVASet1 ################################################### projectDir = system.file("extdata", "AVASet", package = "R453Plus1Toolbox") avaSet = AVASet(dirname=projectDir) ################################################### ### code chunk number 4: createRocheAVASet2 (eval = FALSE) ################################################### ## projectDir = "My/AVA/Project" ## avaSet = AVASet(dirname=projectDir, avaBin="/home/User/AVA/bin") ################################################### ### code chunk number 5: createRocheAVASet3 (eval = FALSE) ################################################### ## projectDir = system.file("extdata", "AVASet_doAmplicon", package="R453Plus1Toolbox") ## avaSetExample = AVASet(dirname=projectDir, file_sample="sample.csv", ## file_amp="amp.csv", file_reference="reference.csv", file_variant="variant.csv", ## file_variantHits="variantHits.csv") ################################################### ### code chunk number 6: showAVASet ################################################### avaSet ################################################### ### code chunk number 7: assayDataAVA ################################################### assayData(avaSet)$totalForwCount[1:3, ] ################################################### ### code chunk number 8: fDataAVA ################################################### fData(avaSet)[1:3, ] ################################################### ### code chunk number 9: pDataAVA ################################################### pData(avaSet) ################################################### ### code chunk number 10: assayDataAmpAVA ################################################### assayDataAmp(avaSet)$forwCount ################################################### ### code chunk number 11: fDataAmpAVA ################################################### fDataAmp(avaSet) ################################################### ### code chunk number 12: referenceSequences ################################################### library(ShortRead) referenceSequences(avaSet) sread(referenceSequences(avaSet)) ################################################### ### code chunk number 13: AVASubSet1 ################################################### avaSubSet = avaSet[1:10, "Sample_1"] ################################################### ### code chunk number 14: AVASubSet2 ################################################### avaSubSet = subset(avaSet, subset=1:10, dimension="variants") ################################################### ### code chunk number 15: AVASubSet3 ################################################### avaSubSet = subset(subset(avaSet, subset=1:10, dimension="variants"), subset="Sample_1", dimension="samples") ################################################### ### code chunk number 16: AVASubSet4 ################################################### avaSubSet = subset(avaSet, subset=c("TET2_E11.04", "TET2_E06"), dimension="amplicons") ################################################### ### code chunk number 17: filterAVA1 ################################################### avaSetFiltered1 = setVariantFilter(avaSet, filter=0.05) ################################################### ### code chunk number 18: filterAVA2 ################################################### avaSetFiltered2 = setVariantFilter(avaSet, filter=c(0.1, 0.05)) ################################################### ### code chunk number 19: filterAVA4 ################################################### avaSet = setVariantFilter(avaSetFiltered1, filter=0) ################################################### ### code chunk number 20: filterAVA5 ################################################### avaSet = setVariantFilter(avaSetFiltered2) ################################################### ### code chunk number 21: varPercentages1 ################################################### getVariantPercentages(avaSet, direction="both")[20:25, 1:4] ################################################### ### code chunk number 22: varPercentages2 ################################################### (assayData(avaSet)[[1]] + assayData(avaSet)[[3]]) / (assayData(avaSet)[[2]] + assayData(avaSet)[[4]]) ################################################### ### code chunk number 23: alignShortReads (eval = FALSE) ################################################### ## library(BSgenome.Hsapiens.UCSC.hg19) ## seqNames = names(Hsapiens)[1:24] ## avaSet = alignShortReads(avaSet, bsGenome=Hsapiens, ## seqNames=seqNames, ensemblNotation=TRUE) ################################################### ### code chunk number 24: annotateVariants (eval = FALSE) ################################################### ## avaSet = setVariantFilter(avaSet, filter=0.05) ## avaAnnot = annotateVariants(avaSet) ################################################### ### code chunk number 25: htmlReport (eval = FALSE) ################################################### ## 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) ################################################### ### code chunk number 26: plotAmpCov1 ################################################### plotAmpliconCoverage(avaSet[, 2], type="amplicon") ################################################### ### code chunk number 27: plotAmpCov2 ################################################### plotAmpliconCoverage(avaSet, bothDirections=TRUE, type="amplicon") ################################################### ### code chunk number 28: loadXLS ################################################### file = system.file("extdata", "AVAVarFreqExport", "AVAVarFreqExport.xls", package="R453Plus1Toolbox") ################################################### ### code chunk number 29: plotVarFreq ################################################### plotVariationFrequency(file, plotRange=c(50, 150)) ################################################### ### code chunk number 30: loadXLS ################################################### data(plotVariantsExample) ################################################### ### code chunk number 31: plotVariants ################################################### geneInfo = plotVariants(data=variants, gene="ENSG00000168769", transcript="ENST00000513237", regions=regions, mutationInfo=mutationInfo, horiz=TRUE, cex=0.8) ################################################### ### code chunk number 32: vcfexport (eval = FALSE) ################################################### ## ava2vcf(avaSet, filename="variants.vcf", annot=avaAnnot) ################################################### ### code chunk number 33: gsmDir1 ################################################### dir_sample01 = system.file("extdata", "MapperSet", "N01", package = "R453Plus1Toolbox") ################################################### ### code chunk number 34: gsmDir2 ################################################### dir_sample03 = system.file("extdata", "MapperSet", "N03", package = "R453Plus1Toolbox") ################################################### ### code chunk number 35: gsmDir3 ################################################### dir_sample04 = system.file("extdata", "MapperSet", "N04", package = "R453Plus1Toolbox") ################################################### ### code chunk number 36: gsmDir4 ################################################### dirs = c(dir_sample01, dir_sample03, dir_sample04) ################################################### ### code chunk number 37: createRocheGSMSet ################################################### mapperSet = MapperSet(dirs=dirs, samplenames=c("N01", "N03", "N04")) ################################################### ### code chunk number 38: showMapperSet ################################################### mapperSet ################################################### ### code chunk number 39: assayDataMapper1 ################################################### assayData(mapperSet)$variantForwCount[1:4, ] ################################################### ### code chunk number 40: assayDataMapper2 ################################################### assayData(mapperSet)$totalForwCount[1:4, ] ################################################### ### code chunk number 41: fDataMapper ################################################### fData(mapperSet)[1:4, ] ################################################### ### code chunk number 42: pDataMapper ################################################### pData(mapperSet) ################################################### ### code chunk number 43: annotateVarMapper (eval = FALSE) ################################################### ## mapperAnnot = annotateVariants(mapperSet) ################################################### ### code chunk number 44: htmlReportMapper (eval = FALSE) ################################################### ## htmlReport(mapperSet, annot=mapperAnnot, dir="htmlReportExampleMapper", ## title="htmlReport Example", minMut=3) ################################################### ### code chunk number 45: demultiplex ################################################### fnaFile = system.file("extdata", "SVDetection", "R_2009_07_30", "D_2009_07_31", "1.TCA.454Reads.fna", package="R453Plus1Toolbox") seqs = readDNAStringSet(fnaFile, format="fasta") MIDSeqs = genomeSequencerMIDs(c("MID1", "MID2", "MID3")) dmReads = demultiplexReads(seqs, MIDSeqs, numMismatches=2, trim=TRUE) length(seqs) sum(sapply(dmReads, length)) ################################################### ### code chunk number 46: removeLinker ################################################### 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) }) ################################################### ### code chunk number 47: writeFASTA (eval = FALSE) ################################################### ## write.XStringSet(trimReads[["MID1"]], file="/tmp/N01.fasta", format="fasta") ################################################### ### code chunk number 48: readBam ################################################### library("Rsamtools") bamFile = system.file("extdata", "SVDetection", "bam", "N01.bam", package="R453Plus1Toolbox") parameters = ScanBamParam(what=scanBamWhat()) bam = scanBam(bamFile, param=parameters) ################################################### ### code chunk number 49: filterReads ################################################### library("rtracklayer") bedFile = system.file("extdata", "SVDetection", "chip", "CaptureArray_hg19.bed", package="R453Plus1Toolbox") chip = import.ucsc(bedFile, subformat="bed") chip = split(ranges(chip[[1]]), seqnames(chip[[1]])) names(chip) = gsub("chr", "", names(chip)) linker = sequenceCaptureLinkers("gSel3")[[1]] filterReads = filterChimericReads(bam, targetRegion=chip, linkerSeq=linker) filterReads$log ################################################### ### code chunk number 50: detectBreakpoints ################################################### bp = detectBreakpoints(filterReads, minClusterSize=1) bp table(bp) mbp = mergeBreakpoints(bp) summary(mbp) ################################################### ### code chunk number 51: plotCR1 ################################################### plotChimericReads(mbp[1], legend=TRUE) ################################################### ### code chunk number 52: plotCR2 ################################################### plotChimericReads(mbp[1], plotBasePairs=TRUE, maxBasePairs=30) ################################################### ### code chunk number 53: importSFFfiles ################################################### file <- system.file("extdata", "SFF", "example.sff", package="R453Plus1Toolbox") sffContainer <- readSFF(file) ################################################### ### code chunk number 54: SFFcontainer ################################################### showClass("SFFContainer") ################################################### ### code chunk number 55: sffreads ################################################### reads(sffContainer) ################################################### ### code chunk number 56: sffsubsetting ################################################### subSffContainer <- sffContainer[1:5] ################################################### ### code chunk number 57: sffsubsetting (eval = FALSE) ################################################### ## writeSFF(subSffContainer, subSffFile.sff) ################################################### ### code chunk number 58: positionQualityBoxplot ################################################### positionQualityBoxplot(sffContainer) ################################################### ### code chunk number 59: dinucleotideOddsRatio ################################################### dinucleotideOddsRatio(sffContainer)