TarSeqQC-vignette.R
In TarSeqQC: TARgeted SEQuencing Experiment Quality Control

### R code from vignette source 'TarSeqQC-vignette.Rnw'

###################################################
### code chunk number 1: General options for R
###################################################
options(prompt="> ", continue="+  ", width=90, useFancyQuotes=FALSE, digits=3)


###################################################
### code chunk number 2: TarSeqQC-vignette.Rnw:171-174 (eval = FALSE)
###################################################
## if (!requireNamespace("BiocManager", quietly=TRUE))
##     install.packages("BiocManager")
## BiocManager::install("TarSeqQC")


###################################################
### code chunk number 3: TarSeqQC-vignette.Rnw:185-186
###################################################
citation("TarSeqQC")


###################################################
### code chunk number 4: TarSeqQC-vignette.Rnw:286-287
###################################################
bedFile<-system.file("extdata", "mybed.bed", package="TarSeqQC", mustWork=TRUE)


###################################################
### code chunk number 5: TarSeqQC-vignette.Rnw:310-311
###################################################
bamFile<-system.file("extdata", "mybam.bam", package="TarSeqQC", mustWork=TRUE)


###################################################
### code chunk number 6: TarSeqQC-vignette.Rnw:324-326
###################################################
fastaFile<-system.file("extdata", "myfasta.fa", package="TarSeqQC", 
                        mustWork=TRUE)


###################################################
### code chunk number 7: library call
###################################################
suppressMessages(library("TarSeqQC"))
suppressMessages(library("BiocParallel"))


###################################################
### code chunk number 8: bedFasta control
###################################################
checkBedFasta(bedFile , fastaFile)


###################################################
### code chunk number 9: constructor (eval = FALSE)
###################################################
## library("TarSeqQC")
## library("BiocParallel")
## BPPARAM<-bpparam()
## myPanel<-TargetExperiment(bedFile, bamFile, fastaFile, feature="amplicon", 
##                         attribute="coverage", BPPARAM=BPPARAM)


###################################################
### code chunk number 10: loading TargetExperiment object
###################################################
BPPARAM<-bpparam()
data(ampliPanel, package="TarSeqQC")
myPanel<-ampliPanel
rm(ampliPanel)
setBamFile(myPanel)<-system.file("extdata", "mybam.bam", package="TarSeqQC",
                                    mustWork=TRUE)
setFastaFile(myPanel)<-system.file("extdata", "myfasta.fa", 
                                    package="TarSeqQC", mustWork=TRUE)



###################################################
### code chunk number 11: setters
###################################################
# set feature slot value
setFeature(myPanel)<-"amplicon"
# set attribute slot value
setAttribute(myPanel)<-"coverage"


###################################################
### code chunk number 12: setters (eval = FALSE)
###################################################
## # set scanBamP slot value
## scanBamP<-ScanBamParam()
## #set scanBamP which slot
## bamWhich(scanBamP)<-getBedFile(myPanel)
## setScanBamP(myPanel)<-scanBamP
## # set pileupP slot value
## setPileupP(myPanel)<-PileupParam(max_depth=1000)
## # build the featurePanel again
## setFeaturePanel(myPanel)<-buildFeaturePanel(myPanel, BPPARAM)
## # build the genePanel again
## setGenePanel(myPanel)<-summarizePanel(myPanel, BPPARAM)


###################################################
### code chunk number 13: loading TargetExperiment object
###################################################
data(ampliPanel, package="TarSeqQC")


###################################################
### code chunk number 14: re-defininf file paths
###################################################
# Defining bam file and fasta file names and paths
setBamFile(ampliPanel)<-system.file("extdata", "mybam.bam", 
    package="TarSeqQC", mustWork=TRUE)
setFastaFile(ampliPanel)<-system.file("extdata", "myfasta.fa", 
    package="TarSeqQC", mustWork=TRUE)


###################################################
### code chunk number 15: exploration
###################################################
# show/print
myPanel
# summary
summary(myPanel)
#summary at feature level
summaryFeatureLev(myPanel)
#summary at gene level
summaryGeneLev(myPanel)


###################################################
### code chunk number 16: boxplot and density plot (eval = FALSE)
###################################################
## g<-plotAttrExpl(myPanel,level="feature",join=TRUE, log=FALSE, color="blue")
## x11(type="cairo");
## g


###################################################
### code chunk number 17: length exploration (eval = FALSE)
###################################################
## # explore amplicon length distribution
## plotMetaDataExpl(myPanel, "length", log=FALSE, join=FALSE, color=
##                     "blueviolet")


