How to Normalize Your Own Proteomics Data

By the end of this guide you will have a normalized protein abundance matrix for your own samples. It will have proteins in rows and draws in columns on a log2 scale, with plate and run-date effects removed. You will also have a written record of which normalization you chose and the evidence you chose it on. You will also have a per-protein estimate of your own technical and biological variation, which is what makes a second draw interpretable relative to a first. You need your platform’s output file, the sample metadata that came with it, and R 4.3 or later. The output file is a DIA-NN report.pg_matrix.tsv, a MaxQuant proteinGroups.txt, a Spectronaut protein report, an Olink NPX table, or a SomaScan .adat. The metadata covers plate, well, run date, acquisition order, and any bridging or control samples. You also need the following packages.

• limma
• vsn
• preprocessCore
• NormalyzerDE
• sva
• data.table

Budget an afternoon. Deciding what to trust takes the time. The analysis itself is not hard.

1. Find out what normalization your platform already applied

Normalization exists because the number you get for a protein is a product of the biology you want and a stack of things you do not want. The stack includes how much peptide was loaded, how the autosampler behaved that morning, where the sample sat on the plate and how the mass spectrometer’s response drifted over a 30-hour acquisition. The point is to make abundances comparable across samples while preserving real differences between them. Applying a second global correction on top of one the vendor already applied usually removes signal rather than noise, so start by reading the pipeline you were given.

For Olink Explore, the delivered NPX (Normalized Protein eXpression) value is already log2 scale and already plate-normalized: extension controls correct within-sample variation and plate controls or a set of bridging samples correct between-plate shifts. If you received an NPX file, your global normalization is largely done and your remaining job is bridging across shipments and quality control flags. For SomaScan, the .adat has passed hybridization control normalization, median signal normalization and plate scaling. It has usually also passed ANML, which is adaptive normalization by maximum likelihood, against a reference population. The relevant header fields are RowCheck, PlateId, and the per-sample scale factors. Samples whose scale factors fall outside roughly 0.4 to 2.5 are conventionally dropped.

Mass spectrometry output is the case where you do the work yourself. DIA-NN’s report.pg_matrix.tsv contains MaxLFQ intensities, which are already cross-run normalized in a relative sense but still carry loading and drift effects, particularly if you ran with --no-norm or across separated batches. MaxQuant’s proteinGroups.txt gives you raw intensity, LFQ intensity, and iBAQ columns. Use LFQ for label-free comparison across your draws and treat raw intensity as the diagnostic for loading differences. Isobaric labeling (TMT) has a distinct problem: each plex is internally ratio-comparable but each plex has its own reference, so plex-to-plex comparison requires a bridging channel, which is step 6.

2. Build the matrix and record the design

Everything downstream depends on having a clean matrix and a metadata table with the technical variables you might later need to model. Load the protein report, drop contaminants and decoys, and keep only the quantitative columns.

library(data.table)

pg <- fread("report.pg_matrix.tsv")
# DIA-NN: first five columns are identifiers
ids <- pg[, 1:5]
mat <- as.matrix(pg[, -(1:5)])
rownames(mat) <- ids$Protein.Group

# MaxQuant equivalent:
# pg <- fread("proteinGroups.txt")
# pg <- pg[Reverse != "+" & `Potential contaminant` != "+" & `Only identified by site` != "+"]
# mat <- as.matrix(pg[, grep("^LFQ intensity ", names(pg)), with = FALSE])

mat[mat == 0] <- NA

meta <- data.table(
  sample    = colnames(mat),
  draw_date = as.Date(c("2025-01-14", "2025-04-02", "2025-07-09", "2025-10-15")),
  plate     = c("P1", "P1", "P2", "P2"),
  well      = c("A01", "B01", "A01", "B01"),
  run_order = 1:4,
  digest_batch = c("D1", "D1", "D2", "D2")
)

Write down which technical variables are confounded with time. In a personal longitudinal series they almost always are, because your January draw was run in January. This is the central difficulty of self-quantification with proteomics: a drift between plates and a real change in your physiology produce the same pattern in the data. Bridging samples, described below, are the only clean way out. You have to plan them before the fact.

3. Filter on missingness before you normalize, and do not leave it until after

Global normalization methods estimate a per-sample statistic such as a median, a set of quantiles or a variance-stabilizing transform. Proteins present in only some samples bias that statistic. Filter first. We keep proteins quantified in at least 70% of samples for exploratory work and require complete cases for any protein we intend to track over time.

