## 1 Introduction

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.

## 2 Pipeline overview

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 Section 4; enrichment (Section 9) is detailed separately in Figure 2.

Figure 1: 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.

## 3 Example data

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.

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:

Code

    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)
    #> [1] 405   8
    head(counts_mat[, 1:4])
    #>                  S01  S02  S03  S04
    #> ENSG00000000001 3589 2202 3057 2013
    #> ENSG00000000002 4531 2728 2530 2586
    #> ENSG00000000003 2679 1827 1491 1229
    #> ENSG00000000004  188  206  110  290
    #> ENSG00000000005 8871 9745 6369 5543
    #> ENSG00000000006 3571 4441 3389 3936
    samp_info_raw
    #>   SampleName   Group Subject
    #> 1        S01 control      P1
    #> 2        S02 control      P2
    #> 3        S03 control      P3
    #> 4        S04 control      P4
    #> 5        S05 treated      P1
    #> 6        S06 treated      P2
    #> 7        S07 treated      P3
    #> 8        S08 treated      P4

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:

Code

    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
    #> [1] 405   8
    head(counts[, 1:4])
    #>                  S01  S02  S03  S04
    #> ENSG00000000001 3589 2202 3057 2013
    #> ENSG00000000002 4531 2728 2530 2586
    #> ENSG00000000003 2679 1827 1491 1229
    #> ENSG00000000004  188  206  110  290
    #> ENSG00000000005 8871 9745 6369 5543
    #> ENSG00000000006 3571 4441 3389 3936

## 4 Choosing an analysis mode

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

## 5 Running the whole pipeline

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.

Code

    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)

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.

### 5.1 Paired and blocked analyses

Code

    # '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)

## 6 Running the steps manually

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.

Code

    counts_f <- diff_expr_filter_counts(counts, samp.info, strict = TRUE)
    dim(counts_f)   # low-expression genes and the __-rows have gone
    #> [1] 367   8

Code

    groups <- stats::relevel(samp.info$Groups, ref = "control")

    d <- edgeR::DGEList(counts = counts_f, group = groups)
    d <- edgeR::calcNormFactors(d)
    #> calcNormFactors has been renamed to normLibSizes

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

    design
    #>   control treated
    #> 1       1       0
    #> 2       1       0
    #> 3       1       0
    #> 4       1       0
    #> 5       0       1
    #> 6       0       1
    #> 7       0       1
    #> 8       0       1
    #> attr(,"assign")
    #> [1] 1 1
    #> attr(,"contrasts")
    #> attr(,"contrasts")$groups
    #> [1] "contr.treatment"
    contrasts
    #>          Contrasts
    #> Levels    treated-control
    #>   control              -1
    #>   treated               1

Code

    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))
    #> Coefficient:  -1*control 1*treated
    #>                     logFC    logCPM        F       PValue          FDR
    #> ENSG00000000007  3.582967  9.386085 64.45008 1.595131e-15 5.854132e-13
    #> ENSG00000000053  3.507002 14.836965 62.13925 4.989322e-15 9.155405e-13
    #> ENSG00000000018 -3.581821  7.194548 59.54555 1.798209e-14 2.199809e-12
    #> ENSG00000000027  3.227979 10.799305 56.32489 8.864424e-14 8.133109e-12
    #> ENSG00000000038 -3.117599 13.557293 51.49895 9.746956e-13 6.591007e-11
    #> ENSG00000000011  3.093844 12.804865 51.29741 1.077549e-12 6.591007e-11

Because the ground truth is known, the result can be checked directly: the truly differentially expressed genes are the first sixty, `ENSG00000000001` to `ENSG00000000060`.

Code

    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")
    #> Of the 50 top-ranked genes, 50 are truly differentially expressed

## 7 Reshaping outputs without re-running

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"]]`.

Code

    ## 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

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:

Code

    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"

Code

    invisible(grDevices::dev.off())

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

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:

Code

    v <- diff_expr_volcano_plot(de_tab, id = "ID", base.size = 11)   # 16 suits the 15-inch PDF
    v$FDR

Every threshold, palette and size argument behaves the same way: change it, re-draw, and the fit is never recomputed.

## 8 Output, logging and verbosity

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.

Code

    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

