RNA-Parallel: identical results, measured speedup

Parallel companions for ComBat-seq, limma and edgeR

rnaparallel runs the original function itself and parallelises its hot paths. ComBat_seq_parallel corrects the counts; calcNormFactors_parallel, lmFit_parallel, and duplicateCorrelation_parallel run the differential expression that follows. None of them reimplements a method from its description: the vendor is called unchanged on blocks of genes or samples, with at most a symbol rebound in a child of the vendor’s own environment.

The question is narrow on purpose: does each companion return exactly what the original returns, and how much time does it save. Not close, not correlated at 0.999. Identical, or it is not usable.

Two datasets do two jobs. Simulated counts, generated with the ComBat-seq paper’s parameters, establish that the correction behaves correctly against known truth. TCGA supplies the scale at which parallelism is worth measuring.

citation("rnaparallel") prints the methods to cite alongside this package.

section	what it shows
B. The claim in ten lines	one pipeline, run both ways, identical()
C. Simulated counts	ComBat-seq then limma and edgeR, each against its original, at 2, 4 and 8 workers
D. TCGA	the cohort, then the same two comparisons at cohort scale
E. Equivalence summary	every quantity, computed both ways and compared
F. Speed summary	every speed measurement in one figure

A. Setup

if (!requireNamespace("remotes", quietly = TRUE)) install.packages("remotes")
remotes::install_github("GenomeRx/RNA-Parallel")

if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("sva", "edgeR", "limma", "statmod", "SummarizedExperiment", "TCGAbiolinks",
                       "ComplexHeatmap"))

cran <- c("ggplot2", "ggrepel", "patchwork", "callr", "gtsummary", "circlize", "RColorBrewer")
need <- cran[!vapply(cran, requireNamespace, logical(1), quietly = TRUE)]
if (length(need)) install.packages(need)

library(sva)

## Warning: package 'BiocParallel' was built under R version 4.4.3

library(edgeR)
library(limma)
library(rnaparallel)
library(SummarizedExperiment)
library(ggplot2)

theme_set(
  theme_minimal(base_size = 15) +
    theme(plot.title    = element_text(face = "bold", size = 18, margin = margin(b = 4)),
          plot.subtitle = element_text(size = 13, colour = "grey30", margin = margin(b = 10)),
          axis.title    = element_text(size = 14),
          axis.text     = element_text(size = 12),
          strip.text    = element_text(face = "bold", size = 14),
          legend.title  = element_text(size = 13),
          legend.text   = element_text(size = 12),
          plot.margin   = margin(12, 16, 10, 12)))

# workers is run as asked: 8 means 8 forked processes, bounded only by detectCores(). Added
# workers return less throughput each, so whether 8 beats 4 is what the sweep measures rather
# than assumes. Not a core-type effect: forks migrate across the performance and efficiency
# clusters, so no worker is stranded on a slow one.
# TCGA is corrected at every count; the last arm is the one the figures below use.
tcga_workers <- c(2L, 4L, 8L)
n_workers    <- max(tcga_workers)

B. The claim in ten lines

Before the two datasets, the smallest possible version of the whole argument. The counts are the ones the ComBat-seq README itself uses for its example, so the call below is the upstream call with one argument added.

set.seed(1)
demo_counts <- matrix(rnbinom(400, size = 10, prob = 0.1), nrow = 50, ncol = 8)
demo_batch  <- c(rep(1, 4), rep(2, 4))
demo_group  <- rep(c(0, 1), 4)

demo_ref <- sva::ComBat_seq(demo_counts, batch = demo_batch, group = demo_group)

## Found 2 batches
## Using full model in ComBat-seq.
## Adjusting for 1 covariate(s) or covariate level(s)
## Estimating dispersions
## Fitting the GLM model
## Shrinkage off - using GLM estimates for parameters
## Adjusting the data

demo_par <- ComBat_seq_parallel(demo_counts, batch = demo_batch, group = demo_group,
                                workers = 2L)

## Found 2 batches
## Using full model in ComBat-seq.
## Adjusting for 1 covariate(s) or covariate level(s)
## Estimating dispersions
## Fitting the GLM model
## Shrinkage off - using GLM estimates for parameters
## Adjusting the data

data.frame(genes = nrow(demo_counts), samples = ncol(demo_counts),
           batches = length(unique(demo_batch)),
           identical_to_original = identical(demo_ref, demo_par),
           differing_cells = sum(demo_ref != demo_par))

genes	samples	batches	identical_to_original	differing_cells
50	8	2	TRUE	0

stopifnot(identical(demo_ref, demo_par))

Fifty genes by eight samples sits far below combat.min.cells, so every row-split dispatch runs serially: a fork would cost more than the work. Only the common dispersion is handed out, one whole-matrix estimate per batch, which carries no size gate because it is always worth dispatching. Either way the output is the same matrix, which is the only thing the gate is allowed to leave alone.

C. Simulated counts

Counts keep the ComBat-seq paper’s negative binomial draw and its batch mean and dispersion parameters. That pipeline called polyester::simulate_experiment, whose count step is rnbinom(mu = reads_per_transcript * fold_change, size = size) with zero draws replaced by 1. The draw is reproduced directly here because polyester is no longer on current Bioconductor. The mean structure is this notebook’s own: the paper’s two-batch design cannot show what ten batches look like in a projection.

set.seed(123)

G          <- 10000L    # genes
n_batch    <- 10L       # sequencing batches
n_per_cell <- 50L       # per batch-condition cell, so 1,000 patients in total
bio_fold   <- 1.8       # biological signal
batch_fold <- 3         # paper: mean of batch 2 is 1.5, 2 or 3 times batch 1
disp_fold  <- 4         # paper: dispersion of batch 2 is 2, 3 or 4 fold batch 1
disp_1     <- 0.15      # paper: dispersion in batch 1

batch <- factor(rep(sprintf("B%02d", seq_len(n_batch)), each = 2 * n_per_cell),
                levels = sprintf("B%02d", seq_len(n_batch)))
group <- factor(rep(rep(c("control", "treated"), each = n_per_cell), n_batch))
n     <- length(batch)

base <- pmax(5, round(rlnorm(G, log(120), 0.9)))   # lognormal, the shape a transcriptome has

# Technical batch effects concentrate on a few axes in real data: library depth, degradation,
# GC bias. Two dominant technical directions, with the ten batches scattered across
# that plane, give ten separated clouds. An independent shift per gene per batch instead
# spreads the effect over nine dimensions, where no two-component projection can separate it.
t_load <- matrix(rnorm(G * 2L, 0, 0.22), G, 2L)

# Batch positions are scattered rather than placed on a ring: evenly spaced angles at one
# shared radius produce a visibly drawn circle. Candidates are drawn at random and rejected
# when they land on an existing batch, which keeps them separable without regimenting them.
b_pos  <- matrix(NA_real_, n_batch, 2L)
placed <- 0L
while (placed < n_batch) {
  cand <- rnorm(2L, 0, 2.3)
  far  <- placed == 0L ||
    min(sqrt(rowSums((b_pos[seq_len(placed), , drop = FALSE] -
                        rep(cand, each = placed))^2))) > 1.9
  if (far) { placed <- placed + 1L; b_pos[placed, ] <- cand }
}

# batches also differ in how tightly they cluster, so the clouds are not ten copies of one
# another
b_spread <- runif(n_batch, 0.32, 0.62)
jit <- matrix(rnorm(n * 2L), n, 2L) * b_spread[rep(seq_len(n_batch), each = 2L * n_per_cell)]

# a smaller per-gene per-batch component, so the effect is not exactly rank two
b_idio <- matrix(rnorm(G * n_batch, 0, 0.20), G, n_batch)
b_disp <- runif(n_batch, 1, disp_fold)

# Continuous patient heterogeneity, three latent programmes. Without it every sample in a cell
# shares one mean vector, the only spread left is sampling noise, and that noise averages away
# across thousands of genes: batches collapse to points and the condition becomes a straight
# line once batch is removed.
n_lat     <- 3L
lat_load  <- matrix(rnorm(G * n_lat, 0, 0.18), G, n_lat)
lat_score <- matrix(rnorm(n * n_lat), n, n_lat)
lib       <- rlnorm(n, 0, 0.25)                    # sequencing depth varies per sample

n_de     <- G %/% 20L                              # 5% of genes truly differential
de_idx   <- sample.int(G, n_de)
ups      <- de_idx[seq_len(n_de %/% 2)]; downs <- setdiff(de_idx, ups)
de_shift <- numeric(G); de_shift[ups] <- log(bio_fold); de_shift[downs] <- -log(bio_fold)

# one independent draw per sample, not one per cell
bi  <- as.integer(batch); gi <- as.integer(group)
sim <- vapply(seq_len(n), function(j) {
  lmu <- log(base) + as.vector(t_load %*% (b_pos[bi[j], ] + jit[j, ])) + b_idio[, bi[j]] +
         as.vector(lat_load %*% lat_score[j, ]) + (gi[j] == 2L) * de_shift + log(lib[j])
  x <- rnbinom(G, mu = exp(lmu), size = 1 / (disp_1 * b_disp[bi[j]]))
  x[x == 0] <- 1L                                  # polyester replaces zero draws with 1
  x
}, numeric(G))
storage.mode(sim) <- "integer"
rownames(sim) <- sprintf("gene%05d", seq_len(G))
colnames(sim) <- sprintf("s%04d", seq_len(n))

truth <- rep("null", G); truth[ups] <- "up"; truth[downs] <- "down"
data.frame(genes = G, patients = n, batches = nlevels(batch), truly_DE = n_de,
           median_count = median(sim), max_count = max(sim),
           dispersion_ratio_range = sprintf("%.1f-%.1f", min(b_disp), max(b_disp)))

genes	patients	batches	truly_DE	median_count	max_count	dispersion_ratio_range
10000	1000	10	500	99	602239	1.2-3.8

Batch effects act mainly along two technical directions, the way library depth, degradation and GC bias do, with the 10 batches scattered across that plane and a smaller idiosyncratic shift per gene per batch on top. Each batch has its own dispersion. Samples also carry continuous patient heterogeneity, without which every sample in a cell would share one mean vector and the batches would collapse to points. 500 of 10,000 genes are truly differential at 1.8-fold, half up and half down.

# scoring helpers, defined before the equivalence checks that use them
pc_sim <- function(mat) {
  p <- prcomp(t(cpm(mat, log = TRUE, prior.count = 1)), scale. = FALSE)
  # centre and loadings are kept so another matrix can be projected onto this basis, which is
  # what makes a component index mean the same thing across two arms
  list(x = p$x, ve = p$sdev^2 / sum(p$sdev^2), center = p$center, rotation = p$rotation)
}
sim_de <- function(mat) {
  d   <- calcNormFactors(DGEList(mat))
  des <- model.matrix(~ group)
  topTable(eBayes(lmFit(voom(d, des), des)), coef = 2, number = Inf, sort.by = "none")
}

1. ComBat-seq against the original

One correction, computed twice. Everything in this sub-section reads these two objects.

sim_ref <- sva::ComBat_seq(sim, batch = batch, group = group)

## Found 10 batches
## Using full model in ComBat-seq.
## Adjusting for 1 covariate(s) or covariate level(s)
## Estimating dispersions
## Fitting the GLM model
## Shrinkage off - using GLM estimates for parameters
## Adjusting the data

sim_par <- ComBat_seq_parallel(sim, batch = batch, group = group, workers = n_workers)

