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


## -----------------------------------------------------------------------------
#| label: load-packages
#| include: false

library(diffwrap)


## -----------------------------------------------------------------------------
#| label: example-data
set.seed(20240721)

n_genes <- 400L; n_true_de <- 60L
samples <- sprintf("S%02d", 1:8)
group   <- rep(c("control", "treated"), each = 4)
subject <- rep(c("P1", "P2", "P3", "P4"), times = 2)   # each subject: one of each group

## baseline mean expression per gene, with a low-expression tail
base_mu <- exp(rnorm(n_genes, mean = 7.5, sd = 1.5))
base_mu[(n_genes - 39):n_genes] <- runif(40, 0, 14)

## true log2 fold changes: first 60 genes differentially expressed, half up half down
lfc <- numeric(n_genes)
lfc[seq_len(n_true_de)] <- rep(c(1, -1), length.out = n_true_de) * runif(n_true_de, 1.2, 3)

## per-subject offset, so the paired/blocked analyses have a real effect to remove
subj_off <- stats::setNames(runif(4, 0.85, 1.18), c("P1", "P2", "P3", "P4"))

counts_mat <- sapply(seq_along(samples), function(j) {
  mu <- base_mu * subj_off[subject[j]]
  if (group[j] == "treated") mu <- mu * 2^lfc
  rnbinom(n_genes, mu = mu, size = 1 / 0.15)     # negative binomial, dispersion 0.15
})
dimnames(counts_mat) <- list(sprintf("ENSG%011d", seq_len(n_genes)), samples)

## the htseq-count summary rows that filtering is expected to strip
special <- c("__no_feature", "__ambiguous", "__too_low_aQual",
             "__not_aligned", "__alignment_not_unique")
counts_mat <- rbind(counts_mat,
                    matrix(sample(20000:90000, length(special) * 8, replace = TRUE),
                           nrow = length(special), dimnames = list(special, samples)))

samp_info_raw <- data.frame(SampleName = samples, Group = group,
                            Subject = subject, stringsAsFactors = FALSE)

dim(counts_mat)
head(counts_mat[, 1:4])
samp_info_raw


## -----------------------------------------------------------------------------
#| label: read-from-file
counts_file <- file.path(tempdir(), "example_counts.tsv")
write.table(counts_mat, counts_file, sep = "\t", quote = FALSE, col.names = NA)

samp.info <- diff_expr_get_samp_info(samp_info_raw,
                                     samples = "SampleName",
                                     groups  = "Group")
counts <- diff_expr_read_counts(counts_file, samp.info)
dim(counts)          # the htseq-count summary rows are still present here
head(counts[, 1:4])


## -----------------------------------------------------------------------------
#| label: full-run
#| eval: false
# out.dir <- file.path(tempdir(), "diffwrap_example")
# dir.create(out.dir, showWarnings = FALSE)
# 
# res <- diffExpr(expr.dat  = counts_file,
#                 samp.info = samp_info_raw,
#                 samples   = "SampleName",
#                 groups    = "Group",
#                 control   = "control",
#                 analysis.name = "demo",
#                 out.dir   = out.dir,
#                 enr.do = FALSE)


## -----------------------------------------------------------------------------
#| label: paired-run
#| eval: false
# # 'Subject' pairs each control sample with a treated sample from the same subject
# res_paired <- diffExpr(expr.dat  = counts_file,
#                        samp.info = samp_info_raw,
#                        samples   = "SampleName",
#                        groups    = "Group",
#                        pairs     = "Subject",     # fixed effect, intercept design
#                        control   = "control",
#                        analysis.name = "demo_paired",
#                        out.dir   = out.dir,
#                        enr.do = FALSE)
# 
# # the same column used as a correlation block instead (forces voom)
# res_blocked <- diffExpr(expr.dat  = counts_file,
#                         samp.info = samp_info_raw,
#                         samples   = "SampleName",
#                         groups    = "Group",
#                         pairs     = "Subject",    # 'pairs' is the blocking variable
#                         block     = TRUE,         # duplicateCorrelation
#                         control   = "control",
#                         analysis.name = "demo_blocked",
#                         out.dir   = out.dir,
#                         enr.do = FALSE)


