Bio Experimental Design Multiple Testing

Name: Bio Experimental Design Multiple Testing
Author: gptomics

gptomics/bioskills

2 installs
1.1k repo stars
Updated July 25, 2026
gptomics/bioskills

Apply multiple-testing correction (FDR, Bonferroni, q-value) to genomics p-values using R p.adjust/qvalue or Python statsmodels multipletests.

About

Applies multiple-testing correction methods to adjust raw p-values from thousands of simultaneous genomics tests, controlling false discovery or family-wise error rate. A developer uses it when filtering differential expression results or setting significance thresholds.

BH/Bonferroni/q-value methods across R and Python
Guidance on choosing correction per study design

Bio Experimental Design Multiple Testing by the numbers

2 all-time installs (skills.sh)
Ranked #1,759 of 2,064 Data Science & ML skills by installs in the Skillselion catalog
Data as of Aug 5, 2026 (Skillselion catalog sync)

npx skills add https://github.com/gptomics/bioskills --skill bio-experimental-design-multiple-testing

Add your badge

Show developers this skill is listed on Skillselion. Paste this into your README.

[![Listed on Skillselion](https://skillselion.com/badge/skills/gptomics/bioskills/bio-experimental-design-multiple-testing.svg)](https://skillselion.com/skills/gptomics/bioskills/bio-experimental-design-multiple-testing)

Installs	2
repo stars	★ 1.1k
Last updated	July 25, 2026
Repository	gptomics/bioskills ↗

What it does

Apply multiple-testing correction (FDR, Bonferroni, q-value) to genomics p-values using R p.adjust/qvalue or Python statsmodels multipletests.

Files

SKILL.mdMarkdownGitHub ↗

Version Compatibility

Reference examples tested with: R stats (base), statsmodels 0.14+

Before using code patterns, verify installed versions match. If versions differ:

Python: pip show <package> then help(module.function) to check signatures
R: packageVersion('<pkg>') then ?function_name to verify parameters

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Multiple Testing Correction

"Correct p-values for multiple testing" → Adjust raw p-values from thousands of simultaneous tests to control false discovery rate or family-wise error rate.

R: p.adjust(pvalues, method = 'BH'), qvalue::qvalue()
Python: statsmodels.stats.multitest.multipletests()

The Problem

Testing 20,000 genes at p < 0.05 yields ~1,000 false positives by chance. Correction is essential.

Common Methods

Bonferroni (Most Conservative)

# Strict family-wise error rate control
p_adj <- p.adjust(pvalues, method = 'bonferroni')
# Threshold: alpha / n_tests
# Use for: small gene sets, confirmatory studies

Benjamini-Hochberg FDR (Standard)

# Controls false discovery rate
p_adj <- p.adjust(pvalues, method = 'BH')
# Most common for genomics
# FDR 0.05 = expect 5% of significant results to be false

q-value (Recommended for Large-Scale)

Goal: Estimate the false discovery rate for each gene in a genome-wide test while maximizing detection power by estimating the proportion of true nulls.

Approach: Fit the q-value model to the p-value distribution, which estimates pi0 (fraction of true null hypotheses) and converts each p-value to a q-value representing the minimum FDR at which that gene would be called significant.

library(qvalue)
qobj <- qvalue(pvalues)
qvalues <- qobj$qvalues
pi0 <- qobj$pi0  # Estimated proportion of true nulls

# q-value directly estimates FDR for each gene
# More powerful than BH when many true positives exist

Method Selection Guide

Scenario	Recommended Method	Threshold
Genome-wide DE	BH or q-value	FDR < 0.05
Candidate genes	Bonferroni	p < 0.05/n
Exploratory	BH	FDR < 0.10
Validation study	Bonferroni	p < 0.05/n
GWAS	Bonferroni	p < 5e-8

Python Equivalent

from statsmodels.stats.multitest import multipletests

# Benjamini-Hochberg
rejected, pvals_corrected, _, _ = multipletests(pvalues, method='fdr_bh')

# Bonferroni
rejected, pvals_corrected, _, _ = multipletests(pvalues, method='bonferroni')

Interpreting Results

FDR 0.05: Among genes called significant, ~5% are false positives
FDR 0.01: More stringent, fewer false positives but more false negatives
padj vs qvalue: Both estimate FDR; q-value is slightly more powerful

Related Skills

differential-expression/de-results - Applying corrections to DE output
population-genetics/association-testing - GWAS significance thresholds
pathway-analysis/go-enrichment - Correcting enrichment p-values

# Reference: R stats (base), statsmodels 0.14+ | Verify API if version differs
# Multiple Testing Correction Examples
# Demonstrates different correction methods and when to use them

# =============================================================================
# Generate Example P-values
# =============================================================================

set.seed(42)
n_genes <- 10000
n_true_de <- 500  # 5% truly DE

# Simulate p-values
# True nulls: uniform(0,1)
# True positives: beta distribution skewed toward 0
pvalues <- c(
  rbeta(n_true_de, 0.3, 5),      # True DE genes (small p-values)
  runif(n_genes - n_true_de)     # True nulls (uniform)
)
is_true_de <- c(rep(TRUE, n_true_de), rep(FALSE, n_genes - n_true_de))

cat('Simulated', n_genes, 'genes with', n_true_de, 'truly DE\n')
cat('Significant at p < 0.05:', sum(pvalues < 0.05), '\n\n')

# =============================================================================
# Bonferroni Correction
# =============================================================================

# Most conservative - controls family-wise error rate (FWER)
# FWER = P(at least one false positive)
# Use for: small targeted studies, confirmatory analyses

p_bonf <- p.adjust(pvalues, method = 'bonferroni')
sig_bonf <- sum(p_bonf < 0.05)

cat('=== Bonferroni Correction ===\n')
cat('Threshold: p <', 0.05/n_genes, '\n')
cat('Significant genes:', sig_bonf, '\n')
cat('True positives:', sum(p_bonf < 0.05 & is_true_de), '\n')
cat('False positives:', sum(p_bonf < 0.05 & !is_true_de), '\n\n')

# =============================================================================
# Benjamini-Hochberg FDR
# =============================================================================

# Standard for genomics - controls false discovery rate
# FDR = E[false positives / total positives]
# More powerful than Bonferroni when many true positives exist

p_bh <- p.adjust(pvalues, method = 'BH')
sig_bh <- sum(p_bh < 0.05)

cat('=== Benjamini-Hochberg FDR ===\n')
cat('Significant at FDR < 0.05:', sig_bh, '\n')
cat('True positives:', sum(p_bh < 0.05 & is_true_de), '\n')
cat('False positives:', sum(p_bh < 0.05 & !is_true_de), '\n')
cat('Observed FDR:', round(sum(p_bh < 0.05 & !is_true_de) / sig_bh, 3), '\n\n')

# =============================================================================
# q-value Method
# =============================================================================

# More powerful than BH by estimating pi0 (proportion true nulls)
# Directly estimates FDR for each feature

library(qvalue)
qobj <- qvalue(pvalues)
qvalues <- qobj$qvalues
sig_qval <- sum(qvalues < 0.05)

cat('=== q-value Method ===\n')
cat('Estimated pi0 (true null proportion):', round(qobj$pi0, 3), '\n')
cat('Significant at q < 0.05:', sig_qval, '\n')
cat('True positives:', sum(qvalues < 0.05 & is_true_de), '\n')
cat('False positives:', sum(qvalues < 0.05 & !is_true_de), '\n\n')

# =============================================================================
# Method Comparison
# =============================================================================

cat('=== Method Comparison ===\n')
comparison <- data.frame(
  Method = c('None (p < 0.05)', 'Bonferroni', 'BH FDR', 'q-value'),
  Significant = c(sum(pvalues < 0.05), sig_bonf, sig_bh, sig_qval),
  TruePos = c(sum(pvalues < 0.05 & is_true_de),
              sum(p_bonf < 0.05 & is_true_de),
              sum(p_bh < 0.05 & is_true_de),
              sum(qvalues < 0.05 & is_true_de)),
  FalsePos = c(sum(pvalues < 0.05 & !is_true_de),
               sum(p_bonf < 0.05 & !is_true_de),
               sum(p_bh < 0.05 & !is_true_de),
               sum(qvalues < 0.05 & !is_true_de))
)
comparison$Sensitivity <- round(comparison$TruePos / n_true_de, 3)
comparison$FDR <- round(comparison$FalsePos / comparison$Significant, 3)

print(comparison, row.names = FALSE)

# =============================================================================
# Method Selection Guide
# =============================================================================

cat('\n=== When to Use Each Method ===\n')
cat('Bonferroni:\n')
cat('  - Small, targeted gene panels (<100 genes)\n')
cat('  - Confirmatory/validation studies\n')
cat('  - When ANY false positive is unacceptable\n\n')

cat('Benjamini-Hochberg:\n')
cat('  - Standard genome-wide DE analysis\n')
cat('  - Exploratory studies\n')
cat('  - When some false positives are acceptable\n\n')

cat('q-value:\n')
cat('  - Large-scale studies with many true positives\n')
cat('  - Maximum power needed\n')
cat('  - When pi0 < 0.9 (many true effects expected)\n\n')

cat('GWAS:\n')
cat('  - Use genome-wide threshold p < 5e-8\n')
cat('  - Bonferroni for ~1 million independent tests\n')

Multiple Testing Correction Usage Guide

Overview

This guide covers applying multiple testing corrections including FDR, Bonferroni, and q-value methods.

Prerequisites

# R/Bioconductor
install.packages('BiocManager')
BiocManager::install('qvalue')

# Python
pip install statsmodels scipy

Quick Start

Tell your AI agent what you want to do:

"Apply FDR correction to my differential expression p-values"
"Which multiple testing correction should I use for my GWAS results?"
"Calculate q-values for my DE results"
"Filter my results at FDR < 0.05"

Example Prompts

Differential Expression

"I have p-values from DESeq2 for 20,000 genes. Apply Benjamini-Hochberg correction and filter at FDR 0.05"

"Compare the number of significant genes using Bonferroni vs BH correction"

Method Selection

"I'm doing an exploratory analysis. Should I use FDR 0.05 or 0.10?"

"What's the difference between adjusted p-value and q-value?"

GWAS

"Apply genome-wide significance threshold to my GWAS results"

What the Agent Will Do

1. Identify the analysis context 2. Select appropriate correction method 3. Apply correction to p-values 4. Report number of significant results 5. Explain interpretation of corrected values