## Found 10 batches
## Using full model in ComBat-seq.
## Adjusting for 1 covariate(s) or covariate level(s)
## Estimating dispersions

## Fitting the GLM model
## Shrinkage off - using GLM estimates for parameters
## Adjusting the data

identical(sim_ref, sim_par)

## [1] TRUE

stopifnot(identical(sim_ref, sim_par))

Every stage below is computed from both corrected matrices and compared, not assumed equal because the inputs were.

# decomposed here and reused by the metrics table and the figure below, rather than repeated
q_rf       <- pc_sim(sim_ref); q_pr <- pc_sim(sim_par)
sim_tt_ref <- sim_de(sim_ref); sim_tt_par <- sim_de(sim_par)

data.frame(
  stage = c("corrected counts", "storage mode", "PCA scores", "PCA variance explained",
            "log fold changes", "t statistics", "adjusted p-values", "significant gene set"),
  identical = c(
    identical(sim_ref, sim_par),
    identical(storage.mode(sim_ref), storage.mode(sim_par)),
    identical(q_rf$x, q_pr$x),
    identical(q_rf$ve, q_pr$ve),
    identical(sim_tt_ref$logFC, sim_tt_par$logFC),
    identical(sim_tt_ref$t, sim_tt_par$t),
    identical(sim_tt_ref$adj.P.Val, sim_tt_par$adj.P.Val),
    identical(rownames(sim_tt_ref)[sim_tt_ref$adj.P.Val < 0.05],
              rownames(sim_tt_par)[sim_tt_par$adj.P.Val < 0.05])))

stage	identical
corrected counts	TRUE
storage mode	TRUE
PCA scores	TRUE
PCA variance explained	TRUE
log fold changes	TRUE
t statistics	TRUE
adjusted p-values	TRUE
significant gene set	TRUE

stopifnot(identical(sim_ref, sim_par),
          identical(q_rf$x, q_pr$x),
          identical(q_rf$ve, q_pr$ve),
          identical(sim_tt_ref$logFC, sim_tt_par$logFC),
          identical(sim_tt_ref$t, sim_tt_par$t),
          identical(sim_tt_ref$adj.P.Val, sim_tt_par$adj.P.Val))

q_un <- pc_sim(sim)   # the corrected arms were decomposed in the identity check above

# how far apart the batch centroids sit, in units of the spread inside a batch, so the number
# does not move when the components are rescaled
batch_sep <- function(q) {
  ctr <- t(vapply(levels(batch), function(l)
    colMeans(q$x[batch == l, 1:2, drop = FALSE]), numeric(2)))
  spread <- mean(vapply(levels(batch), function(l)
    sqrt(mean(apply(q$x[batch == l, 1:2, drop = FALSE], 2, var))), numeric(1)))
  mean(dist(ctr)) / spread
}
# the condition is not expected on PC1: it is a smaller effect than patient heterogeneity, and
# reporting only PC1 would score a preserved variable as destroyed
best_r2 <- function(q, f, k = 10L)
  max(vapply(seq_len(k), function(i) summary(lm(q$x[, i] ~ f))$r.squared, numeric(1)))

# The null for a maximum over ten components is not the null for one, so it is measured by
# permuting the labels rather than borrowing the single-component floor.
set.seed(7)
null_best <- mean(replicate(200, {
  g <- sample(group)
  max(vapply(1:10, function(i) summary(lm(q_un$x[, i] ~ g))$r.squared, numeric(1)))
}))

data.frame(quantity = c("batch R2 on PC1", "batch separation on PC1-PC2",
                        "condition R2, best of PC1-PC10"),
           uncorrected = round(c(summary(lm(q_un$x[, 1] ~ batch))$r.squared,
                                 batch_sep(q_un), best_r2(q_un, group)), 4),
           corrected   = round(c(summary(lm(q_rf$x[, 1] ~ batch))$r.squared,
                                 batch_sep(q_rf), best_r2(q_rf, group)), 4),
           chance_floor = c(round((nlevels(batch) - 1) / (ncol(sim) - 1), 4), NA,
                            round(null_best, 4)))

quantity	uncorrected	corrected	chance_floor
batch R2 on PC1	0.9830	0.0003	0.0090
batch separation on PC1-PC2	12.8284	0.0271	NA
condition R2, best of PC1-PC10	0.6370	0.8998	0.0039

ve_sim <- function(q, lab) sprintf("%s\nPC1 %.1f%%, PC2 %.1f%%", lab, 100 * q$ve[1], 100 * q$ve[2])
arm_of <- function(q, lab) data.frame(PC1 = q$x[, 1], PC2 = q$x[, 2], batch = batch, arm = lab)

sd_df <- rbind(arm_of(q_un, "Uncorrected"), arm_of(q_rf, "sva::ComBat_seq"),
               arm_of(q_pr, "ComBat_seq_parallel"))
sd_df$arm <- factor(sd_df$arm, levels = c("Uncorrected", "sva::ComBat_seq", "ComBat_seq_parallel"))
sim_lab <- setNames(c(ve_sim(q_un, "Uncorrected"), ve_sim(q_rf, "sva::ComBat_seq"),
                      ve_sim(q_pr, "ComBat_seq_parallel")), levels(sd_df$arm))

ggplot(sd_df, aes(PC1, PC2, colour = batch)) +
  geom_point(size = 1.5, alpha = 0.8) +
  stat_ellipse(aes(group = batch), linewidth = 0.6, type = "norm", level = 0.68) +
  scale_colour_manual(values = setNames(RColorBrewer::brewer.pal(10, "Paired"), levels(batch)),
                      name = "batch") +
  guides(colour = guide_legend(override.aes = list(size = 3, alpha = 1), ncol = 1)) +
  facet_wrap(~ arm, nrow = 1, scales = "free", labeller = labeller(arm = sim_lab)) +
  labs(title = "Simulated counts, coloured by batch",
       subtitle = sprintf("%s genes, %s patients, %d batches. Ten clouds become one.",
                           format(G, big.mark = ","), format(ncol(sim), big.mark = ","),
                           nlevels(batch)),
       x = "PC1", y = "PC2")

Ground truth is known by construction, so the correction can be scored rather than described.

score <- function(tt, lab) {
  called <- tt$adj.P.Val < 0.05
  data.frame(analysis = lab,
             called = sum(called),
             true_positives = sum(called & truth != "null"),
             false_positives = sum(called & truth == "null"),
             sensitivity = round(sum(called & truth != "null") / sum(truth != "null"), 3),
             FDR_observed = round(ifelse(sum(called) > 0,
                                         sum(called & truth == "null") / sum(called), 0), 3))
}
# the corrected arms were fitted in the identity check above
rbind(score(sim_de(sim), "uncorrected"),
      score(sim_tt_ref, "sva::ComBat_seq"),
      score(sim_tt_par, "ComBat_seq_parallel"))

analysis	called	true_positives	false_positives	sensitivity	FDR_observed
uncorrected	502	500	2	1	0.004
sva::ComBat_seq	516	500	16	1	0.031
ComBat_seq_parallel	516	500	16	1	0.031

stopifnot(identical(sim_tt_ref$logFC, sim_tt_par$logFC),
          identical(sim_tt_ref$adj.P.Val, sim_tt_par$adj.P.Val))

a. Speed

The sweep below runs the same simulated matrix in fresh sessions through callr::r, so no arm inherits a warm cache, a loaded package set, or a cluster from the one before it. Each arm is compared against the reference matrix itself, not against a summary of it.

The baseline arm is sva::ComBat_seq itself, not this package held to one worker.

sim_file <- file.path(tempdir(), "sim-bench.rds")
saveRDS(list(counts = sim, batch = batch, group = group, ref = sim_ref), sim_file)

sim_timings <- do.call(rbind, lapply(sort(unique(c(1L, 2L, 4L, n_workers))), function(w) {
  r <- callr::r(function(f, w) {
    suppressMessages({library(sva); library(rnaparallel)})
    d  <- readRDS(f)
    t0 <- Sys.time()
    out <- if (w == 1L) sva::ComBat_seq(d$counts, batch = d$batch, group = d$group)
           else ComBat_seq_parallel(d$counts, batch = d$batch, group = d$group, workers = w)
    # compared against the reference matrix itself. A checksum over 1e7 cells lands near
    # 1e16, past 2^53, where one unit in the last place is about 2 and a single-cell corruption
    # rounds away: the summary that was meant to catch it cannot represent it.
    list(secs = as.numeric(difftime(Sys.time(), t0, units = "secs")),
         same = identical(out, d$ref))
  }, args = list(f = sim_file, w = w))
  data.frame(workers = w, seconds = r$secs, same = r$same)
}))
unlink(sim_file)

sim_timings$implementation <- ifelse(sim_timings$workers == 1L, "sva::ComBat_seq",
                                     sprintf("ComBat_seq_parallel, %d workers", sim_timings$workers))
sim_timings$speedup <- sim_timings$seconds[1] / sim_timings$seconds
# the baseline arm is the reference, so it has nothing to compare against
sim_timings$identical_to_original <- ifelse(sim_timings$workers == 1L, NA, sim_timings$same)
sim_timings[, c("implementation", "seconds", "speedup", "identical_to_original")]

implementation	seconds	speedup	identical_to_original
sva::ComBat_seq	77.24844	1.000000	NA
ComBat_seq_parallel, 2 workers	33.70378	2.291982	TRUE
ComBat_seq_parallel, 4 workers	20.45763	3.776021	TRUE
ComBat_seq_parallel, 8 workers	16.21880	4.762895	TRUE

sim_best <- sim_timings[which.min(sim_timings$seconds), ]
sim_timings$arm <- factor(
  ifelse(sim_timings$workers == 1L, "sva::ComBat_seq\n(original)",
         sprintf("parallel\n%d workers", sim_timings$workers)),
  levels = ifelse(sim_timings$workers == 1L, "sva::ComBat_seq\n(original)",
                  sprintf("parallel\n%d workers", sim_timings$workers)))

ggplot(sim_timings, aes(arm, speedup)) +
  geom_col(width = 0.65,
           fill = ifelse(sim_timings$workers == 1L, "grey55",
                         ifelse(sim_timings$workers == sim_best$workers, "#b2182b", "#22506e"))) +
  geom_text(aes(label = sprintf("%.2fx\n%.1f s", speedup, seconds)), vjust = -0.3,
            size = 4.6, fontface = "bold", lineheight = 0.95) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.22))) +
  labs(title = "Simulated counts: original against the parallel companion",
       subtitle = sprintf("%s genes by %s patients, %d batches. Every arm returns a matrix identical to the original.",
                          format(G, big.mark = ","), format(ncol(sim), big.mark = ","), nlevels(batch)),
       x = NULL, y = "speedup vs sva::ComBat_seq")

par_arms <- sim_timings$identical_to_original[sim_timings$workers != 1L]
stopifnot(length(par_arms) >= 1L, !anyNA(par_arms), all(par_arms))

b. Preserving the biological variable

ComBat_seq is told which variation to keep. Whether that matters depends on how the batches were filled. The same simulation is rerun with the condition partly aligned to batch, which is what a real cohort looks like when one arm concentrates in a few cancer types.