## -----------------------------------------------------------------------------
#| label: manual-filter
counts_f <- diff_expr_filter_counts(counts, samp.info, strict = TRUE)
dim(counts_f)   # low-expression genes and the __-rows have gone


## -----------------------------------------------------------------------------
#| label: manual-model
groups <- stats::relevel(samp.info$Groups, ref = "control")

d <- edgeR::DGEList(counts = counts_f, group = groups)
d <- edgeR::calcNormFactors(d)

design    <- diff_expr_make_design(samp.info = samp.info, groups = groups)
contrasts <- diff_expr_make_contrasts(design = design, groups = groups)

design
contrasts


## -----------------------------------------------------------------------------
#| label: manual-fit
fit.l <- diff_expr_fit(counts = counts_f, d = d, design = design,
                       do.voom = FALSE, quasi.likelihood = TRUE)

de <- edgeR::glmQLFTest(fit.l$fit, contrast = contrasts[, 1])
head(edgeR::topTags(de, n = 10))


## -----------------------------------------------------------------------------
#| label: check-truth
tt      <- as.data.frame(edgeR::topTags(de, n = Inf))
top50   <- rownames(tt)[seq_len(50)]
n_true  <- sum(as.integer(sub("^ENSG0*", "", top50)) <= 60)
cat("Of the 50 top-ranked genes,", n_true, "are truly differentially expressed\n")


## -----------------------------------------------------------------------------
#| label: reshape-table
## identical in structure to res$contrasts[["treated-control"]] from a full run
de_tab <- merge(edgeR::cpm(fit.l$d2, log = TRUE), tt, by = "row.names")
names(de_tab)[1]   <- "ID"
de_tab$gene_symbol <- de_tab$ID       # offline demo; a real run already carries symbols


## -----------------------------------------------------------------------------
#| label: reshape-heatmap
#| fig-width: 7
#| fig-height: 8
grDevices::pdf(tempfile(fileext = ".pdf"))          # swallow the companion correlograms
hm <- pheatmap_plots(de_tab, id = "ID",
                     samp.info = samp.info, samples = "SampleNames", groups = groups,
                     fdr.thr = 0.01, topn = 30,
                     color.blind.pal = "RdBu")       # default palette is "PuOr"
invisible(grDevices::dev.off())

hm$fdr$regular                                       # the FDR-filtered heatmap, new settings


## -----------------------------------------------------------------------------
#| label: reshape-volcano
#| eval: false
# v <- diff_expr_volcano_plot(de_tab, id = "ID", base.size = 11)   # 16 suits the 15-inch PDF
# v$FDR


## -----------------------------------------------------------------------------
#| label: logging
#| eval: false
# res <- diffExpr(..., out.dir = out.dir, verbose = FALSE)          # quiet console
# res <- diffExpr(..., out.dir = out.dir, verbose = "all")          # everything echoed
# res <- diffExpr(..., out.dir = out.dir,
#                 log.file = file.path(out.dir, "my_run.log"))      # custom log location


## -----------------------------------------------------------------------------
#| label: enrichment
#| eval: false
# res <- diffExpr(expr.dat  = counts_file,
#                 samp.info = samp_info_raw,
#                 samples   = "SampleName",
#                 groups    = "Group",
#                 control   = "control",
#                 analysis.name = "demo_annotated",
#                 out.dir   = out.dir,
#                 biom.use  = TRUE,
#                 biom.data.set = "hsapiens_gene_ensembl",
#                 enr.do = TRUE,
#                 enr.methods = c("clusterProfilerGO", "gProfileR"))


## -----------------------------------------------------------------------------
#| label: session-info
sessionInfo()