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

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


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


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


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


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


###################################################
### code chunk number 6: 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 7: setters
###################################################
# set feature slot value
setFeature(myPanel)<-"amplicon"
# set attribute slot value
setAttribute(myPanel)<-"coverage"


###################################################
### code chunk number 8: setters
###################################################
# set scanBamP slot value
scanBamP<-ScanBamParam()
#set which slot
bamWhich(scanBamP)<-getBedFile(myPanel)
setScanBamP(myPanel)<-scanBamP
# set attribute 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 9: loading TargetExperiment object
###################################################
data(ampliPanel, package="TarSeqQC")


###################################################
### code chunk number 10: 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 11: exploration
###################################################
# show/print
myPanel

# summary
summary(myPanel)

#summary at feature level
summaryFeatureLev(myPanel)

#summary at gene level
summaryGeneLev(myPanel)


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


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


###################################################
### code chunk number 14: plot (eval = FALSE)
###################################################
## 
## # plot panel overview
## g<-plot(myPanel, attributeThres, chrLabels =TRUE)
## g


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


###################################################
### code chunk number 16: frequency table
###################################################

# summaryIntervals
summaryIntervals(myPanel, attributeThres)


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


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


###################################################
### code chunk number 19: 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 20: pileupCounts
###################################################
# 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)
head(myCounts)


###################################################
### code chunk number 21: 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 22: plotRegion (eval = FALSE)
###################################################
## g<-plotRegion(myPanel, region=c(4500,6800), seqname="chr10", SNPs=TRUE,  
##     xlab="", title="gene7 amplicons",size=0.5)
## x11(type="cairo")
## g


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


###################################################
### code chunk number 24: plotNtd (eval = FALSE)
###################################################
## g<-plotNtdPercentage(myPanel, featureID="AMPL20")
## g


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


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


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


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


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


###################################################
### code chunk number 30: 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 31: buildFeaturePanel using TarSeqQC example data (eval = FALSE)
###################################################
## setFeaturePanel<-buildFeaturePanel(ampliPanel)


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