unbal <- local({
  set.seed(123)
  n1c <- 85L; n1t <- 15L; n2c <- 15L; n2t <- 85L        # batch 1 mostly control, batch 2 mostly treated
  b <- factor(c(rep("1", n1c + n1t), rep("2", n2c + n2t)))
  g <- factor(c(rep("control", n1c), rep("treated", n1t),
                rep("control", n2c), rep("treated", n2t)))
  lv <- c("1.control", "2.control", "1.treated", "2.treated")
  cl <- factor(paste(b, g, sep = "."), levels = lv)
  bs <- round(runif(G, 20, 600))
  di <- sample.int(G, 100); u <- di[1:50]; d <- di[51:100]
  bu <- sample.int(G, G / 2)
  f <- matrix(1, G, 4, dimnames = list(NULL, lv))
  f[bu, c(2, 4)] <- batch_fold; f[-bu, c(2, 4)] <- 1 / batch_fold
  f[u, c(3, 4)] <- f[u, c(3, 4)] * bio_fold
  f[d, c(1, 2)] <- f[d, c(1, 2)] * bio_fold
  sz <- matrix(1 / disp_1, G, 4); sz[, c(2, 4)] <- 1 / (disp_1 * disp_fold)
  kk <- as.integer(cl)
  m <- vapply(seq_along(kk), function(j) {
    x <- rnbinom(G, mu = bs * f[, kk[j]], size = sz[, kk[j]]); x[x == 0] <- 1L; x
  }, numeric(G))
  storage.mode(m) <- "integer"; rownames(m) <- sprintf("gene%04d", seq_len(G))
  list(counts = m, batch = b, group = g, de = di)
})

score_un <- function(mat, lab) {
  des <- model.matrix(~ unbal$group)
  tt  <- topTable(eBayes(lmFit(voom(calcNormFactors(DGEList(mat)), des), des)),
                  coef = 2, number = Inf, sort.by = "none")
  sig <- tt$adj.P.Val < 0.05; is_de <- seq_len(G) %in% unbal$de
  data.frame(correction = lab, called = sum(sig), true_positives = sum(sig & is_de),
             false_positives = sum(sig & !is_de),
             sensitivity = round(sum(sig & is_de) / length(unbal$de), 3),
             FDR_observed = round(ifelse(sum(sig) > 0, sum(sig & !is_de) / sum(sig), 0), 3))
}
rbind(score_un(unbal$counts, "none"),
      score_un(ComBat_seq(unbal$counts, batch = unbal$batch, group = NULL),
               "corrected, nothing preserved"),
      score_un(ComBat_seq(unbal$counts, batch = unbal$batch, group = unbal$group),
               "corrected, condition preserved"))

## Found 2 batches
## Using null model in ComBat-seq.
## Adjusting for 0 covariate(s) or covariate level(s)
## Estimating dispersions
## Fitting the GLM model
## Shrinkage off - using GLM estimates for parameters
## Adjusting the data
## Found 2 batches
## Using full model in ComBat-seq.
## Adjusting for 1 covariate(s) or covariate level(s)
## Estimating dispersions
## Fitting the GLM model
## Shrinkage off - using GLM estimates for parameters
## Adjusting the data

correction	called	true_positives	false_positives	sensitivity	FDR_observed
none	9975	75	9900	0.75	0.992
corrected, nothing preserved	40	40	0	0.40	0.000
corrected, condition preserved	114	99	15	0.99	0.132

Without correction the batch effect produces mostly false positives. Correcting without naming the biology is worse in the other direction: where the condition partly aligns with batch, its signal is absorbed into the batch estimate and removed. Naming it recovers the genes.

Preserving the variable that is later tested makes the reported p-values slightly optimistic, so the observed false discovery rate above tends to sit over the nominal 0.05.

2. limma and edgeR against the originals

The correction is half a pipeline; the differential expression that follows is the other half. calcNormFactors_parallel, lmFit_parallel and duplicateCorrelation_parallel call the original edgeR and limma functions unchanged on blocks of genes or samples, and each stage below is timed twice on the same input, once through the original and once through the companion.

The contrast is the planted one, treated against control, on the corrected counts from C.1 at the same worker counts. Because the truth is known here, a divergence would show up as a different recovered gene set, not only as a failed identical().

sim_des <- model.matrix(~ group)

sim_nf_v <- calcNormFactors(DGEList(sim_ref))
sim_nf_p <- calcNormFactors_parallel(DGEList(sim_ref), workers = n_workers)

sim_vv    <- voom(sim_nf_v, sim_des)
sim_fit_v <- lmFit(sim_vv, sim_des)
sim_fit_p <- lmFit_parallel(sim_vv, sim_des, workers = n_workers)

sim_tt_v <- topTable(eBayes(sim_fit_v), coef = 2, number = Inf, sort.by = "none")
sim_tt_p <- topTable(eBayes(sim_fit_p), coef = 2, number = Inf, sort.by = "none")

# batch is the repeated unit here. A subset keeps the serial arm to minutes, and the per-gene
# work does not change with gene count, so the ratio is the same either way.
set.seed(20260820)
sim_dc_genes <- sort(sample.int(nrow(sim_vv), min(3000L, nrow(sim_vv))))
sim_dc_v <- duplicateCorrelation(sim_vv[sim_dc_genes, ], sim_des, block = batch)
sim_dc_p <- duplicateCorrelation_parallel(sim_vv[sim_dc_genes, ], sim_des, block = batch,
                                          workers = n_workers)

sim_limma_checks <- data.frame(
  stage = c("calcNormFactors (TMM)", "lmFit", "duplicateCorrelation",
            "topTable, final gene list"),
  identical = c(identical(sim_nf_v, sim_nf_p), identical(sim_fit_v, sim_fit_p),
                identical(sim_dc_v, sim_dc_p), identical(sim_tt_v, sim_tt_p)))
sim_limma_checks

stage	identical
calcNormFactors (TMM)	TRUE
lmFit	TRUE
duplicateCorrelation	TRUE
topTable, final gene list	TRUE

stopifnot(all(sim_limma_checks$identical))

Scored against the planted truth, so the two implementations are compared on what they recover rather than only on whether their numbers match.

rbind(score(sim_tt_v, "limma"),
      score(sim_tt_p, "rnaparallel companions"))

analysis	called	true_positives	false_positives	sensitivity	FDR_observed
limma	516	500	16	1	0.031
rnaparallel companions	516	500	16	1	0.031

# Seconds only: holding every fitted object would keep a second copy of the whole pipeline
# resident. Each original is timed on both sides of the companion arms and averaged, so drift
# moves both halves of a ratio together instead of inflating whichever arms ran last.
sim_secs <- function(e) { t0 <- Sys.time(); force(e)
                          as.numeric(difftime(Sys.time(), t0, units = "secs")) }

sim_stages <- c("calcNormFactors (TMM)", "lmFit", "duplicateCorrelation")
sim_grid   <- sort(unique(c(2L, 4L, n_workers)))
sim_pre <- c(sim_secs(calcNormFactors(DGEList(sim_ref))),
             sim_secs(lmFit(sim_vv, sim_des)),
             sim_secs(duplicateCorrelation(sim_vv[sim_dc_genes, ], sim_des, block = batch)))

sim_limma_speed <- do.call(rbind, lapply(sim_grid, function(w) {
  p <- c(sim_secs(calcNormFactors_parallel(DGEList(sim_ref), workers = w)),
         sim_secs(lmFit_parallel(sim_vv, sim_des, workers = w)),
         sim_secs(duplicateCorrelation_parallel(sim_vv[sim_dc_genes, ], sim_des,
                                                block = batch, workers = w)))
  invisible(gc(verbose = FALSE))
  data.frame(stage = sim_stages, workers = w, companion = p)
}))

sim_post <- c(sim_secs(calcNormFactors(DGEList(sim_ref))),
              sim_secs(lmFit(sim_vv, sim_des)),
              sim_secs(duplicateCorrelation(sim_vv[sim_dc_genes, ], sim_des, block = batch)))

sim_limma_speed$original <- rep((sim_pre + sim_post) / 2, length(sim_grid))
sim_limma_speed$speedup  <- sim_limma_speed$original / sim_limma_speed$companion
sim_drift <- max(abs(sim_post - sim_pre) / sim_pre)
sim_limma_speed[, c("stage", "workers", "original", "companion", "speedup")]

stage	workers	original	companion	speedup
calcNormFactors (TMM)	2	3.8752520	1.3183370	2.939500
lmFit	2	0.9672674	0.5578258	1.733995
duplicateCorrelation	2	46.2719774	23.3297691	1.983388
calcNormFactors (TMM)	4	3.8752520	0.7822599	4.953919
lmFit	4	0.9672674	0.3707290	2.609096
duplicateCorrelation	4	46.2719774	12.8427708	3.602959
calcNormFactors (TMM)	8	3.8752520	0.6043172	6.412613
lmFit	8	0.9672674	0.3329549	2.905101
duplicateCorrelation	8	46.2719774	9.8068671	4.718324

Each original was timed twice, once before the companion arms and once after, and the mean of the two is the denominator. The two readings differ by at most 5.9%, which is the measurement noise every speedup here carries.

ggplot(sim_limma_speed, aes(factor(workers), speedup, fill = stage)) +
  geom_col(position = position_dodge(0.8), width = 0.72) +
  geom_hline(yintercept = 1, linetype = 2, colour = "grey50") +
  geom_text(aes(label = sprintf("%.2fx", speedup)),
            position = position_dodge(0.8), vjust = -0.35, size = 3.1) +
  scale_fill_manual(values = c("calcNormFactors (TMM)" = "#4C72B0", "lmFit" = "#DD8452",
                               "duplicateCorrelation" = "#55A868")) +
  labs(x = "workers", y = "speedup against the original", fill = NULL,
       title = "limma and edgeR companions, simulated counts",
       subtitle = sprintf("%s genes by %d samples. duplicateCorrelation on a %s gene subset. Identity checked at %d workers.",
                           format(nrow(sim_ref), big.mark = ","), ncol(sim_ref),
                           format(length(sim_dc_genes), big.mark = ","), n_workers)) +
  expand_limits(y = 0)

D. TCGA

1. Cohort

TCGA cohorts recording tobacco smoking status. Sequencing plate is the batch variable. Both smoking and cancer type are preserved, passed together through covar_mod, so the correction removes plate variation while holding those two columns. It cannot promise to remove nothing else: variation confounded with plate that is not in the model goes with it.

Preserving cancer type requires that plates cross cancer types. Where a cohort shares no plate with any other, its cancer-type column is an exact sum of its plate columns. The design loses rank and ComBat-seq refuses, rather than return an answer that cannot be attributed to either variable.

The first run downloads the cohorts once. Later runs read the cached per-project files.

# The three largest cohorts with smoking status. Runtime is driven by the number of plates,
# because the GLM design carries one column per plate and its cost grows with the square of the
# design width, so cohort count is the lever that decides whether this finishes in a usable
# time. Three keeps the design narrow enough, and ComBat-seq still adjusts every gene that
# survives filtering rather than skipping any.
projects <- c("TCGA-HNSC", "TCGA-LUAD", "TCGA-LUSC")

# Two locations only, both chosen by the user: RNAPARALLEL_TCGA_DIR if set, otherwise the
# standard per-user cache. An earlier version also searched parent directories, which meant
# anyone able to write to an ancestor path on a shared machine could plant a directory the
# notebook would silently read from. The pre-rename cache name is still read, so an
# existing download is not fetched a second time.
cache_new <- tools::R_user_dir("rnaparallel", "cache")
cache_old <- tools::R_user_dir("combatseqparallel", "cache")
cache_dir <- if (!dir.exists(file.path(cache_new, "smoking")) &&
                 dir.exists(file.path(cache_old, "smoking"))) cache_old else cache_new
