ProteoForge: What Happens When You Stop Averaging Your Peptides

I spent the better part of my PhD staring at peptide intensity matrices. Thousands of rows, dozens of columns, lots of missing values. The standard workflow says: roll those peptides up into protein-level quantities, run your differential expression, make a volcano plot, write the paper ◆ If you have never seen a peptide intensity matrix, imagine a spreadsheet where most of the cells in columns 5 through 12 are empty, and the ones that are not empty disagree with each other about what the protein is doing. That is bottom-up proteomics. . I did that. Multiple times. It always felt like we were throwing away information.

This post is about ProteoForge, the framework that came out of that frustration. It is not a summary of the paper (you can read the preprint for that). It is about the problem as I experienced it, the decisions we made while building the tool, and what surprised us when we finally applied it to real data.

The averaging problem

The human genome codes for roughly 20,000 proteins. But the actual number of distinct protein forms in a cell is orders of magnitude larger¹. Alternative splicing, post-translational modifications, proteolytic processing: these create proteoforms², and proteoforms from the same gene can have opposing biological functions.

Bottom-up proteomics measures peptides, not intact proteins. The field has lived with this trade-off for years: you get broad coverage but you lose proteoform identity. The standard fix is to group peptides by protein accession, average their intensities, and call it a protein quantity.

That averaging step is where the information dies.

If three peptides from a protein go up under treatment and two go down, the protein-level average might show no change at all. A flat line on your volcano plot. You move on. You never see that the protein was actually doing something complicated at the proteoform level.

What already existed (and where it broke)

The core issue with both tools, and really with every method we tested: missing data ◆ A typical DIA dataset from an Orbitrap will have about 25% missing values at the precursor level. By the time you filter for valid peptides per protein, that number often exceeds 40%. Most statistical methods assume complete data and crash or bias badly when confronted with this. . In a typical DIA proteomics experiment, 20 to 40 percent of peptide measurements are missing. Some are missing at random (instrument did not pick them up). Some are missing because the peptide is genuinely absent in that condition. The distinction matters enormously, and most tools either ignore it or require you to delete incomplete cases.

That was the gap.

The idea behind ProteoForge

ProteoForge does four things, in order. I will walk through each one, not with equations (the paper has those), but with the reasoning.

Figure 1: The four-module ProteoForge pipeline. (A) Framework schematic. (B) Effect of normalization steps on peptide distributions. (C) Three example proteins showing peptide profiles, distance heatmaps, and final dPF assignments.

Module 1: Normalize relative to a control

Raw peptide intensities are noisy and have different baselines across samples. We log-transform, z-score, and then subtract the mean control intensity for each peptide. After this step, every value represents the deviation of that peptide from its own control baseline. Sounds simple. It is. But it makes the downstream models much cleaner because you are comparing changes, not absolute values.

Module 2: Find discordant peptides

This is the statistical core. For each protein, we fit a linear model that asks: does this peptide behave differently from its siblings across conditions?

Intensity ~ Condition * Peptide

The interaction term (Condition * Peptide) captures the discordance. If a peptide responds to the condition differently than the other peptides from the same protein, the interaction term will be significant.

We tested several model types. Ordinary least squares, weighted least squares with custom imputation-aware weights, generalized linear models, and robust linear models (RLM) with Huber M-estimation³.

We settled on RLM as the default because of a property that turned out to be exactly what we needed: it automatically down-weights outliers. Imputed values, especially the down-shifted low values used for condition-complete missingness, look like outliers to the model. RLM treats them accordingly. No manual weight specification needed. The user does not have to think about which values were imputed and which were real.

The p-values from the interaction terms go through a two-step Benjamini-Hochberg FDR correction⁴: first within each protein, then globally ◆ The first correction accounts for multiple peptides per protein. The second accounts for multiple proteins genome-wide. A miscleaved peptide or technical artifact would have to survive both rounds, which is why we use this approach despite it being conservative. . Peptides that pass the threshold (we used adjusted p < 0.001) are flagged as significantly discordant.

Module 3: Cluster the discordant peptides

Finding discordant peptides is not enough. You need to know which discordant peptides move together, because a proteoform is defined by a group of co-varying peptides, not a single one.

We compute the Euclidean distance between median adjusted intensity profiles and cluster with Ward linkage⁵. The default cut method (hybrid_outlier_cut) determines the number of clusters automatically.

This is the part that separates ProteoForge from PeCorA. PeCorA stops at "this peptide is different." ProteoForge goes further: "these three peptides are different in the same way, and they likely belong to the same proteoform."

Module 4: Build proteoforms

The naming convention tells the story:

The result is a peptide-to-dPF mapping table. Each protein can have zero, one, or multiple differential proteoforms, plus a canonical form.

Benchmarking: what we actually learned

Benchmarking a proteoform discovery tool is tricky because there is no ground truth for proteoforms in most real datasets. We used two approaches: real data with in-silico perturbations (the COPF benchmark approach, using a SWATH-MS interlab HEK293 dataset⁶⁷), and controlled simulations.

