---
title: "Differential expression analysis with diffwrap"
output: html
editor: visual
author: Vidal Fey (vidal.fey@gmail.com)
date: 'Last updated: `r format(Sys.Date(), "%B %d, %Y")`'
format:
  html:
    code-fold: true
    toc: true
    html-math-method: katex
    toc-title: Contents
    number-sections: true
    number-depth: 3
    code-copy: true
vignette: >
  %\VignetteIndexEntry{Differential expression analysis with diffwrap}
  %\VignetteEngine{quarto::html}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r}
#| label: load-packages
#| include: false

library(diffwrap)
```

## Introduction

::: {style="text-align: justify;"}
The ***diffwrap*** package provides functionality for differential expression analysis of read counts from mRNA sequencing data or miRNA expression values generated by the CAP-miRSeq *expression_reports.sh* script. The workflow follows the **edgeR**-**limma** expression data analysis pipeline providing options for different approaches, e.g., using only **edgeR** functions or tools available in the **limma** package such as `voom`. The functions in the package generate text files with differential expression lists, optionally annotated with information from *biomart*, expression summary plots as well as several QC plots.

The main entry point is `diffExpr()`, a convenience wrapper that runs the whole pipeline with sensible defaults. Every step it performs is also available as an exported function, so the workflow can equally well be run manually, or individual steps reused in a different context. This vignette shows both routes.
:::

## Pipeline overview

::: {style="text-align: justify;"}
Pre-processing runs across the top of the diagram below. The analysis diverges into five possible routes at a *Mode resolution* hub \[five colour-coded lanes, one per analysis mode: each lane shows the flags that select it (`pairs`, `block`, `do.voom`) and the resulting recipe (design · contrasts · fit)\]. The routes re-converge at contrast extraction, and the required pipeline ends at the DE tables; QC, per-contrast plots, Venn diagrams and enrichment analyses are optional outputs. The five modes are described in @sec-modes; enrichment (@sec-enrichment) is detailed separately in @fig-enrichment.
:::

![The diffwrap `diffExpr()` workflow. Pre-processing feeds a Mode resolution hub that fans into the five analysis modes (role × engine), each labelled with its selecting flags and its design/contrasts/fit recipe. The modes reconverge at contrast extraction; the required pipeline ends at the DE tables, with QC, plots, Venn and enrichment as optional branches.](diffwrap-pipeline.svg){#fig-pipeline}

## Example data

::: {style="text-align: justify;"}
The package includes a small **simulated** data set so that the workflow can be demonstrated without any external data or network access. Note that it is *not* a real experiment: counts were drawn from a negative binomial model with a known ground truth, which makes it convenient for illustration and testing.

The design is eight samples in two groups of four, `control` and `treated`. Each of the four subjects (`P1`--`P4`) contributed one control and one treated sample, so the same data can be analysed as an unpaired comparison, as a paired design, or as a blocked design. A subject-specific offset is built into the simulation, giving the paired analyses a real effect to remove.

Of the 400 genes, 60 are truly differentially expressed (30 up, 30 down) with absolute log2 fold changes between 1.2 and 3. A further 40 genes are given deliberately low expression so that they trigger the default filtering threshold, and the five `htseq-count` summary rows (`__no_feature` and similar) are retained so that the filtering step has something to remove.
:::

::: {style="text-align: justify;"}
The package ships this data set as the lazy-loaded objects `diffwrap_counts` and `diffwrap_samp_info` (see `?diffwrap_counts`), and as the plain tab-separated files `example_counts.tsv` and `example_samp_info.tsv` under the package's `extdata` directory.

To ensure this vignette is completely self-contained and reproducible, we regenerate an equivalent data set inline from the same simulation model rather than loading the shipped copy:
:::

```{r}
#| 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
```

::: {style="text-align: justify;"}
Real analyses usually begin from a counts file rather than an in-memory matrix. To demonstrate that entry point we write the matrix to a temporary file and read it back with `diff_expr_read_counts()`. `diff_expr_get_samp_info()` first standardises the sample sheet to the conventions the package expects:
:::

```{r}
#| 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])
```

## Choosing an analysis mode {#sec-modes}

::: {style="text-align: justify;"}
Three arguments together determine how an analysis is carried out: `pairs`, `block` and `do.voom`. Rather than acting independently, they resolve into one of a small number of well-defined modes. The distinction that matters most is **how the `pairs` column enters the model**:

- **not used at all** --- a simple means model `~0 + groups`, with contrasts formed explicitly by `makeContrasts()`.
- **as a fixed effect** (`pairs` given, `block = FALSE`) --- an additive model `~pairs + groups` with an intercept. Group coefficients are then already expressed relative to the control group and are extracted directly by coefficient name.
- **as a correlation block** (`block = TRUE`) --- the means model again, but with the subject effect modelled through `limma::duplicateCorrelation()`.

The second axis determines the model used for testing differential expression: **edgeR**'s generalised linear models, or **limma** after a `voom` transformation, selected with `do.voom`.

These two axes give five valid combinations. The only constraint is that a blocked design forces the **limma** engine, because **edgeR** has no equivalent of `duplicateCorrelation()`; if `block = TRUE` is combined with `do.voom = FALSE`, `voom` is enabled automatically and the reason is written to the log.
:::

| `pairs` | `block` | `do.voom`     | Design            | Engine     | Role of `pairs`   |
|------------|------------|------------|------------|------------|------------|
| ---     | `FALSE` | `FALSE`       | `~0 + groups`     | edgeR GLM  | ---               |
| ---     | `FALSE` | `TRUE`        | `~0 + groups`     | limma/voom | ---               |
| given   | `FALSE` | `FALSE`       | `~pairs + groups` | edgeR GLM  | fixed effect      |
| given   | `FALSE` | `TRUE`        | `~pairs + groups` | limma/voom | fixed effect      |
| given   | `TRUE`  | forced `TRUE` | `~0 + groups`     | limma/voom | correlation block |

## Running the whole pipeline

::: {style="text-align: justify;"}
In the simplest case a single call is enough. Note that `out.dir` is **required**: the package deliberately has no default output location, so that nothing is ever written to the working directory unintentionally. Here the results go to a temporary directory, which is also what users should do in their own examples and tests.
:::

```{r}
#| 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)
```

::: {style="text-align: justify;"}
The chunk above and the two concrete full-pipeline examples below are not evaluated when the vignette is built, because a full run writes a fair number of files and is too time-consuming for a package build; it also needs no network access with these settings. Run it interactively to see the complete output. The manual walk-through in the next section *is* evaluated, and shows real differential expression results.

`diffExpr()` returns a list of the objects created along the way --- the `DGEList` or `voom` object, the model fits, and one annotated result table per contrast --- but it is first and foremost called for its side effects. The `out.dir` will contain, per contrast, a full and a filtered differential expression table, a PDF of the M-A, volcano and p-value plots, heatmaps, and a Venn diagram plus intersection tables across contrasts.
:::

### Paired and blocked analyses

```{r}
#| 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)
```

## Running the steps manually

::: {style="text-align: justify;"}
The wrapper is convenient, but each step is exported and can be called on its own. This is useful when you want to inspect intermediate objects, or plug a step into a different workflow.
:::

```{r}
#| 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
```

```{r}
#| 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
```

```{r}
#| 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))
```

::: {style="text-align: justify;"}
Because the ground truth is known, the result can be checked directly: the truly differentially expressed genes are the first sixty, `ENSG00000000001` to `ENSG00000000060`.
:::

```{r}
#| 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")
```

## Reshaping outputs without re-running {#sec-reshape}

::: {style="text-align: justify;"}
Because `diffExpr()` returns every intermediate object, plots and tables can be regenerated with different thresholds, palettes or gene sets **without repeating the analysis**. The returned list `res` holds the model fit(s), one annotated table per contrast in `res$contrasts`, and the plot objects in `res$MAplots`, `res$volcanoPlots` and `res$heatmapPlots`. Each `res$contrasts[[<contrast>]]` is a data frame with the per-sample normalised expression columns, the gene symbol and the `logFC`/`PValue`/`FDR` statistics --- everything the standalone functions need.

Here we assemble that table from the manual objects so the vignette stays self-contained and offline; after a full run you would simply take `de_tab <- res$contrasts[["treated-control"]]`.
:::

```{r}
#| 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
```

::: {style="text-align: justify;"}
Regenerating the heatmap is then a single call to `pheatmap_plots()` --- the same function the wrapper uses --- at a **stricter** FDR cut-off (0.01 vs. the default 0.05) and in a **different** palette, reusing the numbers already computed. `pheatmap_plots()` draws the heatmap together with its gene and sample correlograms (as it does into the pipeline's PDF); here we send those companion panels to a scratch device and keep the heatmap itself:
:::

```{r}
#| 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
```

::: {style="text-align: justify;"}
The same pattern covers the other outputs: `diff_expr_volcano_plot()` and `diff_expr_ma_plot()` take the same `res$contrasts[[<contrast>]]` table with their own `p.thr`/`fdr.thr`/`logfc.thr` and colour arguments, and `diffr_expr_generate_cleaned_de_table_output()` re-writes a trimmed, significance-filtered result table at the cut-offs you choose and without touching the model fit.

Both plot functions also take `base.size`, a single parameter that scales legend, axes, titles, point labels and points together. The pipeline writes its per-contrast plots to a large (15 × 15 inch) PDF, so it passes a correspondingly large default via the `de.plot.base.size` argument of `diffExpr()`. When re-drawing a plot for a smaller figure --- a slide, or a vignette chunk like this one --- lower it:
:::

```{r}
#| 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
```

::: {style="text-align: justify;"}
Every threshold, palette and size argument behaves the same way: change it, re-draw, and the fit is never recomputed.
:::

## Output, logging and verbosity

::: {style="text-align: justify;"}
A pipeline of this size produces a lot of progress information. Rather than writing it all to the console, `diffExpr()` always writes a complete log file to `out.dir` and prints only the major steps to the console. Console output is written to *stdout*, and `verbose` is the only switch that controls it.

- `verbose = TRUE` (default) --- major workflow steps on the console, full detail in the log file.
- `verbose = FALSE` --- nothing on the console; the log file is still written in full.
- `verbose = "all"` --- mirror the complete log to the console as well, which is handy while debugging.

Note that `suppressMessages()` has no effect on this output: stdout is used deliberately, because packages loaded during a run (**VennDiagram** via **futile.logger**, for instance) can disturb the *stderr* stream that `message()` writes to. Use `verbose = FALSE` to silence the run, or capture the output with `capture.output()` or `sink()`. The log file is named after the analysis and is written next to the results; pass `log.file` to save it somewhere else.
:::

```{r}
#| 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
```

## Annotation and enrichment {#sec-enrichment}

::: {style="text-align: justify;"}
Setting `biom.use = TRUE` annotates the result tables with information retrieved from Ensembl, and `enr.do = TRUE` runs over-representation and gene set enrichment analyses through **clusterProfiler**, **gprofiler2** and **topGO**.

Both require network access, and enrichment additionally requires the relevant organism annotation package (`org.Hs.eg.db` for human, `org.Mm.eg.db` for mouse).

Because a run can take a long time, `diffExpr()` does not simply assume that Ensembl is reachable. When `biom.use = TRUE`, it probes the configured host once at the start of the run, before any modelling. Each attempt is given `biom.timeout` seconds (10 by default), so an unreachable server costs seconds rather than a full network timeout. If the probe fails, the run does not stop: it warns, annotates the result tables with offline gene symbols from `convertid::convertId2()` instead, and continues. The same fallback applies if a query fails part-way through a run.

KEGG pathway enrichment is the one part that cannot degrade this way, since it needs Entrez identifiers that only biomart provides. If biomart is off or unreachable, KEGG is therefore dropped from `enr.methods` automatically and the remaining methods run as usual.

Enrichment runs as a separate stage that consumes the differential-expression results: over-representation (ORA) on the significant genes and gene-set enrichment (GSEA) on the ranked gene list, across the methods selected in `enr.methods` (@fig-enrichment).
:::

![Enrichment sub-pipeline. The per-contrast DE results are used in over-representation (ORA) and gene-set enrichment (GSEA), both branches run through the methods chosen in enr.methods. KEGG needs Entrez IDs and therefore biomart annotation. Results are written as tables and optional network plots.](diffwrap-enrichment.svg){#fig-enrichment}

::: {style="text-align: justify;"}
Neither step will produce anything meaningful for the simulated example data, whose gene identifiers are well-formed but fictitious. The chunk below is therefore shown for reference only.
:::

```{r}
#| 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"))
```

## Session information

```{r}
#| label: session-info
sessionInfo()
```