env_dir   <- Sys.getenv("RNAPARALLEL_TCGA_DIR", Sys.getenv("COMBATSEQ_TCGA_DIR", ""))
root      <- if (nzchar(env_dir)) env_dir else cache_dir
data_dir  <- normalizePath(file.path(root, "smoking"), mustWork = FALSE)

dir.create(data_dir, recursive = TRUE, showWarnings = FALSE)
absent <- projects[!file.exists(file.path(data_dir, paste0(projects, "-counts.rds")))]

if (length(absent)) {
  suppressMessages(library(TCGAbiolinks))
  for (p in absent) {
    message("downloading ", p)
    q <- GDCquery(project = p, data.category = "Transcriptome Profiling",
                  data.type = "Gene Expression Quantification",
                  workflow.type = "STAR - Counts")
    GDCdownload(q, method = "api", files.per.chunk = 40,
                directory = file.path(data_dir, "GDCdata"))
    se  <- GDCprepare(q, directory = file.path(data_dir, "GDCdata"), summarizedExperiment = TRUE)
    bc  <- colnames(se)
    idx <- which(substr(bc, 14, 15) == "01")          # primary tumour
    idx <- idx[order(bc[idx])]                        # deterministic aliquot choice
    idx <- idx[!duplicated(substr(bc[idx], 1, 12))]   # one per patient
    pc  <- rowData(se)$gene_type == "protein_coding"
    cts <- assay(se, "unstranded")[pc, idx, drop = FALSE]
    storage.mode(cts) <- "integer"
    saveRDS(list(counts = cts, barcode = colnames(se)[idx], project = p,
                 symbol = rowData(se)$gene_name[pc]),
            file.path(data_dir, paste0(p, "-counts.rds")))
    rm(se, cts); invisible(gc())
  }
}
parts <- lapply(projects, function(p) readRDS(file.path(data_dir, paste0(p, "-counts.rds"))))

suppressMessages(library(TCGAbiolinks))
bad  <- c("not reported", "unknown", "not available", "")
clin <- do.call(rbind, lapply(projects, function(p) {
  d <- suppressMessages(GDCquery_clinic(p, type = "clinical"))
  if (!"tobacco_smoking_status" %in% names(d)) return(NULL)
  data.frame(patient = d$submitter_id, status = as.character(d$tobacco_smoking_status),
             stringsAsFactors = FALSE)
}))
clin <- clin[!is.na(clin$status) & !(tolower(clin$status) %in% bad), ]
# The GDC label as recorded, not a collapsed grouping. Lifelong Non-Smoker is the reference,
# so the differential expression coefficient below is current smokers against never smokers and
# the three reformed categories keep their own terms rather than being pooled into "ever".
clin$smoking <- clin$status

genes   <- Reduce(intersect, lapply(parts, function(x) rownames(x$counts)))
counts  <- do.call(cbind, lapply(parts, function(x) x$counts[genes, , drop = FALSE]))
cancer  <- factor(unlist(lapply(parts, function(x) rep(sub("TCGA-", "", x$project), ncol(x$counts)))))
barcode <- unlist(lapply(parts, function(x) x$barcode))

sym <- parts[[1]]$symbol; names(sym) <- rownames(parts[[1]]$counts); sym <- sym[genes]
rm(parts); invisible(gc())

smoking <- clin$smoking[match(substr(barcode, 1, 12), clin$patient)]
has <- !is.na(smoking)
counts <- counts[, has]; cancer <- droplevels(cancer[has])
barcode <- barcode[has]

# Lifelong Non-Smoker first, so it is the reference level every other category is compared to
gdc_lv  <- c("Lifelong Non-Smoker",
             sort(setdiff(unique(smoking[has]), "Lifelong Non-Smoker")))
smoking <- droplevels(factor(smoking[has], levels = gdc_lv))

# A plate holding one sample carries no within-batch information, so ComBat-seq cannot estimate
# anything from it. The singletons are pooled into one SMALL level rather than dropped, which
# keeps their samples in the correction but means that one level is not a real plate: it is a
# bag of unrelated singletons sharing one estimated effect. It is excluded from the plate
# rankings below for that reason.
plate <- substr(barcode, 22, 25)
tab   <- table(plate)
plate <- factor(ifelse(plate %in% names(tab)[tab >= 2], plate, "SMALL"))

keep_g <- filterByExpr(counts, group = smoking)
counts <- counts[keep_g, ]; sym <- sym[keep_g]
storage.mode(counts) <- "integer"

stopifnot(nlevels(smoking) >= 2L, sum(smoking == "Current Smoker") >= 30L,
          levels(smoking)[1] == "Lifelong Non-Smoker")   # a GDC label change would
                                                            # otherwise surface an hour later
data.frame(genes = nrow(counts), samples = ncol(counts), plates = nlevels(plate),
           cancer_types = nlevels(cancer),
           never = sum(smoking == "Lifelong Non-Smoker"),
           current = sum(smoking == "Current Smoker"),
           other_categories = nlevels(smoking) - 2L)

genes	samples	plates	cancer_types	never	current	other_categories
18270	1500	54	3	210	430	3

a. Cohort characteristics

library(gtsummary)

## Warning: package 'gtsummary' was built under R version 4.4.3

theme_gtsummary_journal(journal = "jama")
theme_gtsummary_compact()

tbl1_data <- droplevels(data.frame(
  cancer      = cancer,
  lib_size    = colSums(counts) / 1e6,
  genes_det   = colSums(counts > 0),
  plate_n     = as.integer(table(plate)[as.character(plate)]),
  smoking     = smoking
))

tbl1 <- tbl_summary(
  tbl1_data,
  by = smoking,
  missing = "ifany", missing_text = "Unknown",
  missing_stat = "{N_miss} ({p_miss}%)",
  label = list(cancer     ~ "Cancer type",
               lib_size  ~ "Library size, millions of reads",
               genes_det ~ "Genes detected",
               plate_n   ~ "Samples sharing the plate"),
  statistic = list(all_continuous()  ~ "{mean}",
                   all_categorical() ~ "{n} ({p}%)")) |>
  add_overall(last = FALSE) |>
  add_ci(include = c(lib_size, genes_det, plate_n),
         method = list(all_continuous() ~ "t.test"), pattern = "{stat} ({ci})") |>
  add_p(test = list(cancer     ~ "chisq.test",
                    lib_size  ~ "kruskal.test",
                    genes_det ~ "kruskal.test",
                    plate_n   ~ "kruskal.test")) |>
  modify_header(all_stat_cols() ~ "**{level}**  \nN = {n}") |>
  bold_labels()
tbl1

Characteristic	Overall N = 1500	Lifelong Non-Smoker N = 210	Current Reformed Smoker for < or = 15 yrs N = 558	Current Reformed Smoker for > 15 yrs N = 291	Current Reformed Smoker, Duration Not Specified N = 11	Current Smoker N = 430	p-value1
Cancer type, n (%)							<0.001
HNSC	508 (34%)	117 (56%)	139 (25%)	73 (25%)	2 (18%)	177 (41%)
LUAD	503 (34%)	75 (36%)	169 (30%)	135 (46%)	4 (36%)	120 (28%)
LUSC	489 (33%)	18 (8.6%)	250 (45%)	83 (29%)	5 (45%)	133 (31%)
Library size, millions of reads, Mean	50 (49, 51)	51 (49, 53)	50 (49, 52)	50 (48, 52)	51 (35, 66)	50 (49, 52)	0.76
Genes detected, Mean	17,154 (17,137, 17,170)	17,091 (17,046, 17,136)	17,184 (17,159, 17,209)	17,105 (17,066, 17,145)	17,137 (16,870, 17,405)	17,178 (17,149, 17,208)	<0.001
Samples sharing the plate, Mean	35 (35, 36)	37 (35, 39)	35 (34, 36)	35 (34, 37)	29 (20, 39)	36 (34, 37)	0.77
1 Pearson’s Chi-squared test; Kruskal-Wallis rank sum test
Abbreviation: CI = Confidence Interval

Smoking status is the GDC label as recorded, not a collapsed grouping, so each reformed category keeps its own term and its own count. Lifelong Non-Smoker is the reference level. The differential expression in D.2 contrasts Current Smoker against it, so someone who stopped fifteen years ago is neither counted as exposed nor discarded from the correction.

The correction is applied once across the whole table above.

# built once here; the correction, the design-rank check and the DE model all read this
covar <- model.matrix(~ smoking + cancer)

data.frame(
  quantity = c("samples corrected", "genes", "batch variable",
               "batch levels", "singleton plates pooled",
               "biology preserved", "covariate columns"),
  value = c(format(ncol(counts), big.mark = ","),
            format(nrow(counts), big.mark = ","),
            "sequencing plate, barcode characters 22-25",
            nlevels(plate),
            sum(table(substr(barcode, 22, 25)) < 2),
            "smoking status and cancer type",
            ncol(covar) - 1))

quantity	value
samples corrected	1,500
genes	18,270
batch variable	sequencing plate, barcode characters 22-25
batch levels	54
singleton plates pooled	2
biology preserved	smoking status and cancer type
covariate columns	6

# both biological variables are preserved, so both go in covar_mod and group stays NULL
design <- cbind(model.matrix(~ -1 + plate), covar[, -1, drop = FALSE])
design_rank <- qr(design)$rank
data.frame(design_columns = ncol(design), rank = design_rank,
           estimable = design_rank >= ncol(design))

design_columns	rank	estimable
60	60	TRUE

stopifnot(design_rank >= ncol(design))

b. Argument parity

a <- formals(sva::ComBat_seq)
b <- formals(rnaparallel::ComBat_seq_parallel)
data.frame(argument = names(a),
           position_matches = match(names(a), names(b)) == seq_along(a),
           default_matches = vapply(names(a), function(n)
             identical(deparse(a[[n]]), deparse(b[[n]])), logical(1)),
           row.names = NULL)

argument	position_matches	default_matches
counts	TRUE	TRUE
batch	TRUE	TRUE
group	TRUE	TRUE
covar_mod	TRUE	TRUE
full_mod	TRUE	TRUE
shrink	TRUE	TRUE
shrink.disp	TRUE	TRUE
gene.subset.n	TRUE	TRUE

stopifnot(identical(names(a), names(b)[seq_along(a)]),
          all(vapply(names(a), function(n)
            identical(deparse(a[[n]]), deparse(b[[n]])), logical(1))))
data.frame(parallel_only_arguments = setdiff(names(b), names(a)))

parallel_only_arguments
workers
chunks
parallel_backend
backend

2. ComBat-seq against the original

One correction of the whole cohort per worker count, each compared against sva::ComBat_seq on the same call. Only the last arm is kept for the figures below.

t0  <- Sys.time()
ref <- sva::ComBat_seq(counts, batch = plate, group = NULL, covar_mod = covar)

## Found 54 batches
## Using null model in ComBat-seq.
## Adjusting for 6 covariate(s) or covariate level(s)
## Estimating dispersions
## Fitting the GLM model
## Shrinkage off - using GLM estimates for parameters
## Adjusting the data

secs_ref <- as.numeric(difftime(Sys.time(), t0, units = "secs"))