keep <- rowSums(!is.na(mat)) >= 0.7 * ncol(mat)
mat  <- mat[keep, ]
nrow(mat)  # typical plasma DIA run: 400-1200 protein groups after filtering

Do not impute at this stage. In label-free MS, missingness is a mixture of missing-not-at-random (below detection) and missing-at-random (stochastic sampling), and imputing before normalization propagates a guess into the scale factors. If you need a complete matrix for PCA, impute after normalization and keep the imputed version separate from the one you use for statistics. Tools such as PANDA-view implement several imputation and normalization combinations with paired visualization, which is useful for seeing how much of your structure is an artifact of the imputation choice.¹

4. Log transform, then apply the candidate normalizations

Protein intensities are right-skewed with variance that scales with the mean, so nearly every method assumes you work on the log scale. The formulas people search for are simple, and it is worth writing them out since the names get used loosely.

Log2 transformation: y_ij = log2(x_ij), for protein i in sample j.

Median centering, sometimes called global median or global normalization: subtract each sample’s median from every value in that sample, y'_ij = y_ij - median_i(y_ij). This assumes total protein load is the same across samples and that most proteins do not change. For plasma from one person over months, that assumption is reasonable.

Mean centering with scaling, sometimes called Z-score or standardization, is y'_ij = (y_ij - mean_j) / sd_j across samples for each protein. Use it for clustering and heatmaps. Do not use it as your quantitative normalization, because it discards the magnitude information you want when comparing draws.

Quantile normalization forces every sample to share one empirical distribution: rank the values within each sample, average the values at each rank across samples, and replace each value by the average for its rank. It is the most aggressive of the common options and correspondingly the most likely to erase a real global shift.

Variance stabilizing normalization (VSN) fits an arcsinh-family transform whose parameters are estimated from the data, producing values whose variance is roughly independent of the mean.

library(limma); library(vsn); library(preprocessCore)

log2mat <- log2(mat)

norm_median <- sweep(log2mat, 2, apply(log2mat, 2, median, na.rm = TRUE), "-")

norm_quant <- normalize.quantiles(log2mat, copy = TRUE)
dimnames(norm_quant) <- dimnames(log2mat)

norm_vsn <- justvsn(mat)          # takes raw intensities, returns ~log2 scale

norm_loess <- normalizeCyclicLoess(log2mat, method = "fast", iterations = 3)

Our default for a personal plasma MS series is median centering on log2 values, with cyclic loess as the alternative we check against. Median centering is the least destructive method that fixes the dominant problem (loading and total-signal differences), and it leaves a genuine global change in plasma protein concentration visible instead of normalizing it away. Quantile normalization and VSN are stronger tools and both make sense for large cohorts. On four to twelve samples from one person they force agreement between distributions that carry information.

5. Choose between methods with evidence, and do not go by preference

The literature is consistent that no single normalization wins across datasets, and that the choice materially changes which proteins come out as different. In a SWATH-MS biomarker study, the normalization method altered both the number and the identity of the differentiating proteins reported, which is exactly the failure mode to guard against.² Broader evaluations across mass-spectrometry omics data reach the same conclusion. Method performance depends on the platform, the amount of missingness and the size of the true effect.³ Treat normalization as a parameter you tune with metrics. It is not a default you inherit.

The practical approach is to run the candidates and compare them on quantitative criteria. NormalyzerDE implements this directly. It applies a panel of normalizations and reports pooled coefficient of variation (PCV), pooled median absolute deviation (PMAD), pooled estimate of variance (PEV) and intragroup correlation. It also produces MA and RLE plots.⁴

library(NormalyzerDE)

# designMatrix needs 'sample' and 'group' columns; use plate or subject-state as group
normalyzer(jobName = "pmp_norm",
           designPath = "design.tsv",
           dataPath   = "raw_matrix.tsv",
           outputDir  = "norm_out")

Read the output this way. Lower PCV and PMAD within a set of samples you believe are similar (technical replicates, or two draws a week apart with no intervening change) indicate the method is removing technical variation. But a method can achieve low within-group variance by flattening everything, so pair the variance metrics with a check that known differences survive. Spike a positive control into your design if you can: two aliquots of the same draw run on different plates should converge after normalization, and two draws separated by a real intervention should not. Supervised and unsupervised machine learning has been used as a formal version of this two-sided check, evaluating whether a normalization improves class separation without inflating it.⁵

