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 original is called unchanged on blocks of genes or
samples, with at most a symbol rebound in a child of the original’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, 6 and 8 workers, across every backend |
| 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 |
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", "BiocParallel"))
# future/future.apply, foreach/doParallel and BiocParallel are needed here and are not
# needed by the macOS or Linux reports: those sweep worker counts on one backend, this
# one sweeps the backends themselves.
cran <- c("ggplot2", "ggrepel", "patchwork", "callr", "gtsummary", "broom", "circlize",
"RColorBrewer", "future", "future.apply", "foreach", "doParallel")
need <- cran[!vapply(cran, requireNamespace, logical(1), quietly = TRUE)]
if (length(need)) install.packages(need)library(sva)
library(edgeR)
library(limma)
library(rnaparallel)
library(SummarizedExperiment)
library(ggplot2)
# Refuse a stale install before anything expensive runs. `git pull` updates this file and the
# package SOURCE; it does not reinstall the package, and nothing downstream would notice. The
# report stamps its own sessionInfo and keys every timing cache on packageVersion("rnaparallel"),
# so rendering against a previous install produces a report labelled with the wrong version that
# measured the wrong code, silently, after an hour. Skipped when the working tree is not to hand,
# which is the case for a reader rendering from an installed copy of the package.
.desc <- file.path("..", "..", "DESCRIPTION")
if (file.exists(.desc)) {
.want <- as.character(read.dcf(.desc, fields = "Version")[[1L]])
.have <- as.character(packageVersion("rnaparallel"))
if (!identical(.want, .have)) {
stop("this working tree is rnaparallel ", .want, " but the installed package is ", .have,
". A pull does not reinstall. From the repository root: R CMD INSTALL . ",
"(or devtools::install()), restart R, and render again.", call. = FALSE)
}
}
edgeR_norm <- if (exists("normLibSizes.default", envir = asNamespace("edgeR"), inherits = FALSE)) {
edgeR::normLibSizes
} else {
edgeR::calcNormFactors
}
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)))
# one colour per companion, shared by every bar figure below. Three figures each carried their
# own copy and removeBatchEffect was in none of them, so it drew as a grey NA series
pal_stage <- c("ComBat_seq" = "#8172B3",
"calcNormFactors (TMM)" = "#4C72B0",
"lmFit" = "#DD8452",
"duplicateCorrelation" = "#55A868",
"removeBatchEffect" = "#937860")
# workers is run as asked: 8 means 8 worker processes, bounded only by detectCores(). Added
# workers return less throughput each, so whether 8 beats 6 is what the sweep measures rather
# than assumes.
#
# WINDOWS SWEEP. 2, 4, 6, 8. Every platform's sweep brackets its PERFORMANCE core count, which
# is not its core count. This is an Ultra 9 185H: 16 physical cores, but only 6 of them are
# performance cores. The other ten are efficiency cores -- eight E plus two low-power E -- and
# they carry no second thread, which is how the split is readable at all: 22 logical threads
# against 16 physical cores means six cores with SMT, and on this family only P-cores have it.
#
# So a 16-worker arm is not sixteen equal workers. It is six fast ones and ten slow ones, and
# every dispatch waits on the slowest, which measures the efficiency cores rather than the
# scaling. Measured here: with the per-chunk payload sliced before dispatch, the curve stops
# improving between 4 and 8 workers and flattens, exactly where the P-cores run out.
#
# The macOS run used 2, 4, 6, 8 on an M3 with 4 performance cores, and the Linux run used
# 4, 8, 16, 32 on a dual Xeon whose 16 cores all carry SMT. Same rule, three answers.
# TCGA is corrected at every count; the last arm is the one the figures below use.
tcga_workers <- c(2L, 4L, 6L, 8L)
n_workers <- max(tcga_workers)
# BACKEND. This is the one thing that genuinely differs on Windows, and it is not a detail.
# There is no fork() here, so the default backend cannot do what it does on the other two
# platforms. The package handles that rather than failing, and each backend handles it
# differently:
#
# mclapply cannot fork, so it runs SERIALLY and says so once. Correct output, no speedup.
# BiocParallel the package asks for MulticoreParam, which BiocParallel itself substitutes
# with a serial param here. Also correct, also no speedup.
# future real parallelism over PSOCK, but only if the CALLER sets the plan. The package
# deliberately never touches a session's plan, so an unset plan resolves every
# future in this process. This report sets one, here and per arm below.
# foreach real parallelism over PSOCK, and it needs nothing from the caller: the package
# registers its own cluster when no backend is registered.
#
# Two of the four are serial on this platform and two are not, and the only thing that tells
# them apart is the clock. That is why C.1.a sweeps backends and not worker counts alone.
#
# The sweep in C.1.a is what chose the backend on the line below, and it did not choose the one
# the package's own documentation recommends here. `foreach` degrades as workers are added,
# rather than merely failing to improve: measured across a 2-to-16 sweep it fell from 1.18x to
# 0.28x, which at the wide end is three and a half times SLOWER than not parallelising at all.
# ComBat-seq dispatches 3 + 2 * n_batch times per correction, and the cached cluster is rebuilt
# whenever the requested width changes, which it does on every alternation between the
# batch-level and row-level dispatches. `future` holds one `multisession` set of workers across
# all of them and shows the ordinary shape instead. Both return identical() output; only the
# clock separates them, which is the whole reason this section exists. The table below is the
# current measurement, not this note -- read it there.
win_backends <- c("mclapply", "BiocParallel", "future", "foreach")
# Every backend in the sweep needs its framework present. The package does raise a clear
# error for a missing one, but it raises it in section C, after the simulation has already
# been built and cached. Cheaper to find out here.
backend_pkgs <- c("future", "future.apply", "BiocParallel", "foreach", "doParallel")
absent_backend <- backend_pkgs[!vapply(backend_pkgs, requireNamespace, logical(1),
quietly = TRUE)]
if (length(absent_backend))
stop("the backend sweep in C.1.a needs these packages and they are not installed: ",
paste(absent_backend, collapse = ", "), ". See the install chunk in section A.",
call. = FALSE)
options(combat.backend = "future")
options(future.globals.maxSize = 8 * 1024^3)
# `future` is the one backend whose parallelism the caller owns: the package will not touch a
# session's plan, and without one every future resolves in this process and the run is serial
# while the tables still say "workers". Set once here so the calls outside the sweeps are
# parallel, and reset per arm inside each sweep so `workers` means what the column says.
#
# Sizing also controls memory, which is the real ceiling here rather than cores. A forked
# worker shares the parent's pages; a multisession worker is a whole R process holding its own
# copy, measured at about 1.2 GB each on this cohort. Sixteen of them is 20 GB of 31 GB, so a
# plan is torn down and rebuilt at each arm's width rather than left at the maximum.
future::plan(future::multisession, workers = n_workers)
# Reported, not assumed. Every backend note above is about the absence of fork(), so a run of
# this file on a platform that has one would be describing something that is not happening.
on_windows <- identical(.Platform$OS.type, "windows")
if (!on_windows)
warning("This is the Windows report and this is not Windows. Section A describes the absence ",
"of fork(), which this platform has. Render RNA_Parallel.Rmd or ",
"RNA_Parallel_linux.Rmd instead.")
# Windows R ships reference BLAS, which is single-threaded, so there is no OPENBLAS_NUM_THREADS
# to pin as there is on Linux. A user-installed OpenBLAS or MKL build is multi-threaded, and
# then every PSOCK worker opens its own pool on top of the cluster. Thread counts are read when
# the library loads, so render_windows.ps1 sets them before R starts. This reports what happened.
blas_raw <- sessionInfo()$BLAS
blas_name <- if (is.null(blas_raw) || !nzchar(blas_raw)) "R internal (reference)" else basename(blas_raw)
blas_thread <- Sys.getenv("OPENBLAS_NUM_THREADS", Sys.getenv("MKL_NUM_THREADS", ""))
blas_pinned <- identical(blas_thread, "1")
if (grepl("openblas|mkl", blas_name, ignore.case = TRUE) && !blas_pinned)
warning("A multi-threaded BLAS (", blas_name, ") is loaded and its thread count is not ",
"pinned to 1. Each PSOCK worker will open its own pool, oversubscribing the machine ",
"and deflating every speedup below. Render with .\\render_windows.ps1")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
## 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 |
Fifty genes by eight samples sits far below
combat.min.cells, so every row-split dispatch runs
serially: starting a PSOCK worker 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.
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 <- edgeR_norm(DGEList(mat))
des <- model.matrix(~ group)
topTable(eBayes(lmFit(voom(d, des), des)), coef = 2, number = Inf, sort.by = "none")
}One correction, computed twice. Everything in this sub-section reads these two objects.
## 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
## 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
## [1] TRUE
# Release the pool before the next section. This matters far more here than on the other two
# platforms: a forked worker shares the parent's pages, so an idle pool costs almost nothing,
# while a PSOCK worker is a whole R process holding its own copy. Measured on this machine,
# sixteen idle workers held 20 GB of 31 GB and never released it, and free memory reached
# 1.4 GB during the section below -- which starts fresh sessions of its own on top of them.
combat_cluster_stop()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))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.
Every backend is swept here, not just the one the rest of the report
uses, because on Windows the backend decides whether anything runs in
parallel at all. mclapply and BiocParallel
both fall back to serial and are expected to land on the baseline.
future and foreach run over PSOCK and are the
two arms with something to measure. The fallbacks are run at the top
worker count only: repeating a serial run at four sizes measures this
machine’s noise, not the backend.
Every arm is checked against the reference whatever its speed. A backend that quietly ran serially still has to return the right answer, and that is the claim this report exists to support.
sim_file <- file.path(tempdir(), "sim-bench.rds")
saveRDS(list(counts = sim, batch = batch, group = group, ref = sim_ref), sim_file)
# The two that cannot parallelise here are run once, at the top worker count. Sweeping them
# would time the same serial call four times over.
serial_here <- c("mclapply", "BiocParallel")
sim_grid <- rbind(
data.frame(backend = "serial", workers = 1L, stringsAsFactors = FALSE),
data.frame(backend = serial_here, workers = n_workers, stringsAsFactors = FALSE),
expand.grid(backend = setdiff(win_backends, serial_here), workers = tcga_workers,
stringsAsFactors = FALSE, KEEP.OUT.ATTRS = FALSE)[, c("backend", "workers")])
sim_timings <- do.call(rbind, lapply(seq_len(nrow(sim_grid)), function(i) {
b <- sim_grid$backend[i]; w <- sim_grid$workers[i]
r <- callr::r(function(f, b, w) {
suppressMessages({library(sva); library(rnaparallel)})
if (b == "future") suppressMessages(library(future))
d <- readRDS(f)
# Both backends pay for their own workers, on the clock. `future` is the one whose
# parallelism the caller owns -- the package refuses to touch a session's plan, so without
# this the arm resolves every future in this process and times as serial while the table
# calls it parallel. Starting that plan BEFORE the timer, which is where it used to sit,
# handed `future` its worker startup for free while `foreach` built its cluster inside the
# timed call. That is a real cost on a platform with no fork(), and charging it to one
# backend and not the other is the difference between measuring a backend and flattering it.
t0 <- Sys.time()
if (b == "future") {
future::plan(future::multisession, workers = w)
on.exit(future::plan(future::sequential), add = TRUE)
}
out <- if (b == "serial") 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, parallel_backend = b)
# 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, b = b, w = w))
data.frame(backend = b, workers = w, seconds = r$secs, same = r$same,
stringsAsFactors = FALSE)
}))
unlink(sim_file)
# Hand SIGCHLD back to `parallel`. processx, which callr runs on, calls sigaction(SIGCHLD, ...)
# on every process it starts and restores SIG_DFL on teardown rather than the handler it
# displaced; `parallel` installs its own once, latched, and never reinstalls. So from the first
# callr call onward `parallel` can no longer observe its forked workers exiting, and every
# dispatch leaves one unreaped zombie per worker for the rest of the session. They hold a
# process-table slot each and nothing else -- no memory, no CPU -- but slots are the binding
# limit here, not memory: measured 493 zombies from a single eight-worker correction at this
# plate count, against a 4,000 per-user ceiling. This call clears parallel's latch so the next
# mcfork reinstalls the handler, and it returns immediately because nothing is stuck yet.
#
# Deliberately NOT inside a companion: with shutdown = TRUE it acts on every child past every
# cleanup mark, so it would kill a caller's own concurrent workers, which is exactly what
# combat_reap()'s spare snapshot exists to prevent. The report is the layer that mixes the two
# packages, so the report is where this belongs. Unix only; parallel::cleanup is in the fork
# half of the package and does not exist on Windows, where mclapply is serial and there are no
# forked children to lose in the first place.
if (.Platform$OS.type != "windows") {
ok <- tryCatch({ parallel:::cleanup(kill = TRUE, detach = TRUE, shutdown = TRUE); TRUE },
error = function(e) FALSE)
if (!ok) warning("could not hand SIGCHLD back to parallel after the callr sweep; forked ",
"workers will accumulate as zombies for the rest of this render")
}
sim_timings$implementation <- ifelse(
sim_timings$backend == "serial", "sva::ComBat_seq",
sprintf("ComBat_seq_parallel, %s, %d workers", sim_timings$backend, sim_timings$workers))
sim_base <- sim_timings$seconds[sim_timings$backend == "serial"][1]
sim_timings$speedup <- sim_base / sim_timings$seconds
# the baseline arm is the reference, so it has nothing to compare against
sim_timings$identical_to_original <- ifelse(sim_timings$backend == "serial", NA, sim_timings$same)
sim_timings[, c("implementation", "seconds", "speedup", "identical_to_original")]| implementation | seconds | speedup | identical_to_original |
|---|---|---|---|
| sva::ComBat_seq | 396.1835 | 1.0000000 | NA |
| ComBat_seq_parallel, mclapply, 8 workers | 323.4414 | 1.2249005 | TRUE |
| ComBat_seq_parallel, BiocParallel, 8 workers | 315.7031 | 1.2549245 | TRUE |
| ComBat_seq_parallel, future, 2 workers | 158.8915 | 2.4934217 | TRUE |
| ComBat_seq_parallel, foreach, 2 workers | 239.9220 | 1.6513017 | TRUE |
| ComBat_seq_parallel, future, 4 workers | 124.3928 | 3.1849401 | TRUE |
| ComBat_seq_parallel, foreach, 4 workers | 290.1784 | 1.3653101 | TRUE |
| ComBat_seq_parallel, future, 6 workers | 113.2315 | 3.4988808 | TRUE |
| ComBat_seq_parallel, foreach, 6 workers | 380.6239 | 1.0408792 | TRUE |
| ComBat_seq_parallel, future, 8 workers | 122.5176 | 3.2336863 | TRUE |
| ComBat_seq_parallel, foreach, 8 workers | 471.9253 | 0.8395048 | TRUE |
# sim_arms, not sim_par: `sim_par` is the corrected count matrix from the sim-correct chunk,
# and the end-to-end assertion in section E still compares it against sim_ref a thousand lines
# below. Reusing the name here overwrote the matrix with this data frame, so that assertion was
# comparing a matrix against a table of timings and failed with "identical(sim_ref, sim_par) is
# not TRUE" -- which reads like a correctness failure in the companion and is not one.
sim_arms <- sim_timings[sim_timings$backend != "serial", ]
sim_arms$backend <- factor(sim_arms$backend, levels = win_backends)
sim_best <- sim_arms[which.max(sim_arms$speedup), ]
ggplot(sim_arms, aes(factor(workers), speedup)) +
geom_col(aes(fill = speedup >= 1.05), width = 0.65) +
geom_hline(yintercept = 1, linetype = 2, colour = "grey40") +
geom_text(aes(label = sprintf("%.2fx\n%.0f s", speedup, seconds)), vjust = -0.3,
size = 3.5, fontface = "bold", lineheight = 0.95) +
facet_wrap(~ backend, nrow = 1, scales = "free_x") +
scale_fill_manual(values = c(`FALSE` = "grey60", `TRUE` = "#b2182b"), guide = "none") +
scale_y_continuous(expand = expansion(mult = c(0, 0.3))) +
labs(title = "Simulated counts: every backend against the original",
subtitle = sprintf(paste("%s genes by %s patients, %d batches. Dashed line is",
"sva::ComBat_seq at %.0f s. Every arm returns a matrix identical",
"to it, including the two that ran serially."),
format(G, big.mark = ","), format(ncol(sim), big.mark = ","),
nlevels(batch), sim_base),
x = "workers", y = "speedup vs sva::ComBat_seq")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(edgeR_norm(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.
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 <- edgeR_norm(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_rbe_v <- removeBatchEffect(sim_vv$E, batch = batch, design = sim_des)
sim_rbe_p <- removeBatchEffect_parallel(sim_vv$E, batch = batch, design = sim_des,
workers = n_workers)
sim_limma_checks <- data.frame(
stage = c("calcNormFactors (TMM)", "lmFit", "duplicateCorrelation",
"removeBatchEffect", "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_rbe_v, sim_rbe_p),
identical(sim_tt_v, sim_tt_p)))
sim_limma_checks| stage | identical |
|---|---|
| calcNormFactors (TMM) | TRUE |
| lmFit | TRUE |
| duplicateCorrelation | TRUE |
| removeBatchEffect | TRUE |
| topTable, final gene list | TRUE |
Scored against the planted truth, so the two implementations are compared on what they recover rather than only on whether their numbers match.
| 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_best <- function(f, reps = 3L) min(vapply(seq_len(reps), function(i) sim_secs(f()), numeric(1)))
sim_stages <- c("calcNormFactors (TMM)", "lmFit", "duplicateCorrelation", "removeBatchEffect")
sim_grid <- sort(unique(tcga_workers))
sim_pre <- c(sim_best(function() edgeR_norm(DGEList(sim_ref))),
sim_best(function() lmFit(sim_vv, sim_des)),
sim_secs(duplicateCorrelation(sim_vv[sim_dc_genes, ], sim_des, block = batch)),
sim_best(function() removeBatchEffect(sim_vv$E, batch = batch, design = sim_des)))
sim_limma_speed <- do.call(rbind, lapply(sim_grid, function(w) {
# sized to this arm, for the same reason as the correction sweep: under `future` the plan
# starts the processes and `workers` only caps the futures in flight
future::plan(future::multisession, workers = w)
p <- c(sim_best(function() calcNormFactors_parallel(DGEList(sim_ref), workers = w)),
sim_best(function() lmFit_parallel(sim_vv, sim_des, workers = w)),
sim_secs(duplicateCorrelation_parallel(sim_vv[sim_dc_genes, ], sim_des,
block = batch, workers = w)),
sim_best(function() removeBatchEffect_parallel(sim_vv$E, batch = batch,
design = sim_des, workers = w)))
invisible(gc(verbose = FALSE))
future::plan(future::sequential); combat_cluster_stop() # release; see the correction sweep
data.frame(stage = sim_stages, workers = w, companion = p)
}))
future::plan(future::multisession, workers = n_workers)
sim_post <- c(sim_best(function() edgeR_norm(DGEList(sim_ref))),
sim_best(function() lmFit(sim_vv, sim_des)),
sim_secs(duplicateCorrelation(sim_vv[sim_dc_genes, ], sim_des, block = batch)),
sim_best(function() removeBatchEffect(sim_vv$E, batch = batch, design = sim_des)))
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 | 9.164310 | 6.2498329 | 1.4663289 |
| lmFit | 2 | 1.782099 | 1.7296410 | 1.0303287 |
| duplicateCorrelation | 2 | 52.764947 | 34.2665510 | 1.5398383 |
| removeBatchEffect | 2 | 0.508590 | 0.5520260 | 0.9213152 |
| calcNormFactors (TMM) | 4 | 9.164310 | 8.3227260 | 1.1011188 |
| lmFit | 4 | 1.782099 | 1.9855082 | 0.8975529 |
| duplicateCorrelation | 4 | 52.764947 | 27.4164920 | 1.9245696 |
| removeBatchEffect | 4 | 0.508590 | 0.5199332 | 0.9781833 |
| calcNormFactors (TMM) | 6 | 9.164310 | 10.7566619 | 0.8519660 |
| lmFit | 6 | 1.782099 | 1.7076731 | 1.0435830 |
| duplicateCorrelation | 6 | 52.764947 | 30.8527460 | 1.7102188 |
| removeBatchEffect | 6 | 0.508590 | 0.5297902 | 0.9599838 |
| calcNormFactors (TMM) | 8 | 9.164310 | 12.8471289 | 0.7133353 |
| lmFit | 8 | 1.782099 | 1.9032190 | 0.9363603 |
| duplicateCorrelation | 8 | 52.764947 | 36.7351840 | 1.4363599 |
| removeBatchEffect | 8 | 0.508590 | 0.5505791 | 0.9237365 |
Every short stage is the best of three runs rather than one, because a stage lasting a second or two is defined by whatever else the machine did during it rather than by its own cost. Each original is also timed on both sides of the companion arms and the mean is the denominator; those two readings differ by at most 13.2% here, which is the drift the bracketing absorbs.
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 = pal_stage) +
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)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$statusgenes <- 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 |
library(gtsummary)
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 |
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 |
| label |
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.
## 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()
# The plan is rebuilt at this arm's width, inside the timer. Under `future` the `workers`
# argument does not start processes -- the caller's plan does, and the argument only caps how
# many futures are in flight -- so a plan left at one size would run all four arms at that
# size and the column would be decoration. Rebuilding also bounds memory: each multisession
# worker is a whole R process holding its own copy of the cohort, so the 16 arm must not be
# paid for while the 2 arm runs. Startup sits inside the clock because the caller really does
# pay it here, and it is cheap enough to include honestly: measured 0.3 to 0.7 s on this
# machine, flat in the worker count, against arms of several hundred seconds. What is NOT
# flat is the per-dispatch shipping cost, which is why the curve turns over -- a PSOCK worker
# receives a copy where a forked one would have shared the page.
future::plan(future::multisession, workers = tcga_workers[i])
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))
# Every arm starts from no pool, and no arm leaves one behind. On a forking platform this is
# tidiness; here it is what keeps the widest arm reachable at all, because each worker holds
# its own copy of the cohort rather than sharing the parent's pages. Both are released:
# `plan(sequential)` shuts down the multisession workers this arm started, and
# `combat_cluster_stop()` clears any cached cluster from a run on another backend. Timing is
# unaffected -- the clock above closes before either runs.
future::plan(future::sequential)
combat_cluster_stop()
}## 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
## 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
# leave a working plan behind for the sections that follow this sweep
future::plan(future::multisession, workers = n_workers)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)## [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")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 <- edgeR_norm(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 | 201 | 60 | 141 |
| ~ smoking + cancer | 76 | 60 | 16 |
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)
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,
# columns are split by the two groups being contrasted, so the panel shows the
# contrast rather than testing whether the groups separate on their own. Samples
# still cluster inside each group; genes cluster across the whole matrix.
column_split = droplevels(smoking[hm_keep]),
cluster_column_slices = FALSE,
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, columns split by smoking group",
top_n, top_n + hm_more)
else
sprintf("All %d significant genes, columns split by smoking group", 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 | 76 | 0 | -3.56 to 4.87 |
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.
# The macOS run read this with sysctl and the Linux run read sysfs and /proc. Neither exists
# here. CIM is the authoritative source on Windows, and every call is wrapped because a
# locked-down machine can refuse it and a missing CPU name is not worth failing a render over.
win_cim <- function(cls, prop) {
out <- tryCatch(suppressWarnings(system2("powershell",
c("-NoProfile", "-NonInteractive", "-Command",
sprintf("(Get-CimInstance %s | Select-Object -First 1 -ExpandProperty %s)", cls, prop)),
stdout = TRUE, stderr = FALSE)), error = function(e) character())
out <- trimws(out[nzchar(trimws(out))])
if (length(out)) out[1] else NA_character_
}
win_num <- function(cls, prop) suppressWarnings(as.numeric(win_cim(cls, prop)))
n_log <- parallel::detectCores()
n_phys <- win_num("Win32_Processor", "NumberOfCores")
if (is.na(n_phys)) n_phys <- suppressWarnings(parallel::detectCores(logical = FALSE))
if (is.na(n_phys) || n_phys < 1L) n_phys <- n_log
mem_gb <- win_num("Win32_ComputerSystem", "TotalPhysicalMemory") / 2^30
# The same thirteen rows as the macOS and Linux reports, read from this platform's own
# sources. Keeping the schema identical is what makes the three tables comparable side by
# side; only where each value comes from changes.
#
# performance cores are read from the topology rather than an original table, the same way the
# package's own default does it: on this family only P-cores carry SMT, so the number of
# logical processors above the physical count is the number of cores with a second thread.
# 22 logical against 16 physical gives 6, which is right, and the guard degrades to the
# physical count on a uniformly threaded or unthreaded machine.
smt <- n_log - n_phys
n_perf <- if (!is.na(smt) && smt > 0 && smt < n_phys) smt else n_phys
# Windows exposes no single NUMA-node count that is present on every edition, so this asks
# the OS and falls back to one rather than erroring or printing NA. A laptop is one node; the
# row exists so the three tables line up, and it is the Linux report that has more than one.
win_numa <- function() {
n <- tryCatch(suppressWarnings(as.numeric(win_cim("Win32_ComputerSystem",
"NumberOfProcessors"))),
error = function(e) NA_real_)
if (is.na(n) || n < 1) 1L else as.integer(n)
}
data.frame(
quantity = c("CPU", "performance cores", "physical cores", "logical CPUs",
"threads per core", "NUMA nodes", "memory", "OS",
"BLAS", "BLAS pinned to 1 thread", "fork available", "backend",
"worker counts"),
value = as.character(c(
win_cim("Win32_Processor", "Name"),
n_perf,
n_phys,
n_log,
sprintf("%.0f", n_log / n_phys),
win_numa(),
if (is.na(mem_gb)) NA_character_ else sprintf("%.0f GB", mem_gb),
paste(Sys.info()[["sysname"]], Sys.info()[["release"]]),
blas_name,
blas_pinned,
!identical(.Platform$OS.type, "windows"),
getOption("combat.backend"),
paste(tcga_workers, collapse = ", "))))| quantity | value |
|---|---|
| CPU | Intel(R) Core(TM) Ultra 9 185H |
| performance cores | 6 |
| physical cores | 16 |
| logical CPUs | 22 |
| threads per core | 1 |
| NUMA nodes | 1 |
| memory | 31 GB |
| OS | Windows 10 x64 |
| BLAS | R internal (reference) |
| BLAS pinned to 1 thread | TRUE |
| fork available | FALSE |
| backend | future |
| worker counts | 2, 4, 6, 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 | 2573.3 | 1.00x | NA |
| ComBat_seq_parallel, 2 workers | 1215.5 | 2.12x | TRUE |
| ComBat_seq_parallel, 4 workers | 770.0 | 3.34x | TRUE |
| ComBat_seq_parallel, 6 workers | 720.1 | 3.57x | TRUE |
| ComBat_seq_parallel, 8 workers | 727.6 | 3.54x | 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")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.
nf_v <- edgeR_norm(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")
rbe_v <- removeBatchEffect(vv$E, batch = plate, design = de_design)
rbe_p <- removeBatchEffect_parallel(vv$E, batch = plate, design = de_design,
workers = n_workers)
limma_checks <- data.frame(
stage = c("calcNormFactors (TMM)", "lmFit", "removeBatchEffect",
"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(rbe_v, rbe_p),
identical(tt_limma_v, tt_limma_p), identical(tt_limma_v, tt_ref)))
limma_checks| stage | identical |
|---|---|
| calcNormFactors (TMM) | TRUE |
| lmFit | TRUE |
| removeBatchEffect | TRUE |
| topTable, final gene list | TRUE |
| matches the gene list D.2 already drew | TRUE |
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.
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.003267 |
| identical() | TRUE |
Timing for this stage is in D.3.d, where the baseline is bracketed rather than read once.
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 PSOCK worker 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.
# 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 dispatches to the workers.
tick <- function(e) { t0 <- Sys.time(); force(e)
as.numeric(difftime(Sys.time(), t0, units = "secs")) }
# The BEST of three, not one reading and not their mean. A stage that takes a second or two is
# not robust to whatever else the machine decides to do during it: one render timed lmFit at
# four workers at 14.19 s against 1.70 s for the identical call minutes earlier, and reported
# it as 0.39x, slower than serial. The minimum is the reading least contaminated by load,
# which is the quantity this table is actually claiming. Long stages are left at one reading:
# duplicateCorrelation runs eleven minutes serial, long enough to absorb a transient spike
# rather than be defined by one, and three of them would cost an hour to say the same thing.
best <- function(f, reps = 3L) min(vapply(seq_len(reps), function(i) tick(f()), numeric(1)))
# 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", "removeBatchEffect")
pre <- c(best(function() edgeR_norm(DGEList(ref))),
best(function() lmFit(vv, de_design)),
tick(duplicateCorrelation(v_dc, de_design, block = plate)),
best(function() removeBatchEffect(vv$E, batch = plate, design = de_design)))
limma_speed <- do.call(rbind, lapply(tcga_workers, function(w) {
future::plan(future::multisession, workers = w) # sized to this arm; see the correction sweep
p <- c(best(function() calcNormFactors_parallel(DGEList(ref), workers = w)),
best(function() lmFit_parallel(vv, de_design, workers = w)),
tick(duplicateCorrelation_parallel(v_dc, de_design, block = plate, workers = w)),
best(function() removeBatchEffect_parallel(vv$E, batch = plate,
design = de_design, workers = w)))
invisible(gc(verbose = FALSE))
# same reason as the correction sweep above: a pool left resident is a whole R process per
# worker, and the `post` timings below would otherwise be measured against a machine still
# holding the widest arm's workers
future::plan(future::sequential)
combat_cluster_stop()
data.frame(stage = stages, workers = w, companion = p)
}))
future::plan(future::multisession, workers = n_workers)
post <- c(best(function() edgeR_norm(DGEList(ref))),
best(function() lmFit(vv, de_design)),
tick(duplicateCorrelation(v_dc, de_design, block = plate)),
best(function() removeBatchEffect(vv$E, batch = plate, design = de_design)))
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 | 26.530530 | 15.248802 | 1.7398435 |
| lmFit | 2 | 10.125268 | 9.561065 | 1.0590105 |
| duplicateCorrelation | 2 | 784.751489 | 433.858971 | 1.8087709 |
| removeBatchEffect | 2 | 4.399733 | 4.895648 | 0.8987029 |
| calcNormFactors (TMM) | 4 | 26.530530 | 17.229926 | 1.5397936 |
| lmFit | 4 | 10.125268 | 9.196090 | 1.1010406 |
| duplicateCorrelation | 4 | 784.751489 | 379.302855 | 2.0689311 |
| removeBatchEffect | 4 | 4.399733 | 4.942493 | 0.8901850 |
| calcNormFactors (TMM) | 6 | 26.530530 | 24.384648 | 1.0880013 |
| lmFit | 6 | 10.125268 | 9.695114 | 1.0443681 |
| duplicateCorrelation | 6 | 784.751489 | 475.624918 | 1.6499377 |
| removeBatchEffect | 6 | 4.399733 | 5.543204 | 0.7937166 |
| calcNormFactors (TMM) | 8 | 26.530530 | 30.323008 | 0.8749307 |
| lmFit | 8 | 10.125268 | 10.684164 | 0.9476893 |
| duplicateCorrelation | 8 | 784.751489 | 560.474884 | 1.4001546 |
| removeBatchEffect | 8 | 4.399733 | 5.924219 | 0.7426688 |
All four worker counts are shown, matching every other sweep in the
report. These arms run on foreach, the backend section A
selected, so they measure a PSOCK cluster rather than forked children:
every worker is a fresh R process that has to be handed its data over a
socket before it can start. That cost is paid per chunk and does not
shrink as workers are added, which is why this curve flattens earlier
than the fork curves in the macOS and Linux reports, and can turn over
entirely. A row below 1.00x means serialising the input cost more than
the work it saved. It does not mean the companion misbehaved: the
identity checks above passed on these same arms.
Method as in C.2: calcNormFactors and lmFit
are the best of three runs, and originals are bracketed before and after
the companion arms, those two readings differing by at most 10.4%.
duplicateCorrelation is a single reading at each point,
because eleven minutes serial is long enough to absorb a transient spike
rather than be defined by one.
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 = pal_stage) +
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 dispatching, 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 original
returns before a PSOCK worker could finish starting.
eBayes, fitFDist and topTable
pool across every gene by construction. None of them ships as a
workers argument that does nothing.
# 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:
# Both columns are read from the checks already computed, never recomputed: one identity, one
# place. They are also indexed BY NAME rather than by position, because adding a stage to either
# table silently shifted the positional indices and bound the wrong TRUE to the wrong row.
pick <- function(tbl, nm) tbl$identical[match(nm, tbl$stage)]
sum_stages <- c("calcNormFactors (TMM)", "lmFit", "duplicateCorrelation",
"removeBatchEffect", "topTable, final gene list")
data.frame(
stage = sum_stages,
simulated = pick(sim_limma_checks, sum_stages),
TCGA = c(pick(limma_checks, "calcNormFactors (TMM)"), pick(limma_checks, "lmFit"), dc_same,
pick(limma_checks, "removeBatchEffect"),
pick(limma_checks, "topTable, final gene list")))| stage | simulated | TCGA |
|---|---|---|
| calcNormFactors (TMM) | TRUE | TRUE |
| lmFit | TRUE | TRUE |
| duplicateCorrelation | TRUE | TRUE |
| removeBatchEffect | 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))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 simulated ComBat-seq bars are the foreach arms
of the C.1.a backend sweep, so that every bar in this figure is the same
backend at cohort and simulated scale. 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$backend == "foreach"],
speedup = sim_timings$speedup[sim_timings$backend == "foreach"]),
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 = names(pal_stage))
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 = pal_stage) +
labs(x = "workers", y = "speedup against the original", fill = NULL,
title = "Every companion against its original",
# read from tcga_workers rather than spelled out: this line said "2, 4, 8 and 16" for a
# sweep that no longer runs, and a hardcoded caption goes stale silently
subtitle = sprintf("Dashed line is the original. Every arm at %s workers.",
paste(tcga_workers, collapse = ", "))) +
expand_limits(y = 0) +
theme(legend.position = "bottom")## R version 4.6.1 (2026-06-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## 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.28.0 patchwork_1.3.2
## [3] ggrepel_0.9.8 gtsummary_2.6.0
## [5] TCGAbiolinks_2.40.0 ggplot2_4.0.3
## [7] SummarizedExperiment_1.42.0 Biobase_2.72.0
## [9] GenomicRanges_1.64.0 Seqinfo_1.2.0
## [11] IRanges_2.46.0 S4Vectors_0.50.1
## [13] BiocGenerics_0.58.1 generics_0.1.4
## [15] MatrixGenerics_1.24.0 matrixStats_1.5.0
## [17] rnaparallel_0.5.0 edgeR_4.10.1
## [19] limma_3.68.4 sva_3.60.0
## [21] BiocParallel_1.46.0 genefilter_1.94.0
## [23] mgcv_1.9-4 nlme_3.1-169
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 shape_1.4.6.1
## [3] jsonlite_2.0.0 magrittr_2.0.5
## [5] farver_2.1.2 cardx_0.3.4
## [7] rmarkdown_2.31 GlobalOptions_0.1.4
## [9] fs_2.1.0 vctrs_0.7.3
## [11] memoise_2.0.1 htmltools_0.5.9
## [13] S4Arrays_1.12.0 progress_1.2.3
## [15] curl_7.1.0 broom_1.0.13
## [17] SparseArray_1.12.2 sass_0.4.10
## [19] parallelly_1.48.0 bslib_0.12.0
## [21] plyr_1.8.9 httr2_1.3.0
## [23] cachem_1.1.0 gt_1.3.0
## [25] commonmark_2.0.0 lifecycle_1.0.5
## [27] iterators_1.0.14 pkgconfig_2.0.3
## [29] Matrix_1.7-5 R6_2.6.1
## [31] fastmap_1.2.0 clue_0.3-68
## [33] future_1.75.0 digest_0.6.39
## [35] colorspace_2.1-3 AnnotationDbi_1.74.0
## [37] RSQLite_3.53.3 filelock_1.0.3
## [39] labeling_0.4.3 httr_1.4.8
## [41] abind_1.4-8 compiler_4.6.1
## [43] bit64_4.8.4 withr_3.0.3
## [45] doParallel_1.0.17 downloader_0.4.1
## [47] S7_0.2.2 backports_1.5.1
## [49] DBI_1.3.0 biomaRt_2.68.0
## [51] DelayedArray_0.38.2 rjson_0.2.23
## [53] tools_4.6.1 otel_0.2.0
## [55] future.apply_1.20.2 glue_1.8.1
## [57] callr_3.8.0 cluster_2.1.8.2
## [59] gtable_0.3.6 tzdb_0.5.0
## [61] tidyr_1.3.2 data.table_1.18.4
## [63] hms_1.1.4 xml2_1.6.0
## [65] XVector_0.52.0 foreach_1.5.2
## [67] pillar_1.11.1 markdown_2.0
## [69] stringr_1.6.0 circlize_0.4.18
## [71] splines_4.6.1 dplyr_1.2.1
## [73] BiocFileCache_3.2.0 lattice_0.22-9
## [75] survival_3.8-6 bit_4.6.0
## [77] annotate_1.90.0 tidyselect_1.2.1
## [79] locfit_1.5-9.12 Biostrings_2.80.1
## [81] knitr_1.51 litedown_0.11
## [83] xfun_0.60 statmod_1.5.2
## [85] stringi_1.8.9 yaml_2.3.12
## [87] TCGAbiolinksGUI.data_1.32.0 evaluate_1.0.5
## [89] codetools_0.2-20 tibble_3.3.1
## [91] cli_3.6.6 xtable_1.8-8
## [93] processx_3.9.0 jquerylib_0.1.4
## [95] Rcpp_1.1.2 globals_0.19.1
## [97] dbplyr_2.6.0 png_0.1-9
## [99] XML_3.99-0.24 parallel_4.6.1
## [101] readr_2.2.0 blob_1.3.0
## [103] prettyunits_1.2.0 listenv_1.0.0
## [105] scales_1.4.0 purrr_1.2.2
## [107] crayon_1.5.3 GetoptLong_1.1.1
## [109] rlang_1.3.0 KEGGREST_1.52.2
## [111] rvest_1.0.5 cards_0.9.0