# Each worker count corrects the whole matrix again. Only the last arm is kept: a corrected
# copy of this cohort is about 110 MB, and holding every arm buys nothing once its output has
# been compared against the reference.
secs_w <- numeric(length(tcga_workers))
same_w <- logical(length(tcga_workers))
for (i in seq_along(tcga_workers)) {
  t1  <- Sys.time()
  out <- ComBat_seq_parallel(counts, batch = plate, group = NULL, covar_mod = covar,
                             workers = tcga_workers[i])
  secs_w[i] <- as.numeric(difftime(Sys.time(), t1, units = "secs"))
  same_w[i] <- identical(out, ref)
  if (i == length(tcga_workers)) par <- out
  rm(out); invisible(gc(verbose = FALSE))
}

## Found 54 batches
## Using null model in ComBat-seq.
## Adjusting for 6 covariate(s) or covariate level(s)
## Estimating dispersions
## Fitting the GLM model
## Shrinkage off - using GLM estimates for parameters
## Adjusting the data
## Found 54 batches
## Using null model in ComBat-seq.
## Adjusting for 6 covariate(s) or covariate level(s)
## Estimating dispersions
## Fitting the GLM model
## Shrinkage off - using GLM estimates for parameters
## Adjusting the data
## Found 54 batches
## Using null model in ComBat-seq.
## Adjusting for 6 covariate(s) or covariate level(s)
## Estimating dispersions
## Fitting the GLM model
## Shrinkage off - using GLM estimates for parameters
## Adjusting the data

stopifnot(all(same_w), identical(ref, par))

a. Principal components

pca_of <- pc_sim                      # same decomposition the simulated section uses
p_un <- pca_of(counts); p_rf <- pca_of(ref); p_pr <- pca_of(par)

identical(p_rf$x, p_pr$x) && identical(p_rf$ve, p_pr$ve)

## [1] TRUE

pdf_of <- function(p, lab) data.frame(PC1 = p$x[, 1], PC2 = p$x[, 2],
                                      cancer = cancer, arm = lab)
pd <- rbind(pdf_of(p_un, "Uncorrected"), pdf_of(p_rf, "sva::ComBat_seq"),
            pdf_of(p_pr, "ComBat_seq_parallel"))
pd$arm <- factor(pd$arm, levels = c("Uncorrected", "sva::ComBat_seq", "ComBat_seq_parallel"))

# Each arm is decomposed separately, so its components explain their own share of a different
# total variance. One shared axis label would report the uncorrected figure under all three.
ve_strip <- function(p, lab) sprintf("%s\nPC1 %.1f%%, PC2 %.1f%%", lab, 100 * p$ve[1], 100 * p$ve[2])
arm_lab <- setNames(c(ve_strip(p_un, "Uncorrected"), ve_strip(p_rf, "sva::ComBat_seq"),
                      ve_strip(p_pr, "ComBat_seq_parallel")), levels(pd$arm))

ggplot(pd, aes(PC1, PC2, colour = cancer)) +
  geom_point(size = 0.9, alpha = 0.8) +
  scale_colour_brewer(palette = "Paired", name = "cancer type") +
  guides(colour = guide_legend(override.aes = list(size = 3.4, alpha = 1), ncol = 1)) +
  facet_wrap(~ arm, nrow = 1, scales = "free", labeller = labeller(arm = arm_lab)) +
  labs(title = "Principal components, coloured by cancer type",
       subtitle = sprintf("%s samples, %d cancer types, %d plates. Right two panels: same numbers, two ways.",
                           format(ncol(counts), big.mark = ","), nlevels(cancer), nlevels(plate)),
       x = "PC1", y = "PC2")

# Which plate actually deviates most? Rank them by how far the plate centroid sits from the
# cohort centroid on PC1 and PC2, in units of the within-plate spread, so a large plate with a
# small offset does not outrank a small plate that is genuinely displaced. The worst plate is
# selected on the uncorrected arm, so its own shrinkage regresses toward the mean whether or not
# the correction works; the median across all plates is reported beside it for that reason.
plate_dev <- function(x, i1 = pcx, i2 = pcy) {
  vapply(levels(plate), function(l) {
    i <- plate == l
    if (sum(i) < 3) return(NA_real_)
    ctr <- c(mean(x[i, i1]), mean(x[i, i2]))
    sqrt(sum(ctr^2)) / sqrt(mean(c(var(x[i, i1]), var(x[i, i2]))))
  }, numeric(1))
}
dev_un <- plate_dev(p_un$x); dev_rf <- plate_dev(p_rf$x)

# PC1 and PC2 are mostly cancer type, and covar_mod is told to keep cancer type. Because each
# cohort has its own plates, a plate centroid on those axes carries the biology the correction
# must not remove, so plate displacement there is not expected to fall and can rise as the
# within-plate spread tightens. The components where plate acts on its own are the ones that
# say whether the correction worked, so both numbers are reported.
plate_gap <- vapply(1:8, function(k)
  summary(lm(p_un$x[, k] ~ plate))$r.squared - summary(lm(p_un$x[, k] ~ cancer))$r.squared,
  numeric(1))
pdom <- order(plate_gap, decreasing = TRUE)[1:2]

# Component k of one decomposition is not component k of another: separately fitted PCAs can
# rotate and reorder, so a reduction read off the same index compares two different directions.
# The corrected arm is projected onto the uncorrected loadings instead, where the index means
# one thing and the before/after difference is the correction rather than a change of basis.
onto_un <- function(mat) {
  lc <- t(cpm(mat, log = TRUE, prior.count = 1))
  sweep(lc, 2, p_un$center, "-") %*% p_un$rotation
}
dev_un_p <- plate_dev(p_un$x, pdom[1], pdom[2])
dev_rf_p <- plate_dev(onto_un(ref), pdom[1], pdom[2])
keep_pl  <- setdiff(names(dev_un), "SMALL")
med_red_p <- median(1 - dev_rf_p[keep_pl] / dev_un_p[keep_pl], na.rm = TRUE)
worst  <- setdiff(names(sort(dev_un, decreasing = TRUE)), "SMALL")[1:10]

med_red <- median(1 - dev_rf[keep_pl] / dev_un[keep_pl], na.rm = TRUE)

data.frame(plate = c(worst, "median of all plates"),
           n = c(as.integer(table(plate)[worst]), NA_integer_),
           deviation_uncorrected = round(c(dev_un[worst], NA), 3),
           deviation_corrected = round(c(dev_rf[worst], NA), 3),
           reduction_PC1_PC2 = sprintf("%.0f%%", 100 * c(1 - dev_rf[worst] / dev_un[worst], med_red)),
           reduction_plate_PCs = sprintf("%.0f%%",
             100 * c(1 - dev_rf_p[worst] / dev_un_p[worst], med_red_p)),
           row.names = NULL)

plate	n	deviation_uncorrected	deviation_corrected	reduction_PC1_PC2	reduction_plate_PCs
2287	17	5.271	4.807	9%	39%
A24X	34	4.970	4.166	16%	65%
2241	35	4.784	4.823	-1%	-12%
1107	11	4.709	4.169	11%	59%
2170	22	4.186	4.300	-3%	90%
2039	33	3.933	4.015	-2%	81%
A39D	13	3.923	3.892	1%	58%
0946	21	3.488	3.935	-13%	61%
1755	36	3.110	3.509	-13%	83%
1858	36	3.070	3.477	-13%	62%
median of all plates	NA	NA	NA	-13%	60%

# the ten most displaced plates in colour, every other plate grey behind them
hl <- factor(ifelse(plate %in% worst, as.character(plate), "other plates"),
             levels = c(worst, "other plates"))

pdp <- rbind(data.frame(X = p_un$x[, pcx], Y = p_un$x[, pcy], arm = "Uncorrected"),
             data.frame(X = p_rf$x[, pcx], Y = p_rf$x[, pcy], arm = "sva::ComBat_seq"),
             data.frame(X = p_pr$x[, pcx], Y = p_pr$x[, pcy], arm = "ComBat_seq_parallel"))
pdp$hl  <- rep(hl, 3)
pdp$arm <- factor(pdp$arm, levels = c("Uncorrected", "sva::ComBat_seq", "ComBat_seq_parallel"))

pal_plate <- setNames(c(RColorBrewer::brewer.pal(10, "Paired"), "grey86"), levels(hl))

arm_lab_p <- setNames(c(ve_strip(p_un, "Uncorrected"), ve_strip(p_rf, "sva::ComBat_seq"),
                        ve_strip(p_pr, "ComBat_seq_parallel")), levels(pdp$arm))

ggplot(pdp, aes(X, Y, colour = hl)) +
  geom_point(data = subset(pdp, hl == "other plates"), size = 0.5, alpha = 0.25) +
  geom_point(data = subset(pdp, hl != "other plates"), size = 1.3, alpha = 0.9) +
  stat_ellipse(data = subset(pdp, hl != "other plates"), aes(group = hl),
               linewidth = 0.9, type = "norm", level = 0.68) +
  scale_colour_manual(values = pal_plate, name = "plate",
                      guide = guide_legend(override.aes = list(size = 3.4, alpha = 1))) +
  facet_wrap(~ arm, nrow = 1, scales = "free", labeller = labeller(arm = arm_lab_p)) +
  labs(title = sprintf("The %d most displaced plates of %d, and what the correction does to them",
                       length(worst), nlevels(plate)),
       subtitle = sprintf("PC1-PC2 carry cancer type, which covar_mod preserves: median plate moves %.0f%%. On PC%d and PC%d, where plate acts alone, it falls %.0f%%.",
                           100 * med_red, pdom[1], pdom[2], 100 * med_red_p),
       x = "PC1", y = "PC2")

b. Differential expression

The contrast is current smokers against lifelong non-smokers. Former smokers stay in the cohort and in the correction, but are not part of this comparison. Cancer type sits in the design as an adjustment, not as a second contrast, because smoking is heavily imbalanced across cohorts: LUSC has 18 lifelong non-smokers against 471 patients with a smoking history, and LUAD has 75 against 428. Both models are reported below, so the size of that confounding is visible rather than asserted.

de_design <- model.matrix(~ smoking + cancer)          # smoking contrast, adjusted for cancer
de_design_unadj <- model.matrix(~ smoking)             # smoking alone, for comparison

# the current-vs-never column, matched rather than hardcoded: the GDC labels contain spaces and
# punctuation, so the coefficient name is not something to type from memory
de_coef <- grep("Current Smoker$", colnames(de_design), value = TRUE)
stopifnot(length(de_coef) == 1L)

run_de <- function(mat, des = de_design) {
  d <- calcNormFactors(DGEList(mat))
  topTable(eBayes(lmFit(voom(d, des), des)), coef = de_coef,
           number = Inf, sort.by = "none")
}
tt_ref <- run_de(ref); tt_par <- run_de(par)
tt_unadj <- run_de(ref, de_design_unadj)

data.frame(quantity = c("log fold change", "t statistic", "adjusted p", "significant set"),
           identical = c(identical(tt_ref$logFC, tt_par$logFC),
                         identical(tt_ref$t, tt_par$t),
                         identical(tt_ref$adj.P.Val, tt_par$adj.P.Val),
                         identical(rownames(tt_ref)[tt_ref$adj.P.Val < 0.05],
                                   rownames(tt_par)[tt_par$adj.P.Val < 0.05])))

quantity	identical
log fold change	TRUE
t statistic	TRUE
adjusted p	TRUE
significant set	TRUE

