Bio Proteomics Quantification

Version Compatibility

Reference examples tested with: MSstats 4.10+, MSnbase 2.28+, iq 1.9+, numpy 1.26+, pandas 2.2+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures
• R: packageVersion('<pkg>') then ?function_name to verify parameters

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Protein Quantification -- Reconstructing Protein Abundance from Ions Whose Physical Origin Dictates the Irreducible Error

"Quantify proteins from my mass spec data" -> Reconstruct a protein-by-sample abundance matrix from peptide/reporter ion signals, choosing a summarizer and normalizer that match where the signal physically came from -- because the measurement's physical origin sets an error that no normalization can remove.

• R: MSstats::dataProcess() (Tukey median polish summarization) for label-free feature-to-protein
• R: iq::maxLFQ() for the real MaxLFQ algorithm (delayed normalization + maximal peptide-ratio least-squares)
• R: MSnbase::quantify(reporters=TMT10) + purityCorrect() for isobaric reporter extraction
• R/Python: sample-loading + IRS scaling to bridge TMT plexes; per-sample median centering for LFQ

Scope: this skill OWNS converting peptide/PSM/reporter signal into a normalized protein abundance matrix (LFQ, TMT/iTRAQ, SILAC; summarization; normalization; IRS). Statistical testing of the matrix -> differential-abundance. DIA quant mechanics and DIA-NN execution -> dia-analysis. Parsing MaxQuant/DIA-NN outputs -> data-import. Razor/shared-peptide group assignment -> protein-inference. OUT OF SCOPE: missing-value imputation and the downshift false-positive trap (modeled in differential-abundance), and absolute copy-number calibration beyond a one-line pointer.

The Single Most Important Modern Insight

1. Every quant method answers "where does the signal physically come from?" differently, and that physical origin dictates the error structure no normalization can remove. LFQ measures MS1 precursor area (or DIA fragment area) in SEPARATE runs -> the irreducible error is run-to-run variation plus stochastic, left-censored (MNAR) missingness. Isobaric TMT/iTRAQ measures low-m/z reporter ions from CO-ISOLATED, co-eluting peptides in ONE spectrum -> the irreducible error is RATIO COMPRESSION toward 1:1, a PHYSICAL co-isolation effect (interloper reporters add roughly equally to every channel), attacked at the instrument by SPS-MS3 and never fully undone in software. SILAC measures a heavy/light MS1 pair in the SAME scan -> lowest per-ratio variance, but its irreducible vulnerabilities are incomplete labeling and Arg->Pro label scrambling, which bias every ratio and cannot be corrected post hoc because channels are combined before any MS (a correctly mixed sample carries no inherent mixing error). Teach the signal origin and every threshold and failure mode below follows from it.

2. The peptide-to-protein SUMMARIZATION choice is the highest-leverage decision in the pipeline, and it is invisible in the output. In log space a peptide intensity is protein abundance + a peptide effect (ionization efficiency, flyability, missed cleavages, modifications) that spans orders of magnitude and is partly context-dependent, plus a run effect. Sum is dominated by the highest-flying peptide (dropout collapses it -> a fold change driven by detectability, not biology); mean is unbiased only if the detected peptide SET is identical across runs (it is not under MNAR); median discards relative-intensity information. Benchmarks confirm the quantification method is a dominant driver of which proteins are called differential (Lin 2022). Report the summarizer as prominently as the test, and run a sensitivity analysis across >=2 summarizers -- that is where the answer is most likely to move.

Tool Taxonomy

Decision Tree by Scenario

Default when uncertain: label-free DDA -> MSstats::dataProcess with summaryMethod='TMP', normalization='equalizeMedians'; report the summarizer alongside results and sanity-check against iq::maxLFQ().

Label-Free Summarization and Normalization

Summarize peptides to proteins with MSstats

Goal: Turn MaxQuant feature-level evidence into a normalized protein-level abundance matrix.

Approach: Reformat to MSstats input, then dataProcess applies median equalization and Tukey median polish (robust to outlier peptides, 50% breakdown) with censored-value handling for label-free missingness.

library(MSstats)

maxquant_input <- MaxQtoMSstatsFormat(
    evidence = read.table('evidence.txt', sep = '\t', header = TRUE),
    proteinGroups = read.table('proteinGroups.txt', sep = '\t', header = TRUE),
    annotation = read.csv('annotation.csv')
)

# TMP = Tukey median polish; censoredInt='NA' treats missing intensities as left-censored
processed <- dataProcess(maxquant_input, normalization = 'equalizeMedians',
                         summaryMethod = 'TMP', censoredInt = 'NA', MBimpute = FALSE)

protein_abundance <- processed$ProteinLevelData

Run the real MaxLFQ (not median centering)

Goal: Produce MaxLFQ protein intensities from a peptide quant matrix.

Approach: Call iq::maxLFQ(), which implements Cox 2014 delayed normalization and maximal peptide-ratio least-squares. Per-sample median centering shares only the name and silently gives a different answer.

library(iq)

# rows = peptide ions, columns = samples, values = log2 intensities for ONE protein group
result <- maxLFQ(peptide_log2_matrix)
protein_estimate <- result$estimate    # one MaxLFQ value per sample

Median-center label-free intensities (a normalizer, not a summarizer)

Goal: Correct per-sample loading differences before testing.

