How to Run Enrichment Analysis on Your Own Proteomics Data

Oak
proteomics enrichment gsea clusterprofiler bioinformatics
A glowing vertical column of colored beads inside a lab instrument, with bright clustered bands lit by a scanning light.

By the end of this guide you will have a ranked list of proteins from your own plasma or tissue proteomics, a set of enrichment results from both over-representation analysis and gene set enrichment analysis, and a collapsed network of non-redundant pathway terms you can read. You need three things to start: a protein abundance matrix (Olink NPX, SomaScan RFU, or MaxQuant/DIA-NN/Spectronaut protein-level intensities), the identifier column that came with it (UniProt accessions are best), and R 4.3 or later with clusterProfiler, fgsea, limma, msigdbr, and org.Hs.eg.db installed. Enrichment analysis asks a narrow question: given a list or ranking of proteins, are members of some predefined biological set concentrated near the top or bottom more than chance would produce? That is the whole idea, and almost every way it goes wrong is a violation of the word “chance” in that sentence.

1. Load the matrix and fix identifiers first

Identifier mapping is the step that silently destroys the most analyses, so do it before anything else and count what you lose. Affinity platforms report one row per assay target, and Olink gives you a UniProt accession per assay; SomaScan gives you an aptamer ID (SeqId) plus a UniProt mapping in the ADAT metadata. Mass spectrometry search engines report protein groups, which frequently contain several semicolon-separated accessions where the peptides could not distinguish isoforms or family members.

library(tidyverse); library(org.Hs.eg.db); library(AnnotationDbi)

prot <- readr::read_tsv("report.pg_matrix.tsv")   # DIA-NN protein groups

prot <- prot |>
  mutate(uniprot = str_split_i(Protein.Group, ";", 1)) |>   # keep leading accession
  filter(!str_detect(Protein.Group, "^(CON__|REV__)"))      # contaminants, decoys

map <- AnnotationDbi::select(org.Hs.eg.db,
                             keys    = unique(prot$uniprot),
                             keytype = "UNIPROT",
                             columns = c("ENTREZID", "SYMBOL"))

# one Entrez ID per protein group; drop ambiguous multimaps
map <- map |> group_by(UNIPROT) |> filter(n() == 1) |> ungroup()
message(nrow(map), " of ", length(unique(prot$uniprot)), " accessions mapped")

Expect to lose somewhere between two and ten percent of accessions to obsolete entries, non-human contaminants, and one-to-many mappings. If you lose more than that, your file has a version problem: UniProt accessions get merged and demerged across releases, and an old MaxQuant FASTA against a current org.Hs.eg.db will drop a visible chunk. Record the mapped set. It becomes your statistical background in step 3.

2. Build a ranking statistic that reflects your design

Enrichment consumes either a list (over-representation) or a full ranking (GSEA), and the quality of the ranking determines everything downstream. For a conventional two-group comparison, use limma on log2-transformed abundances. Olink NPX is already on a log2 scale; SomaScan RFU and MS intensities are not, so log-transform them and, for MS data, normalize by median or by limma::normalizeBetweenArrays(method = "quantile") before modeling.

library(limma)
X <- as.matrix(log2(intensities))          # proteins x samples
design <- model.matrix(~ group)
fit <- eBayes(lmFit(X, design), trend = TRUE, robust = TRUE)
res <- topTable(fit, coef = 2, number = Inf)

For a single person measured repeatedly, there is no group to contrast, and this changes the statistic rather than the method. Two designs work. If you have a perturbation with before and after timepoints, block on the subject and fit a paired model with duplicateCorrelation(), which is the correct way to handle repeated measures in limma. If you have a longitudinal series with no discrete perturbation, compute a within-protein z-score at each timepoint against that protein’s own historical mean and standard deviation, then rank by the z at the timepoint of interest. The second approach is weaker statistically because a single person gives you one draw per protein per timepoint, but it has the advantage of controlling for the enormous between-person variation in plasma protein concentrations that dominates any cross-sectional reference range.

Rank by the moderated t statistic, or by sign(logFC) * -log10(P.Value). Do not rank by fold change alone: low-abundance proteins near the detection floor produce large ratios with no evidence behind them, and they will sit at the extremes of your list where GSEA weights them most heavily.

3. Choose the background, and be strict about it

This is the single most consequential decision in a proteomics enrichment analysis, and it is the one most often made by accident. Over-representation analysis computes a hypergeometric p-value: given N proteins in the universe, K of which belong to a pathway, and n proteins in your significant list, how surprising is it to find k pathway members? If you set N to all ~20,000 human protein-coding genes when your assay only measured 1,200 of them, every term whose members happen to be well-represented in the measurable plasma proteome becomes significant. Plasma mass spectrometry without depletion is heavily biased toward secreted, complement, coagulation, and apolipoprotein families, so you will discover, with tiny p-values, that your data are enriched for complement activation. They are, and it tells you nothing about your biology.

Set the universe to the mapped, quantified proteins that passed your filters.

universe <- map$ENTREZID          # everything you could have detected
sig      <- res |> filter(adj.P.Val < 0.05, abs(logFC) > 0.3) |> pull(entrez)