sig_adj <- rownames(tt_ref)[tt_ref$adj.P.Val < 0.05 & abs(tt_ref$logFC) > 1]
sig_una <- rownames(tt_unadj)[tt_unadj$adj.P.Val < 0.05 & abs(tt_unadj$logFC) > 1]

data.frame(model = c("~ smoking", "~ smoking + cancer"),
           significant = c(length(sig_una), length(sig_adj)),
           shared = length(intersect(sig_una, sig_adj)),
           only_this_model = c(length(setdiff(sig_una, sig_adj)),
                               length(setdiff(sig_adj, sig_una))))

model	significant	shared	only_this_model
~ smoking	200	59	141
~ smoking + cancer	71	59	12

Genes called by the unadjusted model alone are the cost of ignoring cancer type. Smoking is not randomised across these cohorts, so an unadjusted contrast partly measures which cancers the smokers have.

stopifnot(identical(tt_ref$logFC, tt_par$logFC),
          identical(tt_ref$t, tt_par$t),
          identical(tt_ref$adj.P.Val, tt_par$adj.P.Val))

# label helper: GDC ships versioned Ensembl identifiers, the figures want symbols
lab_of <- function(ids) ifelse(is.na(sym[ids]) | sym[ids] == "", ids, sym[ids])

library(ggrepel); library(patchwork)

## Warning: package 'ggrepel' was built under R version 4.4.3

vd <- rbind(data.frame(gene = rownames(tt_ref), logFC = tt_ref$logFC, p = tt_ref$adj.P.Val,
                       arm = "sva::ComBat_seq"),
            data.frame(gene = rownames(tt_par), logFC = tt_par$logFC, p = tt_par$adj.P.Val,
                       arm = "ComBat_seq_parallel"))
vd$arm <- factor(vd$arm, levels = c("sva::ComBat_seq", "ComBat_seq_parallel"))
vd$symbol <- lab_of(vd$gene)
vd$state <- "not significant"
vd$state[vd$p < 0.05 & vd$logFC >  1] <- "up in current smokers"
vd$state[vd$p < 0.05 & vd$logFC < -1] <- "up in never smokers"
vd$state <- factor(vd$state, levels = c("up in current smokers", "up in never smokers", "not significant"))
pal_state <- c(`up in current smokers` = "#B2182B", `up in never smokers` = "#2166AC",
               `not significant` = "grey82")

# counts per arm, split left and right so each sits over the side it describes
n_lab <- do.call(rbind, lapply(split(vd, vd$arm), function(d) data.frame(
  arm  = d$arm[1],
  up   = sum(d$p < 0.05 & d$logFC >  1, na.rm = TRUE),
  down = sum(d$p < 0.05 & d$logFC < -1, na.rm = TRUE))))
n_lab$x_up   <- max(vd$logFC, na.rm = TRUE) * 0.78
n_lab$x_down <- min(vd$logFC, na.rm = TRUE) * 0.78
# a strip above the tallest point, so labels and counts never share space
n_lab$y_lab  <- max(-log10(vd$p), na.rm = TRUE) * 1.06

# selected by adjusted p, then ordered by fold change so the side table reads as two blocks,
# everything raised in current smokers above everything lowered
top10 <- head(tt_par[order(tt_par$adj.P.Val), ], 10)
top10 <- top10[order(top10$logFC, decreasing = TRUE), ]
top10$symbol <- lab_of(rownames(top10))
lab <- vd[vd$gene %in% rownames(top10), ]

p_vol <- ggplot(vd, aes(logFC, -log10(p), colour = state)) +
  geom_point(size = 0.7, alpha = 0.6) +
  geom_point(data = lab, size = 1.7, colour = "black") +
  geom_text_repel(data = lab, aes(label = symbol), size = 3.5, fontface = "bold",
                  colour = "black", min.segment.length = 0, max.overlaps = Inf,
                  box.padding = 0.45, segment.colour = "grey45", seed = 1) +
  # the counts sit on their own strip above the cloud, so they cannot collide with a
  # gene label placed near the top of the panel
  geom_text(data = n_lab, aes(x = x_up, y = y_lab, label = paste0("Up: ", up)),
            colour = "#B2182B", fontface = "bold", size = 5, hjust = 0.5, vjust = 0,
            inherit.aes = FALSE) +
  geom_text(data = n_lab, aes(x = x_down, y = y_lab, label = paste0("Down: ", down)),
            colour = "#2166AC", fontface = "bold", size = 5, hjust = 0.5, vjust = 0,
            inherit.aes = FALSE) +
  scale_y_continuous(expand = expansion(mult = c(0.02, 0.16))) +
  scale_colour_manual(values = pal_state, name = NULL) +
  guides(colour = guide_legend(override.aes = list(size = 3.2, alpha = 1))) +
  facet_wrap(~ arm, nrow = 1) +
  labs(x = "log2 fold change", y = "-log10 adjusted p") +
  theme(legend.position = "bottom")

tab10 <- data.frame(y = rev(seq_len(nrow(top10))), symbol = top10$symbol,
                    fc = sprintf("%+.1f", top10$logFC),
                    padj = format(top10$adj.P.Val, digits = 2, scientific = TRUE))
p_tab <- ggplot(tab10, aes(y = y)) +
  geom_text(aes(x = 0.00, label = symbol), hjust = 0, size = 3.9, fontface = "bold") +
  geom_text(aes(x = 1.02, label = fc), hjust = 1, size = 3.7,
            colour = ifelse(top10$logFC > 0, "#B2182B", "#2166AC")) +
  geom_text(aes(x = 1.92, label = padj), hjust = 1, size = 3.3, colour = "grey35") +
  annotate("text", x = c(0, 1.02, 1.92), y = nrow(top10) + 1, hjust = c(0, 1, 1),
           label = c("gene", "logFC", "FDR"), size = 3.5, fontface = "bold", colour = "grey20") +
  scale_x_continuous(limits = c(-0.04, 1.96)) +
  scale_y_continuous(limits = c(0.4, nrow(top10) + 1.6)) +
  labs(title = "Top 10 by adjusted p") +
  theme_void(base_size = 13) +
  theme(plot.title = element_text(face = "bold", size = 13, margin = margin(b = 8)),
        plot.margin = margin(6, 6, 6, 14))

p_vol + p_tab + plot_layout(widths = c(4, 1.5)) +
  plot_annotation(
    title = "Differential expression, current versus never smokers",
    subtitle = sprintf("%s genes at FDR < 0.05 and |logFC| > 1 in both arms. Several top genes are Y-linked: part of this contrast is sex.",
                        format(sum(tt_par$adj.P.Val < 0.05 & abs(tt_par$logFC) > 1), big.mark = ",")),
    theme = theme(plot.title = element_text(face = "bold", size = 18, margin = margin(b = 4)),
                  plot.subtitle = element_text(size = 13, colour = "grey30", margin = margin(b = 10))))

library(ComplexHeatmap)

# the heatmap shows the two groups being contrasted; former smokers are corrected and kept in
# the cohort but are not part of this comparison
hm_keep  <- smoking %in% c("Lifelong Non-Smoker", "Current Smoker")

# every gene the volcano calls significant, on the same two thresholds, rather than a round
# number off the top of the table
hm_sig  <- which(tt_ref$adj.P.Val < 0.05 & abs(tt_ref$logFC) > 1)
hm_sig  <- hm_sig[order(tt_ref$adj.P.Val[hm_sig])]
hm_cap  <- 120L                       # a stated ceiling, reported in the title when it bites
hm_more <- max(0L, length(hm_sig) - hm_cap)
top_g   <- rownames(tt_ref)[head(hm_sig, hm_cap)]
top_n   <- length(top_g)
top_lab <- lab_of(top_g)
hm_fontsize <- max(4.5, min(9, 633 / top_n))   # 633pt of panel once annotations are drawn

# Scaled before plotting: each gene is centred and divided by its own standard deviation across
# the displayed samples, so the colour reads as contrast rather than expression level.
# ComplexHeatmap does no scaling of its own, so this is the only place it happens.
zmat <- function(mat) {
  # library sizes must come from the whole matrix; deriving them from the displayed genes shifts
  # each sample by its own subset total, and that shift does not cancel in a per-gene z-score
  m <- cpm(mat, log = TRUE, prior.count = 1)[top_g, hm_keep, drop = FALSE]
  rownames(m) <- top_lab
  z <- t(scale(t(m)))            # scale() works down columns, so transpose to scale per gene
  # a gene with no variance across the displayed samples divides by zero and comes back NaN,
  # which draws as a blank row and reads as missing data rather than as no contrast
  z[!is.finite(z)] <- 0
  z
}
z_ref <- zmat(ref); z_par <- zmat(par)
hm_flat <- sum(!is.finite(rowSums(t(scale(t(
  cpm(ref, log = TRUE, prior.count = 1)[top_g, hm_keep, drop = FALSE]))))))

col_fun <- circlize::colorRamp2(c(-2, 0, 2), c("#2166AC", "grey96", "#B2182B"))

# gene-level annotation: each row carries its own effect size and direction from the DE result,
# so the heatmap is readable without cross-referencing the volcano
top_fc  <- tt_ref$logFC[match(top_g, rownames(tt_ref))]
row_ann <- rowAnnotation(
  logFC = anno_barplot(top_fc, baseline = 0, bar_width = 0.8,
                       gp = gpar(fill = ifelse(top_fc > 0, "#B2182B", "#2166AC"), col = NA),
                       axis_param = list(gp = gpar(fontsize = 9)), width = unit(2.2, "cm")),
  Direction = ifelse(top_fc > 0, "up in current", "up in never"),
  col = list(Direction = c(`up in current` = "#B2182B", `up in never` = "#2166AC")),
  annotation_name_gp = gpar(fontsize = 11, fontface = "bold"),
  annotation_name_rot = 90,
  annotation_legend_param = list(Direction = list(title_gp = gpar(fontsize = 12, fontface = "bold"),
                                                  labels_gp = gpar(fontsize = 11))))

# the annotation name is drawn into the neighbouring heatmap and clipped, so only the
# right-hand panel carries it
mk_ann <- function(show_name) HeatmapAnnotation(
  Smoking = smoking[hm_keep],
  Cancer  = cancer[hm_keep],
  col = list(Smoking = setNames(
               RColorBrewer::brewer.pal(max(3L, nlevels(droplevels(smoking[hm_keep]))), "Set2")[
                 seq_len(nlevels(droplevels(smoking[hm_keep])))],
               levels(droplevels(smoking[hm_keep]))),
             Cancer  = setNames(RColorBrewer::brewer.pal(nlevels(cancer), "Paired"),
                                levels(cancer))),
  show_annotation_name = show_name,
  annotation_name_side = "right",
  annotation_name_gp = gpar(fontsize = 12, fontface = "bold"),
  annotation_legend_param = list(
    Smoking = list(title_gp = gpar(fontsize = 12, fontface = "bold"),
                   labels_gp = gpar(fontsize = 11)),
    Cancer  = list(title_gp = gpar(fontsize = 12, fontface = "bold"),
                   labels_gp = gpar(fontsize = 11), ncol = 2)))