The diagnostic we look at first, before any metric, is the relative log expression (RLE) plot.

rle <- norm_median - rowMedians(norm_median, na.rm = TRUE)
boxplot(rle, las = 2, ylim = c(-1.5, 1.5),
        main = "RLE after median centering", ylab = "log2 ratio to protein median")
abline(h = 0, col = "red")

Every box should sit on zero with similar width. A sample whose box is shifted has residual loading error. A sample whose box is much wider has a degradation or digestion problem that no global normalization will fix, and the right move is to exclude it rather than to normalize harder.

6. Handle batch effects separately from normalization

Normalization aligns sample-level scale and distribution. Batch correction removes systematic differences between groups of samples processed together, and it is a distinct operation that comes after. Conflating the two is the most common mistake we see. Batch effects in mass spectrometry arise from plate, digestion day, column change, calibration state, and instrument, and they are large enough to dominate the first principal component if uncorrected.⁶ A multicentric high-sensitivity proteomics study showed that protocol and instrument contributions can be separated and modeled explicitly, which is the right conceptual frame even for a single-person series run across different sites or instrument generations.⁷

Start by quantifying the problem. Run PCA on the normalized matrix and check whether PC1 tracks plate or run order.

pc <- prcomp(t(na.omit(norm_median)), scale. = FALSE)
summary(pc)$importance[, 1:3]
plot(pc$x[, 1], pc$x[, 2], col = factor(meta$plate), pch = 19)

If plate separates the samples, you have three options in increasing order of assumption. First, include batch as a covariate in your statistical model rather than modifying the data: design <- model.matrix(~ state + plate) in limma, which is the statistically cleanest route when batch and biology are not confounded. Second, limma::removeBatchEffect, for visualization only. Third, ComBat from sva, an empirical Bayes method that shrinks per-protein batch means and variances.

library(sva)
mod <- model.matrix(~ 1, data = meta)
corrected <- ComBat(dat = as.matrix(na.omit(norm_median)),
                    batch = meta$plate, mod = mod, par.prior = TRUE)

ComBat needs several samples per batch and will overcorrect when batch is confounded with the variable you care about, which in a longitudinal personal series it frequently is. When you have samples that appear in more than one batch, use them. Internal reference scaling for TMT and the tunable median polish of ratio implemented as TAMPOR work the same way. For each protein they compute the ratio of each sample to a within-batch reference, which is either dedicated bridging samples or the batch’s own trimmed mean. They then iteratively polish rows and columns until convergence. This approach harmonizes datasets without requiring that the biological groups be balanced across batches, and the tunable form lets you use all samples as the reference when dedicated bridges are unavailable.⁸

The operational recommendation: on every plate, include one aliquot of a single frozen reference draw of your own plasma. That one well per plate costs almost nothing and converts an unidentifiable batch effect into a measurable one. Without it, a change you observe between a January and a July draw cannot be separated from the plate they ran on.

7. Build your own reference interval and re-normalize the whole series

Once the matrix is normalized and batch-corrected, the useful quantity for a personal profile is not the absolute concentration but your own variation around your own baseline. Compute the mean and standard deviation across your draws for each protein. Express each new draw as a deviation in your own units.

prot_mu <- rowMeans(corrected, na.rm = TRUE)
prot_sd <- apply(corrected, 1, sd, na.rm = TRUE)
z <- sweep(sweep(corrected, 1, prot_mu, "-"), 1, prot_sd, "/")

cv_intra <- prot_sd / abs(prot_mu)
head(sort(cv_intra))   # your most stable analytes
head(sort(cv_intra, decreasing = TRUE))  # noisiest, interpret with caution

With four to six draws these estimates are crude, and a protein whose intra-individual coefficient of variation exceeds the platform’s technical CV (typically 5 to 15% for Olink and SomaScan, 10 to 20% for label-free DIA) is the only kind where a single change is worth a second look. Anything at or below technical CV is noise.

Two rules follow. First, re-normalize the entire series whenever a new draw arrives rather than normalizing the new sample against the old parameters, because quantile normalization, VSN, loess, and ComBat all estimate parameters from the full sample set and their output changes when the set changes. Keep the raw matrices and rerun the pipeline from a script so this is one command. Second, freeze and version everything.

• tool versions
• the FASTA and spectral library
• the filtering thresholds
• the normalization method
• the batch model

Record them in the same directory as the output. A change in the normalization method between two analyses will produce apparent biological change, and you will not be able to tell unless you logged it.

