RNA-Parallel

Parallel companions for RNA-seq tools. Each one runs the original function unmodified and returns output identical() to it, bit for bit. Same arguments, same defaults, same result, faster.

Read the rendered analysis

companion	complements	measured
ComBat_seq_parallel()	sva::ComBat_seq	5.40x
calcNormFactors_parallel()	edgeR::calcNormFactors	8.39x
lmFit_parallel()	limma::lmFit	4.52x
duplicateCorrelation_parallel()	limma::duplicateCorrelation	4.34x

Measured on the rendered TCGA analysis, 18,270 genes by 1,500 tumours, each arm identical() to the original it is timed against. Best of 2, 4 and 8 workers: every companion above peaks at 8 except lmFit, which peaks at 4 and is slower at 8. The full grid is in section 5. calcNormFactors clears the worker count because part of its gain is a serial fix to TMM's double rank() rather than parallelism.

The method is never reimplemented from a description of it. The original function is called with its hot paths rebound to row-parallel versions in a child of that function's own environment, so every other symbol still resolves to the original code and no ungated second copy exists to diverge from it. A companion earns its place only if it returns exactly what the function it replaces returns, so the test suite asserts identical() rather than a correlation or a tolerance.

Everything below documents the ComBat-seq companion, the one that ships today.

Equivalence	identical() at every pipeline stage, across every argument path, five execution backends, and chunk layouts from one to more than there are genes
Benchmark	5.40x on the TCGA correction at 8 workers, 8.39x on TMM normalisation
Interface	every sva::ComBat_seq argument, unchanged, plus four parallel controls
Method	ComBat-seq, run unmodified except for one pinned, byte-gated transcription of its hottest inner loop

ComBat_seq_parallel() accepts every sva::ComBat_seq() argument in the same order with the same defaults, and adds four parallel controls: workers, chunks, parallel_backend and backend.

This repository does not reimplement ComBat-seq, and holds one gated excerpt of it. The original function from zhangyuqing/ComBat-seq, distributed in Bioconductor's sva, is the one executed. It is called with six symbols rebound in a child of its own environment: glmFit, glmFit.default, match_quantiles and estimateGLMTagwiseDisp resolve to row-parallel implementations, and sapply dispatches the per-batch common-dispersion estimate across batches and lapply the per-batch tagwise estimate. Every other symbol resolves to ComBat-seq's own code.

One of the six is more than a split. match_quantiles runs a row-vectorised transcription of the sva 3.54.0 body rather than its cell loop, because that path is 66.5% of serial time. The vendor body is pinned as text and compared byte for byte on every call, and one changed character sends the work back to the vendor's own function. duplicateCorrelation_parallel lifts statmod's design-invariant QR and SVD out of the per-gene loop behind the same kind of gate. Those two excerpts are the only vendor code here, both are marked inline, and both stand down on any upstream change.

The parallel implementation is by Nguyen N (GenomeRx). The method is Zhang, Parmigiani and Johnson's, unmodified. Cite both.

1. Install

if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install(c("sva", "edgeR"))

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

2. Use

library(rnaparallel)
library(edgeR)   # DGEList, calcNormFactors
library(limma)   # voom, eBayes, topTable

# batch correction
adjusted <- ComBat_seq_parallel(counts, batch = batch, group = NULL, workers = 4L)

# a limma-voom differential expression pipeline
dge <- calcNormFactors_parallel(DGEList(counts), workers = 8L)
v   <- voom(dge, design)
fit <- lmFit_parallel(v, design, workers = 8L)
tt  <- topTable(eBayes(fit), coef = 2, number = Inf)

# a blocked design
cor <- duplicateCorrelation_parallel(v, design, block = subject, workers = 8L)
fit <- lmFit_parallel(v, design, block = subject,
                      correlation = cor$consensus.correlation, workers = 8L)

Each companion takes its vendor's arguments in the same order with the same defaults and adds the same four controls. Swapping lmFit for lmFit_parallel changes runtime and nothing else. The one exception is duplicateCorrelation_parallel, which requires block where the vendor defaults it to NULL; section 5a says why.

group, covar_mod, full_mod, shrink, shrink.disp and gene.subset.n behave exactly as they do in sva::ComBat_seq. chunks sets how many row blocks the genes are cut into, which trades fork overhead against peak memory per worker; backend wraps a sourced copy of upstream ComBat-seq instead of sva.

a. The ten-line check

set.seed(1)
counts <- matrix(rnbinom(1600, mu = 50, size = 5), nrow = 200)
batch  <- rep(1:2, each = 4)

