inst/doc/vignette.R

### 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)

Try the R453Plus1Toolbox package in your browser

Any scripts or data that you put into this service are public.

R453Plus1Toolbox documentation built on Nov. 8, 2020, 5:59 p.m.