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 <-60Lsamples <-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 tailbase_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 downlfc <-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 removesubj_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^lfcrnbinom(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 stripspecial <-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 8head(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 3936samp_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:
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.
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 subjectres_paired <-diffExpr(expr.dat = counts_file,samp.info = samp_info_raw,samples ="SampleName",groups ="Group",pairs ="Subject", # fixed effect, intercept designcontrol ="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 variableblock =TRUE, # duplicateCorrelationcontrol ="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 normLibSizesdesign <-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
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 runde_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 correlogramshm <-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 PDFv$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.
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.