The same reasoning applies to affinity panels with even more force, because the panel content is a deliberate curatorial choice. An Olink inflammation panel is enriched for inflammation gene sets by construction. Reviews of pathway analysis practice repeatedly flag background selection as a primary source of false enrichment, and it is worth treating as a hard requirement rather than an option.12

4. Run over-representation analysis

With a defensible list and universe, over-representation analysis is a few lines. We use clusterProfiler because it handles GO, KEGG, Reactome, and arbitrary MSigDB collections through one interface and returns a consistent result object.3

library(clusterProfiler); library(msigdbr)

ego <- enrichGO(gene          = sig,
                universe      = universe,
                OrgDb         = org.Hs.eg.db,
                keyType       = "ENTREZID",
                ont           = "BP",
                pAdjustMethod = "BH",
                pvalueCutoff  = 0.05,
                qvalueCutoff  = 0.2,
                minGSSize     = 15,
                maxGSSize     = 500,
                readable      = TRUE)

head(as.data.frame(ego)[, c("Description","GeneRatio","BgRatio","p.adjust","Count")])

The minGSSize and maxGSSize bounds matter. Sets smaller than about 15 members produce unstable p-values driven by two or three proteins, and sets larger than 500 are so general (“metabolic process”) that significance carries no information. Fifteen to five hundred is the convention used in the GSEA literature and we see no reason to deviate.2

Read GeneRatio and BgRatio together, never p.adjust alone. A term with GeneRatio = 8/120 and BgRatio = 30/1100 has a fold enrichment of about 2.4. A term with three of your proteins out of a five-member set may have a beautiful p-value and describe an artifact of a single protein complex being co-quantified.

5. Run rank-based GSEA on the full list

Over-representation throws away the ranking and imposes an arbitrary significance cutoff. Rank-based gene set enrichment analysis uses every protein, which matters in proteomics where coverage is limited and coordinated small shifts across a pathway are common and individually non-significant.4

library(fgsea)

stats <- res |>
  mutate(score = sign(logFC) * -log10(P.Value)) |>
  filter(!is.na(entrez)) |>
  arrange(desc(score))
ranks <- setNames(stats$score, stats$entrez)

msig <- msigdbr(species = "Homo sapiens", collection = "H")   # hallmark
pathways <- split(as.character(msig$ncbi_gene), msig$gs_name)

set.seed(42)
fg <- fgseaMultilevel(pathways = pathways, stats = ranks,
                      minSize = 15, maxSize = 500, eps = 0)
fg <- fg[order(fg$padj), ]

Two things about this ranking. Ties are bad: if many proteins share an identical score, the enrichment score walk becomes order-dependent, and fgsea will warn you. Use a continuous statistic rather than a rounded one. And restrict pathways to genes present in ranks, which fgsea does internally when it computes set sizes, but check the reported size column to confirm that a “significant” 400-member hallmark set was not represented by 18 measured proteins.

Start with the MSigDB hallmark collection (50 sets, deliberately non-redundant) before touching GO biological process (over 12,000 terms, deeply nested, heavily overlapping). If hallmark gives you nothing, GO will give you a thousand correlated terms that also mean nothing, at greater cost in multiple testing.

6. Read the sign and the magnitude correctly

The normalized enrichment score (NES) is the enrichment score for a set, divided by the mean enrichment score of permutations of the same set size, which makes scores comparable across sets. Its sign is defined entirely by where the set’s members sit in the ranking you supplied. A positive NES means the proteins in that set are concentrated at the top of your ranked list. A negative NES means they are concentrated at the bottom.

What “top” and “bottom” mean is your choice, not the algorithm’s, and this is the origin of most confusion about negative scores. If you ranked by sign(logFC) * -log10(p) for a post-versus-pre contrast, a negative NES means that set’s proteins are lower after than before. It does not mean the pathway is inhibited, and it does not mean a biological process was switched off. Plasma protein abundance reflects secretion, clearance, and dilution as much as it reflects transcriptional activity of the pathway named in the set. Nothing in the enrichment statistic distinguishes those causes.

Magnitude is easier to over-read than the sign. NES above roughly 2 in absolute value with an adjusted p below 0.05 is a strong signal in a well-powered comparison. Focus on the leading edge, the subset of proteins that drove the score, because that is what you can go back and inspect at the level of individual measurements.

le <- fg[fg$pathway == "HALLMARK_COMPLEMENT", ]$leadingEdge[[1]]
AnnotationDbi::mapIds(org.Hs.eg.db, le, "SYMBOL", "ENTREZID")

If three proteins carry a 2.1 NES across a 60-member set, you have a three-protein finding wearing a pathway’s name.

7. Collapse redundancy before you interpret anything

GO enrichment returns parent and child terms that share most of their genes, and a table of 200 significant terms is usually 12 distinct signals repeated. Two tools solve this. clusterProfiler::simplify() removes GO terms whose semantic similarity exceeds a cutoff, keeping the more significant member of each pair. For a broader view, export to Cytoscape and build an EnrichmentMap, where nodes are gene sets, edges are gene overlap, and clusters of nodes become the units you interpret. The protocol by Reimand and colleagues walks through this end to end and is the reference we follow.2