identical(ComBat_seq_parallel(counts, batch, group = NULL, workers = 4L),
          sva::ComBat_seq(counts, batch, group = NULL))
#> TRUE

citation("rnaparallel") prints how to cite this and the method it runs.

3. Architecture

ComBat-seq is not rewritten. ComBat_seq_parallel() creates a child of the backend's own environment, binds six names in it, and calls the vendor function unchanged. Unbound symbols resolve to ComBat-seq's own code, so no ungated second copy of the algorithm exists to diverge.

Shares are of serial time on a 10,000 by 500 run with group = NULL.

ComBat-seq calls	share	execution	basis for exactness
match_quantiles	66.5%	by gene row	each output cell reads only its own gene
estimateGLMTagwiseDisp	14.0%	by gene row	valid at prior.df = 0, which collapses the term coupling genes. Other values go to edgeR unsplit
estimateGLMCommonDisp	13.1%	across batches	it sums over all genes, so a row split changes accumulation order. Batches are independent and RNG-free
glmFit, glmFit.default	remainder	by gene row	offset, dispersion, weights and start arrive explicitly and slice with the rows
monte_carlo_int_NB	small	serial	draws depend on the previous batch's state

Assembly is what makes the split safe rather than merely fast. Chunks are interleaved, so an expression-sorted matrix does not leave workers idle; each carries a tag and is reassembled by it, so a backend returning results out of order cannot bind genes to the wrong rows; and a dead worker, duplicate chunk, wrong chunk height or short result count halts the run instead of being patched. The backend loader refuses to start when a rebind target is no longer a bare call, since an unreachable rebind returns correct output at serial speed and no equivalence test can see it.

workers, chunks, parallel_backend and backend are added. Every sva::ComBat_seq() argument passes through unchanged.

ComBat-seq's speedup depends on whether a covariate is supplied, and the difference is large enough to state up front. With covar_mod the design has covariate columns and every hot path splits by gene row. With group = NULL and no covariates the design is batch-only, which is the one-group layout edgeR sends to a kernel that is not a pure function of the gene it fits, so neither glmFit nor the tagwise dispersion may be split that way.

They are split a different way instead. ComBat-seq computes the tagwise dispersion once per batch, and batches are independent of each other with no shared state and no RNG, so the work is dispatched across batches rather than across gene rows. That never reaches the kernel that makes a row split unsafe, and on the profile it is the stage that matters: at 6,000 genes by 400 samples across 10 batches the tagwise dispersion was 60% of the corrected run while the quantile match was 7%. Measured on 10,000 genes by 1,000 samples across 10 batches with group = NULL: 3.49x at 2 workers, 5.33x at 4, 5.97x at 8.

4. Benchmark

Two measurements against sva::ComBat_seq itself, not against this package held to one worker. In the simulated sweep each arm runs in a fresh R session through callr::r, so none inherits a warm cache or a cluster from the one before it. The TCGA arms share one session and run after the serial reference, so they run warm; the simulated sweep is the isolated measurement. Every arm returns a matrix identical() to the original.

Machine: Apple M3, 4 performance and 4 efficiency cores, 24 GB. workers is run as asked. Eight workers do not return twice the throughput of four, and the reason is not that half of them sit stranded on efficiency cores. Given byte-identical CPU-bound work, eight forked children finish within 1.08x of each other, where stranding half of them on a slower core would show up as roughly threefold: macOS migrates the forks across both clusters. What costs is concurrency itself. On repeated QR decompositions of a 500 by 500 matrix under mclapply, each of eight concurrent children ran 1.34x slower than one running alone, and eight workers returned 1.57x the aggregate throughput of four rather than 2x. The size of that penalty depends on the work, so treat it as the shape of the ceiling rather than a constant.

a. Simulated counts, no download

10,000 genes by 1,000 patients across 10 batches, each arm in a fresh R session through callr::r, from the rendered report.

implementation	seconds	speedup
sva::ComBat_seq	77.2	1.00x
ComBat_seq_parallel, 2 workers	33.7	2.29x
ComBat_seq_parallel, 4 workers	20.5	3.78x
ComBat_seq_parallel, 8 workers	16.2	4.76x

This arm passes group, so the planted biology is preserved and the design carries a condition column. A batch-only design, group = NULL with no covariates, routes its GLM fit and tagwise dispersion differently; section 3 describes that path and its measured curve.

b. TCGA

18,270 protein-coding genes, 1,500 primary tumours across 3 cancer types, 54 sequencing plates as the batch variable, smoking status and cancer type both preserved through covar_mod. The differential expression contrast is 430 current smokers against 210 lifelong non-smokers.

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

