Running DEP on Your Own Plasma Proteomics Data

By the end of this guide you will have a SummarizedExperiment holding your own protein abundances across several blood draws, normalized and quality-checked, plus a table of proteins whose abundance differs between the states you defined (say, three fasted mornings against three post-prandial afternoons, or pre- and post-intervention blocks), with log2 fold changes, moderated t-statistics, and Benjamini-Hochberg adjusted p-values. You need R 4.4 or later with Bioconductor 3.20, the packages DEP, SummarizedExperiment, limma, vsn, and dplyr, roughly 8 GB of RAM (protein-level matrices are small, a few thousand rows), and the search-engine output from your proteomics run: a MaxQuant proteinGroups.txt, a DIA-NN report.pg_matrix.tsv, or a Spectronaut protein-level report. You also need to know how many samples you have per state, because that number sets the ceiling on what any of this can tell you.

DEP was written for designed experiments with treatment groups and replicates, and applying it to one person’s longitudinal series means being precise about what plays the role of “condition” and what plays the role of “replicate.” We will get to that in step 6, and it is the step that decides whether your output is a hypothesis list or noise with a p-value column.

1. Locate the quantitative columns in your search output

DEP does not parse raw spectra. It expects a wide matrix of quantitative values with one row per protein group and one column per sample, plus identifier columns. Where those live depends on the upstream pipeline, and the naming matters because DEP selects columns by regular expression.

From MaxQuant, the file is txt/proteinGroups.txt, tab-separated, with identifier columns Protein.IDs, Majority.protein.IDs, Gene.names, and Protein.names, and quantitative columns named LFQ.intensity.<experiment> after R’s read.delim converts spaces to dots. There are also Reverse and Potential.contaminant columns carrying + for decoy hits and contaminant matches, which must be dropped before anything else.

library(DEP); library(dplyr); library(SummarizedExperiment)

pg <- read.delim("txt/proteinGroups.txt", stringsAsFactors = FALSE)
pg <- filter(pg, Reverse != "+", Potential.contaminant != "+")
grep("^LFQ\\.intensity", colnames(pg), value = TRUE)

From DIA-NN, use report.pg_matrix.tsv, which gives Protein.Group, Genes, and one column per raw file path. Rename those columns to short sample labels before import, because the file paths contain dots and slashes that will wreck your regexes. DIA-NN’s MaxLFQ-normalized protein groups are the right input; do not feed DEP the precursor-level report.tsv unless you intend to use DEP2’s peptide workflow (step 8). If your data came out of a cloud pipeline rather than a local search, the same rule applies: you want protein-group-level quantities with a clear record of which normalization the platform already applied, since one-stop processing platforms differ in whether they hand back raw or normalized intensities ¹.

2. Build the experimental design table

This is a three-column data frame, and DEP is strict about it. label must match the quantitative column names exactly once the LFQ.intensity. prefix is stripped. condition is the grouping variable you will contrast. replicate is an integer distinguishing samples within a condition.

expdesign <- data.frame(
  label     = c("D1_fast","D2_fast","D3_fast","D1_fed","D2_fed","D3_fed"),
  condition = c("fasted","fasted","fasted","fed","fed","fed"),
  replicate = c(1,2,3,1,2,3),
  stringsAsFactors = FALSE
)

For a personal profile, condition is whatever state you can define crisply and sample repeatedly: fasting versus fed, on-plan versus off-plan training blocks, before and after a documented change in sleep or altitude. Vague conditions produce uninterpretable contrasts. Keep a fourth column recording the draw date and the processing batch even though DEP ignores it; you will need it in step 4.

3. Assemble the SummarizedExperiment

DEP needs unique row names, and gene symbols are not unique in proteinGroups.txt (isoform groups share symbols, and some rows have empty Gene.names). make_unique resolves this by falling back to protein IDs.

data_unique <- make_unique(pg, names = "Gene.names", ids = "Protein.IDs", delim = ";")
lfq_cols    <- grep("^LFQ\\.intensity", colnames(data_unique))
se          <- make_se(data_unique, lfq_cols, expdesign)