Figure 2: Benchmark results across real and simulated data. (A) SWATH-MS interlab test. (B-E) Simulations 1-4: imputation impact, missingness tolerance, signal strength, multi-condition complexity.

Four simulation scenarios tested the things that matter in practice:

For peptide grouping (which PeCorA cannot do), ProteoForge consistently hit MCC around 0.45 to 0.61, while COPF hovered around 0.27. At a p-value cutoff of 0.001, COPF's grouping MCC was near zero in most simulated conditions.

One result surprised us. In the imputation simulation, ProteoForge's grouping performance actually improved slightly after imputation (MCC from 0.60 to 0.61). Why? When a peptide is completely missing in one condition and gets imputed with a down-shifted low value, it creates a strong, unambiguous signal. The clustering picks it up correctly.

The hypoxia application

Benchmarks are one thing. Real biology is another.

We applied ProteoForge to a published dataset from Tomin et al. (2025)⁸, who measured the proteome of H358 lung cancer cells under normoxia (21% O2) and hypoxia (1% O2) at 48 and 72 hours. The data came from DIA-NN⁹ as precursor-level quantities. After filtering and imputation, we had 95,369 peptides across 7,161 proteins.

The numbers alone tell a story. At 48 hours, about 15% of proteins had multiple discordant peptides. At 72 hours, that jumped to 36%. Hypoxia is not a subtle perturbation, and the proteoform-level response grows over time in a way that protein-level analysis completely misses. Figure 3A shows this distribution shift from 48 to 72 hours. The fraction of proteins with multiple significantly discordant peptides more than doubles.

The original study by Tomin et al.⁸ reported increased levels of S-lactoylglutathione (a GLO1 product) under hypoxia, but no change in GLO1 protein abundance. That is a contradiction. If the product goes up, something about the enzyme should be different.

ProteoForge resolved it.

At 48 hours, one peptide in GLO1 (harboring a known phosphorylation site at T35 and a mutagenesis site at Q34) showed elevated intensity under hypoxia while the other peptides stayed flat. By 72 hours, this expanded into a multi-peptide dPF (three peptides, including one with ubiquitination/acetylation at K148). The protein average saw nothing because the canonical peptides diluted the signal.

This is what I mean by "the averaging problem." GLO1 was regulated. The original analysis just could not see it.

FPGT: catching isoform-driven variance

Fucose-1-phosphate guanylyltransferase (FPGT) provided a different kind of example. At 48 hours, all peptides behaved the same. At 72 hours, two peptides spanning amino acids 60 to 95 dropped under hypoxia while the rest stayed stable.

Figures 3F and 3G show the FPGT peptide profiles and protein schematic, with the two discordant peptides highlighted in the N-terminal region and isoform boundaries marked.

What happens to your pathway analysis

This was, for me, the most striking result.

We ran four different protein-level quantification strategies through QuEStVar¹⁰ (our equivalence testing framework from the earlier paper) and then into pathway enrichment with g:Profiler¹¹:

That is not a marginal improvement. That is the difference between "hypoxia did almost nothing to this proteome" and "hypoxia restructured a quarter of the proteome at the proteoform level."

The pathway results were telling. The original analysis found broad equivalence across GO categories at 72 hours. The dPF analysis found canonical hypoxia response pathways correctly classified as elevated. Cell cycle terms were reduced. Metabolic terms showed complex, time-dependent shifts. Signaling terms appeared elevated.

The "Original" analysis missed all of this. Not because the data was bad. Because the quantification averaged it away.

Things I wish I had known earlier

A few practical notes for anyone thinking about using ProteoForge.

It is not a package yet. The analysis repository on GitHub contains the scripts and notebooks used for the manuscript. A standalone Python package is under development but not released. If you use the current code, expect some rough edges.

Why this matters (to me)

I started this work because of pediatric cancer. My PhD focused on computational proteomics for childhood leukemia, and the proteoform question kept coming up. In our leukemia data, we would see proteins where the overall level looked stable between diagnosis and relapse, but specific peptides shifted dramatically. We could not make sense of it with protein-level tools.

ProteoForge came from that frustration. The hypoxia application in the paper is a clean, well-controlled dataset that made for a good proof of concept. But the motivation was always translational: can we find proteoform-level changes in disease that protein-level analysis misses?

The GLO1 result gives me confidence that the answer is yes.

What comes next

The preprint is currently under review at the Journal of Proteome Research. The standalone Python package is in development.

If you want to try it on your data, the analysis repository has everything. The Zenodo snapshot (doi:10.5281/zenodo.17795845) freezes the exact version used in the manuscript.

Paper: Ergin, E. K., Conrrero, A., Ferguson, K. M., Lange, P. F. (2025). ProteoForge: An Imputation-Aware Framework for Differential Proteoform Discovery in Bottom-Up Proteomics. bioRxiv. doi:10.64898/2025.12.12.694008

Code: github.com/LangeLab/ProteoForge_Analysis

Data: PRIDE repository, PXD062503 (Tomin et al. hypoxia dataset)