mk_hm <- function(z, title, show_name, nm = "z-score", right = NULL) {
  Heatmap(z, name = nm, col = col_fun, top_annotation = mk_ann(show_name),
          right_annotation = right,
          show_heatmap_legend = identical(nm, "z-score"),
          column_title = title,
          column_title_gp = gpar(fontsize = 15, fontface = "bold"),
          # the panel is a fixed height, so the label size has to follow the row count
          show_row_names = TRUE, row_names_gp = gpar(fontsize = hm_fontsize),
          show_column_names = FALSE,
          # both axes cluster on the z-scores alone. Splitting the columns by smoking status
          # would draw the separation the figure is supposed to be evidence for.
          cluster_columns = TRUE, cluster_rows = TRUE,
          show_column_dend = TRUE, column_dend_height = unit(9, "mm"),
          heatmap_legend_param = list(title_gp = gpar(fontsize = 12, fontface = "bold"),
                                      labels_gp = gpar(fontsize = 11)))
}

draw(mk_hm(z_ref, "sva::ComBat_seq", FALSE) +
       mk_hm(z_par, "ComBat_seq_parallel", TRUE, "z2", right = row_ann),
     column_title = if (hm_more > 0L)
       sprintf("Top %d of %d significant genes, clustered unsupervised", top_n, top_n + hm_more)
     else
       sprintf("All %d significant genes, clustered unsupervised", top_n),
     column_title_gp = gpar(fontsize = 18, fontface = "bold"),
     merge_legends = TRUE)

data.frame(z_scores_identical = identical(z_ref, z_par),
           genes_shown = top_n,
           genes_with_no_variance = hm_flat,
           z_range = sprintf("%.2f to %.2f", min(z_ref), max(z_ref)))

z_scores_identical	genes_shown	genes_with_no_variance	z_range
TRUE	71	0	-3.56 to 4.87

stopifnot(identical(z_ref, z_par),
          all(is.finite(z_ref)),
          # a correctly scaled gene has mean 0 and unit variance across the displayed samples
          all(abs(rowMeans(z_ref)) < 1e-8))

c. Runtime

Every arm corrects the same matrix: 18,270 genes by 1,500 samples across 54 plates, the same call each time. These arms share one session and run in a fixed order, sva::ComBat_seq first, so the parallel arms run warm. The C.1 sweep is the isolated fresh-session measurement.

Speed depends on the machine, so the machine is reported rather than left to be inferred.

sysctl <- function(k) tryCatch(system2("sysctl", c("-n", k), stdout = TRUE, stderr = FALSE),
                               error = function(e) NA_character_)
mac <- identical(Sys.info()[["sysname"]], "Darwin")

data.frame(
  quantity = c("CPU", "performance cores", "logical cores", "memory", "OS", "worker counts"),
  value = as.character(c(
    if (mac) sysctl("machdep.cpu.brand_string") else NA_character_,
    if (mac) sysctl("hw.perflevel0.physicalcpu") else NA_character_,
    parallel::detectCores(),
    if (mac) sprintf("%.0f GB", as.numeric(sysctl("hw.memsize")) / 2^30) else NA_character_,
    paste(Sys.info()[["sysname"]], Sys.info()[["release"]]),
    paste(tcga_workers, collapse = ", "))))

quantity	value
CPU	Apple M3
performance cores	4
logical cores	8
memory	24 GB
OS	Darwin 25.6.0
worker counts	2, 4, 8

data.frame(implementation = c("sva::ComBat_seq",
                              sprintf("ComBat_seq_parallel, %d workers", tcga_workers)),
           seconds = round(c(secs_ref, secs_w), 1),
           speedup = sprintf("%.2fx", c(1, secs_ref / secs_w)),
           identical_to_original = c(NA, same_w))

implementation	seconds	speedup	identical_to_original
sva::ComBat_seq	1587.6	1.00x	NA
ComBat_seq_parallel, 2 workers	848.0	1.87x	TRUE
ComBat_seq_parallel, 4 workers	450.0	3.53x	TRUE
ComBat_seq_parallel, 8 workers	294.2	5.40x	TRUE

rt <- data.frame(
  arm = c("sva::ComBat_seq\n(original)", sprintf("parallel\n%d workers", tcga_workers)),
  seconds = c(secs_ref, secs_w),
  speedup = c(1, secs_ref / secs_w))
rt$arm <- factor(rt$arm, levels = rt$arm)

ggplot(rt, aes(arm, speedup)) +
  geom_col(width = 0.55, fill = c("grey55",
    colorRampPalette(c("#EF8A62", "#b2182b"))(length(tcga_workers)))) +
  geom_text(aes(label = sprintf("%.2fx\n%.0f s", speedup, seconds)), vjust = -0.3,
            size = 5, fontface = "bold", lineheight = 0.95) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.24))) +
  labs(title = "Original ComBat-seq against the parallel companion",
       subtitle = sprintf("%s genes by %s samples, %d plates, %d cancer types. Output identical to the original.",
                          format(nrow(counts), big.mark = ","),
                          format(ncol(counts), big.mark = ","), nlevels(plate), nlevels(cancer)),
       x = NULL, y = "speedup vs sva::ComBat_seq")

3. limma and edgeR against the originals

The same comparison as C.2, at cohort scale, on the corrected matrix from D.2.

Part a runs the exact differential expression the volcano and heatmap were drawn from, a second way, and the table asserts that rather than claiming it. Part b adds a blocking factor, which is where the expensive per-gene REML fits live. Part c records what was measured and deliberately left serial.

a. The differential expression pipeline, both ways

nf_v <- calcNormFactors(DGEList(ref))
nf_p <- calcNormFactors_parallel(DGEList(ref), workers = n_workers)

vv <- voom(nf_v, de_design)
fit_v <- lmFit(vv, de_design)
fit_p <- lmFit_parallel(vv, de_design, workers = n_workers)

tt_limma_v <- topTable(eBayes(fit_v), coef = de_coef, number = Inf, sort.by = "none")
tt_limma_p <- topTable(eBayes(fit_p), coef = de_coef, number = Inf, sort.by = "none")

limma_checks <- data.frame(
  stage = c("calcNormFactors (TMM)", "lmFit", "topTable, final gene list",
            "matches the gene list D.2 already drew"),
  identical = c(identical(nf_v, nf_p), identical(fit_v, fit_p),
                identical(tt_limma_v, tt_limma_p), identical(tt_limma_v, tt_ref)))
limma_checks

stage	identical
calcNormFactors (TMM)	TRUE
lmFit	TRUE
topTable, final gene list	TRUE
matches the gene list D.2 already drew	TRUE

stopifnot(all(limma_checks$identical))

The last row is the one that matters. A difference of one unit in the last place of sigma is invisible in the fit and becomes a different gene list once Benjamini-Hochberg has run, so the assertion is made on topTable rather than on the object it came from.

b. Blocked designs

duplicateCorrelation is the expensive step in a repeated-measures limma analysis: one REML fit and one LAPACK SVD per gene. Plate is the repeated unit here.

# The serial arm has to run in full for there to be anything to compare against, and the
# per-gene work does not change with gene count, so the ratio is the same on a subset while the
# baseline costs minutes instead of hours. Fixed random draw, and the identical() claim below is
# about these genes.
set.seed(20260820)
dc_genes <- sort(sample.int(nrow(vv), min(4000L, nrow(vv))))
v_dc <- vv[dc_genes, ]

dc_v <- duplicateCorrelation(v_dc, de_design, block = plate)
dc_p <- duplicateCorrelation_parallel(v_dc, de_design, block = plate, workers = n_workers)
dc_same <- identical(dc_v, dc_p)

data.frame(quantity = c("genes used", "consensus correlation", "identical()"),
           value = c(format(length(dc_genes), big.mark = ","),
                     sprintf("%.6f", dc_v$consensus.correlation), dc_same))

quantity	value
genes used	4,000
consensus correlation	-0.003276
identical()	TRUE

Timing for this stage is in D.3.d, where the baseline is bracketed rather than read once.

stopifnot(dc_same)

c. What is deliberately not parallelised

Not everything survives the exactness bar, and the list of rejections is as much a part of this package as the list of companions. voom measured 0.99x, because its lowess trend takes its span as a fraction of the gene count. contrasts.fit returns before a fork could finish starting. eBayes, fitFDist and topTable pool across every gene by construction.

Two edgeR companions were built, measured, and then deleted, both for the same reason: an edgeR kernel that is not a pure function of the gene it fits. estimateDisp diverged once one library was heavily over-sequenced. glmQLFit measured 1.6x and passed every small fixture, then failed on this cohort: mglmLevenberg records a deviance and an iteration count whose values depend on which genes share the block when a fit does not converge cleanly, 22 of 18,270 genes here, with no flag set on any of them. The assert in this report is what caught it, which is the reason every stage in this document carries one.

d. Speed

# Seconds only. Returning the fitted object too held a second copy of four results nothing
# reads, about 400 MB, resident at the moment the next line forks.
tick <- function(e) { t0 <- Sys.time(); force(e)
                      as.numeric(difftime(Sys.time(), t0, units = "secs")) }

# This sweep runs long enough to outlive whatever else the machine is doing, and an original
# timed once at the start would hand every later arm a denominator measured on a quieter machine.
# Each original is timed on both sides of the companion arms and the mean is the denominator, so
# drift moves both halves of a ratio together instead of inflating whichever arms ran last.
stages <- c("calcNormFactors (TMM)", "lmFit", "duplicateCorrelation")
pre <- c(tick(calcNormFactors(DGEList(ref))),
         tick(lmFit(vv, de_design)),
         tick(duplicateCorrelation(v_dc, de_design, block = plate)))

limma_speed <- do.call(rbind, lapply(tcga_workers, function(w) {
  p <- c(tick(calcNormFactors_parallel(DGEList(ref), workers = w)),
         tick(lmFit_parallel(vv, de_design, workers = w)),
         tick(duplicateCorrelation_parallel(v_dc, de_design, block = plate, workers = w)))
  invisible(gc(verbose = FALSE))
  data.frame(stage = stages, workers = w, companion = p)
}))

post <- c(tick(calcNormFactors(DGEList(ref))),
          tick(lmFit(vv, de_design)),
          tick(duplicateCorrelation(v_dc, de_design, block = plate)))

limma_speed$original <- rep((pre + post) / 2, length(tcga_workers))
limma_speed$speedup  <- limma_speed$original / limma_speed$companion
limma_drift <- max(abs(post - pre) / pre)
limma_speed[, c("stage", "workers", "original", "companion", "speedup")]

stage	workers	original	companion	speedup
calcNormFactors (TMM)	2	13.520948	3.521182	3.839889
lmFit	2	6.837098	2.393211	2.856872
duplicateCorrelation	2	657.336409	334.083638	1.967580
calcNormFactors (TMM)	4	13.520948	2.101084	6.435225
lmFit	4	6.837098	1.513995	4.515932
duplicateCorrelation	4	657.336409	209.517805	3.137377
calcNormFactors (TMM)	8	13.520948	1.611812	8.388664
lmFit	8	6.837098	2.050905	3.333698
duplicateCorrelation	8	657.336409	151.448302	4.340335

All three worker counts are shown, matching every other sweep in the report. Eight workers do not return twice what four return, and not because half of them are stranded on efficiency cores: given identical work, eight forked children here finish within 1.08x of each other, where stranding would show as roughly threefold. Concurrency is what costs, so each added worker returns less throughput than the one before it.

Originals were bracketed as in C.2, once before the companion arms and once after, and the two readings differ by at most 11.1%.

