bio-experimental-design-power-analysis

What this skill does

## Version Compatibility

Reference examples tested with: RNASeqPower 1.42+, PROPER 1.34+, powsimR 1.2+ (GitHub), DESeq2 1.42+, edgeR 4.0+, pwr 1.3+.

Before using code patterns, verify installed versions match. If versions differ:
- R: `packageVersion('<pkg>')` then `?function_name` to verify parameters

If code throws an error, introspect the installed package and adapt to the actual API. Notes: `RNASeqPower::rnapower()` solves for whichever of `n` or `power` is omitted; PROPER is a multi-step pipeline (`RNAseq.SimOptions.2grp` -> `simRNAseq` -> `runSims` -> `comparePower`); powsimR is GitHub-only and its `estimateParam`/`Setup`/`simulateDE` signatures drift — pin a commit SHA for reproducible work. Verify each against the installed help before relying on argument names.

# Power Analysis for Genomics Experiments

**"How many replicates does my sequencing experiment need?"** -> Compute the probability of detecting a biologically meaningful effect given replicate number, sequencing depth, and biological variability — modeling counts as negative-binomial and recognizing that power is a per-gene quantity, not one number for the whole transcriptome.
- R: `RNASeqPower::rnapower()` — closed-form NB power/sample size; `PROPER`, `powsimR` — simulation from the mean-dispersion trend

## The Single Most Important Modern Insight -- Genomics Power Is Per-Gene; Simulate, and Never Report Observed Power

Power in a sequencing experiment is not a single number. It is a per-gene quantity that depends on that gene's mean expression and dispersion, so the honest summary is the **marginal (average) power** across the expression distribution at a target FDR — the expected discovery rate. A single coefficient of variation plugged into a closed-form formula mis-states power for low- and high-expressed genes alike, because dispersion varies systematically with the mean; the defensible default for count data is **simulation from the empirical mean-dispersion trend** (PROPER, Wu 2015 *Bioinformatics* 31:233; powsimR, Vieth 2017 *Bioinformatics* 33:3486). The second rule is negative: **observed (post-hoc) power is information-free.** Computed from the effect a study actually estimated, it is a one-to-one function of the p-value and cannot explain a null result (Hoenig & Heisey 2001 *Am Stat* 55:19). Power is a design-stage quantity, computed for hypothesized effects before data exist. Underpowering does not merely miss true effects — it makes the significant ones overstate magnitude (Type-M) and sometimes reverse sign (Type-S), lowering the chance a significant call is real (Button 2013 *Nat Rev Neurosci* 14:365; Gelman & Carlin 2014 *Perspect Psychol Sci* 9:641).

## Algorithmic Taxonomy

| Approach | Model | Tool | Strength | Fails / costs when |
|----------|-------|------|----------|--------------------|
| NB closed-form | negative-binomial, single CV/dispersion | `RNASeqPower::rnapower` | fast; transparent; grant-ready | one CV cannot represent the mean-dispersion trend |
| Simulation, parametric | NB with mean-dispersion relationship | `PROPER` | honest marginal power + EDR at target FDR | needs a dispersion model / pilot |
| Simulation, empirical | resampled from pilot (incl. dropout) | `powsimR` | bulk AND scRNA-seq; realistic | GitHub-only; heavier; version drift |
| Gaussian closed-form | t-test / Cohen's d | `pwr::pwr.t.test` | per-feature ATAC/proteomics after transform | wrong for raw counts; ignores overdispersion |
| Effect-inflation design analysis | retrodesign for Type-S/Type-M | `retrodesign` (Gelman) | exposes exaggeration in noisy small-n | needs a plausible true effect |

## Decision Tree by Scenario