se
assay(se)[1:5, ]

make_se log2-transforms the intensities and converts MaxQuant’s zeros to NA, which is the behavior you want: a zero in an LFQ column means “not quantified in this run,” not “absent from the plasma.” Confirm the transform took by checking that assay(se) values sit in the 15 to 35 range rather than the millions.

If your input is already log-scale (some Spectronaut exports are), build the object with make_se_parse or construct it manually and skip the transform, otherwise you will take the log of a log and flatten every fold change toward zero.

4. Inspect missingness and batch structure before touching the numbers

Proteomics missingness is not random, and this is the single most consequential fact about the data. Low-abundance proteins fall below the detection limit in some runs and not others, so their missing values are left-censored (missing not at random, MNAR), while a peptide lost to a bad injection is missing at random (MAR). Imputation methods for the two cases are different, so look at which you have.

plot_frequency(se)                       # how many proteins in how many samples
filt <- filter_missval(se, thr = 0)      # keep proteins with no NAs in >= 1 condition
plot_numbers(filt); plot_coverage(filt)
plot_missval(filt)                       # heatmap of NA positions
plot_detect(filt)                        # intensity distribution: all vs. NA-containing

plot_detect is the diagnostic that matters. If the density for proteins with missing values sits clearly to the left of the density for complete proteins, your missingness is dominated by the detection limit and an MNAR-appropriate imputation is defensible. If the two densities overlap, something else is going on, usually run-to-run instrument variation, and imputing with a low-value distribution will manufacture fold changes.

plot_missval also reveals batch structure without any statistics: if the NA pattern clusters by processing date rather than by condition, you have a batch effect that will dominate everything downstream. For a longitudinal personal profile this is the normal case, because your samples were collected weeks apart and may have been prepared and injected in different runs. Record the batch, and plan to model it rather than remove it from the assay. Cross-cohort work on proteomic scores shows how poorly signals transfer when technical structure is not handled at the modeling stage ².

filter_missval(se, thr = 0) is aggressive and correct for small designs: it retains only proteins quantified in every sample of at least one condition. With six samples you will typically go from 4,000-6,000 plasma protein groups down to 2,000-3,500. Accept the loss. Testing a protein observed in two of six runs is not a test.

5. Normalize, then decide on imputation deliberately

DEP’s default normalization is variance-stabilizing normalization, which fits a per-intensity variance model rather than assuming constant variance on the log scale.

norm <- normalize_vsn(filt)
plot_normalization(filt, norm)
meanSdPlot(norm)

meanSdPlot should show a roughly flat red trend line. A strong upward slope at high intensity means VSN did not converge, usually because too few proteins remain or because the input was already normalized and log-transformed twice.

For imputation, we use MinProb when plot_detect shows left-censoring, and we impute only after filtering, never before:

imp <- impute(norm, fun = "MinProb", q = 0.01)
plot_imputation(norm, imp)

MinProb draws from a Gaussian centered on a low quantile of each sample’s observed distribution, with q = 0.01 setting that quantile. plot_imputation overlays the pre- and post-imputation densities; you want a modest left shoulder, not a second mode that rivals the main peak. If the imputed values form a large separate hump, you filtered too loosely.

The alternative worth considering is mixed imputation: MinProb for proteins that are missing in an entire condition (plausibly below detection there) and k-nearest-neighbors for scattered single missing values. DEP supports this through impute(se, fun = "mixed", randna = ..., mar = "knn", mnar = "MinProb"), where randna is a logical vector marking MAR proteins. It is more work and it is the more defensible choice when your plot_detect densities overlap.

Our opinion for personal data: prefer filtering over imputing. With six to twelve samples, every imputed value is a meaningful fraction of a group mean, and a protein that is present in three draws and absent in three is more interesting as a presence/absence observation than as an imputed fold change.

6. Define the contrast and the design formula

