Differential Expression Analysis in Proteomics: A Complete R Pipeline (limma, t-test, ANOVA)

Differential expression (DE) analysis is the core step in any quantitative proteomics experiment: which proteins are significantly different between conditions? Having built a complete DE pipeline from scratch (for BioAI Market), here's everything I learned about doing it properly.

Overview: Which Statistical Test to Use?

The right test depends on your data:

Scenario	Recommended Method	Why
2 groups, ≥3 replicates	limma (R)	Standard choice: empirical Bayes moderated t-test
2 groups, ≥3 replicates, no R	Welch's t-test	Robust to unequal variance
2 groups, small n or non-normal	Mann-Whitney U	Non-parametric, no normality assumption
3+ groups	limma (moderated F-test + contrasts)	Same variance moderation as the 2-group case; ANOVA + Tukey HSD only if you can't use R
2 groups, equal variance confirmed	Student's t-test	Only if Levene's test confirms homoscedasticity

Rule of thumb: If you have R available, use limma, for two groups or more. It borrows information across proteins to stabilize variance estimates, which is especially important in proteomics where replicates are often limited.

Step 1: Data Preprocessing

Before any statistical test, preprocess your expression matrix:

Revised October 2026. The order of the steps below changed (log-transform before imputing), the per-protein normality "diagnostics" were removed, and every R block in this post was run on a real MaxQuant 2.8.1 protein table (the HYE benchmark from this post).

Log2 transform first

# data: proteins x samples, raw intensities (e.g. MaxQuant LFQ), 0 = not quantified
data[data == 0] <- NA
data_log <- log2(as.matrix(data))

Log-transform before anything else. Imputation methods such as MinProb assume log-scale data; imputing raw intensities and then taking log2 can produce values ≤ 0 and therefore -Inf/NaN.

Filtering (valid values per group)

group <- factor(c(rep("Control", 3), rep("Treatment", 3)))
min_per_group <- 2    # e.g. n - 1 for n = 3
keep <- sapply(levels(group), function(g) rowSums(!is.na(data_log[, group == g])) >= min_per_group)
data_filtered <- data_log[apply(keep, 1, all), ]
cat(sprintf("Kept %d / %d proteins\n", nrow(data_filtered), nrow(data_log)))

Filter per group, not over all samples: a 70%-of-all-samples rule can keep proteins that were never seen in one of the groups.

Normalization

# Median normalization: a common default, not the only valid choice
col_medians <- apply(data_filtered, 2, median, na.rm = TRUE)
data_norm <- sweep(data_filtered, 2, col_medians - median(col_medians))

Imputation (optional, on the log scale)

# limma can work with the remaining NAs, so imputation is a choice, not a requirement.
# If you impute, do it on log2 values:
library(imputeLCMD)
data_imputed <- impute.MinProb(data_norm, q = 0.01)   # left-censored (MNAR) assumption

The trade-offs between imputation methods are covered in missing values in proteomics.

Step 2: Look at the Samples, Not Per-Protein Tests

An earlier version of this guide ran a Shapiro-Wilk and a Levene test on every protein and used the results to pick a method. With three or four values per group, those tests have almost no power, so they say little about your data. Look at sample-level plots instead:

library(limma)
boxplot(data_norm, las = 2, main = "After median normalization")   # are the distributions aligned?
plotMDS(data_norm[complete.cases(data_norm), ], col = as.integer(group))   # do replicates cluster by group?

After fitting limma (next step), plotSA(fit) shows the mean-variance trend that the empirical Bayes step relies on.

Step 3A: limma (Recommended for 2-Group Comparison)

library(limma)

# Create design matrix
group <- factor(c(rep("Control", 3), rep("Treatment", 3)))
design <- model.matrix(~ 0 + group)
colnames(design) <- levels(group)

# Create contrast
contrast <- makeContrasts(Treatment - Control, levels = design)

# Fit model
fit <- lmFit(data_norm, design)
fit2 <- contrasts.fit(fit, contrast)
fit3 <- eBayes(fit2)

# Extract results
results <- topTable(fit3, number = Inf, sort.by = "none")
results$protein_id <- rownames(results)   # e.g. UniProt accessions, if those are your row names

# Add significance annotation
results$significant <- abs(results$logFC) > 1 & results$adj.P.Val < 0.05
results$regulation <- ifelse(!results$significant, "NS",
                      ifelse(results$logFC > 0, "Up", "Down"))

