## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
)

## ----load_package, include = FALSE--------------------------------------------
# Load the package
library(STADyUM)

## ----install, eval = FALSE----------------------------------------------------
# if (!requireNamespace("BiocManager", quietly=TRUE))
#     install.packages("BiocManager")
# BiocManager::install("STADyUM")

## ----estimate_experiment_rates------------------------------------------------
load(system.file("extdata", "granges_for_read_counting_DLD1_chr21.RData", 
package = "STADyUM"))

controlRates <- estimateTranscriptionRates(system.file("extdata",
"PROseq-DLD1-aoi-NELFC_Auxin_Ctrl-SE_plus_chr21.bw", package = "STADyUM"),
system.file("extdata", "PROseq-DLD1-aoi-NELFC_Auxin_Ctrl-SE_minus_chr21.bw",
package = "STADyUM"),
    bw_pause_filtered, bw_gb_filtered, "Control")

show(controlRates)

treatedRates <- estimateTranscriptionRates(system.file("extdata",
"PROseq-DLD1-aoi-NELFC_Auxin-SE_plus_chr21.bw", package = "STADyUM"),
system.file("extdata", "PROseq-DLD1-aoi-NELFC_Auxin-SE_minus_chr21.bw", 
package = "STADyUM"),
    bw_pause_filtered, bw_gb_filtered, "Auxin Treated")

show(treatedRates)

controlRatesTbl <- rates(controlRates)

head(controlRatesTbl)

## ----plot_exp_rates, fig.width=7, fig.height=5, out.width="75%"---------------
plotChiDistrib(controlRates, file="ChiDistributionPlot.png")

plotExpectedVsActualPauseSiteCounts(treatedRates, 
                                    file="treatedExpActPauseSites.png")

## ----estimate_experiment_rates_steric_hindrance, eval=FALSE-------------------
# controlRatesSH <- estimateTranscriptionRates(system.file("extdata",
# "PROseq-DLD1-aoi-NELFC_Auxin_Ctrl-SE_plus_chr21.bw", package = "STADyUM"),
#     system.file("extdata", "PROseq-DLD1-aoi-NELFC_Auxin_Ctrl-SE_minus_chr21.bw",
#     package = "STADyUM"),
#     bw_pause_filtered, bw_gb_filtered, "Control", stericHindrance=TRUE,
#     omegaScale=12.3768278981277)
# 
# show(controlRatesSH)
# 
# treatedRatesSH <- estimateTranscriptionRates(
#     system.file("extdata", "PROseq-DLD1-aoi-NELFC_Auxin-SE_plus_chr21.bw",
#     package = "STADyUM"),
#     system.file("extdata", "PROseq-DLD1-aoi-NELFC_Auxin-SE_minus_chr21.bw",
#     package = "STADyUM"),
#     bw_pause_filtered, bw_gb_filtered, "Auxin Treated", stericHindrance=TRUE,
#     omegaScale=11.0318379571379)
# 
# show(treatedRatesSH)

## ----likelihood_ratio_tests---------------------------------------------------
lrt <- likelihoodRatioTest(controlRates, treatedRates)

show(lrt)

## ----plot_lrt_results,fig.width=7,fig.height=5,out.width="75%"----------------

plotMeanPauseDistrib(controlRates, file="meanPauseSitesControl.png")
plotMeanPauseDistrib(treatedRates, file="meanPauseSitesTreated.png")

BetaViolinPlot(lrt, file="lrtBetaViolinPlot.png")


## ----create_simpol, fig.width=7,fig.height=5,out.width="70%",out.height="50%"----
simpol <- simulatePolymerase(k=50, ksd=25, kMin=17, kMax=200, geneLen=1950,
alpha=2, beta=0.5, zeta=2000, zetaSd=1000, zetaMin=1500, zetaMax=2500,
cellNum=100, polSize=33, addSpace=17, time=39, timesToRecord=NULL)

show(simpol)

readCounts <- readCounts(simpol)

plotCombinedCells(simpol, file="combinedCellsLollipopPlot.png")

## ----estimate_simulated_rates, fig.width=7, fig.height=5, out.width="75%"-----
simRates <- estimateTranscriptionRates(simpol, name="simRatesb0.5k50",
stericHindrance=TRUE)

show(simRates)

## ----simulated_lrts, fig.width=7, fig.height=5, out.width="75%", eval=FALSE----
# 
# simpol2 <- simulatePolymerase(k=100, ksd=25, kMin=75, kMax=200, geneLen=1950,
# alpha=2, beta=10, zeta=2000, zetaSd=1000, zetaMin=1500, zetaMax=2500,
# cellNum=100, polSize=33, addSpace=17, time=30, timesToRecord=NULL)
# 
# simRates2 <- estimateTranscriptionRates(simpol2, name="simRatesb10k100",
# stericHindrance=TRUE)
# 
# simLRT <- likelihoodRatioTest(simRates, simRates2)
# 
# plotPauseSiteContourMapTwoConditions(simLRT, file="simPauseSitesContourMap.png")
# BetaViolinPlot(simLRT, file="simBetaViolinPlot.png")
# 

## ----session_info-------------------------------------------------------------
sessionInfo()