How to Analyze Your Own Plasma Proteomics Data in R

This guide walks through a complete analysis of plasma proteomics data for a single person measured repeatedly over time, from the vendor’s raw deliverable to a fitted model of how your proteins move. By the end you will have your data in a SummarizedExperiment (or QFeatures) object with sample metadata attached. You will also have a QC report that tells you which proteins and which draws you can trust, a normalized log2-scale matrix, and per-protein estimates of your own within-person variation. Finally, you will have a fitted model of protein trajectories over time with false discovery rate control.

The software requirements are modest. You need R 4.4 or newer and Bioconductor 3.20 or newer, plus the vendor deliverable itself. That deliverable will be an Olink NPX file (CSV or parquet), a SomaScan .adat, or a DIA-NN / Spectronaut / MaxQuant report from a mass spectrometry lab. You also need a sample manifest recording draw date, time of day, and fasting status. It should also record the plate or batch identifier and anything else that plausibly moves protein levels. The manifest matters more than the software. Almost every analysis failure we see in single-person proteomics traces back to unrecorded pre-analytical variation rather than to a modeling mistake.

1. Install the packages you will use

R’s proteomics tooling is mature and lives mostly in Bioconductor, where the design goal has been to keep quantitative proteomics data in the same container classes used for microarray and RNA-seq so that the downstream statistical machinery transfers ¹. Install the set below once and pin the versions with renv so that a rerun a year from now produces the same numbers.

install.packages(c("renv", "BiocManager", "data.table", "tidyverse", "lme4", "broom.mixed"))
BiocManager::install(c(
  "QFeatures", "SummarizedExperiment", "limma", "vsn",
  "msqrob2", "MsCoreUtils", "imputeLCMD",
  "clusterProfiler", "ReactomePA", "org.Hs.eg.db", "fgsea", "mixOmics"
))
# vendor readers
install.packages("OlinkAnalyze")        # Olink NPX
# SomaDataIO from Bioconductor or GitHub, depending on version
BiocManager::install("SomaDataIO")

MSnbase still works, but QFeatures is the object model we would build on now. It links precursor, peptide, and protein assays in one object and keeps the aggregation provenance. For affinity platforms such as Olink and SomaScan you only ever have one assay level, so a plain SummarizedExperiment is enough.

2. Read the vendor file and build a SummarizedExperiment

Each platform hands you a differently shaped file, so the first job is to get everything into a common form: a proteins-by-samples matrix on a log2 scale, accompanied by a colData table of sample annotations. Everything downstream assumes that structure.

For Olink, the NPX file arrives in long format with one row per sample-assay pair. NPX is already a log2-scale relative quantity, so do not log it again.

library(OlinkAnalyze); library(data.table); library(SummarizedExperiment)

npx <- read_NPX("Q-12345_NPX_2026-03-02.parquet") |> as.data.table()

# NPX columns: SampleID, OlinkID, UniProt, Assay, Panel, NPX, LOD,
#              QC_Warning, Assay_Warning, PlateID, Normalization
mat <- dcast(npx, OlinkID + Assay + UniProt ~ SampleID, value.var = "NPX")
rd  <- mat[, .(OlinkID, Assay, UniProt)]
m   <- as.matrix(mat[, -(1:3)]); rownames(m) <- rd$OlinkID

cd <- unique(npx[, .(SampleID, PlateID, QC_Warning)])
manifest <- fread("draws.csv")   # SampleID, draw_date, fasting, time_of_day, ...
cd <- merge(cd, manifest, by = "SampleID")
setkey(cd, SampleID); cd <- cd[colnames(m)]

se <- SummarizedExperiment(
  assays  = list(npx = m),
  rowData = rd,
  colData = DataFrame(cd, row.names = cd$SampleID)
)

For SomaScan, SomaDataIO::read_adat() returns a wide tibble in which the RFU (relative fluorescence unit) columns are named by SeqId and the sample metadata sits in the leading columns. Take log2() of the RFU values and keep the RowCheck and ColCheck flags. Record which normalization the lab applied, since ANML against a reference population is the common default and changes what your values mean.

For mass spectrometry, work from the protein-group matrix rather than the long precursor report, unless you intend to redo the quantification yourself. DIA-NN’s report.pg_matrix.tsv is already MaxLFQ-summarized protein groups by run. If you start instead from report.tsv, filter before summarizing. Require Q.Value < 0.01, PG.Q.Value < 0.01, and Lib.Q.Value < 0.01. Then drop precursors with fewer than two identifications across your runs. Log2 the intensities once you have them.