ego2 <- simplify(ego, cutoff = 0.7, by = "p.adjust", select_fun = min)
# fgsea results: collapse leading-edge-redundant pathways
keep <- collapsePathways(fg[fg$padj < 0.05, ], pathways, ranks)
fg_main <- fg[fg$pathway %in% keep$mainPathways, ]

Name each resulting cluster in plain language and treat the cluster, not the individual term, as the finding.

8. Run the negative controls

Before you believe anything, break the analysis on purpose. Permute your sample labels 10 times, rerun the full pipeline, and count significant terms at your chosen threshold. You should see close to zero. If permuted labels reliably produce significant enrichment, your ranking statistic is confounded, most often by batch, plate, or sample-processing date, all of which affect plasma proteomics strongly. Second, replace your gene sets with size-matched random sets drawn from your universe and confirm the NES distribution centers near zero.

For the longitudinal single-person case, the equivalent control is to compare two timepoints where nothing was different and see what your pipeline claims. Within-person plasma proteomics has meaningful technical and diurnal variation, and a baseline-versus-baseline comparison tells you the floor of what counts as a signal in your own data.

Common problems

Coverage bias masquerading as biology. Undepleted plasma mass spectrometry typically quantifies a few hundred to roughly a thousand protein groups, dominated by high-abundance secreted proteins; affinity platforms such as Olink Explore and SomaScan reach several thousand targets but only the targets they were designed for. Either way, your universe is a biased sample of the proteome, and enrichment can only speak about that sample. Say so when you write results down. Reviews of MS-based multi-omics discuss this dynamic-range constraint directly.5

Pathway databases disagree. KEGG, Reactome, GO, and WikiPathways annotate the same gene differently, and the same query can be significant in one and absent in another. Report the database and version. The reproducibility problem here is well documented, and a term that only appears in one source deserves less weight.16

Mixing over-representation results with GSEA results in one table. They answer different questions with different nulls. Present them separately and note where they agree.

Aptamer and antibody specificity. A SomaScan or Olink signal is a binding measurement, not a direct concentration, and cross-reactivity or a coding variant in the epitope can shift a single target substantially. If one protein drives your leading edge, verify it on an orthogonal platform before building interpretation on it.

Treating enrichment as mechanism. Enrichment is a statement about the coordinates of your gene list within an annotation database. It is a hypothesis generator, and pathway-level results from proteomics, metabolomics, and transcriptomics of the same samples frequently converge only partially.789 Where a result touches a clinical question, take the underlying protein measurements and the analysis to a physician. Enrichment output is not a diagnosis and cannot support one.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

  1. Davide Chicco, Giuseppe Agapito. Nine quick tips for pathway enrichment analysis. PLOS Computational Biology, 2022. https://doi.org/10.1371/journal.pcbi.1010348 ↩ ↩2

  2. Jüri Reimand, Ruth Isserlin, Veronique Voisin, et al. Pathway enrichment analysis and visualization of omics data using g:Profiler, GSEA, Cytoscape and EnrichmentMap. Nature Protocols, 2019. https://doi.org/10.1038/s41596-018-0103-9 ↩ ↩2 ↩3

  3. Tianzhi Wu, Erqiang Hu, Shuangbin Xu, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation, 2021. https://doi.org/10.1016/j.xinn.2021.100141 ↩

  4. Kangmei Zhao, Seung Yon Rhee. Interpreting omics data with pathway enrichment analysis. Trends in Genetics, 2023. https://doi.org/10.1016/j.tig.2023.01.003 ↩

  5. Andrew T. Rajczewski, Pratik D. Jagtap, Timothy J. Griffin. An overview of technologies for MS-based proteomics-centric multi-omics. Expert Review of Proteomics, 2022. https://doi.org/10.1080/14789450.2022.2070476 ↩

  6. William G. Ryan V., Smita Sahay, John Vergis, et al. Pathway Analysis Interpretation in the Multi-Omic Era. BioTech, 2025. https://doi.org/10.3390/biotech14030058 ↩

  7. Tien-Chueh Kuo, Tze-Feng Tian, Yufeng Jane Tseng. 3Omics: a web-based systems biology tool for analysis, integration and visualization of human transcriptomic, proteomic and metabolomic data. BMC Systems Biology, 2013. https://doi.org/10.1186/1752-0509-7-64 ↩

  8. Maurizio Bruschi, Simona Granata, Silvia Lai, et al. Multi-Omics Analysis of PBMCs Revealed Distinct Biological Differences Between Hemodialysis and Peritoneal Dialysis. International Journal of Molecular Sciences, 2026. https://doi.org/10.3390/ijms27167470 ↩

  9. Shijiani Li, Yilan Tang, Jie Fang, et al. Multi-omics profiling of migraine links plasma proteomic signatures with lipid metabolomic measures. The Journal of Headache and Pain, 2026. https://doi.org/10.1186/s10194-026-02493-x ↩