cat(sprintf("Significant: %d (%d up, %d down)\n",
    sum(results$significant),
    sum(results$regulation == "Up"),
    sum(results$regulation == "Down")))

Why limma Is the Standard Choice for Proteomics

1. Empirical Bayes moderation — Shrinks protein-level variance toward a common prior, solving the "small n" problem
2. Works with missing values — lmFit fits each protein on the values it has. It does not model why values are missing (e.g. below detection), so very incomplete proteins should still be filtered
3. Multiple testing built in — adj.P.Val uses BH correction by default
4. Widely used — published by Smyth in 2004 and maintained in Bioconductor ever since

Step 3B: Welch's t-test (Python Alternative)

If R is not available:

from scipy import stats
from statsmodels.stats.multitest import multipletests
import numpy as np

def welch_ttest_de(data, protein_ids, control_idx, treatment_idx, log2fc_threshold=1.0, p_threshold=0.05):
    results = []
    p_values = []

    for i in range(data.shape[0]):
        control = data[i, control_idx]
        treatment = data[i, treatment_idx]

        # Remove NaN
        control = control[~np.isnan(control)]
        treatment = treatment[~np.isnan(treatment)]

        if len(control) < 2 or len(treatment) < 2:
            continue

        log2fc = np.mean(treatment) - np.mean(control)  # Already log2
        _, pval = stats.ttest_ind(control, treatment, equal_var=False)

        results.append({
            'protein': protein_ids[i],
            'log2FC': log2fc,
            'pvalue': pval,
        })
        p_values.append(pval)

    # BH correction (CRITICAL — never skip this!)
    _, adj_pvals, _, _ = multipletests(p_values, method='fdr_bh')
    for i, r in enumerate(results):
        r['adj_pvalue'] = adj_pvals[i]
        r['significant'] = abs(r['log2FC']) >= log2fc_threshold and r['adj_pvalue'] < p_threshold   # threshold on the log2 scale

    return results

⚠️ CRITICAL: Always apply multiple testing correction (BH/FDR). With 5000+ proteins tested, you'll get hundreds of false positives at p < 0.05 without it.

Step 3C: Three or More Groups (limma)

For more than two groups, stay in limma: a moderated F-test asks "does this protein differ anywhere?", and contrasts give the pairwise comparisons, all with the same variance moderation.

library(limma)
group3 <- factor(c("A","A","A","B","B","B","C","C","C"))
design <- model.matrix(~ 0 + group3); colnames(design) <- levels(group3)
contr <- makeContrasts(BvsA = B - A, CvsA = C - A, CvsB = C - B, levels = design)
fit3g <- eBayes(contrasts.fit(lmFit(data_norm3, design), contr))   # data_norm3: proteins x 9 samples

anyDiff <- topTable(fit3g, number = Inf, sort.by = "none")       # moderated F-test across all contrasts
BvsA    <- topTable(fit3g, coef = "BvsA", number = Inf, sort.by = "none")
summary(decideTests(fit3g, method = "separate", p.value = 0.05, lfc = 1))

[An earlier version used a hand-written Python ANOVA with Tukey HSD and reported a 215 s → 8 s speed-up. The code had bugs (it returned an undefined variable) and the timing had no data behind it, so both are gone.]

Without R, scipy.stats.f_oneway per protein, BH correction with statsmodels, and Tukey HSD (statsmodels.stats.multicomp.pairwise_tukeyhsd) on the significant proteins is a workable fallback, but you lose the variance moderation.

Step 4: Volcano Plot

library(ggplot2)

ggplot(results, aes(x = logFC, y = -log10(adj.P.Val), color = regulation)) +
  geom_point(alpha = 0.6, size = 1.5) +
  scale_color_manual(values = c("Up" = "#e74c3c", "Down" = "#3498db", "NS" = "#bdc3c7")) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "gray40") +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "gray40") +
  labs(
    title = "Volcano Plot — Differential Expression",
    x = "Log2 Fold Change",
    y = "-Log10 Adjusted P-value"
  ) +
  theme_minimal() +
  annotate("text", x = max(results$logFC) * 0.8, y = max(-log10(results$adj.P.Val)) * 0.9,
           label = sprintf("%d Up", sum(results$regulation == "Up")), color = "#e74c3c") +
  annotate("text", x = min(results$logFC) * 0.8, y = max(-log10(results$adj.P.Val)) * 0.9,
           label = sprintf("%d Down", sum(results$regulation == "Down")), color = "#3498db")