As a sanity check on yield, expect 400 to 800 protein groups from neat plasma without depletion, and a few thousand with depletion or deep fractionation. Albumin and immunoglobulins will dominate the intensity range and should be treated as a normalization hazard rather than as signal.

3. Run QC before anything else

Quality control in single-person longitudinal proteomics carries a different emphasis than in a case-control study. Outlier subjects are not the concern here. What you are looking for is draws and plates that will masquerade as biology.

Start with the vendor’s own flags, then examine per-sample distributions and the overall structure of the matrix.

library(limma)

# 1. vendor QC
table(colData(se)$QC_Warning)

# 2. per-sample central tendency and spread
x <- assay(se, "npx")
qc <- data.frame(
  sample = colnames(x),
  median = apply(x, 2, median, na.rm = TRUE),
  iqr    = apply(x, 2, IQR, na.rm = TRUE),
  n_na   = colSums(is.na(x))
)

# 3. structure: does PC1 track plate or draw date?
pc <- prcomp(t(x[complete.cases(x), ]), scale. = FALSE)
summary(pc)$importance[, 1:4]
plot(pc$x[, 1], pc$x[, 2], col = factor(colData(se)$PlateID), pch = 19)

If the first principal component separates plates, you have a batch effect that will be confounded with time whenever draws were batched by collection date, which they usually are. This is the single most consequential design decision in a personal profile. Ask the lab to randomize your samples across plates, or bank samples and run them together. If neither is possible, at minimum keep a bridging sample on every plate so that plate offsets remain estimable. A bridging sample is an aliquot of one pooled draw.

For mass spectrometry data, add two further checks. Track the number of quantified protein groups per run, where a drop of more than roughly 20 percent below the median usually means a bad injection or a column problem, and track the median coefficient of variation across technical replicates if you have them.

4. Handle missingness deliberately, and know which kind you have

Missing values need different treatment depending on how they arose, and affinity and mass spectrometry platforms fail in different ways. Olink reports a limit of detection (LOD) per assay, so low-abundance proteins produce values that are reported but uninformative. Mass spectrometry instead produces explicit NA values that are mostly left-censored, meaning the protein was present below the detection threshold rather than missing at random.

Our default rule follows from that distinction: filter first, impute as little as possible, and never impute a protein you intend to interpret individually.

# Olink: keep assays above LOD in at least 75% of draws
lod <- dcast(npx, OlinkID ~ SampleID, value.var = "LOD")
above <- assay(se, "npx") > as.matrix(lod[match(rownames(se), lod$OlinkID), -1])
keep  <- rowMeans(above, na.rm = TRUE) >= 0.75
se    <- se[keep, ]

# MS: require quantification in 70% of runs, then impute left-censored
library(MsCoreUtils)
keep <- rowMeans(!is.na(assay(se))) >= 0.70
se   <- se[keep, ]
assay(se, "imp") <- impute_matrix(assay(se), method = "MinDet", q = 0.01)

MinDet and QRILC from imputeLCMD assume left-censoring and are the right family for mass spectrometry missingness. k-nearest-neighbor imputation, by contrast, assumes values are missing at random and will pull censored values upward. That inflates apparent stability in exactly the low-abundance proteins you were curious about. If a protein is missing in a structured way, present in every summer draw and absent in every winter draw, stop and check whether the runs were batched by season before you interpret anything.

5. Normalize, then decide what to do about batch

Normalization removes sample-level technical differences in scale, and the appropriate method depends on the platform. Log2 plus median centering is adequate for Olink, whose intra-plate normalization is already applied upstream. For mass spectrometry we prefer variance-stabilizing normalization, which handles the intensity-dependent variance of MS data better than quantile normalization does at the low end.

# affinity platforms: median-center columns
assay(se, "norm") <- limma::normalizeBetweenArrays(assay(se, "npx"), method = "scale")

# mass spec: VSN on the linear-scale intensity matrix
library(vsn)
fit  <- vsn2(2^assay(se, "imp"))
assay(se, "norm") <- exprs(fit)   # returns glog2-scale values

Batch is a separate problem, and the right place to handle plate or batch is as a model term rather than by overwriting your data. limma::removeBatchEffect() is appropriate for making a PCA plot interpretable and inappropriate as an input to a hypothesis test, because the downstream model then underestimates residual variance. Fit ~ batch + time instead. If your design confounds batch with time completely, no method fixes it, and the correct conclusion is that plate-level differences are not separable from temporal change.

6. Quantify your own within-person variation

