Experimental Design & Stats

Power, sample size, batch balance — before you run.

Overview

Problem. How many samples; how to avoid batch confounding?

Learning goals

• Do the stats before, not after
• Batch–group confounding is the top killer

Figures

Tutorial

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

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

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

Detailed input requirements: references/experimental_design_best_practices.md#input-requirements

Outputs

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()

Without pilot data (alternative):

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

Decision: Pilot data provides more accurate estimates. See power_analysis_guidelines.md#pilot-vs-literature

Step 2 - Calculate Design

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
)

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
)

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")
)

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

CRITICAL RULE: Batch must NEVER confound with condition. See batch_effect_mitigation.md#cardinal-rule

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"
)

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

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")

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

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

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

References

Code preview

scripts/batch_assignment.R

# Batch Assignment Functions
# Create balanced batch designs preventing confounding
# Based on: Leek et al. (2010) Nat Rev Genet and osat package

#' Assign samples to batches with optimal balance
#'
#' @param metadata Data frame with sample metadata (must include sample IDs)
#' @param batch_size Number of samples per batch
#' @param balance_vars Vector of column names to balance across batches
#' @param sample_id_col Name of sample ID column (default: "sample_id")
#' @return Data frame with batch assignments added
#' @export
#' @examples
#' # Create sample metadata
#' metadata <- data.frame(
#'   sample_id = paste0("S", 1:24),
#'   condition = rep(c("Control", "Treatment"), each = 12),
#'   sex = rep(c("M", "F"), 12),
#'   age_group = rep(c("Young", "Old"), each = 12)
#' )
#'
#' # Generate balanced batch assignment
#' batch_design <- assign_samples_to_batches(
#'   metadata = metadata,
#'   batch_size = 8,
#'   balance_vars = c("condition", "sex", "age_group")
#' )
assign_samples_to_batches <- function(metadata,
                                      batch_size,
                                      balance_vars,
                                      sample_id_col = "sample_id") {

  if (!requireNamespace("OSAT", quietly = TRUE)) {
    stop("Package 'OSAT' is required. Install with BiocManager::install('OSAT')")
  }

  # Validate inputs
  if (!sample_id_col %in% colnames(metadata)) {
    stop(paste0("Column '", sample_id_col, "' not found in metadata"))
  }

  missing_vars <- balance_vars[!balance_vars %in% colnames(metadata)]
  if (length(missing_vars) > 0) {
    stop(paste0("Balance variables not found in metadata: ", paste(missing_vars, collapse = ", ")))
  }

  if (batch_size <= 0 || batch_size > nrow(metadata)) {
    stop("Invalid batch_size: must be positive and <= number of samples")
  }

  # Calculate number of batches
  n_batches <- ceiling(nrow(metadata) / batch_size)

  # Prepare data for osat
  sample_ids <- metadata[[sample_id_col]]

  # Create factor variables for balancing
  balance_data <- metadata[, balance_vars, drop = FALSE]

  # Convert all balance variables to factors
  balance_data <- as.data.frame(lapply(balance_data, as.factor))

  # Use osat to generate optimal assignment
  message("Generating optimal batch assignment...")

  tryCatch({
    # Run osat optimization
    osat_result <- OSAT::osat(
      x = balance_data,
      batch.size = batch_size,
      n.batch = n_batches
    )

    # Extract batch assignments
    batch_assignments <- osat_result$BatchLabel

    # Add batch column to metadata
    metadata$batch <- batch_assignments

    # Reorder by batch for convenience

scripts/batch_validation.R

# Batch Design Validation Functions
# Validate batch assignments to prevent confounding and check balance

#' Internal helper to save plots in both PNG and SVG formats
#' @param plot ggplot object
#' @param base_path Base file path (extension will be replaced)
#' @param width Plot width in inches
#' @param height Plot height in inches
#' @param dpi Resolution for PNG
.save_plot <- function(plot, base_path, width = 8, height = 6, dpi = 300) {
  # Always save PNG
  png_path <- sub("\\.(svg|png)$", ".png", base_path)
  ggplot2::ggsave(png_path, plot = plot, width = width, height = height, dpi = dpi, device = "png")
  cat("   Saved:", png_path, "\n")

  # Always try SVG - try ggsave first, fall back to svg() device
  svg_path <- sub("\\.(svg|png)$", ".svg", base_path)
  tryCatch({
    ggplot2::ggsave(svg_path, plot = plot, width = width, height = height, device = "svg")
    cat("   Saved:", svg_path, "\n")
  }, error = function(e) {
    # If ggsave fails, try base R svg() device directly
    tryCatch({
      svg(svg_path, width = width, height = height)
      print(plot)
      dev.off()
      cat("   Saved:", svg_path, "\n")
    }, error = function(e2) {
      cat("   (SVG export failed)\n")
    })
  })
}

#' Check for confounding between batch and condition
#'
#' @param batch_assignment Data frame with batch assignments
#' @param condition_var Name of condition variable column
#' @param batch_var Name of batch variable column (default: "batch")
#' @return List with test results and interpretation
#' @export
#' @examples
#' # Check if batch is confounded with condition
#' # confounding_check <- check_confounding(batch_design, "condition")
check_confounding <- function(batch_assignment,
                              condition_var,
                              batch_var = "batch") {

  # Validate inputs
  if (!condition_var %in% colnames(batch_assignment)) {
    stop(paste0("Condition variable '", condition_var, "' not found in data"))
  }

  if (!batch_var %in% colnames(batch_assignment)) {
    stop(paste0("Batch variable '", batch_var, "' not found in data"))
  }

  # Create contingency table
  cont_table <- table(batch_assignment[[batch_var]], batch_assignment[[condition_var]])

  # Perform chi-square test
  chi_test <- chisq.test(cont_table)

  # Interpretation
  is_confounded <- chi_test$p.value < 0.05

  if (is_confounded) {
    result_message <- paste0(
      "WARNING: Batch is CONFOUNDED with ", condition_var, "!\n",
      "Chi-square test p-value: ", format(chi_test$p.value, digits = 4), " < 0.05\n",
      "This design will not allow separation of batch effects from biological effects.\n",
      "RECOMMENDATION: Regenerate batch assignment."
    )
    status <- "FAILED"
  } else {
    result_message <- paste0(
      "PASS: No confounding detected between batch and ", condition_var, "\n",
      "Chi-square test p-value: ", format(chi_test$p.value, digits = 4), " > 0.05\n",
      "Batch effects can be separated from biological effects."
    )
    status <- "PASSED"

scripts/export_design.R

# Export Experimental Design Functions
# Generate lab-ready files and documentation

#' Export batch assignment layout for lab use
#'
#' @param batch_assignment Data frame with batch assignments
#' @param output_file Output CSV file path
#' @export
export_batch_layout <- function(batch_assignment, output_file = "batch_layout_for_lab.csv") {

  # Reorder columns for lab convenience
  priority_cols <- c("sample_id", "batch", "condition", "processing_order",
                     "plate", "well", "overall_sequence")

  # Get columns that exist
  existing_priority <- priority_cols[priority_cols %in% colnames(batch_assignment)]
  other_cols <- setdiff(colnames(batch_assignment), existing_priority)

  # Reorder
  export_data <- batch_assignment[, c(existing_priority, other_cols)]

  # Sort by batch and processing order (if available)
  if ("batch" %in% colnames(export_data)) {
    if ("processing_order" %in% colnames(export_data)) {
      export_data <- export_data[order(export_data$batch, export_data$processing_order), ]
    } else {
      export_data <- export_data[order(export_data$batch), ]
    }
  }

  # Write CSV
  write.csv(export_data, output_file, row.names = FALSE)

  message(paste0("Batch layout exported to: ", output_file))
  message(paste0("Samples: ", nrow(export_data)))
  message(paste0("Batches: ", length(unique(export_data$batch))))

  # Return invisibly
  invisible(export_data)
}


#' Export statistical analysis plan
#'
#' @param design_params List with experimental design parameters
#' @param output_file Output markdown file path
#' @export
export_statistical_plan <- function(design_params, output_file = "statistical_analysis_plan.md") {

  # Create markdown document
  plan_text <- paste0(
    "# Statistical Analysis Plan\n\n",
    "**Generated:** ", format(Sys.time(), "%Y-%m-%d %H:%M:%S"), "\n\n",
    "---\n\n",
    "## Experimental Design\n\n",
    "**Assay Type:** ", design_params$assay, "\n\n",
    "**Experimental Groups:**\n",
    "- ", paste(design_params$conditions, collapse = "\n- "), "\n\n",
    "**Sample Size:** ", design_params$n_per_group, " biological replicates per group (",
    design_params$n_per_group * length(design_params$conditions), " total samples)\n\n"
  )

  # Add batch information if present
  if (!is.null(design_params$batches)) {
    plan_text <- paste0(
      plan_text,
      "**Batch Structure:** ", design_params$batches, " processing batches\n",
      "- Each batch contains all experimental conditions to prevent confounding\n\n"
    )
  }

  # Add power analysis results
  if (!is.null(design_params$power)) {
    plan_text <- paste0(
      plan_text,
      "## Power Analysis\n\n",
      "**Statistical Power:** ", round(design_params$power, 3), "\n",
      "**Minimum Detectable Effect:** ", design_params$effect_size, "-fold change\n",
      "**Significance Threshold:** α = ", design_params$alpha, "\n\n"
    )

Companion files