###################################################
### code chunk number 18: gene exploration (eval = FALSE)
###################################################
## # explore gene's relative frequencies 
## plotMetaDataExpl(myPanel, "gene", abs=FALSE)


###################################################
### code chunk number 19: read percentages and plot (eval = FALSE)
###################################################
## readFrequencies(myPanel)
## plotInOutFeatures(readFrequencies(myPanel))


###################################################
### code chunk number 20: attributeThres
###################################################
# definition of the interval extreme values
attributeThres<-c(0,1,50,200,500, Inf)


###################################################
### code chunk number 21: plot (eval = FALSE)
###################################################
## # plot panel overview
## plot(myPanel, attributeThres=attributeThres, chrLabels =TRUE)


###################################################
### code chunk number 22: plotFeatPerform (eval = FALSE)
###################################################
## # plot panel overview
## g<-plotFeatPerform(myPanel, attributeThres, complete=TRUE, log=FALSE, 
##     featureLabs=TRUE, sepChr=TRUE, legend=TRUE)
## g


###################################################
### code chunk number 23: gc exploration (eval = FALSE)
###################################################
## x11(type="cairo")
## biasExploration(myPanel, source="gc", dens=TRUE)


###################################################
### code chunk number 24: frequency table
###################################################
# summaryIntervals
summaryIntervals(myPanel, attributeThres)


###################################################
### code chunk number 25: coveragePerform (eval = FALSE)
###################################################
## plotAttrPerform(myPanel, attributeThres)


###################################################
### code chunk number 26: low counts features
###################################################
getLowCtsFeatures(myPanel, level="gene", threshold=50)


###################################################
### code chunk number 27: low counts features
###################################################
getLowCtsFeatures(myPanel, level="feature", threshold=50)


###################################################
### code chunk number 28: geneAttrPerFeat (eval = FALSE)
###################################################
## g<-plotGeneAttrPerFeat(myPanel, geneID="gene4")
## # adjust text size
## g<-g+theme(title=element_text(size=16), axis.title=element_text(size=16),
##     legend.text=element_text(size=14))
## g


###################################################
### code chunk number 29: geneAttrPerFeat (eval = FALSE)
###################################################
## g<-plotGeneAttrPerFeat(myPanel, geneID="gene4", overlap = TRUE, level="both")
## g


###################################################
### code chunk number 30: getOverlappedRegions
###################################################
dfRegions<-getOverlappedRegions(myPanel,collapse=TRUE)
head(dfRegions)
#change the region_id field name to consistency of the new bed file
names(dfRegions)[5]<-"name"


###################################################
### code chunk number 31: getOverlappedRegions (eval = FALSE)
###################################################
## #save the new bed file
## write.table(dfRegions[,1:5], file="myRegions.bed", sep="\t", col.names=T,
## row.names=F,quote=F)


###################################################
### code chunk number 32: pileupCounts (eval = FALSE)
###################################################
## # define function parameters
## bed<-getBedFile(myPanel)
## bamFile<-system.file("extdata", "mybam.bam", package="TarSeqQC", mustWork=TRUE)
## fastaFile<-system.file("extdata", "myfasta.fa", package="TarSeqQC", 
##                         mustWork=TRUE)
## scanBamP<-getScanBamP(myPanel)
## pileupP<-getPileupP(myPanel)
## #call pileupCounts function
## myCounts<-pileupCounts(bed=bed, bamFile=bamFile, fastaFile=fastaFile, 
##                         scanBamP=scanBamP, pileupP=pileupP, BPPARAM=BPPARAM)


###################################################
### code chunk number 33: pileupCounts
###################################################
data("myCounts", package="TarSeqQC")


###################################################
### code chunk number 34: pileupCounts
###################################################
head(myCounts)


###################################################
### code chunk number 35: getRegion
###################################################
#complete information for gene7
getRegion(myPanel, level="gene", ID="gene7", collapse=FALSE)
#summarized information for gene7
getRegion(myPanel, level="gene", ID="gene7", collapse=TRUE)


###################################################
### code chunk number 36: plotRegion (eval = FALSE)
###################################################
## g<-plotRegion(myPanel, region=c(4500,6800), seqname="chr10", SNPs=TRUE, 
##     minAAF = 0, minRD = 0, xlab="", title="gene7 amplicons",size=0.5)
## x11(type="cairo")
## g


###################################################
### code chunk number 37: plotFeature (eval = FALSE)
###################################################
## g<-plotFeature(myPanel, featureID="AMPL20")
## x11(type="cairo")
## g


