## ----setup, message=FALSE, warning=FALSE--------------------------------------
has_connectivity <- requireNamespace("BiocFileCache", quietly = TRUE) &&
    curl::has_internet()
knitr::opts_chunk$set(dpi = 72, fig.width = 5, fig.height = 4,
                      eval = has_connectivity)
if (!has_connectivity)
    message("No internet connection; vignette code not evaluated.")

library(STADyUM)           # For transcription rate estimation
library(GenomicRanges)     # For genomic data manipulation
library(dplyr)
library(plyranges)
#library(tidyverse)         # For data manipulation and visualization
library(ggplot2)           # For plotting
library(pracma)
library(BiocFileCache)

## ----file-paths---------------------------------------------------------------

bfc <- BiocFileCache(ask = FALSE)
zip_path <- bfcrpath(bfc, "https://zenodo.org/records/20618059/files/lrt_vignette_data.zip")
unzip(zip_path, exdir = tempdir())
data_dir <- file.path(tempdir(), "vignettes", "lrt_test_data")

# STADyUM `ExperimentTranscriptionRates` objects containing transcription rates estimated from human CD4+ T cell and human CD14+ monocyte PRO-seq data.
cd4_rate <- readRDS(file.path(data_dir, 'cd4_rate.RDS'))
cd14_rate <- readRDS(file.path(data_dir, 'cd14_rate.RDS'))
scale_factor <- 38359083 / 33294862


## ----lrt_calculation----------------------------------------------------------

lrt <- likelihoodRatioTest(cd4_rate, cd14_rate, scale_factor)

# The merged per-gene LRT summary table is built automatically and stored in
# the lrtTbl slot, ready to be passed straight to the plotting methods below.
head(lrtTbl(lrt))

## ----beta_heatmap-------------------------------------------------------------
# Each plotting method takes the TranscriptionRatesLRT object directly and
# reads the merged table from its lrtTbl slot internally.
p0 <- plotBetaQuantileHeatmap(lrt)

print(p0)

## ----beta_scatter-------------------------------------------------------------
p1 <- plotBetaScatter(lrt,
                      label_x = -6,
                      label_y = -1,
                      x_lim   = c(-6, -1),
                      y_lim   = c(-6, -1))

print(p1)

## ----delta_beta_sigma---------------------------------------------------------
p2 <- plotDeltaBetaSigma(lrt, filter_lfc = TRUE)

print(p2)

## ----fksd_density-------------------------------------------------------------
p3 <- plotFksdDensity(lrt)

print(p3)

## ----session-info-------------------------------------------------------------
sessionInfo()