# Experimental Design and Statistical Planning

Comprehensive workflow for statistical experimental design in genomics, from power analysis and sample size determination to batch-balanced experimental layouts and multiple testing strategy.

## When to Use This Skill

Use this skill when you need to:
- ✅ **Plan new experiments** - Design from scratch with statistical rigor
- ✅ **Justify sample sizes** - Calculate required replicates for grant proposals
- ✅ **Perform power analysis** - Determine statistical power for proposed designs
- ✅ **Design batch layouts** - Create balanced assignments preventing confounding
- ✅ **Optimize budgets** - Balance sequencing depth vs. number of replicates
- ✅ **Select correction methods** - Choose appropriate multiple testing approaches

**Don't use this skill for:**
- ❌ Post-experiment analysis → Use appropriate DE analysis skills
- ❌ Simple two-sample comparisons with fixed n → Use power calculators directly

**Key Concept:** Proper experimental design is the foundation of reproducible science. Good design prevents confounding, maximizes statistical power, and ensures results are interpretable.

## Installation

### Required Software

| Software | Version | License | Commercial Use | Installation |
| --- | --- | --- | --- | --- |
| DESeq2 | ≥1.30.0 | LGPL (≥3) | ✅ Permitted | `BiocManager::install('DESeq2')` |
| RNASeqPower | ≥1.30.0 | LGPL | ✅ Permitted | `BiocManager::install('RNASeqPower')` |
| ssizeRNA | ≥1.3.2 | GPL-2 | ✅ Permitted | `BiocManager::install('ssizeRNA')` |
| powsimR | ≥1.2.3 | GPL-3 | ✅ Permitted | `BiocManager::install('powsimR')` |
| IHW | ≥1.18.0 | Artistic-2.0 | ✅ Permitted | `BiocManager::install('IHW')` |
| osat | ≥1.38.0 | Artistic-2.0 | ✅ Permitted | `install.packages('osat')` |
| ggplot2 | ≥3.3.0 | MIT | ✅ Permitted | `install.packages('ggplot2')` |
| ggprism | ≥1.0.3 | GPL-3 | ✅ Permitted | `install.packages('ggprism')` |
| jsonlite | ≥1.7.0 | MIT | ✅ Permitted | `install.packages('jsonlite')` |
| pasilla | ≥1.18.0 | Artistic-2.0 | ✅ Permitted | `BiocManager::install('pasilla')` |

**Optional (for example data):**
- pasilla - Example RNA-seq pilot data for testing

**Quick install:**

```
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install(c("DESeq2", "RNASeqPower", "ssizeRNA", "powsimR", "IHW", "pasilla"))
install.packages(c("osat", "ggplot2", "ggprism", "jsonlite"))
```

**License Compliance:** All packages use LGPL, GPL, Artistic-2.0, or MIT licenses that permit commercial use in AI agent applications.