This is the step most guides skip, and it is the one that makes a personal proteome usable. Before you can call a protein “changed,” you need to know how much it moves in you when nothing is happening. Compute the within-person standard deviation on the log2 scale from your baseline draws, then convert it to an approximate coefficient of variation.

library(tidyverse)

baseline <- se[, colData(se)$phase == "baseline"]
wp <- data.frame(
  protein = rownames(baseline),
  mean_l2 = rowMeans(assay(baseline, "norm")),
  sd_l2   = apply(assay(baseline, "norm"), 1, sd)
) |>
  mutate(cv_pct = 100 * (2^sd_l2 - 1),           # approx CV on linear scale
         rcv_pct = 100 * (2^(2.77 * sd_l2) - 1)) # 95% reference change value

The reference change value is the fold change you need to see between two draws before it exceeds analytical plus biological noise at roughly 95 percent confidence. For a protein with a within-person log2 SD of 0.15, that works out to about a 34 percent change on the linear scale. Proteins with a within-person CV above 30 or 40 percent are not useful for single-timepoint interpretation no matter how precise the assay is. Publish this table to yourself and consult it before you get excited about any individual number.

7. Model trajectories over time

With four or more draws you can fit time as a covariate and ask whether a protein is trending. For a small number of proteins measured at many timepoints, a mixed model per protein is the transparent choice. For the full matrix, limma with empirical Bayes moderation of the variance estimates is faster and better calibrated when timepoints are few. The limma machinery carries over directly to quantitative proteomics matrices, which is much of the reason R became the default environment for this work ¹. Earlier R-based platforms such as DanteR packaged the same sequence of normalization, rollup, imputation, and linear modeling for label-free quantitative data ².

library(limma); library(splines)

cd  <- as.data.frame(colData(se))
cd$t <- as.numeric(cd$draw_date - min(cd$draw_date)) / 365.25

# linear trend, adjusting for plate and fasting status
design <- model.matrix(~ t + PlateID + fasting, data = cd)
fit <- lmFit(assay(se, "norm"), design) |> eBayes(trend = TRUE, robust = TRUE)
res <- topTable(fit, coef = "t", number = Inf, adjust.method = "BH")

# nonlinear trend, if you have >= 8 draws
design2 <- model.matrix(~ ns(t, df = 3) + PlateID, data = cd)
fit2 <- lmFit(assay(se, "norm"), design2) |> eBayes(trend = TRUE)
res2 <- topTable(fit2, coef = grep("^ns", colnames(design2)), number = Inf)

Two points about interpreting the output are worth keeping in mind. The t coefficient is expressed in log2 units per year, so multiply by log(2) for an approximate proportional rate. And a BH-adjusted p-value below 0.05 across 3,000 assays in a single person means the trend is unlikely to be noise given your measured variance. It says nothing about cause. Intercurrent illness, a change in training load, an altered sleep schedule, and a lab switching reagent lots all produce clean linear trends.

If you have technical replicates or bridging samples, use duplicateCorrelation rather than averaging them. This recovers the within-replicate correlation and carries it into the variance estimate:

corfit <- duplicateCorrelation(assay(se, "norm"), design, block = cd$aliquot_id)
fit <- lmFit(assay(se, "norm"), design, block = cd$aliquot_id,
             correlation = corfit$consensus) |> eBayes()

8. Interpret sets, not single proteins

Movement in a single protein in a single person is fragile evidence. Coordinated movement across a pathway is more informative, because the correlated measurement error across assays is usually smaller than the shared biological signal. When you test for enrichment, use the proteins you measured as the background rather than the whole genome. Otherwise every plasma panel will look enriched for secreted and complement proteins.

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

universe <- bitr(rowData(se)$UniProt, "UNIPROT", "ENTREZID", org.Hs.eg.db)$ENTREZID
hits <- res |> filter(adj.P.Val < 0.05) |> rownames()
hit_entrez <- bitr(rowData(se)[hits, "UniProt"], "UNIPROT", "ENTREZID", org.Hs.eg.db)$ENTREZID

enrichPathway(gene = hit_entrez, universe = universe, pvalueCutoff = 0.05,
              readable = TRUE) |> as.data.frame() |> head(20)

For ranked analysis, running fgsea on the moderated t-statistic from topTable is preferable to a cutoff-based test, since it uses the whole ranking and does not depend on an arbitrary significance threshold.

9. Join the proteome to your other assays