ggplot(limma_speed, aes(factor(workers), speedup, fill = stage)) +
  geom_col(position = position_dodge(0.8), width = 0.72) +
  geom_hline(yintercept = 1, linetype = 2, colour = "grey50") +
  geom_text(aes(label = sprintf("%.2fx", speedup)),
            position = position_dodge(0.8), vjust = -0.35, size = 3.1) +
  scale_fill_manual(values = c("calcNormFactors (TMM)" = "#4C72B0", "lmFit" = "#DD8452",
                               "duplicateCorrelation" = "#55A868")) +
  labs(x = "workers", y = "speedup against the original", fill = NULL,
       title = "limma and edgeR companions, TCGA",
       subtitle = sprintf("%s genes by %s samples. duplicateCorrelation on a %s gene subset. Identity checked at %d workers.",
                           format(nrow(ref), big.mark = ","), format(ncol(ref), big.mark = ","),
                           format(length(dc_genes), big.mark = ","), n_workers)) +
  expand_limits(y = 0)

Not everything in limma is worth forking, and the list of what was left alone is as much a part of the result as the list of what was not. voom measured 0.99x, because its lowess trend takes its span as a fraction of the gene count, so a block of half the genes fits a different curve and only the arithmetic after the trend splits. contrasts.fit is exactly splittable and the vendor returns before a fork could finish starting. eBayes, fitFDist and topTable pool across every gene by construction. None of them ships as a workers argument that does nothing.

E. Equivalence summary

# every pipeline stage below, computed from both implementations and compared. NA marks a
# stage that half does not have: the heatmap is drawn only for TCGA.
same <- function(a, b) if (is.null(a) || is.null(b)) NA else identical(a, b)

data.frame(
  stage = c("corrected counts", "storage mode", "PCA scores", "PCA variance explained",
            "log fold changes", "t statistics", "adjusted p-values",
            "significant gene set", "heatmap z-scores"),
  simulated = c(
    same(sim_ref, sim_par),
    same(storage.mode(sim_ref), storage.mode(sim_par)),
    same(q_rf$x, q_pr$x),
    same(q_rf$ve, q_pr$ve),
    same(sim_tt_ref$logFC, sim_tt_par$logFC),
    same(sim_tt_ref$t, sim_tt_par$t),
    same(sim_tt_ref$adj.P.Val, sim_tt_par$adj.P.Val),
    same(rownames(sim_tt_ref)[sim_tt_ref$adj.P.Val < 0.05],
         rownames(sim_tt_par)[sim_tt_par$adj.P.Val < 0.05]),
    NA),
  TCGA = c(
    same(ref, par),
    same(storage.mode(ref), storage.mode(par)),
    same(p_rf$x, p_pr$x),
    same(p_rf$ve, p_pr$ve),
    same(tt_ref$logFC, tt_par$logFC),
    same(tt_ref$t, tt_par$t),
    same(tt_ref$adj.P.Val, tt_par$adj.P.Val),
    same(rownames(tt_ref)[tt_ref$adj.P.Val < 0.05],
         rownames(tt_par)[tt_par$adj.P.Val < 0.05]),
    same(z_ref, z_par)))

stage	simulated	TCGA
corrected counts	TRUE	TRUE
storage mode	TRUE	TRUE
PCA scores	TRUE	TRUE
PCA variance explained	TRUE	TRUE
log fold changes	TRUE	TRUE
t statistics	TRUE	TRUE
adjusted p-values	TRUE	TRUE
significant gene set	TRUE	TRUE
heatmap z-scores	NA	TRUE

The companions that run after the correction, on the corrected counts from sections C and D:

data.frame(
  stage = c("calcNormFactors (TMM)", "lmFit", "duplicateCorrelation",
            "topTable, final gene list"),
  # read from the checks already computed, not recomputed: one identity, one place
  simulated = sim_limma_checks$identical,
  TCGA = c(limma_checks$identical[1], limma_checks$identical[2], dc_same,
           limma_checks$identical[3]))

stage	simulated	TCGA
calcNormFactors (TMM)	TRUE	TRUE
lmFit	TRUE	TRUE
duplicateCorrelation	TRUE	TRUE
topTable, final gene list	TRUE	TRUE

stopifnot(identical(sim_ref, sim_par),
          identical(q_rf$x, q_pr$x),
          identical(sim_tt_ref$adj.P.Val, sim_tt_par$adj.P.Val),
          identical(ref, par),
          identical(p_rf$x, p_pr$x), identical(p_rf$ve, p_pr$ve),
          identical(tt_ref$adj.P.Val, tt_par$adj.P.Val),
          identical(z_ref, z_par),
          all(limma_checks$identical), dc_same,
          all(sim_limma_checks$identical))

F. Speed summary

Every speed measurement in this report, one figure. Each bar is a companion timed against the original it replaces, on the same input. The ComBat-seq arms were compared against the reference at every worker count. The limma and edgeR arms are timed only; their identity was checked at 8 workers in C.2 and D.3.

all_speed <- rbind(
  data.frame(dataset = "simulated", stage = "ComBat_seq",
             workers = sim_timings$workers[sim_timings$workers > 1L],
             speedup = sim_timings$seconds[sim_timings$workers == 1L] /
                       sim_timings$seconds[sim_timings$workers > 1L]),
  data.frame(dataset = "TCGA", stage = "ComBat_seq",
             workers = tcga_workers, speedup = secs_ref / secs_w),
  data.frame(dataset = "simulated", stage = sim_limma_speed$stage,
             workers = sim_limma_speed$workers, speedup = sim_limma_speed$speedup),
  data.frame(dataset = "TCGA", stage = limma_speed$stage,
             workers = limma_speed$workers, speedup = limma_speed$speedup))
all_speed$stage <- factor(all_speed$stage,
  levels = c("ComBat_seq", "calcNormFactors (TMM)", "lmFit", "duplicateCorrelation"))

ggplot(all_speed, aes(factor(workers), speedup, fill = stage)) +
  geom_col(position = position_dodge2(0.8, preserve = "single"), width = 0.72) +
  geom_hline(yintercept = 1, linetype = 2, colour = "grey50") +
  geom_text(aes(label = sprintf("%.1fx", speedup)),
            position = position_dodge2(0.8, preserve = "single"), vjust = -0.35, size = 2.9) +
  facet_wrap(~ dataset, ncol = 1, scales = "free_x") +
  scale_y_continuous(expand = expansion(mult = c(0, 0.12))) +
  scale_fill_manual(values = c("ComBat_seq" = "#8172B3",
                               "calcNormFactors (TMM)" = "#4C72B0",
                               "lmFit" = "#DD8452",
                               "duplicateCorrelation" = "#55A868")) +
  labs(x = "workers", y = "speedup against the original", fill = NULL,
       title = "Every companion against its original",
       subtitle = "Dashed line is the original. Every arm at 2, 4 and 8 workers.") +
  expand_limits(y = 0) +
  theme(legend.position = "bottom")

G. Session

sessionInfo()

## R version 4.4.2 (2024-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS 26.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] ComplexHeatmap_2.22.0       patchwork_1.3.2            
##  [3] ggrepel_0.9.8               gtsummary_2.5.0            
##  [5] TCGAbiolinks_2.34.1         ggplot2_4.0.3              
##  [7] SummarizedExperiment_1.36.0 Biobase_2.66.0             
##  [9] GenomicRanges_1.58.0        GenomeInfoDb_1.42.3        
## [11] IRanges_2.40.1              S4Vectors_0.44.0           
## [13] BiocGenerics_0.52.0         MatrixGenerics_1.18.1      
## [15] matrixStats_1.5.0           rnaparallel_0.4.3          
## [17] edgeR_4.4.2                 limma_3.62.2               
## [19] sva_3.54.0                  BiocParallel_1.40.2        
## [21] genefilter_1.88.0           mgcv_1.9-4                 
## [23] nlme_3.1-169               
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3          shape_1.4.6.1              
##   [3] rstudioapi_0.18.0           jsonlite_2.0.0             
##   [5] magrittr_2.0.5              magick_2.9.1               
##   [7] farver_2.1.2                cardx_0.3.2                
##   [9] rmarkdown_2.30              GlobalOptions_0.1.4        
##  [11] fs_2.1.0                    zlibbioc_1.52.0            
##  [13] vctrs_0.7.3                 Cairo_1.7-0                
##  [15] memoise_2.0.1               htmltools_0.5.9            
##  [17] S4Arrays_1.6.0              progress_1.2.3             
##  [19] curl_7.1.0                  broom_1.0.12               
##  [21] SparseArray_1.6.2           sass_0.4.10                
##  [23] bslib_0.10.0                plyr_1.8.9                 
##  [25] httr2_1.2.2                 cachem_1.1.0               
##  [27] gt_1.3.0                    commonmark_2.0.0           
##  [29] lifecycle_1.0.5             iterators_1.0.14           
##  [31] pkgconfig_2.0.3             Matrix_1.7-5               
##  [33] R6_2.6.1                    fastmap_1.2.0              
##  [35] clue_0.3-68                 GenomeInfoDbData_1.2.13    
##  [37] digest_0.6.39               colorspace_2.1-2           
##  [39] AnnotationDbi_1.68.0        RSQLite_2.4.6              
##  [41] filelock_1.0.3              labeling_0.4.3             
##  [43] httr_1.4.8                  abind_1.4-8                
##  [45] compiler_4.4.2              bit64_4.8.0                
##  [47] withr_3.0.2                 doParallel_1.0.17          
##  [49] downloader_0.4.1            S7_0.2.2                   
##  [51] backports_1.5.1             DBI_1.3.0                  
##  [53] biomaRt_2.62.1              rappdirs_0.3.4             
##  [55] DelayedArray_0.32.0         rjson_0.2.23               
##  [57] tools_4.4.2                 otel_0.2.0                 
##  [59] glue_1.8.1                  cluster_2.1.8.2            
##  [61] generics_0.1.4              gtable_0.3.6               
##  [63] tzdb_0.5.0                  tidyr_1.3.2                
##  [65] data.table_1.18.4           hms_1.1.4                  
##  [67] xml2_1.5.2                  XVector_0.46.0             
##  [69] foreach_1.5.2               pillar_1.11.1              
##  [71] markdown_2.0                stringr_1.6.0              
##  [73] circlize_0.4.18             splines_4.4.2              
##  [75] dplyr_1.2.1                 BiocFileCache_2.14.0       
##  [77] lattice_0.22-9              survival_3.8-6             
##  [79] bit_4.6.0                   annotate_1.84.0            
##  [81] tidyselect_1.2.1            locfit_1.5-9.12            
##  [83] Biostrings_2.74.1           knitr_1.51                 
##  [85] litedown_0.9                xfun_0.57                  
##  [87] statmod_1.5.1               stringi_1.8.7              
##  [89] UCSC.utils_1.2.0            yaml_2.3.12                
##  [91] TCGAbiolinksGUI.data_1.26.0 evaluate_1.0.5             
##  [93] codetools_0.2-20            tibble_3.3.1               
##  [95] cli_3.6.6                   xtable_1.8-8               
##  [97] jquerylib_0.1.4             dichromat_2.0-0.1          
##  [99] Rcpp_1.1.1-1.1              dbplyr_2.5.2               
## [101] png_0.1-9                   XML_3.99-0.23              
## [103] parallel_4.4.2              readr_2.2.0                
## [105] blob_1.3.0                  prettyunits_1.2.0          
## [107] scales_1.4.0                purrr_1.2.2                
## [109] crayon_1.5.3                GetoptLong_1.1.1           
## [111] rlang_1.2.0                 KEGGREST_1.46.0            
## [113] rvest_1.0.5                 cards_0.7.1