No amount of sophisticated modelling rescues a confounded design. The most important statistics in a mass spectrometry study are decided before a single sample enters the instrument.
Every analysis chapter so far assumed the data already existed. This chapter steps back to the decision that determines whether an experiment can answer its question at all: how it is designed. Get the design right — enough replication, no confounding, randomised run order — and the modelling in the chapters that follow is straightforward. Get it wrong, and no normalization, imputation, or empirical-Bayes moderation can recover a signal that was never separable from the noise.
This chapter is deliberately placed at the start of the statistical-modelling part: it is the bridge between how MS data is generated (Parts I–IV) and how it is analysed (Chapters 20–23).
Confounding your condition with batch or run order. If every case runs on Monday and every control on Tuesday, no statistical method can separate biology from day. Randomize and block before you acquire a single sample.
By the end of this chapter you will be able to:
Every measured intensity is a sum of signal and several kinds of noise. Good design maximises the first relative to the rest.
| Source | Example in MS | Design lever |
|---|---|---|
| Biological | Genuine between-subject differences | This is your signal — sample enough subjects |
| Technical | Sample prep, digestion, injection variability | Technical replicates, internal standards |
| Batch | Reagent lots, columns, instrument recalibration, acquisition day | Randomisation, blocking, QC pools |
| Run-order / drift | Sensitivity decline across a long sequence | Randomise injection order; interleave QCs |
The central design question is not “how do I remove noise?” but “how do I keep noise from aligning with my factor of interest?” Noise that is random can be modelled; noise that is confounded with the biology is fatal.
Treating technical replicates as if they were biological replicates is pseudoreplication — it inflates the apparent sample size and produces anti-conservative p-values. In limma, model technical replicates with duplicateCorrelation() and block on the biological unit (Chapter 21); never let them count as independent samples.
The honest answer is it depends on the effect size and variance — which is exactly what the power analysis below quantifies. As rough field guidance for discovery MS:
| Study type | Typical minimum per group | Notes |
|---|---|---|
| Cell-line proteomics (low variance) | 3–4 | Controlled system; small biological variability |
| Tissue / animal proteomics | 5–8 | Higher biological variance |
| Human clinical cohorts | 20+ per group | Large inter-individual variability; often imbalanced |
| Untargeted metabolomics (human) | 30+ per group | High variance + many features to correct over |
These are starting points, not substitutes for a power calculation with your own pilot variance.
Two samples measured back-to-back are more alike than two measured a day apart. If all your cases run on Monday and all controls on Tuesday, day and condition are confounded and no analysis can separate them.
# Randomise injection order while keeping QC pools at fixed positions
set.seed(2024)
samples <- paste0("S", 1:40)
run_order <- sample(samples) # randomised biological samples
# Insert a QC pool every 10 injections
schedule <- append(run_order,
values = rep("QC_pool", 4),
after = c(0, 10, 20, 30))Before acquisition, tabulate your factor of interest against every technical factor. Any strong association is a red flag.
# Is condition balanced across batch? (want roughly uniform counts)
table(metadata$condition, metadata$batch)
# Is condition correlated with run order? (want ~0)
cor(as.integer(factor(metadata$condition)), metadata$run_order)Statistical power is the probability of detecting a true effect of a given size. Planning for adequate power is what separates a study that can succeed from one that is underpowered before it begins. Power depends on four quantities; fix any three and the fourth follows:
For a single feature, a two-sample t-test power calculation gives a first estimate. The effect size (Cohen’s d) is the log₂ fold-change divided by the standard deviation.
library(pwr)
log2fc <- 1.0 # detect a two-fold change
sd_est <- 0.8 # residual SD on the log2 scale (from pilot data)
d <- log2fc / sd_est
pwr::pwr.t.test(d = d, sig.level = 0.05, power = 0.80,
type = "two.sample") # -> required n per groupA single-feature calculation is optimistic: an MS experiment tests thousands of features and controls the false discovery rate, so the effective per-feature threshold is far stricter than 0.05. A practical adjustment is to target the Bonferroni-style level for the number of features you expect to test.
n_features <- 5000
target_fdr <- 0.05
approx_alpha <- target_fdr / n_features # conservative planning threshold
pwr::pwr.t.test(d = 1.0 / 0.8, sig.level = approx_alpha,
power = 0.80, type = "two.sample")Analytical formulas ignore MS realities — non-normal intensities, missing values, and the limma moderation you will actually use. A simulation that mimics your planned analysis gives a far more honest estimate: generate data under a realistic model, run the real pipeline, and count how often true effects are recovered.
library(limma)
simulate_power <- function(n_per_group, n_features = 5000, n_de = 200,
log2fc = 1, sigma = 0.8, n_sim = 100) {
hits <- numeric(n_sim)
for (s in seq_len(n_sim)) {
group <- factor(rep(c("A", "B"), each = n_per_group))
design <- model.matrix(~ group)
# Null features + a block of truly differential features
mat <- matrix(rnorm(n_features * 2 * n_per_group, sd = sigma),
nrow = n_features)
mat[1:n_de, group == "B"] <- mat[1:n_de, group == "B"] + log2fc
fit <- eBayes(lmFit(mat, design))
padj <- p.adjust(fit$p.value[, 2], method = "BH")
# Power = fraction of true DE features recovered at FDR 5%
hits[s] <- mean(padj[1:n_de] < 0.05)
}
mean(hits)
}
# Power curve across candidate sample sizes
sizes <- c(3, 5, 8, 10, 15, 20)
power_curve <- vapply(sizes, simulate_power, numeric(1))
data.frame(n_per_group = sizes, power = round(power_curve, 2))Plot power against sample size and choose the smallest n that clears your target (commonly 0.80). Because the simulation uses the same limma pipeline as your real analysis, its estimate already accounts for empirical-Bayes moderation.
Before acquisition:
limma pipeline over analytical shortcuts.The next chapter takes a well-designed dataset and builds the statistical model that extracts its differentially abundant features.
limma (Chapter 21)?pwr, find the per-group sample size needed to detect a log₂FC of 1.5 with residual SD 0.9 at 80% power, first at α = 0.05 and then at a Bonferroni threshold for 4,000 features. Compare.simulate_power() function to include 30% random missingness (set values to NA before fitting) and re-draw the power curve. How much does missingness raise the required n?sessionInfo()