### R code from vignette source 'InPAS.Rnw'

###################################################
### code chunk number 1: style
###################################################
BiocStyle::latex()


###################################################
### code chunk number 2: loadLibrary
###################################################
library(InPAS)
library(BSgenome.Mmusculus.UCSC.mm10)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)
path <- file.path(find.package("InPAS"), "extdata")


###################################################
### code chunk number 3: prepareAnno
###################################################
library(org.Hs.eg.db)
samplefile <- system.file("extdata", "hg19_knownGene_sample.sqlite",
            package="GenomicFeatures")
txdb <- loadDb(samplefile)
utr3.sample.anno <- utr3Annotation(txdb=txdb, 
                   orgDbSYMBOL="org.Hs.egSYMBOL")
utr3.sample.anno


###################################################
### code chunk number 4: loadAnno
###################################################
##step1 annotation
# load from dataset
data(utr3.mm10)


###################################################
### code chunk number 5: step2HugeDataF
###################################################
load(file.path(path, "polyA.rds"))
library(cleanUpdTSeq)
data(classifier)
bedgraphs <- c(file.path(path, "Baf3.extract.bedgraph"), 
               file.path(path, "UM15.extract.bedgraph"))
hugeData <- FALSE
##step1 Calculate coverage
coverage <- coverageFromBedGraph(bedgraphs, 
                                 tags=c("Baf3", "UM15"),
                                 genome=BSgenome.Mmusculus.UCSC.mm10,
                                 hugeData=hugeData)

##step2 Predict cleavage sites
CPsites <- CPsites(coverage=coverage, 
                   gp1="Baf3", gp2="UM15", 
                   genome=BSgenome.Mmusculus.UCSC.mm10,
                   utr3=utr3.mm10, 
                   search_point_START=200, 
                   cutEnd=.2, 
                   long_coverage_threshold=3,
                   PolyA_PWM=pwm, 
                   classifier=classifier,
                   shift_range=100)

##step3 Estimate 3UTR usage
res <- utr3UsageEstimation(CPsites=CPsites, 
                           coverage=coverage, 
                           utr3=utr3.mm10,
                           genome=BSgenome.Mmusculus.UCSC.mm10,
                           gp1="Baf3", gp2="UM15", 
                           short_coverage_threshold=15, 
                           long_coverage_threshold=3, 
                           adjusted.P_val.cutoff=0.05, 
                           dPDUI_cutoff=0.3, 
                           PDUI_logFC_cutoff=0.59)
res[1]


###################################################
### code chunk number 6: OneStep
###################################################
if(interactive()){
    res <- inPAS(bedgraphs=bedgraphs, tags=c("Baf3", "UM15"), 
              genome=BSgenome.Mmusculus.UCSC.mm10, 
              utr3=utr3.mm10, gp1="Baf3", gp2="UM15",
              txdb=TxDb.Mmusculus.UCSC.mm10.knownGene,
              search_point_START=200,
              short_coverage_threshold=15,
              long_coverage_threshold=3, 
              cutStart=0, cutEnd=.2, 
              hugeData=FALSE,
              shift_range=50,
              PolyA_PWM=pwm, classifier=classifier)
}


###################################################
### code chunk number 7: SingleGroupData
###################################################
res <- inPAS(bedgraphs=bedgraphs[1], tags=c("Baf3"), 
              genome=BSgenome.Mmusculus.UCSC.mm10, 
              utr3=utr3.mm10, gp1="Baf3", gp2=NULL,
              txdb=TxDb.Mmusculus.UCSC.mm10.knownGene,
              search_point_START=200,
              short_coverage_threshold=15,
              long_coverage_threshold=3, 
              cutStart=0, cutEnd=.2, 
              hugeData=FALSE, 
              PolyA_PWM=pwm, classifier=classifier)
res[1]


###################################################
### code chunk number 8: sessionInfo
###################################################
toLatex(sessionInfo())