Approach: Subtract each sample's median log2 intensity (corrects LOCATION only; it cannot manufacture variance, so it is the safe default). Median centering normalizes; it does NOT summarize peptides to proteins.

import numpy as np
import pandas as pd

log_int = np.log2(intensities.replace(0, np.nan))    # MaxQuant writes 0 for missing; log2(0) = -inf
sample_medians = log_int.median(axis=0)
normalized = log_int - sample_medians + sample_medians.median()

Isobaric (TMT/iTRAQ) Quantification

Extract and impurity-correct reporter ions

Goal: Pull TMT reporter intensities from spectra and correct cross-channel isotope bleed.

Approach: Read spectra on disk, quantify the reporter region, then purityCorrect with a LOT-SPECIFIC impurity matrix from the reagent Certificate of Analysis. readMSnSet reads an already-quantified text matrix and does NOT extract reporters.

library(MSnbase)

raw <- readMSData('experiment.mzML', mode = 'onDisk')
# method='max' for centroided spectra; reporters=TMT10 defines the 126-131 reporter m/z
quant <- quantify(raw, reporters = TMT10, method = 'max')

# makeImpuritiesMatrix has manufacturer-default templates; REPLACE with lot-specific Certificate values
imp <- makeImpuritiesMatrix(x = 10)
quant <- purityCorrect(quant, imp)

Bridge multiple TMT plexes with IRS

Goal: Make reporter intensities comparable across separate TMT runs.

Approach: Absolute reporter intensities for the same protein differ 2-5x between plexes because each plex samples a random point on the elution profile. Sample-loading normalization fixes within-run loading; the Internal Reference Scaling bridge (Plubell 2017) then pins each plex's pooled reference channel to a common per-protein value. Order: SL, then IRS.

import numpy as np
import pandas as pd

# protein_psm_sums: protein x channel, summed PSM reporter ions; one reference channel per plex
def sample_loading_normalize(plex):
    target = plex.sum(axis=0).mean()    # common target = mean column sum within the plex
    return plex * (target / plex.sum(axis=0))

def irs_scale(plexes, ref_cols):
    refs = pd.concat([p[ref] for p, ref in zip(plexes, ref_cols)], axis=1)
    geomean = np.exp(np.log(refs.replace(0, np.nan)).mean(axis=1))    # per-protein geometric mean of references
    out = []
    for p, ref in zip(plexes, ref_cols):
        factor = geomean / p[ref]    # per-protein per-plex scaling factor
        out.append(p.mul(factor, axis=0))
    return out

SILAC Quantification

Goal: Compute heavy/light ratios while preserving on/off biology and flagging label artifacts.

Approach: A protein present only in the heavy channel is the interesting biology, not a NaN to discard. Verify labeling efficiency (>=95%, target 97-98%) on a heavy-only pilot and assess Arg->Pro conversion before trusting any ratio.

import numpy as np

# Arg10/Lys8 is the common pairing (avoids overlap with the +6 isotope envelope)
SILAC_SHIFTS = {'Arg10': 10.008269, 'Lys8': 8.014199, 'Arg6': 6.020129, 'Lys6': 6.020129}

def silac_log2_ratio(heavy, light):
    if heavy > 0 and light > 0:
        return np.log2(heavy / light)
    if heavy > 0 and light == 0:
        return np.inf     # present only in heavy: real on/off biology, do NOT discard as NaN
    if light > 0 and heavy == 0:
        return -np.inf
    return np.nan

Per-Method Failure Modes

MaxLFQ reimplemented as median centering

Trigger: A homebrew function named maxlfq that only subtracts per-sample medians. Mechanism: Real MaxLFQ is delayed normalization plus a maximal peptide-ratio least-squares solve; median centering shares only the name. Symptom: Plausible-looking but systematically different intensities; unbalanced peptide sets handled wrongly. Fix: Call iq::maxLFQ(), DIA-NN, or MaxQuant.

TMT ratio compression

Trigger: MS2-only reporter quant on a complex sample. Mechanism: Co-isolated interloper peptides add reporters roughly equally to every channel, pulling large true ratios toward 1:1; PHYSICAL, not removable by normalization. Symptom: "Nothing is significant"; attenuated fold changes, inflated false negatives. Fix: SPS-MS3 acquisition, narrower isolation windows, FAIMS/ion mobility, or complement-reporter methods; a PIF filter helps but MS1 purity underestimates true interference (Savitski 2013).

TMT plexes concatenated without IRS

Trigger: Stacking reporter intensities from multiple plexes directly. Mechanism: Absolute reporter intensities for one protein differ 2-5x across runs from random elution-profile sampling, unrelated to abundance. Symptom: Plex appears as the dominant axis of variation; spurious cross-plex differences. Fix: Include a pooled reference channel in EVERY plex; apply sample-loading then IRS (Plubell 2017).

Isotopic impurity matrix mis-applied

Trigger: Using the default/example impurity matrix, the wrong lot, or a transposed/mis-ordered (127N vs 127C) matrix. Mechanism: A few percent of each channel bleeds to +/-1 Da neighbors; wrong values mis-subtract, negative corrected intensities get clipped. Symptom: Adjacent channels silently biased; extreme contrasts placed in adjacent channels confounded. Fix: Use the lot-specific Certificate of Analysis values; randomize channel-to-condition assignment.