How to Analyze Mass Spectrometry Proteomics Data

You need the raw files themselves and a protein sequence database in FASTA format. Plan on roughly 32–64 GB of RAM and a few hundred gigabytes of disk, since raw plasma files run 0.5–2 GB each and mzML conversion inflates them two- to four-fold. You also need either Linux with Docker or a Windows machine if you want to run vendor libraries natively. Everything below is research-grade measurement. Plasma protein values from a discovery pipeline are not clinical results, and any value you want to act on needs to be re-measured by a clinical laboratory and interpreted by a physician.

1. Establish what kind of data you have

Before touching a tool, answer four questions: what instrument, what acquisition mode, what labeling scheme, and what sample preparation. The acquisition mode matters most. In data-dependent acquisition (DDA), the instrument takes a survey scan (MS1, a mass spectrum of intact peptide ions). It then picks the most intense precursors one at a time for fragmentation, producing MS/MS or MS2 spectra whose fragment masses encode the peptide sequence. In data-independent acquisition (DIA), the instrument fragments everything within a sliding isolation window, so every peptide is fragmented in every run. The resulting chimeric spectra are deconvolved computationally against a library.

This is also where the common terminology confusion resolves. “MS” is a single stage of mass analysis and gives you intact masses. “MS/MS” (tandem MS) adds a fragmentation step and a second analysis, which is what makes sequence identification possible. “LC-MS” describes the front end: peptides are separated in time by reversed-phase liquid chromatography before entering the source, so each ion has a retention time as well as a mass-to-charge ratio. A modern proteomics run is all three at once, LC-MS/MS, and the raw file is a three-dimensional object of retention time, m/z, and intensity, with MS2 scans hanging off it.

Labeling is the third axis. Label-free quantification (LFQ) compares intensities across separate runs. Isobaric labeling (TMT, tandem mass tags) multiplexes 10–18 samples into one run and reads quantities from low-mass reporter ions in the MS2 spectrum. The pipelines diverge completely after the search step. A tool that handles one well may handle the other poorly, which is a long-standing finding in comparisons of quantitative proteomics software.¹ For personal longitudinal plasma profiling we use DIA with label-free quantification, because it gives the most consistent matrix across samples collected months apart with the fewest missing values.

2. Convert vendor raw files to mzML

Vendor formats such as Thermo .raw, Bruker .d directories and SCIEX .wiff/.wiff2 are readable by open tools only through vendor DLLs. The first step is therefore conversion to mzML, the community XML format. ProteoWizard’s msconvert is the reference implementation and is available as a Docker image that bundles the vendor libraries under Wine.