Proteins become far more interpretable next to transcripts, metabolites, and dense physiological data. This is where the R ecosystem’s multi-omics packages earn their place. mixOmics descends from integrOmics, which introduced regularized canonical correlation analysis and sparse partial least squares for finding correlated structure between two omics matrices measured on the same samples ³. Multiple co-inertia analysis extends the same logic to three or more blocks and is useful for identifying which assays carry shared versus assay-specific variation ⁴. Purpose-built pipelines for paired transcriptomic and proteomic tables exist as well, including web-facing R platforms that handle the annotation joins and visualization for you ⁵.

library(mixOmics)
X <- list(prot = t(assay(se, "norm")[vf, ]),
          rna  = t(vst_counts[vg, ]))
Y <- cd$phase
res_block <- block.splsda(X, Y, ncomp = 2,
                          keepX = list(prot = c(25, 25), rna = c(50, 50)))
plotIndiv(res_block); plotVar(res_block, cutoff = 0.5)

Be realistic about sample size here. With ten to twenty draws you can describe covariation and generate hypotheses, and you cannot validate them. Large-scale work that integrates plasma proteomics with transcriptomic evidence to nominate candidate targets depends on thousands of individuals and genetic instruments to separate correlation from cause ⁶. Your data is a time series of one person, which answers different questions. What is your own baseline, how much does it move, and did anything shift after a specific intervention.

One boundary deserves stating plainly. Nothing in this pipeline is a clinical result. If a protein with an established clinical assay such as troponin, CRP, PSA, NT-proBNP, or ferritin sits far outside its expected range in your data, the right next step is a clinically validated test ordered through a physician rather than further modeling. Research-grade platforms are not calibrated to clinical decision thresholds and can differ from a clinical assay by a large factor even when they correlate well.

Common problems

The failures below account for most of the trouble we see in practice, and each has a straightforward remedy.

• Plate effects confounded with time. If draws were shipped and run in collection order, the plate term absorbs the time term and both become uninterpretable. Randomize across plates, or bank and batch.
• Double-logging Olink data. NPX is already log2. Taking log2() again compresses everything and makes all your fold changes look tiny.
• Imputing then testing, with no record of which values were imputed. Keep a logical matrix of imputed positions alongside the assay and exclude heavily imputed proteins from any per-protein claim.
• Using the whole proteome as the enrichment universe. A targeted panel of 3,000 plasma proteins is enriched for secreted, immune, and inflammatory proteins by construction. Use the measured set as the universe.
• Quantile normalization on plasma mass spectrometry data with variable albumin depletion efficiency. Quantile normalization forces identical distributions and will hide a genuine shift in high-abundance carrier proteins. VSN plus a check on the depletion controls is safer.
• Comparing across platform versions. Olink panel revisions and SomaScan menu changes alter which reagents measure a given protein, and values are not interchangeable across versions. Bridge with a shared sample set, or treat the series as two separate series.
• Reading significance into a single draw. This is the most frequent error. Compute your within-person reference change value first, and treat every single-timepoint number as a measurement with a confidence interval you have estimated yourself.

Woolf Software 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. Laurent Gatto, Andy Christoforou. Using R and Bioconductor for proteomics data analysis. Biochimica et Biophysica Acta (BBA) - Proteins and Proteomics, 2014. https://doi.org/10.1016/j.bbapap.2013.04.032 ↩ ↩²

2. Tom Taverner, Yuliya V. Karpievitch, Ashoka D. Polpitiya, et al. DanteR: an extensible R-based tool for quantitative analysis of -omics data. Bioinformatics, 2012. https://doi.org/10.1093/bioinformatics/bts449 ↩

3. Kim-Anh Lê Cao, Ignacio González, Sébastien Déjean. integrOmics: an R package to unravel relationships between two omics datasets. Bioinformatics, 2009. https://doi.org/10.1093/bioinformatics/btp515 ↩

4. Chen Meng, Bernhard Kuster, Aedín C Culhane, et al. A multivariate approach to the integration of multi-omics datasets. BMC Bioinformatics, 2014. https://doi.org/10.1186/1471-2105-15-162 ↩

5. Punit Tyagi, Mangesh Bhide. Development of a bioinformatics platform for analysis of quantitative transcriptomics and proteomics data: the OMnalysis. PeerJ, 2021. https://doi.org/10.7717/peerj.12415 ↩

6. Shucheng Si, Hongyan Liu, Lu Xu, et al. Identification of novel therapeutic targets for chronic kidney disease and kidney function by integrating multi-omics proteome with transcriptome. Genome Medicine, 2024. https://doi.org/10.1186/s13073-024-01356-x ↩