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

## ----installation, eval = FALSE-----------------------------------------------
# # Install stable version from Bioconductor (once available)
# BiocManager::install("decemedip")

## ----setup, warning=F, message=F----------------------------------------------
library(SummarizedExperiment)
library(dplyr)
library(ggplot2)
library(decemedip)
options(digits = 2)

## ----fig.retina = NULL, fig.align='center', fig.wide = TRUE, echo=FALSE-------
# knitr::include_graphics(system.file("figures/method_main_figure.png", package = "decemedip"), dpi = 300)

## -----------------------------------------------------------------------------
data(hg19.ref.cts.se)
print(hg19.ref.cts.se)
head(granges(hg19.ref.cts.se))

## -----------------------------------------------------------------------------
data(hg19.ref.anc.se)
print(hg19.ref.anc.se)
head(granges(hg19.ref.anc.se))

## ----eval=FALSE---------------------------------------------------------------
# sample_bam_file <- "path/to/bam/files"
# paired <- TRUE # whether the sequencing is paired-end

## ----eval=FALSE---------------------------------------------------------------
# output <- decemedip(
#   sample_bam_file = sample_bam_file,
#   paired = paired,
#   cores = 4
# )

## -----------------------------------------------------------------------------
data(example.hg19.ref.cts.se)
data(example.hg19.ref.anc.se)

data(example.pdx.counts.cts.se)
data(example.pdx.counts.anc.se)

## -----------------------------------------------------------------------------
# read counts of cell type-specific CpGs of the sample 'LuCaP_147CR'
counts_cts <- assay(example.pdx.counts.cts.se)[, "LuCaP_147CR"]
# read counts of anchor CpGs of the sample 'LuCaP_147CR'
counts_anc <- assay(example.pdx.counts.anc.se)[, "LuCaP_147CR"]

## ----eval=TRUE----------------------------------------------------------------
output <- decemedip(
  counts_cts = counts_cts,
  counts_anc = counts_anc,
  ref_cts = example.hg19.ref.cts.se,
  ref_anc = example.hg19.ref.anc.se,
  diagnostics = TRUE,
  cores = 4,
  iter = 500
)

## ----message=FALSE------------------------------------------------------------
library(rstan)

## ----eval=TRUE----------------------------------------------------------------
smr_pi.df <- getSummaryOnPi(output$posterior, cell_type_names = colnames(example.hg19.ref.cts.se))
print(smr_pi.df)

## ----fig.fullwidth=TRUE, warning=FALSE, eval=TRUE-----------------------------
labels <- gsub("_", " ", smr_pi.df$cell_type)
labels <- gsub("(.*) EPIC", "\\1", labels)

smr_pi.df |>
  mutate(cell_type = factor(cell_type, labels = labels)) |>
  ggplot(aes(cell_type, mean)) +
  geom_linerange(aes(ymin = `2.5%`, ymax = `97.5%`),
    position = position_dodge2(width = 0.035),
    linewidth = 7, alpha = 0.3
  ) +
  geom_linerange(aes(ymin = `25%`, ymax = `75%`),
    position = position_dodge2(width = 0.035),
    linewidth = 7, alpha = 1
  ) +
  geom_point(
    position = position_dodge2(width = 0.035),
    fill = "white", shape = 21, size = 8
  ) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

## ----fig.fullwidth=TRUE, warning=FALSE, eval=TRUE-----------------------------
plotDiagnostics(output, plot_type = "y_fit") +
  ylim(0, 300)

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