Equivalence is asserted at every stage: corrected counts, principal components, and the limma result including fold changes, adjusted p-values and the significant gene set. A divergence halts the render rather than producing a report with a quiet inconsistency in it.

Scaling is sublinear, for two reasons rather than one. estimateGLMCommonDisp cannot be split by gene row, so it is dispatched across batches instead, and the Monte Carlo integration and the imbalance between batches stay serial. Removing every serial fraction would still not give linear scaling: added workers return less throughput each, so the machine itself caps what any amount of parallelism can reach.

c. Reproducing the numbers

• Rendered report. The full run: cohort, argument parity, both corrections, PCA before and after, limma differential expression, worker sweep. Download and open it.
• run_example.R. Verifies the whole claim in a few minutes with no download. Simulates counts using the ComBat-seq paper's parameters, runs both implementations, and checks corrected counts, principal components, fold changes, adjusted p-values and the significant gene set are each identical():

# from a clone of this repository
Rscript inst/examples/run_example.R              # genes, samples and workers are optional

# or from the installed package, without cloning
Rscript "$(Rscript -e 'cat(system.file("examples/run_example.R", package="rnaparallel"))')"

• RNA_Parallel.Rmd. The source of that report. It needs no configuration: open it anywhere and knit. The first run downloads the HNSC, LUAD and LUSC cohorts from the GDC into the per-user cache directory and caches each as a per-project counts object; later runs read those in seconds. Two locations are searched, both chosen by the user: RNAPARALLEL_TCGA_DIR if it is set, otherwise the standard per-user cache directory. A download made under the package's previous name is still found, so renaming it did not cost anyone a second fetch.
• tests/. The identical() proof: every argument path, every chunk layout, every backend, worker failure, and dispersion exactness, plus a dispatch count through the public entry point, because identical() alone cannot tell a working parallel layer from a dead one. R CMD check runs fewer than devtools::test(): it skips the two cluster-reuse tests, and one further test runs only where future.apply is absent.

5. limma and edgeR

Same pattern, different vendor. The numbers are the rendered report's, measured on the ComBat-corrected TCGA matrix, 18,270 genes by 1,500 samples. Each original is timed twice, once before the companion arms and once after, and the mean is the denominator, so drift in machine load moves both halves of a ratio together instead of inflating whichever arm ran last. The two readings differed by at most 11.1% across this sweep.

companion	vendor	axis	2 workers	4 workers	8 workers
calcNormFactors_parallel()	edgeR::calcNormFactors	sample columns	3.84x	6.44x	8.39x
lmFit_parallel()	limma::lmFit	gene rows	2.86x	4.52x	3.33x
duplicateCorrelation_parallel()	limma::duplicateCorrelation	gene rows	1.97x	3.14x	4.34x

On simulated counts, 10,000 genes by 1,000 samples at 2, 4 and 8 workers, the same three companions measure 2.94x to 6.41x, 1.73x to 2.91x and 1.98x to 4.72x. duplicateCorrelation's serial arm is the report's longest single measurement, 657 seconds on a 4,000 gene subset, against 151 at 8 workers. Where a 4-worker arm beats the 8-worker one, the second four workers did not pay for themselves: added workers return less throughput each, and on these stages that does not cover the extra dispatch. Both are shown rather than the flattering one.

calcNormFactors_parallel() is two changes, not one. TMM computes rank(logR) and rank(absE) twice each per sample; computing them once is bit-identical and worth 1.58x on its own, with no workers involved. The column split sits on top of that. The axis is samples rather than genes because every normalisation method ranks or takes a median across genes, so the gene dimension is the one that must not be cut.

a. What is not parallelised, and why

This list is not an admission, it is the reason the output is identical. Each of these was measured or read out of the vendor source and then deliberately left alone.