**Full installation instructions and version details:** [references/software\_requirements.md](#)

## Inputs

**Required:**
- **Experimental design info**: Assay type, n conditions, sample relationship, planned n
- **Effect size expectations**: Target fold change, variability (CV or pilot data)
- **Statistical requirements**: Target power (0.80/0.90), α (0.05), multiple testing preference

**Optional:**
- **Practical constraints**: Budget, sample availability, batch structure, sequencing depth, covariates

**Detailed input requirements:** [references/experimental\_design\_best\_practices.md#input-requirements](#)

## Outputs

**Power and sample size:**
- `power_analysis_results.csv` - Power calculations for scenarios
- `sample_size_recommendation.txt` - Required n with justification
- `power_vs_n_curve.png` + `.svg` - Power relationship visualizations

**Batch design:**
- `batch_layout_for_lab.csv` - Lab-ready sample-to-batch assignments
- `batch_design_validation.txt` - Confounding check results
- `batch_design_plot.png` + `.svg` - Visual layout

**Documentation:**
- `statistical_analysis_plan.md` - Complete pre-registration plan
- `lab_protocol_checklist.md` - Step-by-step processing guide
- `design_parameters.json` - All parameters (human-readable)

**Analysis objects (RDS) - For downstream use:**
- `batch_design.rds` - Complete batch design object
- Load with: `batch_design <- readRDS('batch_design.rds')`
- Required for: Batch effect correction workflows
- `design_parameters.rds` - Complete design parameters
- Load with: `design_params <- readRDS('design_parameters.rds')`
- Required for: Analysis validation, replication studies

## Clarification Questions

Before starting, gather the following information:

### 1. **Input Files** (ASK THIS FIRST):

- Do you have pilot data or existing results files to inform the experimental design?
- If uploaded: Are these pilot data files (DESeq2 objects, count matrices) you'd like to use for power calculations?
- Expected formats: RDS (DESeqDataSet), CSV/TSV (count matrices)
- **Or use literature-based estimates?** We can use tissue-specific variability values from published data ([references/cv\_tissue\_database.csv](#)).

### 2. **Assay Type**

- **Bulk RNA-seq**, **scRNA-seq**, **ATAC-seq**, **ChIP-seq**, **Methylation**, **Proteomics**, or **Other**

### 3. **Experimental Structure**

- **Conditions**: 2 (case-control), 3+ (multi-group), or factorial design?
- **Planned n**: Replicates per condition (or "unknown")
- **Sample type**: Independent, paired, or repeated measures
- **Covariates**: Variables to balance (sex, age, batch, site)

### 4. **Effect Size & Variability**

- **Target fold change**: Large (≥2x), moderate (1.5-2x), small (1.2-1.5x)
- **Pilot data**: Available? (Provide DESeq2 object, count matrix, or path)
- **No pilot data**: Will use tissue-specific CV from [references/cv\_tissue\_database.csv](#)

### 5. **Statistical Requirements**

- **Power**: 0.80 (standard), 0.90 (grants), or custom
- **Alpha (α)**: 0.05 (standard), 0.01 (stringent), or custom
- **Multiple testing**: BH-FDR (standard), IHW (more power), Bonferroni (stringent), or need guidance

### 6. **Practical Constraints**

- **Budget**: Maximum samples or runs
- **Sample availability**: Limited? Account for failures (add 10-20% buffer)?
- **Batch structure**: Processed in batches? Batch size?
- **Sequencing depth**: Target reads (RNA-seq: 15-30M, ATAC-seq: 25-50M, scRNA-seq: 50-100K/cell)

### 7. **Primary Objective**

- **Power analysis**, **sample size determination**, **batch design**, **multiple testing guidance**, **complete design**, or **budget optimization**

**Comprehensive clarification guide:** [references/experimental\_design\_best\_practices.md#clarification-questions](#)

## Standard Workflow

🚨 **MANDATORY: USE SCRIPTS EXACTLY AS SHOWN - DO NOT WRITE INLINE CODE** 🚨

The experimental design workflow follows 4 steps: **Load** → **Calculate** → **Visualize** → **Export**

### **Step 1 - Load Parameters**

```
source("scripts/load_example_data.R")
pilot_data <- load_example_data()
```

**With pilot data (preferred):**
- Uses `pilot_data$dds` for power calculations
- Uses `pilot_data$cv$median` for sample size estimation
- Provides realistic variability estimates

**Without pilot data (alternative):**

```
source("scripts/load_example_data.R")
cv_db <- load_cv_database()
# Select appropriate tissue type from cv_db
```

**✅ VERIFICATION:** You MUST see: `"✓ Example pilot data loaded successfully!"`

**Decision:** Pilot data provides more accurate estimates. See [power\_analysis\_guidelines.md#pilot-vs-literature](#)

---

### **Step 2 - Calculate Design**

🚨 **DO NOT write inline calculation code. Use the provided scripts.**

**A. Power Analysis** - Calculate power for your proposed design

```
source("scripts/power_rnaseq.R")
power_result <- calc_power_rnaseq(
  depth = 20,
  n_per_group = 6,
  cv = pilot_data$cv$median,
  fold_change = 2,
  alpha = 0.05
)
```

**DO NOT write inline power calculation code. Just source the script and call the function.**

**B. Sample Size Determination** - Calculate required n from pilot data

```
source("scripts/sample_size_de.R")
required_n <- samplesize_from_pilot(
  pilot_dds = pilot_data$dds,
  fold_change = 1.5,
  power = 0.8,
  fdr = 0.05
)
```

**DO NOT write inline sample size code. Use the function from the script.**

**C. Batch Assignment** - Generate balanced batch layout

```
source("scripts/batch_assignment.R")
batch_design <- assign_samples_to_batches(
  metadata = pilot_data$metadata,
  batch_size = 8,
  balance_vars = c("condition", "sex")
)
```

**DO NOT manually create batch assignments. Use the OSAT-optimized function.**

⚠️ **CRITICAL - DO NOT:**
- ❌ Write inline power calculation code → **STOP: Use calc\_power\_rnaseq()**
- ❌ Write inline plotting code (ggsave, ggplot, etc.) → **STOP: Use visualization scripts**
- ❌ Manually assign samples to batches → **STOP: Use assign\_samples\_to\_batches()**
- ❌ Write custom balancing algorithms → **STOP: Script uses OSAT's optimal algorithms**
- ❌ Try to install svglite → scripts handle SVG fallback automatically

**⚠️ IF SCRIPTS FAIL - Script Failure Hierarchy:**
1. **Fix and Retry (90%)** - Install missing package, re-run script
2. **Modify Script (5%)** - Edit the script file itself, document changes
3. **Use as Reference (4%)** - Read script, adapt approach, cite source
4. **Write from Scratch (1%)** - Only if genuinely impossible, explain why

**NEVER skip directly to writing inline code without trying the script first.**

**✅ VERIFICATION:** You should see:
- After power analysis: `"✓ Power analysis completed successfully!"`
- After sample size: `"✓ Sample size calculation completed successfully!"`
- After batch design: `"✓ Batch design generated successfully!"`

**CRITICAL RULE:** Batch must NEVER confound with condition. See [batch\_effect\_mitigation.md#cardinal-rule](#)

**Decision points:**
- Power ≥0.80? If not, increase n or adjust expectations. See [power\_analysis\_guidelines.md#interpreting-power](#)
- Required n exceed budget? See [budget optimization](#)
- Confounding detected? Regenerate with different constraints. See [batch\_effect\_mitigation.md#troubleshooting](#)

---

### **Step 3 - Visualize Design**

**A. Generate power curves:**

```
source("scripts/plot_power_curves.R")
plot_power_vs_samplesize(
  cv = pilot_data$cv$median,
  fold_changes = c(1.5, 2, 3),
  depth = 20,
  output_file = "design_results/power_vs_n"
)
```

**B. Validate and visualize batch design:**

```
source("scripts/batch_validation.R")
confounding_check <- check_confounding(batch_design, "condition")
visualize_batch_design(
  batch_design,
  condition_var = "condition",
  output_file = "design_results/batch_design"
)
```

🚨 **DO NOT write inline plotting code (ggsave, ggplot, etc.). Use the visualization scripts.** 🚨

**The scripts handle PNG + SVG export with graceful fallback for SVG dependencies.**

**✅ VERIFICATION:** You should see:
- `"Saving power curve plots:"` followed by PNG + SVG file paths
- `"PASS: No confounding detected"` or `"WARNING: Batch is CONFOUNDED"`
- `"Saving batch design plots:"` followed by PNG + SVG file paths

**Output formats:** Always generates both PNG (for presentations) and SVG (for publications) with graceful fallback.

---

### **Step 4 - Export All Results**

```
source("scripts/export_design.R")
export_complete_design(batch_design, design_params, output_dir = "design_results")
```

**DO NOT write custom export code. Use export\_complete\_design().**

**✅ VERIFICATION:** You MUST see: `"=== Export Complete ==="`

This will generate:
1. `batch_layout_for_lab.csv` - Lab-ready sample assignments
2. `statistical_analysis_plan.md` - Pre-registration analysis plan
3. `lab_protocol_checklist.md` - Lab processing checklist
4. `batch_design.rds` - Batch design object (for downstream use)
5. `design_parameters.rds` - Design parameters (for downstream use)
6. `design_parameters.json` - Design parameters (human-readable)

**RDS objects are CRITICAL** for downstream workflows and validation studies.

---

### **Complete Workflow Example**

For a complete experimental design with all steps:

```
# Step 1: Load pilot data
source("scripts/load_example_data.R")
pilot_data <- load_example_data()

# Step 2: Calculate design parameters
source("scripts/power_rnaseq.R")
source("scripts/sample_size_de.R")
source("scripts/batch_assignment.R")

power_result <- calc_power_rnaseq(depth = 20, n_per_group = 6, cv = 0.4, fold_change = 2)
required_n <- samplesize_from_pilot(pilot_data$dds, fold_change = 1.5, power = 0.8)
batch_design <- assign_samples_to_batches(pilot_data$metadata, batch_size = 8,
                                          balance_vars = c("condition"))

# Step 3: Visualize and validate
source("scripts/plot_power_curves.R")
source("scripts/batch_validation.R")

plot_power_vs_samplesize(cv = 0.4, fold_changes = c(1.5, 2, 3),
                         output_file = "design_results/power_vs_n")
check_confounding(batch_design, "condition")
visualize_batch_design(batch_design, "condition", output_file = "design_results/batch_design")

# Step 4: Export all results
source("scripts/export_design.R")
design_params <- list(assay = "RNA-seq", conditions = c("control", "treated"),
                     n_per_group = 6, power = 0.85, alpha = 0.05,
                     effect_size = 2, multiple_testing = "BH-FDR")
export_complete_design(batch_design, design_params, output_dir = "design_results")
```

**Note:** Specific parameters depend on your experimental requirements (see Clarification Questions).

## Decision Guide

**Three critical decisions:**
- **Pilot vs Literature:** Use pilot data if available (more accurate). Literature CV acceptable as fallback.
- **Sample Size vs Depth:** Prioritize more samples over deeper sequencing for DE. 15-20M reads sufficient for RNA-seq.
- **Multiple Testing:** BH-FDR (standard), IHW (more power), Bonferroni (stringent).

**See:** [experimental\_design\_best\_practices.md#decision-guide](#) for comprehensive guidance and common usage patterns (power analysis only, batch design only, budget optimization).

## Common Issues

| Issue | Cause | Solution |
| --- | --- | --- |
| Power <0.80 with max budget | Effect size too small or CV too high | Increase n, increase depth, or revise effect size expectations. See [references/power\_analysis\_guidelines.md#low-power](#) |
| Batch confounding detected | Unequal condition distribution across batches | Regenerate with stricter balance constraints or adjust batch size. See [references/batch\_effect\_mitigation.md#troubleshooting](#) |
| Required n exceeds sample availability | Pilot data shows high variability or small effect | Consider paired design, blocking by major covariates, or revise target fold-change. See [references/experimental\_design\_best\_practices.md#budget-optimization](#) |
| Can't balance all covariates | Too many variables for batch size | Prioritize key covariates (condition > sex > age > others). Some minor imbalance acceptable. See [references/batch\_effect\_mitigation.md#covariate-priority](#) |
| CV estimate varies widely | Pilot data has outliers or low counts | Filter low-count genes (mean <10) before CV calculation. Use median, not mean CV. See [references/power\_analysis\_guidelines.md#cv-estimation](#) |
| Power calculations give n<3 | Very large effect size or low variability | Warning: n<3 too low for valid inference. Plan for minimum n=3-4 even if calculations suggest n=2 |
| Multiple testing correction too stringent | Many tests, low discovery rate | Consider IHW (more powerful than BH-FDR) or independent filtering. See [references/multiple\_testing\_guide.md#choosing](#) |

**Detailed troubleshooting:** [references/troubleshooting\_guide.md](#)

## Suggested Next Steps

After completing experimental design:

1. **Execute Experiment** - Use batch assignment file to guide sample processing
2. **Perform DE Analysis** - Use bulk-rnaseq-counts-to-de-deseq2, scrnaseq-scanpy-core-analysis, or appropriate skill
3. **Apply Multiple Testing** - Use pre-specified correction method from statistical plan
4. **Validate Results** - Check batch effects were controlled, verify power calculations

## Related Skills

**Upstream:** None - this is typically the first step in a project

**Downstream (after data collection):**
- **bulk-rnaseq-counts-to-de-deseq2** - Differential expression analysis
- **functional-enrichment-from-degs** - Pathway analysis
- **de-results-to-plots** - Visualization

**Alternative/complementary:**
- **bulk-omics-clustering** - Discover natural groupings post-hoc
- **batch-correction-combat** - Computational batch correction if needed

## References

**Detailed documentation:**
- [references/experimental\_design\_best\_practices.md](#) - General design principles, decision guide, common patterns
- [references/power\_analysis\_guidelines.md](#) - Detailed power calculation methods, pilot vs literature
- [references/batch\_effect\_mitigation.md](#) - Preventing/controlling batch effects, cardinal rule, troubleshooting
- [references/multiple\_testing\_guide.md](#) - Choosing correction methods
- [references/qc\_guidelines.md](#) - Quality control checkpoints
- [references/troubleshooting\_guide.md](#) - Common problems and solutions
- [references/software\_requirements.md](#) - Installation and licenses
- [references/cv\_tissue\_database.csv](#) - Tissue-specific variability estimates

**Scripts:** See scripts/ directory for all analysis functions:
- Data loading: [load\_example\_data.R](#)
- Power/sample size: [power\_rnaseq.R](#), [power\_atacseq.R](#), [sample\_size\_de.R](#), [sample\_size\_scrna.R](#)
- Batch design: [batch\_assignment.R](#), [batch\_validation.R](#)
- Visualization: [plot\_power\_curves.R](#)
- Export: [export\_design.R](#) (includes RDS saving)

**Key Papers:**
- Hart SN et al. (2013) *J Comput Biol* 20(12):970-978 - RNA-seq sample size
- Schurch NJ et al. (2016) *RNA* 22(6):839-851 - Biological replicates needed
- Leek JT et al. (2010) *Nat Rev Genet* 11(10):733-739 - Batch effects impact
- Benjamini & Hochberg (1995) *J R Stat Soc Series B* 57(1):289-300 - FDR control
- Love MI et al. (2014) *Genome Biol* 15(12):550 - DESeq2 methods