test_diff runs limma’s moderated t-test, which borrows variance information across proteins to stabilize per-protein estimates. That borrowing is what makes small designs workable, and it is also why you still need real replication within each condition.

dep <- test_diff(imp, type = "manual", test = c("fed_vs_fasted"),
                 design_formula = formula(~ 0 + condition + batch))
dep <- add_rejections(dep, alpha = 0.05, lfc = log2(1.5))

Three points about this call. First, type = "manual" with an explicit test string is clearer than type = "control", and the string must be <condition>_vs_<condition> using the exact factor levels. Second, adding + batch to design_formula requires that batch exists in colData(imp) and that it is not collinear with condition; if every fasted sample came from batch 1 and every fed sample from batch 2, the model is unidentifiable and limma will drop a coefficient. Design your collection to cross condition with batch, even if that means running two draws per session. Third, lfc = log2(1.5) imposes a minimum effect size on top of the FDR threshold, which we prefer to reporting statistically significant 8% changes in a plasma proteome measured across weeks.

If you have only one sample per condition, stop here. There is no variance to estimate and no valid p-value. What you can do instead is rank proteins by log2 ratio, report the ranking as a hypothesis list, and treat the spread of ratios among housekeeping-like abundant proteins as your empirical noise floor.

7. Extract results and read them carefully

res <- get_results(dep)
res %>% filter(significant) %>%
  select(name, ID, fed_vs_fasted_ratio, fed_vs_fasted_p.adj) %>%
  arrange(fed_vs_fasted_p.adj) %>% head(25)

plot_pca(dep, x = 1, y = 2, n = 500, point_size = 4)
plot_cor(dep, significant = TRUE, lower = 0, upper = 1)
plot_volcano(dep, contrast = "fed_vs_fasted", label_size = 2, add_names = TRUE)
plot_heatmap(dep, type = "centered", kmeans = TRUE, k = 4, col_limit = 3)
plot_single(dep, proteins = c("APOA1","CRP"), type = "centered")

Look at plot_pca first. If PC1 separates your batches rather than your conditions, the batch term in the design is carrying most of the signal and your effect estimates are correspondingly fragile. plot_cor on the significant subset should show within-condition correlations above between-condition correlations; if it does not, the “significant” set is probably tracking run order.

For annotation, join your hits to protein-level identifier resources rather than guessing from gene symbols, since symbol-to-protein mapping is many-to-many and isoform-specific evidence matters for whether a plasma measurement is even plausible ³. Then apply the sanity filters that matter for blood: is this protein normally detectable in plasma at all, or is it an intracellular protein whose appearance suggests hemolysis or tissue leakage during the draw? Does its direction of change match the obvious physiological driver of your contrast? A fed-versus-fasted contrast that moves apolipoproteins is behaving as expected; one that moves only immunoglobulins is telling you about sample handling.

Nothing in this output is a diagnosis. A differential protein abundance in your own plasma is a measurement about your samples under your conditions, and mapping it to clinical meaning requires a physician who can see your history, your exam, and validated clinical assays. Take the table to them rather than acting on it.

8. Extend to peptides and other omics layers with DEP2

DEP2 is the successor package and is worth adopting when you have peptide-level data or more than one molecular layer. Two capabilities justify the switch. First, it separates peptide import from protein rollup, so you can aggregate with aggregate_pe() using robustSummary or median polish from MsCoreUtils instead of accepting the search engine’s protein inference. This matters when a protein group’s quantity is driven by one shared peptide. Second, it carries the same contrast syntax into transcript-level analysis, which lets you build matched protein and RNA contrasts from a single design.

# sketch; see the DEP2 vignette for current argument names
library(DEP2)
pe  <- import_MaxQuant(peptides, expdesign, quantity_col = "Intensity")
pe  <- aggregate_pe(pe, fun = "robustSummary", reserve = "Gene.names")