###################################################
### code chunk number 38: plotNtd (eval = FALSE)
###################################################
## plotNtdPercentage(myPanel, featureID="AMPL20")


###################################################
### code chunk number 39: featureInfo
###################################################
getFeaturePanel(myPanel)["AMPL20"]


###################################################
### code chunk number 40: readCountsFeat<
###################################################
featureCounts<-myCounts[myCounts[, "seqnames"] =="chr10" & 
                        myCounts[,"pos"] >= 4866 & myCounts[,"pos"] <= 4928,]


###################################################
### code chunk number 41: readCountsFeat<
###################################################
featureCounts[which.min(featureCounts[,"="]),]


###################################################
### code chunk number 42: buildReport (eval = FALSE)
###################################################
## imageFile<-system.file("extdata", "plot.pdf", package="TarSeqQC",
##                         mustWork=TRUE)
## buildReport(myPanel, attributeThres, imageFile ,file="Results.xlsx")


###################################################
### code chunk number 43: TarSeqQC-vignette.Rnw:1139-1140 (eval = FALSE)
###################################################
## panelList<-list(panel1=te1, panel2=te2)  


###################################################
### code chunk number 44: TarSeqQC-vignette.Rnw:1155-1157 (eval = FALSE)
###################################################
## TEList<-TargetExperimentList(TEList=ampliList, feature="amplicon", 
##     attribute="coverage")


###################################################
### code chunk number 45: TarSeqQC-vignette.Rnw:1160-1161
###################################################
data(TEList, package = "TarSeqQC")


###################################################
### code chunk number 46: TarSeqQC-vignette.Rnw:1167-1171
###################################################
# show/print
TEList
# summary
summary(TEList)


###################################################
### code chunk number 47: TarSeqQC-vignette.Rnw:1185-1187 (eval = FALSE)
###################################################
## x11(type="cairo")
## plotAttrExpl(TEList, dens=TRUE, join=FALSE, log=FALSE, pool=TRUE)


###################################################
### code chunk number 48: TarSeqQC-vignette.Rnw:1209-1211 (eval = FALSE)
###################################################
## x11(type="cairo");
## plotPoolPerformance(TEList, dens=TRUE, join=FALSE, log=FALSE)


###################################################
### code chunk number 49: TarSeqQC-vignette.Rnw:1253-1255 (eval = FALSE)
###################################################
## plot(TEList, attributeThres = attributeThres, sampleLabs = TRUE,
##     featureLabs = TRUE)


###################################################
### code chunk number 50: TarSeqQC-vignette.Rnw:1278-1279
###################################################
getLowCtsFeatures(TEList, level="feature", threshold=50)


###################################################
### code chunk number 51: TarSeqQC-vignette.Rnw:1289-1292 (eval = FALSE)
###################################################
## x11(type="cairo")
## plotGlobalAttrExpl(TEList, attributeThres=attributeThres, dens=FALSE,
##     pool=TRUE, featureLabs=FALSE)


###################################################
### code chunk number 52: filtering unmapped reads (eval = FALSE)
###################################################
## bedFile<-system.file("extdata", "mybed.bed", package="TarSeqQC", mustWork=TRUE)
## fastaFile<-system.file("extdata", "myfasta.fa", package="TarSeqQC",
##     mustWork=TRUE)
## bamFile<-system.file("extdata", "mybam.bam", package="TarSeqQC", mustWork=TRUE)
## flag<-scanBamFlag(isUnmappedQuery = FALSE)
## scanBamP<-ScanBamParam(flag=flag)
## myPanel<-TargetExperiment(bedFile, bamFile, fastaFile, feature="amplicon", 
##                         attribute="coverage", scanBamP=scanBamP,BPPARAM=BPPARAM)
## 


###################################################
### code chunk number 53: loading TarSeqQC example data
###################################################
data(ampliPanel, package="TarSeqQC")
ampliPanel


###################################################
### code chunk number 54: using TarSeqQC example data (eval = FALSE)
###################################################
## 
## setBamFile(ampliPanel)<-system.file("extdata", "mybam.bam", package="TarSeqQC",
##                                     mustWork=TRUE)
## setFastaFile(ampliPanel)<-system.file("extdata", "myfasta.fa", 
##                                     package="TarSeqQC", mustWork=TRUE)


###################################################
### code chunk number 55: buildFeaturePanel using TarSeqQC example data (eval = FALSE)
###################################################
## plotFeature(ampliPanel, featureID="AMPL1")


###################################################
### code chunk number 56: Session Info
###################################################
sessionInfo()