## 9 Annotation and enrichment

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` (Figure 2).

Figure 2: 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.

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.

Code

    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"))

## 10 Session information

Code

    sessionInfo()
    #> R version 4.6.0 (2026-04-24)
    #> Platform: aarch64-apple-darwin23
    #> Running under: macOS Tahoe 26.5.2
    #>
    #> Matrix products: default
    #> BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib
    #> LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
    #>
    #> locale:
    #> [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
    #>
    #> time zone: Europe/Helsinki
    #> tzcode source: internal
    #>
    #> attached base packages:
    #> [1] stats     graphics  grDevices utils     datasets  methods   base
    #>
    #> other attached packages:
    #> [1] diffwrap_0.6-3
    #>
    #> loaded via a namespace (and not attached):
    #>   [1] DBI_1.3.0            gridExtra_2.3.1      httr2_1.3.0
    #>   [4] formatR_1.14         biomaRt_2.68.0       rlang_1.3.0
    #>   [7] magrittr_2.0.5       otel_0.2.0           compiler_4.6.0
    #>  [10] RSQLite_3.53.3       venn_1.13            png_0.1-9
    #>  [13] vctrs_0.7.3          stringr_1.6.0        pkgconfig_2.0.3
    #>  [16] crayon_1.5.3         fastmap_1.2.0        backports_1.5.1
    #>  [19] dbplyr_2.6.0         XVector_0.52.0       rmarkdown_2.31
    #>  [22] purrr_1.2.2          bit_4.6.0            xfun_0.60
    #>  [25] cachem_1.1.0         jsonlite_2.0.0       progress_1.2.3
    #>  [28] blob_1.3.0           parallel_4.6.0       prettyunits_1.2.0
    #>  [31] cluster_2.1.8.2      R6_2.6.1             stringi_1.8.7
    #>  [34] RColorBrewer_1.1-3   limma_3.68.4         rpart_4.1.27
    #>  [37] Rcpp_1.1.2           Seqinfo_1.2.0        assertthat_0.2.1
    #>  [40] knitr_1.51           base64enc_0.1-6      IRanges_2.46.0
    #>  [43] VennDiagram_1.8.2    polycor_0.8-2        Matrix_1.7-5
    #>  [46] splines_4.6.0        nnet_7.3-20          tidyselect_1.2.1
    #>  [49] rstudioapi_0.19.0    dichromat_2.0-1      yaml_2.3.12
    #>  [52] curl_7.1.0           admisc_0.40          lattice_0.22-9
    #>  [55] tibble_3.3.1         plyr_1.8.9           Biobase_2.72.0
    #>  [58] KEGGREST_1.52.2      S7_0.2.2             evaluate_1.0.5
    #>  [61] foreign_0.8-91       lambda.r_1.2.4       survival_3.8-9
    #>  [64] futile.logger_1.4.9  BiocFileCache_3.2.0  zip_3.0.1
    #>  [67] xml2_1.6.0           Biostrings_2.80.1    pillar_1.11.1
    #>  [70] filelock_1.0.3       checkmate_2.3.4      stats4_4.6.0
    #>  [73] generics_0.1.4       S4Vectors_0.50.1     hms_1.1.4
    #>  [76] ggplot2_4.0.3        scales_1.4.0         glue_1.8.1
    #>  [79] pheatmap_1.0.13      Hmisc_5.2-6          tools_4.6.0
    #>  [82] data.table_1.18.4    openxlsx_4.2.8.1     locfit_1.5-9.12
    #>  [85] mvtnorm_1.4-2        grid_4.6.0           ltm_1.2-0
    #>  [88] AnnotationDbi_1.74.0 edgeR_4.10.1         colorspace_2.1-3
    #>  [91] htmlTable_2.5.0      Formula_1.2-5        convertid_0.3.4
    #>  [94] cli_3.6.6            rappdirs_0.3.4       futile.options_1.0.1
    #>  [97] expm_1.0-0           dplyr_1.2.1          gtable_0.3.6
    #> [100] digest_0.6.39        BiocGenerics_0.58.1  ggrepel_0.9.8
    #> [103] msm_1.8.2            htmlwidgets_1.6.4    farver_2.1.2
    #> [106] memoise_2.0.1        htmltools_0.5.9      lifecycle_1.0.5
    #> [109] httr_1.4.8           statmod_1.5.2        bit64_4.8.2
    #> [112] MASS_7.3-66