Implementationstooling reviewed 2026-09
27 Statistics and Experimental Design: Implementations
Everything on this page is generated from tools/statistics.json. If something here is wrong or out of date, the fix belongs in that file, not in this prose.
For what these tools are doing and why the decisions are made before the analysis rather than during it, see Section 26.1.
The four sections follow the chapter’s four problems, in the order they bite: controlling the error rate across many tests, separating a batch effect from the biology, sizing the experiment before it is run, and testing a hypothesis that the data suggested. The order in which you will need them is the reverse of the order in which they are listed — the design section is the one that should be read first and is read last.
27.1 Controlling the error rate
Start with the one this page does not survey: p.adjust in base R implements Bonferroni, Holm and Benjamini–Hochberg in one line with no dependency, and for most analyses it is the correct and complete answer. Nothing below replaces it. They are all here because they estimate something p.adjust assumes.
IHW and qvalue each recover power that Benjamini–Hochberg gives away, from opposite directions. qvalue estimates what fraction of your tests were null instead of assuming all of them were; IHW estimates which tests had power in the first place and weights them. They can be combined, and on a large experiment either one alone typically returns noticeably more discoveries at the same false discovery rate.
Weights hypotheses by an independent covariate — mean expression is the usual one — so tests with more power to begin with are held to a looser threshold, and controls the false discovery rate over the whole set anyway. The covariate has to be independent of the p-value under the null while still being informative about power, and a filtering statistic such as mean count is exactly that, which is why the method fits genomics so naturally. (Ignatiadis et al. 2016)
The covariate is a modelling choice and has to be made before you look at results. Trying several and keeping whichever yields the most discoveries is a selection you have not accounted for anywhere.
Reports the quantity a genomics reader actually wants: for a list cut at any threshold, the expected fraction of it that is false. It also estimates the proportion of tests that are truly null rather than assuming that proportion is one, which is where its extra power over Benjamini-Hochberg comes from — on an experiment where a fifth of genes really do move, assuming otherwise discards a fifth of the discoveries you paid for. (Storey and Tibshirani 2003; Storey 2002)
The null-proportion estimate is read off the flat part of the p-value distribution, so it is only as good as that distribution. On a histogram that slopes or spikes near one, the estimate is not trustworthy and neither is anything downstream of it.
The table splits by what it is you doubt. ashr and fdrtool change the reported quantity — a false sign rate and a local false discovery rate respectively — because a rate computed over a whole list is not what you want when one row is about to become an experiment. bacon and swfdr change the model instead: bacon when the null itself has shifted, swfdr when the null fraction varies with a covariate. locfdr and multtest are here because you will meet them, not because a new analysis should start there.
| Tool | Upstream | Languages | Why / why not |
|---|---|---|---|
| ashr | maintainedlast release 2023-08last commit 2026-05 | R |
Shrinks effect sizes toward zero under a unimodal prior fitted from the data, and reports a local false sign rate rather than a false discovery rate. That change of quantity is the reason to use it: for an effect you intend to act on, the useful question is whether the sign is right, not whether the effect is exactly zero — which, at genome scale, it never is. (Stephens 2017) |
| bacon | dormantlast commit 2024-04 | R |
Fits the null distribution from the test statistics themselves and rescales against that, rather than against the standard normal the test assumed, which is Efron’s empirical null carried into association studies (Efron 2004). It is the right response to a QQ plot whose whole distribution has shifted — residual population structure or cell composition — because every method downstream of the p-values would otherwise inherit the shift and report it as discoveries. (Iterson et al. 2017) |
| fdrtool | maintainedlast release 2024-08 | R |
Local and tail-area false discovery rates from p-values, z-scores or correlations, with the null fitted rather than assumed. The local quantity is the reason to reach past qvalue: a q-value describes a whole list above a cut, a local false discovery rate describes the individual test, and when you are deciding which single gene to follow up those are different numbers. (Strimmer 2008) |
| swfdr | dormantlast commit 2020-12 | R |
Estimates the proportion of true nulls as a function of covariates rather than as a single number. It is the same insight as IHW’s — that tests are not exchangeable — arriving through the null fraction instead of through the threshold, and it takes several covariates at once, which IHW does not. Reach for it when the covariate structure is itself of interest. (Boca and Leek 2018) |
| locfdr | dormantlast release 2015-07 | R |
Efron’s own implementation of local false discovery rates and the empirical null, and the reference the later implementations were written against. The ideas are load-bearing throughout this section and worth reading here; the package is not where to start a new analysis, because fdrtool and bacon carry the same machinery with another decade of maintenance behind them. (Efron 2004) |
| multtest | active | R |
You will meet its |
WarningAdjusting the wrong set
Every procedure here takes a set of p-values and returns adjusted values for that set. What counts as the set is your decision and it is not a formality: adjusting within each of six contrasts separately, and then reporting the union of the six, controls nothing at the level of the union. Decide what the family is before you run the tests, and let it be the thing you will report.
27.2 Batch effects and unwanted variation
The most important tool for a batch effect is not on this page: it is the design matrix of whichever differential expression tool you were already using (Section 17.3). A recorded batch is a covariate, and putting it in the model is both the cheapest correction and the only one that leaves the downstream inference intact.
What the two cards add is the case where that is not enough. sva estimates the variation you did not record and hands it back as covariates for that same model — and, in ComBat, also offers the corrected matrix you need for anything that is not a model. variancePartition answers the question that should come first: how much of the variance is batch, and how much is the effect you care about. Running it before deciding on a correction is the difference between treating a batch effect and assuming one.
Two different jobs live in one package and it matters which you are doing. Surrogate variable analysis estimates unmeasured sources of variation and hands you covariates to put in your model; ComBat and ComBat-seq instead hand back a corrected data matrix with the batch removed. Prefer the covariates whenever a model is going to be fitted, because correcting first and testing afterwards produces intervals that are too narrow — the test is never told that the correction consumed degrees of freedom (Nygaard et al. 2016). (Leek et al. 2012; Leek and Storey 2007; Johnson et al. 2007; Zhang et al. 2020)
ComBat-seq is the entry point for counts. The original ComBat assumes continuous, roughly Gaussian data and will return numbers, not warnings, if you hand it raw counts.
Answers the question that belongs before any correction: what fraction of each gene’s variance is attributable to batch, to condition, and to everything else. That turns a batch effect from a worry into a measured quantity, and it is the honest way to discover that your effect of interest is the small term. Its dream model extends the same mixed model to repeated measures, where treating several samples from one donor as independent is the more common and more expensive mistake. (Hoffman and Schadt 2016; Hoffman and Roussos 2021)
A variance fraction is descriptive, not a test. A gene whose condition term explains little variance can still carry a real and reproducible effect, and the plot will not tell you which.
The two in the table are alternatives on different axes. InMoose replaces the language, and exists so that a Python pipeline does not have to shell out to R for one step. RUVSeq replaces the assumption: where sva takes the unwanted variation to be whatever is left after the model, RUVSeq takes it from control genes or replicate samples that should not have varied.
| Tool | Upstream | Languages | Why / why not |
|---|---|---|---|
| InMoose | activelast release 2026-01last commit 2026-01 | py |
The Python route to ComBat and ComBat-seq, ported from the R implementations and validated against them. It absorbed the earlier pyComBat package, and it exists for a real situation: a Python pipeline whose only R dependency was one batch-correction call. If you are already in R there is no reason to prefer it. (Colange et al. 2025; Behdenna et al. 2023) |
| RUVSeq | dormantlast commit 2022-11 | R |
Estimates unwanted variation from genes assumed not to respond to the condition, or from replicate samples that ought to be identical. It is the alternative to surrogate variables when you have such controls — spike-ins being the clean case — and its factors enter the design matrix the same way, so the choice between the two is about what evidence you have rather than about the model. (Risso et al. 2014; Gagnon-Bartsch and Speed 2012) |
WarningCorrected matrices leak
A matrix that has had a batch removed carries less uncertainty than the data it came from, and no downstream test knows that. Save the corrected matrix for plots, clustering and machine learning inputs, and fit models on the original with the batch term included — Section 26.5 is where that argument is made in full.
27.3 Designing the experiment
This is the thinnest section in the part and the honest reason is worth stating plainly: bulk RNA-seq power analysis is a solved calculation with no living software culture around it. Every bulk tool in the table below has been quiet for years, and none of them is quiet because it is broken — the negative binomial power calculation has not changed, so there has been nothing to fix. They are small, they still install, and the risk in using them is that their defaults encode assumptions from the sequencing depths of a decade ago rather than that their arithmetic has rotted.
Single-cell is the exception, and it is why the section has a card at all. The design decision there is genuinely unsolved arithmetic — how to split a fixed budget between donors, cells and depth — and scPower is the tool that treats it as the optimization it is.
The design question for a population-scale single-cell study is not a sample size but a budget split: more donors, more cells per donor, or more reads per cell. scPower models detection power as a function of all three and optimizes the split for a fixed cost, which is the decision an experimenter is actually facing and the one no single-number calculator can answer. (Schmid et al. 2021)
Its power estimates rest on expression and effect-size priors taken from a reference dataset, so the answer is only as transferable as that reference is to your tissue and protocol.
The four viable entries differ in where they get their variance estimate from, which is the only input in a power calculation you cannot simply assert. PROPER simulates from your own pilot data, RNASeqPower asks you to state it, RnaSeqSampleSize ships priors from public tumour cohorts, and ssizeRNA sizes against a false discovery rate rather than a per-test threshold. Pick on which of those you actually have.
| Tool | Upstream | Languages | Why / why not |
|---|---|---|---|
| PROPER | dormantlast commit 2020-11 | R |
Simulates from a real pilot dataset rather than from a parametric guess, then reports power across the whole gene set instead of for a representative gene. That is the right shape for the question whenever a pilot exists, because the distribution of dispersions it draws from is precisely what a closed-form calculation has to assume. (Wu et al. 2015) |
| RNASeqPower | active | R |
A closed-form sample size calculation for RNA-seq on the negative binomial, taking depth, coefficient of variation, effect size and the two error rates. Its virtue is that it is small enough to reason about — you can see which assumption is carrying the answer — and in a design calculation that matters more than sophistication, because every input is a guess you are making before the data exist. (Hart et al. 2013) |
| RnaSeqSampleSize | dormantlast commit 2022-02 | R |
Estimates sample size from a real dataset’s distribution of dispersions, and ships prior data for a long list of TCGA tissue types. The built-in priors are the reason to reach for it rather than for PROPER: they let you do the calculation before you have a pilot of your own, which is when the calculation is most needed and least often done. (Zhao et al. 2018) |
| ssizeRNA | maintainedlast release 2025-04 | R |
Sizes an RNA-seq experiment against a target false discovery rate rather than a per-test error rate, which is the quantity the finished analysis will actually be judged on. That makes it the natural companion to a calculation like RNASeqPower’s, which leaves the multiple testing burden to the caller. (Bi and Liu 2016) |
| powsimR | dormantlast release 2020-11last commit 2023-07 | R |
The most thorough simulation-based power framework for bulk and single-cell RNA-seq, spanning many combinations of normalization and differential expression method. It is legacy for an installation reason rather than a methodological one: it is distributed only from GitHub and imports 81 packages, 32 of which are themselves installed from GitHub rather than from a repository, so the chance that a fresh install builds falls a little further every year. (Vieth et al. 2017) |
NoteThe calculation you should do first
Before any of this: how many independent units — donors, animals, patients — will the experiment have? That number, not the number of cells or reads, sets what is detectable, and if it is three then no tool on this page will make the result generalize. Section 26.6 is the argument; this is the arithmetic that follows from it.
27.4 Inference after selection
The two cards attack the same problem from opposite ends. countsplit makes the selection legitimate: it manufactures two independent datasets from one set of counts, so clustering on the first and testing on the second is an ordinary analysis with ordinary guarantees. scSHC leaves the analysis alone and tests the selection itself, asking whether a proposed split of cells is more than one population’s noise.
Which you want depends on the question. If you need per-gene p-values after clustering, split the counts. If you need to know how many clusters there are, test the clustering. Both are recent, both are small packages, and both are doing something no amount of care with p-values downstream can substitute for.
Splits the counts rather than the samples: thinning a Poisson or negative binomial count yields two independent datasets from one, so you can cluster or select on the first and test on the second with the test’s assumptions intact. Unlike sample splitting, every cell contributes to both halves, which is what makes it usable at the sample sizes single-cell work actually has. (Neufeld et al. 2024)
Thinning relies on the count model being right. The negative binomial version needs an overdispersion estimate, and a bad one reintroduces exactly the dependence between halves that the method exists to remove.
Tests whether a proposed split of cells is more than what a single population would produce by chance, instead of assuming that whatever the clustering returned is real. It answers the question that follows every clustering — how many groups are there — with a hypothesis test rather than with a resolution parameter, and it can either build the clustering itself or evaluate one you already have. (Grabski et al. 2023)
Distributed from GitHub only: installation is a source build from the repository rather than a package-manager step, and there is no released version to pin beyond a commit.
| Tool | Upstream | Languages | Why / why not |
|---|---|---|---|
| selectiveInference | dormantlast release 2019-09last commit 2019-09 | R |
The reference implementation of exact post-selection inference for the lasso and forward stepwise regression, and the clearest place to see what an interval corrected for its own selection looks like. It is listed here for that reason rather than as a working recommendation: it is unmaintained, and its scope is regression model selection rather than the cluster-then-test problem that dominates genomics. (Tibshirani et al. 2016) |
27.5 What is not here
General-purpose statistics libraries — statsmodels, SciPy, pwr — are out of scope, by the same rule that keeps package managers and container runtimes out of the workflows part. They did not grow up in this field and are not dominated by it, so this book has no standing to date them and no way to notice when advice about them goes stale. That is a statement about the registry’s promise, not about the software, and withholding a record should not be a way of withholding a view. So, plainly: statsmodels.stats.multitest.multipletests is the right way to apply Bonferroni, Holm or Benjamini–Hochberg in Python, it is equivalent to p.adjust, and a Python analysis needs nothing else for the common cases. pwr is a faithful implementation of Cohen’s power calculations and is fine for a t-test or a proportion; it simply is not a genomics tool, and the genomics part of the problem — the threshold implied by twenty thousand tests — is yours to supply either way.
PEER is absent for the other reason, the one ANNOVAR established: in the field, and untrackable. It was the standard method for removing hidden confounders from expression data before eQTL mapping, and its factors are in the methods sections of a great many published analyses. Its repository has not moved since 2012, it has never tagged a version, and there is therefore no feed from which a version and a date can be read together — so a record would have to render either “unrefreshed”, which claims nobody looked, or a hand-typed version, which is the practice refresh.py exists to prevent. The book’s view, since an exclusion must not be a way of avoiding a judgment: do not build a new analysis on it. Surrogate variable analysis and RUV do the same job, are maintained, and hand you covariates for a model rather than a black box of factors; and where the point is simply to remove the top axes of unwanted variation, principal components of the expression matrix are a defensible and much simpler substitute.
Gene set enrichment and pathway analysis are not here, and the reason is that they are a different inference rather than another correction. The unit of the test is a set rather than a gene, the null hypothesis has two incompatible common forms — competitive and self-contained — and the dependence between overlapping sets is the hard part rather than a footnote. That is a subject the size of a part, not a section, and three rows here would misrepresent how much of it this page has checked.
27.6 Reading the upstream dates on this page
More records here carry thin evidence than anywhere else in the book, and the pattern is systematic enough to explain once rather than in eleven caveats.
Eleven of the twenty-one release through Bioconductor, which pins a version to its own twice-yearly platform release and publishes no per-package release date. Every such card therefore shows a commit date and no release date. Where even the commit date is missing — multtest and RNASeqPower — there is no development repository outside Bioconductor’s own git at all, and the status word on those two cards is asserted rather than derived. Treat it as the weakest claim on the page.
Bioconductor’s git does record a last-commit date per package, and it cannot be used, for a reason worth writing down: the platform bumps the version of every package at every release cycle, so that date has a floor of a few months for every package it carries, whether or not anyone has touched one. It can prove that work happened and can never show that it stopped, which is the wrong shape for a health signal.
One record is affected specifically. sva’s GitHub repository is a development mirror that stopped moving in 2020, while its Bioconductor history shows real work as recently as October 2025 — so the commit date beside it understates the project, and the dormant derived from that date should be read as a fact about the mirror. This is the branch-work problem in a new house: the work is on a different host rather than a different branch, and the fix is the same, which is prose rather than a heuristic.
27.7 Language coverage: the most one-sided part in the book
Twenty of the twenty-one records on this page are R packages. One is Python. There is no Julia and no Rust, and this is not an oversight in the survey.
Statistical methods for genomics were developed by statisticians, and the statistics community’s working language is R. Bioconductor is where those methods ship, its review process is a real filter, and the packages there inherit an ecosystem of shared data structures that a reimplementation would have to rebuild before it could be useful. The result is that for several things on this page — covariate-weighted multiple testing, surrogate variable analysis, count splitting — the R package is not the best implementation but the only one.
NoteWhat this means for a Python pipeline
The routine work is covered: statsmodels adjusts p-values, InMoose runs ComBat, and neither requires leaving Python. Beyond that, the honest answer is that you will be calling R — through rpy2, through a workflow step, or by writing the matrix out and reading the results back. Structuring for that from the beginning is much easier than discovering it after the analysis is written, and it is the same conclusion the alignment part reaches from the opposite direction, where every tool is C and Python is the caller.