function	reason
voom	The lowess mean-variance trend takes f = span as a fraction of the gene count, so a block of half the genes fits a different curve. Only the post-trend arithmetic splits, which measured 0.99x. Shipped as nothing rather than as a workers argument that does nothing.
eBayes, squeezeVar, fitFDist	fitFDist pools a median and two means across all genes. Also 0.09% of a pipeline, so there would be nothing to win even if it split.
topTable, decideTests	p.adjust is length-dependent for every method. A block would compute a different adjustment.
contrasts.fit	Exactly splittable, and the vendor takes 0.001 s. A fork costs more than the work.
arrayWeights	Reducing across genes is the function's entire purpose. There is no gene axis in its output to reassemble.
normalizeBetweenArrays, normalizeQuantiles	Rank and quantile across genes by construction.
lmFit(method = "robust"), ndups >= 2	unwrapdups reshapes the row axis, so a split pairs different genes. Refused with an error rather than left as a silent path.
duplicateCorrelation(block = NULL)	Same reason: the ndups path pairs rows through unwrapdups, and spacing == "topbottom" branches on the block's own row count. block is required here even though the vendor defaults it to NULL.
one-group designs, inside ComBat_seq_parallel	edgeR's mglmOneGroup kernel is not a pure function of the gene it is fitting when that fit fails to converge: the same gene's coefficient changes with which other genes share the matrix. The layout is detected and the vendor is called whole.
estimateDisp	Built, measured at 1.4x, and then deleted. adjustedProfileLik is 83.6% of it and looks exactly per-gene, but the assembled result was not identical() at stock defaults once one library was heavily over-sequenced: 19,999 of 20,000 tagwise dispersions moved. A modest speedup does not buy a companion that returns different numbers.
glmQLFTest	0.13 s against 6.12 s for estimateDisp. Nothing to win.
glmQLFit	Built, measured at 1.6x, passed every small fixture, and deleted after failing the render's own assert on real TCGA data. 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 there, with no flag set on any of them. No gate can see what the vendor does not mark.

b. Three ways a limma row split goes wrong

Each was reproduced against limma 3.62.2 before being guarded, and each returns a wrong answer quietly rather than failing.

asMatrixWeights dispatches on the block's row count, testing its gene branch before its array branch. A bare per-array weight vector is therefore read as per-gene weights by any block whose row count happens to equal the sample count. Measured: coefficients moved by 0.325. Weights are now expanded to a full matrix once, against the whole matrix, before anything is split.

NoProbeWts is an AND-reduction over every cell that selects between two numerically different algorithms which also return different component sets. An all-finite block flips to the fast path while the full matrix took the slow one. Measured: 114 of 400 sigma differed. The companion now proves the branch cannot flip before splitting, and runs serially when it cannot.

stats::lm.fit drops a one-column response to a vector, so a one-gene block swaps colMeans for mean and loses the gene name. Measured: 4,080 of 16,000 one-gene blocks were not identical(). Every chunk now carries at least two genes.

6. Tuning

workers is the primary parameter. The package default is four; six workers alongside a second forking R session has caused a kernel panic on a 24 GB machine.

parallel_backend selects the framework: "mclapply" (default, forks), "future", "BiocParallel", "foreach", "serial", or any function function(idx, f, workers). All return identical results. Forking is why mclapply is the default: a forked worker reads the count matrix through copy-on-write instead of being sent a copy, while socket-based frameworks ship the whole matrix per task.

Windows cannot fork, so the default mclapply backend and BiocParallel's MulticoreParam fall back to serial there and say so once per session. parallel_backend = "foreach" uses a PSOCK cluster and runs in parallel on Windows out of the box. "future" honours whatever plan the caller has set and warns when that plan resolves in one process; this package never sets a plan itself. options(combat.fork = FALSE) forces serial on any backend, the escape hatch if worker processes upset an IDE. combat_cluster_stop() releases cached clusters.

7. Citation

Cite both. This repository adds parallelism and contributes no statistics.

The method. Zhang Y, Parmigiani G, Johnson WE (2020). ComBat-seq: batch effect adjustment for RNA-seq count data. NAR Genomics and Bioinformatics 2(3), lqaa078. doi:10.1093/nargab/lqaa078

This companion. Nguyen N (2026). rnaparallel: Parallel Companions for RNA-Seq Tools. https://github.com/GenomeRx/RNA-Parallel. citation("rnaparallel") prints the current version and the methods to cite alongside it.

edgeR (Robinson, McCarthy and Smyth) does the GLM fitting underneath either way.

8. License

MIT for this companion, copyright GenomeRx. ComBat-seq itself is by Yuqing Zhang, Giovanni Parmigiani and W. Evan Johnson, distributed in Bioconductor's sva under Artistic-2.0, and this repository executes it rather than forking it.

Two files reproduce vendor source, in both cases so the companion can detect an upstream change and stand down. R/helper_seq_parallel.R carries the deparsed sva::match_quantiles body and a row-vectorised transcription of it, derived from Artistic-2.0 code. R/limma_dupcor_parallel.R carries the deparsed head of statmod::mixedModel2Fit and a transcription of that head, derived from code by Gordon Smyth licensed GPL-2 | GPL-3. Both blocks are marked inline. Nothing else here reproduces vendor code.