docker run --rm -e WINEDEBUG=-all -v "$PWD":/data \
  chambm/pwiz-skyline-i-agree-to-the-vendor-licenses \
  wine msconvert /data/raw/*.raw \
  --mzML --zlib --64 \
  --filter "peakPicking vendor msLevel=1-2" \
  --filter "zeroSamples removeExtra" \
  -o /data/mzML

The peakPicking vendor filter is the one that matters. It applies the manufacturer’s centroiding algorithm to profile-mode data, which is consistently better than ProteoWizard’s own peak picker and is what downstream search engines expect. For Bruker timsTOF data, skip this step entirely: DIA-NN reads .d folders directly and conversion loses ion mobility information. For Thermo files on Linux, ThermoRawFileParser -i run.raw -o mzML/ -f 2 is a lighter alternative that needs no Wine.

Inspect the result before spending hours on a search. A total ion chromatogram that flatlines, spikes, or shows a retention time shift of more than a minute relative to sibling runs is a chromatography problem that no software will fix.

from pyteomics import mzml
import numpy as np

rt, tic, n_ms2 = [], [], 0
with mzml.read("mzML/run1.mzML") as reader:
    for s in reader:
        if s["ms level"] == 1:
            rt.append(s["scanList"]["scan"][0]["scan start time"])
            tic.append(s["total ion current"])
        else:
            n_ms2 += 1
print(len(rt), "MS1 scans,", n_ms2, "MS2 scans, median TIC", np.median(tic))

3. Assemble the sequence database

Identification is a matching problem against a defined search space, so the FASTA file determines what can possibly be found. We use the UniProt human reviewed proteome (Swiss-Prot, roughly 20,400 entries) plus isoforms. We concatenate it with a contaminant set such as cRAP or the MaxQuant contaminants file, which catches keratins, trypsin autolysis peptides, and bovine serum albumin carryover.

cat UP000005640_9606.fasta UP000005640_9606_additional.fasta crap.fasta \
  > human_sp_iso_crap.fasta
grep -c '^>' human_sp_iso_crap.fasta

Two decisions have real consequences. Including the full TrEMBL set (over 180,000 human entries) enlarges the search space several-fold, which inflates the decoy distribution and costs you identifications at fixed FDR without adding much real signal. Including isoforms is reasonable for a personal profile, but understand that shotgun proteomics is peptide-centric and most isoforms differ by peptides you will never observe, so isoform-level calls should be treated as hypotheses. Search engines then generate decoy sequences (usually reversed or pseudo-reversed) from this same file to estimate the error rate.²

4. Run the search

For DIA data, we use DIA-NN in library-free mode, where a deep learning model predicts spectra and retention times for every tryptic peptide in the FASTA and builds an in-silico library that is then refined against your own runs. This removes the need for a separately acquired DDA library and, in our experience, is the most reliable default for plasma.

diann \
  --f /data/mzML/run1.mzML --f /data/mzML/run2.mzML \
  --fasta human_sp_iso_crap.fasta --fasta-search \
  --gen-spec-lib --predictor \
  --cut K*,R* --missed-cleavages 1 \
  --min-pep-len 7 --max-pep-len 30 \
  --min-pr-charge 2 --max-pr-charge 4 \
  --min-fr-mz 200 --max-fr-mz 1800 \
  --unimod4 --var-mods 1 --var-mod UniMod:35,15.994915,M \
  --qvalue 0.01 --matrices --reanalyse --smart-profiling \
  --threads 32 --out /data/out/report.tsv

The --unimod4 flag fixes carbamidomethylation on cysteine, which you should include if the sample prep used iodoacetamide. The --var-mod UniMod:35 flag allows oxidized methionine, the one variable modification worth paying for by default. The --reanalyse flag enables match-between-runs style cross-run reanalysis, which raises the number of quantified precursors substantially and reduces missing values across a longitudinal series. The --smart-profiling flag adapts library refinement to the number of runs.

For DDA data, use FragPipe. It wraps the MSFragger search engine, Philosopher for FDR filtering, and IonQuant for label-free quantification. Run it headless so the analysis is scriptable:

fragpipe --headless \
  --workflow LFQ-MBR.workflow \
  --manifest files.fp-manifest \
  --workdir /data/out_dda \
  --threads 32 --ram 48

MSFragger’s fragment-index search makes open searches (precursor tolerance of ±500 Da, which surfaces unexpected modifications) practical on a workstation, something that was not feasible with earlier engines. The broader ecosystem includes Comet, MS-GF+, MaxQuant, OpenMS, and Skyline for targeted extraction, and each has documented strengths on particular data types.³⁴ If you prefer a graphical, reproducible environment over shell scripts, the Galaxy framework packages many of these tools with provenance tracking.⁵

5. Understand the FDR numbers before trusting the table

Every identification in your output carries a q-value, and you need to know what it controls. Search engines score each spectrum against both real (target) and reversed (decoy) sequences. Because decoy matches are by construction wrong, their score distribution estimates the distribution of false matches among targets, and the q-value is the minimum FDR threshold at which a given identification would be accepted. The standard is 1% at both the precursor and the protein group level, and these are separate thresholds: a set of peptides each at 1% FDR does not give you proteins at 1% FDR. The target-decoy framework and its statistical assumptions are the foundation of how tandem MS results are validated.²⁶

Two further points shape interpretation. First, “protein group” means a set of proteins that share the observed peptides and cannot be distinguished by the data, so a group labeled with one accession may include homologs. Second, FDR is a property of a list. It does not describe any single entry. A protein identified by one peptide in one run at q = 0.009 is much less trustworthy than one identified by eight peptides in every run, even though both pass the same filter. When you filter your matrix, we recommend requiring at least two distinct precursors per protein group and presence in at least 70% of samples.

6. Build the quantification matrix

DIA-NN writes report.pg_matrix.tsv (protein groups by run) and report.pr_matrix.tsv (precursors by run), already filtered at 1% FDR. FragPipe writes combined_protein.tsv and an MSstats.csv for downstream modeling. Load the matrix, log-transform, normalize, and filter in that order.

library(tidyverse); library(limma)

pg <- read_tsv("out/report.pg_matrix.tsv")
mat <- pg %>% select(-(1:5)) %>% as.matrix()
rownames(mat) <- pg$Protein.Group

mat[mat == 0] <- NA
mat <- log2(mat)

keep <- rowMeans(!is.na(mat)) >= 0.7
mat  <- mat[keep, ]

# median-center each run, then quantile-normalize
mat <- sweep(mat, 2, apply(mat, 2, median, na.rm = TRUE), "-")
mat <- normalizeBetweenArrays(mat, method = "quantile")

Log transformation is not optional. Intensities span several orders of magnitude and are roughly log-normal, so variance is proportional to the mean on the raw scale and every downstream test will be dominated by abundant proteins. Median centering corrects for differences in loaded amount between injections. Quantile normalization is more aggressive and assumes the bulk distribution is comparable across samples, which holds for replicate plasma draws from one person and can fail if a run has badly degraded chromatography.

Missing values need a decision rather than a default. Values missing because the protein is below the detection limit in that run are informatively missing and should not be imputed with the row mean, which pulls them toward the middle of the distribution. We prefer to model missingness explicitly (MSstats and msqrob2 both do this) or to impute low-end values from a downshifted normal distribution, and to record which entries were imputed so you can check whether any result depends on them. Different software packages make different assumptions here, and the choices measurably change the list of proteins you call as changing.¹

7. Model the data and set expectations for plasma

For a longitudinal personal profile, the question is usually whether a protein’s trajectory over time differs from its own within-person noise, which is a mixed-model question rather than a two-group comparison. Start by quantifying technical variation from replicate injections of the same sample: the median coefficient of variation for a well-run DIA plasma experiment should be in the range of roughly 5–15%, and a protein whose biological swing is smaller than that is not something you can interpret from this assay.

design <- model.matrix(~ timepoint + batch, data = meta)
fit <- lmFit(mat, design) %>% eBayes(trend = TRUE, robust = TRUE)
topTable(fit, coef = "timepoint", number = 50, sort.by = "P")

Plasma imposes a hard ceiling on depth. Its protein concentrations span about ten orders of magnitude, and albumin alone accounts for roughly half the protein mass, with the twenty or so most abundant proteins accounting for the large majority of the rest. A single-shot DIA run on neat plasma therefore quantifies a few hundred protein groups, concentrated at the high-abundance end. Depletion columns, fractionation, or nanoparticle enrichment push that into the low thousands at the cost of added handling variance. Know which regime your data came from before comparing your protein count to a published number.

8. Reading a spectrum by hand

Automated pipelines answer the question “which peptides are here,” but there are times you want to read a single spectrum, usually when validating a surprising identification. Open the annotated MS2 in a viewer such as Skyline, PDV, or the DIA-NN viewer. Check that consecutive b-ion or y-ion peaks differ by residue masses: 71.037 for alanine, 113.084 for leucine or isoleucine, and 128.059 for glutamine versus 128.095 for lysine. A confident identification has a run of four or more consecutive ions from one series covering the middle of the sequence. Mass accuracy is the other check: on an Orbitrap, fragment errors should sit within a few parts per million, and a match that is consistently 20 ppm off is suspect regardless of score.

Small-molecule spectra are read differently, and you will meet them if your profile includes metabolomics. The classical manual heuristics still apply. The rule of 13 generates a candidate hydrocarbon formula from a nominal mass: divide the mass by 13 to get quotient n and remainder r, and the base formula is C(n)H(n+r). For nominal mass 106, 106/13 gives n = 8 with r = 2, so C8H10. To test formulas containing oxygen, subtract 16 from the mass and remove one CH4 from the formula for each oxygen, giving C7H6O for the same mass. A base peak at m/z 57 in an electron-ionization spectrum is the classic butyl or tert-butyl cation, C4H9+, though the propanoyl acylium ion C3H5O+ is also nominally 57. High-resolution data separates them cleanly at 57.0699 versus 57.0335. Metabolomics pipelines diverge from proteomics here because identification rests on accurate mass, isotope pattern, and spectral library matching rather than sequence databases, with their own tool ecosystems for peak picking, alignment, and annotation.⁷⁸⁹

Common problems

Almost no peptides identified. Check the FASTA first: searching human data against a mouse proteome, or forgetting to specify carbamidomethylation when the sample was alkylated, will each drop identifications by an order of magnitude. Then check the precursor mass tolerance and whether the instrument was properly calibrated.

Keratin dominates the identification list. This is skin and dust contamination introduced during sample handling. It cannot be removed computationally beyond excluding the contaminant entries, and it suppresses real signal by consuming instrument time.

High missing-value rate across a longitudinal series. Enable cross-run reanalysis (--reanalyse in DIA-NN, match-between-runs in FragPipe) and verify retention times are aligned. If missingness correlates with acquisition date, you have a batch effect and should include batch as a covariate rather than correcting it away before modeling.

Everything looks significant. Usually a normalization failure: one run with systematically higher intensity will make half the proteome appear to move. Plot the per-run intensity distributions as boxplots before and after normalization, and confirm they superimpose.

A protein of apparent clinical relevance changes sharply. Discovery proteomics measures relative intensity, not concentration, and is not a validated clinical assay. Confirm with an orthogonal, clinically validated measurement ordered and interpreted by a physician before drawing any conclusion about health.

Footnotes

1. Lukas N. Mueller, Mi-Youn Brusniak, D. R. Mani, et al. An Assessment of Software Solutions for the Analysis of Mass Spectrometry Based Quantitative Proteomics Data. Journal of Proteome Research, 2008. https://doi.org/10.1021/pr700758r ↩ ↩²

2. Eric W. Deutsch, Henry Lam, Ruedi Aebersold. Data analysis and bioinformatics tools for tandem mass spectrometry in proteomics. Physiological Genomics, 2008. https://doi.org/10.1152/physiolgenomics.00298.2007 ↩ ↩²

3. Steve W. Haga, Hui‐Fen Wu. Overview of software options for processing, analysis and interpretation of mass spectrometric proteomic data. Journal of Mass Spectrometry, 2014. https://doi.org/10.1002/jms.3275 ↩

4. Chen Chen, Jie Hou, John J. Tanner, et al. Bioinformatics Methods for Mass Spectrometry-Based Proteomics Data Analysis. International Journal of Molecular Sciences, 2020. https://doi.org/10.3390/ijms21082873 ↩

5. Jorrit Boekel, John M Chilton, Ira R Cooke, et al. Multi-omic data analysis using Galaxy. Nature Biotechnology, 2015. https://doi.org/10.1038/nbt.3134 ↩

6. Matthiesen Rune. Mass Spectrometry Data Analysis in Proteomics. 2006. https://doi.org/10.1385/1597452750 ↩

7. Masahiro Sugimoto, Masato Kawakami, Martin Robert, et al. Bioinformatics Tools for Mass Spectroscopy-Based Metabolomic Data Processing and Analysis. Current Bioinformatics, 2012. https://doi.org/10.2174/157489312799304431 ↩

8. Alexander S. Kirpich, Miguel Ibarra, Oleksandr Moskalenko, et al. SECIMTools: a suite of metabolomics data analysis tools. BMC Bioinformatics, 2018. https://doi.org/10.1186/s12859-018-2134-1 ↩

9. Alice Cambiaghi, Manuela Ferrario, Marco Masseroli. Analysis of metabolomic data: tools, current strategies and future challenges for omics data integration. Briefings in Bioinformatics, 2016. https://doi.org/10.1093/bib/bbw031 ↩