| Scenario | Recommended approach | Why |
|----------|---------------------|-----|
| Bulk RNA-seq, pilot data available | PROPER/powsimR simulation from pilot dispersions | matches the real mean-dispersion trend |
| Bulk RNA-seq, no pilot, quick grant number | `rnapower()` with a literature CV, stated as approximate | transparent; flag as conservative-to-rough |
| scRNA-seq cross-condition DE | powsimR on a pseudobulk model; power scales with samples | population power is set by donors, not cells |
| ATAC/ChIP/methylation per-region | NB simulation (PROPER-style) or pwr after variance-stabilizing | overdispersed counts; per-region power |
| Proteomics (continuous, log-abundance) | `pwr::pwr.t.test` per protein with missingness caveat | Gaussian after transform; MNAR matters |
| Justifying a null result post-hoc | report CI / effect size, NOT observed power | post-hoc power is uninformative (Hoenig-Heisey) |
| Fixed budget: depth vs replicates | favor replicates past ~10-20M mapped reads | biological variance dominates (Liu 2014) |
| Clinical-trial endpoint | -> clinical-biostatistics/power-and-sample-size | regulated regime, different machinery |

## Closed-Form NB Power -- RNASeqPower

**Goal:** Get a fast, transparent power or replicate number for bulk RNA-seq from depth, biological CV, and fold change.

**Approach:** Supply per-gene depth, biological coefficient of variation, the fold change to detect, and alpha; supply `n` to get power, or `power` to get the required `n`. Treat the result as a single-gene approximation and sanity-check against simulation.

```r
library(RNASeqPower)
# depth = reads/gene; cv = biological coefficient of variation; effect = fold change
rnapower(depth = 20, n = 5, cv = 0.4, effect = 2, alpha = 0.05)          # solves for POWER
rnapower(depth = 20, cv = 0.4, effect = 2, alpha = 0.05, power = 0.80)   # solves for n per group
```

## Simulation-Based Power -- the Honest Default for Counts

**Goal:** Estimate marginal power and the true realized FDR across the whole expression distribution, accounting for the mean-dispersion trend.

**Approach:** Build (or fit from pilot) a simulation model of counts with a realistic dispersion-mean relationship and DE-effect distribution, simulate many datasets at each candidate sample size, run the intended DE test, and read the average power at the target FDR.

```r
library(PROPER)
sim_opts <- RNAseq.SimOptions.2grp(ngenes = 20000, p.DE = 0.05,
                                   lOD = 'cheung', lBaselineExpr = 'cheung')  # empirical dispersion/expr priors
sims <- runSims(Nreps = c(3, 5, 8, 12), sim.opts = sim_opts, nsims = 50,
                DEmethod = 'edgeR')
powr <- comparePower(sims, alpha.type = 'fdr', alpha.nominal = 0.05,
                     stratify.by = 'expr', delta = log(1.5))          # delta is NATURAL-log lfc in PROPER; marginal power by expression stratum
summaryPower(powr)
```

## Depth vs Replicates -- the Budget Question

For bulk RNA-seq differential expression, sequencing depth shows diminishing returns once it is adequate — Liu, Zhou & White 2014 (*Bioinformatics* 30:301) found the inflection near **~10 million mapped reads** in MCF7 (commonly generalized to a 10-20M band) — whereas adding biological replicates improves power across the whole range. Under a fixed budget, allocate to more biological units before more depth. ATAC/ChIP have their own depth floors (library complexity, peak detection), but the principle holds: biological variance, not read count, limits discovery once depth is adequate.

## CV / Dispersion Guidelines (estimate from pilot when possible)

| Material | Typical biological CV | Source / note |
|----------|----------------------|---------------|
| Cell lines (technical replicates) | 0.1-0.2 | low biological variability |
| Inbred mice | 0.2-0.3 | moderate |
| Primary cells / donor-derived | 0.3-0.4 | donor-dependent |
| Human population samples | 0.3-0.5 | high; Hart 2013 *J Comput Biol* 20:970 default examples |

These are starting points, not substitutes for a pilot estimate; real dispersion is study-specific and a literature CV can be off by a factor of two (estimate via DESeq2/edgeR `estimateDispersions` — see experimental-design/sample-size).

## Per-Method Failure Modes

### Single CV for the whole transcriptome
- **Trigger:** one `cv` p