Finally, on interpretation. Normalized proteomics gives you relative change in one person over time, with reference intervals you built from a handful of your own draws. The platforms are not clinically validated for individual diagnosis. The result is a measurement and does not amount to a finding. Any result that would change what you do about your health belongs in front of a physician, alongside clinically validated assays ordered and interpreted by them.

Common problems

One sample’s RLE box is shifted well off zero after median centering. Median centering only fixes a multiplicative scale error. A shifted box after centering usually means a differential loss of a subset of proteins, from incomplete digestion or a freeze-thaw cycle. Check the sample’s identification count and the distribution of missing values. If both are outliers, exclude it.

PC1 still separates by plate after ComBat. Either batch and the variable of interest are confounded, or the batch effect is not location-and-scale in form. Try the ratio-based approach instead, computing each sample as a log ratio to a within-plate reference and polishing, which does not assume a shared additive batch shift.⁸

Quantile normalization makes your samples look great and your results vanish. That is the expected behavior when a true global shift exists. Compare against median centering, and if the two disagree on direction for many proteins, the global distribution itself is changing and quantile normalization is the wrong tool.

Olink and SomaScan measurements of the same protein do not agree. The platforms are different assays measuring different epitopes, and no normalization reconciles them. Treat them as separate time series. Cross-platform normalization is a separate problem, and integration is usually better handled at the model level than by forcing the measurements onto a common scale.⁹

Log2 of zero produces -Inf. Set zeros to NA before transforming, as in step 2. A pseudocount (log2(x + 1)) is acceptable for count-like data but distorts low-intensity MS values, where the floor is a detection limit and not a true zero.

Different normalizations give you different differential proteins and you do not know which to report. Report the overlap, and state the methods you tried. This sensitivity is a documented property of proteomics data rather than a mistake on your part.²

Footnotes

1. Cheng Chang, Kaikun Xu, Chaoping Guo, et al. PANDA-view: an easy-to-use tool for statistical analysis and visualization of quantitative proteomics data. Bioinformatics, 2018. https://doi.org/10.1093/bioinformatics/bty408 ↩

2. Mythreyi Narasimhan, Sadhana Kannan, Aakash Chawade, et al. Clinical biomarker discovery by SWATH-MS based label-free quantitative proteomics: impact of criteria for identification of differentiators and data normalization method. Journal of Translational Medicine, 2019. https://doi.org/10.1186/s12967-019-1937-9 ↩ ↩²

3. Chi Yen Tseng, Jessica A. Salguero, Joshua D. Breidenbach, et al. Evaluation of normalization strategies for mass spectrometry-based multi-omics datasets. Metabolomics, 2025. https://doi.org/10.1007/s11306-025-02297-1 ↩

4. Aakash Chawade, Erik Alexandersson, Fredrik Levander. Normalyzer: A Tool for Rapid Evaluation of Normalization Methods for Omics Data Sets. Journal of Proteome Research, 2014. https://doi.org/10.1021/pr401264n ↩

5. Aleesa E. Chua, Leah D. Pfeifer, Emily R. Sekera, et al. Workflow for Evaluating Normalization Tools for Omics Data Using Supervised and Unsupervised Machine Learning. Journal of the American Society for Mass Spectrometry, 2023. https://doi.org/10.1021/jasms.3c00295 ↩

6. Benjamin Balluff, Carsten Hopf, Tiffany Porta Siegel, et al. Batch Effects in MALDI Mass Spectrometry Imaging. Journal of the American Society for Mass Spectrometry, 2021. https://doi.org/10.1021/jasms.0c00393 ↩

7. Gianluca Arauz-Garofalo, Sergio Ciordia, Anne Gonzalez de Peredo, et al. A generalizable normalization framework to decouple protocol and instrument effects: Application to high-sensitivity proteomics multicentric study (PME13). 2026. https://doi.org/10.64898/2026.08.16.744113 ↩

8. Eric B. Dammer, Nicholas T. Seyfried, Erik C. B. Johnson. Batch correction and harmonization of –Omics datasets with a tunable median polish of ratio. Frontiers in Systems Biology, 2023. https://doi.org/10.3389/fsysb.2023.1092341 ↩ ↩²

9. Muhtasim Noor Alif, Khandakar Tanvir Ahmed, Sudipto Baul, et al. Graph designs for deep learning–based multi-omics integration. Briefings in Bioinformatics, 2026. https://doi.org/10.1093/bib/bbag410 ↩