Step 5: Pathway Enrichment

After DE, feed significant proteins into pathway analysis:

library(clusterProfiler)
library(org.Hs.eg.db)

# Convert gene symbols to ENTREZ IDs
sig_genes <- results$gene_symbol[results$significant]

ego <- enrichGO(
  gene = sig_genes,
  OrgDb = org.Hs.eg.db,
  keyType = "SYMBOL",
  ont = "BP",
  pAdjustMethod = "BH",
  pvalueCutoff = 0.05
)

# KEGG pathways
gene_entrez <- bitr(sig_genes, fromType = "SYMBOL",
                     toType = "ENTREZID", OrgDb = org.Hs.eg.db)

ekegg <- enrichKEGG(
  gene = gene_entrez$ENTREZID,
  organism = "hsa",
  pAdjustMethod = "BH"
)

If Your IDs are UniProt Accessions

# UniProt's own accession pattern (a plain [A-Z][0-9A-Z]{5} also matches gene symbols like COL1A1)
uniprot_re <- "^([OPQ][0-9][A-Z0-9]{3}[0-9]|[A-NR-Z][0-9]([A-Z][A-Z0-9]{2}[0-9]){1,2})(-[0-9]+)?$"
sig_ids <- results$protein_id[results$significant]
if (all(grepl(uniprot_re, sig_ids))) {
  ego <- enrichGO(gene = sub("-[0-9]+$", "", sig_ids),   # drop isoform suffixes
                  OrgDb = org.Hs.eg.db, keyType = "UNIPROT", ont = "BP")
}

Common Mistakes to Avoid

1. No Multiple Testing Correction

# WRONG — thousands of false positives
significant <- results$P.Value < 0.05

# RIGHT — BH-corrected
significant <- results$adj.P.Val < 0.05

2. Using Log2FC Threshold Without P-value

# WRONG — large fold changes can be noise
significant <- abs(results$logFC) > 1

# RIGHT — both thresholds
significant <- abs(results$logFC) > 1 & results$adj.P.Val < 0.05

3. Student's t-test Instead of Welch's

Use Welch's t-test (equal_var=False) by default. It costs almost nothing when variances are equal and protects you when they are not.

4. Not Log-Transforming

Raw LFQ intensities span orders of magnitude. Statistical tests assume roughly normal distributions. Log2 transformation is essential.

5. Ignoring Missing Data Patterns

Proteins missing in one condition but present in another may be the most biologically interesting! Consider "missing not at random" (MNAR) imputation methods.

Complete Minimal Pipeline

# 1. Load & preprocess
data <- read.csv("protein_matrix.csv", row.names = 1)
data[data == 0] <- NA
data_log <- log2(as.matrix(data))
col_med <- apply(data_log, 2, median, na.rm = TRUE)
data_norm <- sweep(data_log, 2, col_med - median(col_med))

# 2. limma DE
library(limma)
group <- factor(c(rep("Control", 3), rep("Treatment", 3)))
design <- model.matrix(~ 0 + group)
colnames(design) <- levels(group)
fit <- eBayes(contrasts.fit(lmFit(data_norm, design),
              makeContrasts(Treatment - Control, levels = design)))
res <- topTable(fit, number = Inf, sort.by = "none")
res$sig <- abs(res$logFC) > 1 & res$adj.P.Val < 0.05

# 3. Results
res <- res[!is.na(res$adj.P.Val), ]   # drop proteins with too few values to test
cat(sprintf("%d significant (%d up, %d down) out of %d tested\n",
    sum(res$sig),
    sum(res$sig & res$logFC > 0),
    sum(res$sig & res$logFC < 0),
    nrow(res)))

Tools That Do This For You

If you want a GUI instead of code:

• BioAI Market — our own platform
• Perseus — Max Planck's desktop tool
• Proteome Discoverer — Thermo Fisher's commercial solution
• MSstats — R/Bioconductor package specifically for MS proteomics

---

This guide grew out of building the analysis pipeline for BioAI Market, our own platform.

Code and outputs: the scripts and result tables for this post are on GitHub in sbmlab-benchmarks/de-pipeline-check.