When you put protein and RNA results side by side, expect disagreement rather than confirmation. Transcript and protein abundance decouple through translation rate, secretion, and clearance, and studies that model regulatory layers jointly routinely find discordance between them rather than tight coupling ⁴. Treat concordance as extra evidence and discordance as information about which layer is regulated, not as a failure of one assay. If you plan to build integrated multi-layer classifiers, be aware that machine-learning integration on small cohorts overfits readily and that reported biomarker panels are usually cohort-specific until replicated ⁵. Proteomics is also sensitive enough to register exposure-driven changes across several protein layers at once, which is an argument for logging environment and diet alongside every draw ⁶.

Common problems

make_se throws an error about non-unique row names. You skipped make_unique, or you passed names = "Gene.names" on a file where that column is missing (DIA-NN calls it Genes).

Everything looks significant, hundreds of proteins at FDR below 0.05 with large ratios. Almost always a batch effect aligned with condition, or imputation doing the work. Check plot_pca colored by batch, then rerun with lfc = log2(1.5) and without imputation on a complete-cases-only matrix. If the hit list collapses, believe the collapsed version.

Nothing is significant. With three samples per condition and plasma-level biological variability this is a common and correct answer. Before concluding the effect is absent, check that VSN converged, that filtering left enough proteins, and that your conditions were separated by more than a few hours.

test_diff complains about coefficients or returns NA estimates. Collinearity between condition and batch, or a condition with a single sample. Fix the design, not the code.

Fold changes shrink toward zero across the whole matrix. Double log transform: you fed already-logged values into make_se. Check range(assay(se), na.rm = TRUE).

Gene symbols in your hit list that you cannot place. Map them through protein-level records with isoform evidence before interpreting ³. Many will be albumin fragments, keratins (skin contamination), or serum amyloid components that respond to any inflammatory stimulus.

Hits that do not reproduce in the next draw. This is the expected outcome for most single-contrast findings and is the reason to collect more timepoints rather than more statistics on the same six. Repeated quantification of the same system is what separates a reproducible protein-level difference from run variation ⁷.

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. Jinwen Feng, Chen Ding, Naiqi Qiu, et al. Firmiana: towards a one-stop proteomic cloud platform for data processing and analysis. Nature Biotechnology, 2017. https://doi.org/10.1038/nbt.3825 ↩

2. Iain R. Konigsberg, Luciana B. Vargas, Katherine A. Pratte, et al. Omic risk scores are associated with COPD-related traits across three cohorts. Respiratory Research, 2026. https://doi.org/10.1186/s12931-026-03890-1 ↩

3. Simon Fishilevich, Shahar Zimmerman, Asher Kohn, et al. Genic insights from integrated human proteomics in GeneCards. Database, 2016. https://doi.org/10.1093/database/baw030 ↩ ↩²

4. Alper Bülbül, Özdeyiş Hülya Yılmaz-İşgördü, Meziyet Dilara Reda, et al. Integrative Multi-Omics Analysis of Multiple Sclerosis Reveals Cell-Type-Specific Regulatory Landscapes and Discordant Methylation–Expression Coupling. International Journal of Molecular Sciences, 2026. https://doi.org/10.3390/ijms27167275 ↩

5. Rency S. Varghese, Xinran Zhang, Muhammad S. Sajid, et al. Machine Learning-Based Multi-Omics Integration for Identification of Hepatocellular Carcinoma Biomarkers in an Egyptian Cohort. Journal of Proteome Research, 2025. https://doi.org/10.1021/acs.jproteome.5c00741 ↩

6. Alix Sarah Aldehoff, Isabel Karkossa, Helen Broghammer, et al. Advanced Proteomics Approaches Hold Potential for the Risk Assessment of Metabolism-Disrupting Chemicals as Omics-Based NAM: A Case Study Using the Phthalate Substitute DINCH. Environmental Science & Technology, 2025. https://doi.org/10.1021/acs.est.5c01206 ↩

7. Rui Sun, Weigang Ge, Yi Zhu, et al. Proteomic Dynamics of Breast Cancer Cell Lines Identifies Potential Therapeutic Protein Targets. Molecular & Cellular Proteomics, 2023. https://doi.org/10.1016/j.mcpro.2023.100602 ↩