Imported from huang-sh/DeepScience (packages/science/skills/biology/proteomics/bio-proteomics-proteomics-qc/SKILL.md). Install upstream with npx skills add huang-sh/DeepScience --skill bio-proteomics-proteomics-qc. Copyright stays with the author.

Version Compatibility

Reference examples tested with: pandas 2.2+, numpy 1.26+, scipy 1.12+, matplotlib 3.8+, scikit-learn 1.4+, limma 3.58+

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.

Proteomics Quality Control -- A Three-Level Funnel Where the Matrix Sees the Failure Last

"Check the quality of my proteomics data" -> Read instrument, identification, and quantitative metrics as a descending funnel of silent failures, and inspect raw signal BEFORE normalizing -- because by the time a fault reaches the deliverable matrix, normalization has usually erased the evidence.

• Python: pandas for matrix QC; matplotlib/seaborn for raw boxplots, correlation heatmaps, PCA
• R: PTXQC createReport() for MaxQuant search-table QC; limma::plotMDS()/plotDensities(); MSstatsTMT dataProcessPlotsTMT() for TMT channel balance

Scope: This skill OWNS QC diagnosis across all three levels -- which metric localizes which fault, what threshold means trouble, and the mandatory inspect-before-normalize ordering. Normalization mechanics route to quantification; the differential test routes to differential-abundance; DIA q-value computation routes to dia-analysis. OUT OF SCOPE: running the statistical test, the normalization algorithms themselves, and DIA q-value/FDR internals.

The Single Most Important Modern Insight -- QC Is a Three-Level Funnel and Normalization Hides the Evidence

1. QC is a three-level funnel of silent failures, and the deliverable matrix is the LAST place a problem becomes visible. Faults originate at the instrument (spray, calibration, column) or in identification (digestion, contamination, PTM artifacts), but a protein matrix only shows the downstream symptom -- a low correlation or an outlier sample. A matrix-only QC pass is one-third of the job and blind to where faults actually start. Localize by descending: read every metric together with its co-readouts, never alone.

2. Almost no metric has a universal pass/fail cutoff; the defensible practice is a per-instrument, per-method control chart. Deviation from a lab's own rolling baseline (Levey-Jennings, +/-2 SD warn, +/-3 SD action) detects faults that a constant threshold misses or false-flags (Neely and Palmblad 2024). The numbers below seed a control chart, they are not standards-body limits.

3. Median normalization HIDES loading problems that must be SEEN first. Median/quantile normalization works by forcing a chosen summary statistic of every sample equal. A sample that genuinely loaded 3x low sits visibly shifted down in a RAW boxplot -- an obvious, diagnosable defect. The instant the matrix is median-normalized, the algorithm shifts that sample up by a constant to match everyone's median; boxplots line up perfectly; the evidence is mathematically erased. Worse, the low-loaded sample's noisy low-abundance signal gets stretched up to mid-range and injected into the differential test while QC plots look pristine. MANDATE: inspect raw/un-normalized boxplots plus per-sample ID counts, total signal, and missing fraction BEFORE normalizing; remove loading/injection failures and contaminants; THEN normalize and re-plot on the survivors. The identical principle governs TMT channel-loading balance.

The Three Levels in One Table

The most integrative metric (% MS2 identified / ID count) is the first alarm but the LEAST specific -- it moves whenever anything upstream degrades. The same protein-count drop means spray (erratic TIC + maxed injection time), column (lost RT + broad peaks + rising backpressure), or sample (high contaminant fraction) depending on what co-moves.

Tool Taxonomy

Decision Tree by Scenario

Default when uncertain: plot the RAW per-sample boxplots, ID counts, total signal, and missing fraction first; remove loading/injection failures and contaminants; only then normalize, re-plot, and proceed to correlation/CV/PCA on the survivors.

Inspect Raw Signal and Remove Contaminants Before Normalizing

Goal: Catch loading/injection failures and strip contaminant/decoy rows while they are still visible -- before normalization erases them.

Approach: Load the un-normalized matrix, plot per-sample boxplots plus ID counts and total signal, filter MaxQuant Potential contaminant/Reverse/Only identified by site rows, THEN log-transform and normalize on the survivors.

import pandas as pd
import numpy as np

contaminant_flags = ['Potential contaminant', 'Reverse', 'Only identified by site']

def strip_contaminant_rows(protein_groups):
    keep = pd.Series(True, index=protein_groups.index)
    for col in contaminant_flags:
        match = next((c for c in protein_groups.columns if c.lower() == col.lower()), None)  # MaxQuant casing varies by version -- match case-insensitively
        if match is not None:
            keep &= protein_groups[match].fillna('') != '+'  # MaxQuant marks flagged rows with a literal '+'
    return protein_groups[keep]

def raw_sample_qc(raw_intensities):
    return pd.DataFrame({
        'n_quantified': raw_intensities.notna().sum(),
        'total_signal': raw_intensities.sum(),
        'median_intensity': raw_intensities.median(),
        'missing_pct': 100 * raw_intensities.isna().sum() / len(raw_intensities)})

Read the boxplots before normalizing: a sample shifted >=2-3x below its group median is a loading/injection failure to exclude, not to rescale. The contaminant fraction of summed intensity should be small (PTXQC default flags >1%); keratin and trypsin autolysis dominate LOW-INPUT samples (single-cell, IPs, gel bands) because they are a roughly fixed absolute amount whose fractional share explodes as load shrinks.

Replicate Correlation on log2

Goal: Quantify reproducibility without letting a few abundant proteins fake agreement.

Approach: Correlate on log2 intensities (variance-stabilized, high-abundance tail compressed), report within-group pairs, and flag a sample correlating better with another group as a possible swap.

from itertools import combinations

def replicate_correlation(log2_intensities, sample_groups):
    corr = log2_intensities.corr(method='pearson')  # log2 first: Pearson on raw is a high-abundance artifact
    rows = []
    for group in sample_groups.unique():
        members = sample_groups[sample_groups == group].index
        for s1, s2 in combinations(members, 2):
            rows.append({'group': group, 's1': s1, 's2': s2, 'r': corr.loc[s1, s2]})
    return pd.DataFrame(rows)

Technical replicates r > 0.98 (instrument noise only); biological r ~ 0.90-0.98 (genuine variance, lower is expected and correct); soft floor r > 0.8 to retain a biological replicate. A Spearman check is a robustness aid only -- ranks discard the magnitude that quant QC cares about.

Coefficient of Variation on the Linear Scale

Goal: Summarize per-condition precision with a number that means what it says.

Approach: Compute CV = SD/mean on LINEAR (non-log) intensities; if only logged values exist use the geometric-CV formula. Report the median CV per condition (the per-protein distribution is right-skewed).

def median_cv_linear(linear_intensities, sample_groups):
    rows = []
    for group in sample_groups.unique():
        block = linear_intensities[sample_groups[sample_groups == group].index]
        per_protein_cv = block.std(axis=1) / block.mean(axis=1)  # base CV formula REQUIRES linear scale
        rows.append({'group': group, 'median_cv_pct': 100 * per_protein_cv.median()})
    return pd.DataFrame(rows)

def geometric_cv_from_log(log_intensities):
    sigma = log_intensities.std(axis=1) * np.log(2)  # convert log2 SD to natural-log SD
    return 100 * np.sqrt(np.expm1(sigma ** 2))  # gCV = sqrt(exp(sigma^2) - 1)

Applying the base formula to log-transformed data compresses CV ~14x (most proteins appear to have CV < 1%) -- meaningless (Brenes 2024). State normalization state, transform, and software params or the CV is uninterpretable: DIA-NN "High precision" mode silently median-normalizes, halving median CV vs "High accuracy". Technical median CV < ~10-20%, biological ~20-40%; a LOWER CV is not automatically better (loose FDR or faulty MS1 extraction produce artificially low CVs).

Missingness Mechanism and Completeness

Goal: Decide how to impute by first deciding why values are missing.

Approach: Diagnose the missingness profile -- left-tail concentration means MNAR (left-censored, abundance-dependent), all-abundance scatter means MCAR -- and filter on completeness before imputing only the shallow remainder.

def missingness_profile(log2_intensities, n_bins=20):
    observed = log2_intensities.stack()
    abundance_bins = pd.qcut(observed, n_bins, duplicates='drop')
    present_per_protein = log2_intensities.notna().mean(axis=1)
    mean_abundance = log2_intensities.mean(axis=1)
    return mean_abundance, present_per_protein  # plot present-fraction vs abundance: rising-with-abundance = MNAR

def completeness_filter(log2_intensities, sample_groups, min_valid_frac=0.7):
    keep = pd.Series(False, index=log2_intensities.index)
    for group in sample_groups.unique():
        block = log2_intensities[sample_groups[sample_groups == group].index]
        keep |= block.notna().mean(axis=1) >= min_valid_frac  # valid in >=70% of >=1 condition
    return log2_intensities[keep]

kNN-imputing a genuinely-absent (MNAR) value invents mid-range abundance and KILLS a real present/absent difference; a left-shifted draw (Perseus down-shifted normal, downshift=1.8 SD below the observed mean, width=0.3 of observed SD) on an MCAR gap FABRICATES a false low and inflates a difference. Match the imputer to the mechanism. The imputation mechanics themselves are quantification.

PCA and Batch Detection

Goal: See whether the dominant variance is biology or batch, and flag outlier samples.

Approach: On the normalized survivors, run PCA, color by condition and by batch, and test whether top PCs associate with batch.

from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from scipy.stats import f_oneway

def pca_batch_check(normalized_log2, sample_info, batch_col='batch'):
    imputed = normalized_log2.apply(lambda r: r.fillna(r.median()), axis=1)  # temporary, for PCA only
    pcs = PCA(n_components=5).fit(StandardScaler().fit_transform(imputed.T))
    coords = pd.DataFrame(pcs.transform(StandardScaler().fit_transform(imputed.T)),
                          columns=[f'PC{i+1}' for i in range(5)], index=normalized_log2.columns).join(sample_info)
    for pc in ['PC1', 'PC2', 'PC3']:
        groups = [coords[coords[batch_col] == b][pc] for b in coords[batch_col].unique()]
        _, p = f_oneway(*groups)
        print(f'{pc} ~ {batch_col}: p={p:.4f}')
    return coords, pcs.explained_variance_ratio_

A sample isolated from its group is a removal/re-run candidate. If batch is PC1, correct it explicitly (ComBat, or include batch in the design matrix downstream) and re-inspect; never let batch be the dominant axis going into differential testing. Visualization of the projection routes to data-visualization/dimensionality-reduction-plots.

Per-Method Failure Modes

Median normalization hides loading failures

Trigger: Normalizing the matrix before inspecting raw per-sample signal. Mechanism: median-centering shifts each sample by a constant to equalize the very statistic that was the symptom of a low load. Symptom: flat, clean boxplots that hide a 3x-low sample now stretched into mid-range. Fix: plot RAW boxplots + ID counts + total signal first; exclude failures; then normalize.

Normalizing with contaminants still in the matrix

Trigger: Contaminant/decoy rows left in before log + normalize. Mechanism: keratin/trypsin/albumin inflate the denominator and shift the median; when their load differs across groups the differential gets normalized into the real proteins. Symptom: spurious fold changes; a contaminant fraction that varies by group. Fix: filter Potential contaminant + Reverse + Only identified by site BEFORE log + normalize.

Wrong imputer for the missingness mechanism

Trigger: kNN on MNAR, or left-shift on MCAR. Mechanism: kNN borrows mid-range neighbors for a value that is low because it is absent; left-shift draws a deep low for a value missing at random. Symptom: killed present/absent calls (kNN-on-MNAR) or inflated false lows (left-shift-on-MCAR). Fix: diagnose left-tail vs all-abundance from the histogram first; mechanics route to quantification.

CV computed on log-transformed data

Trigger: Base CV formula applied after log2. Mechanism: SD/mean is defined for linear intensity; logging compresses it ~14x. Symptom: most proteins appear to have CV < 1%. Fix: compute on linear intensity, or use the geometric-CV formula on logged values; always state transform + normalization + software.

Pearson r on raw intensity

Trigger: Correlating un-logged intensities. Mechanism: a few high-abundance proteins dominate the covariance. Symptom: r = 0.99 while the bulk disagrees. Fix: log2 before correlating; Spearman as a robustness check only.

Bimodal mass-error histogram read as calibration

Trigger: A two-peaked ppm-error distribution. Mechanism: almost always monoisotopic mis-assignment or co-isolation (a search/sample problem), NOT calibration drift; a generous tolerance still IDs the mis-assigned precursors so ID rate looks fine. Symptom: bimodal histogram, normal ID rate. Fix: correct isotope-error tolerance/deisotoping, not recalibration.

Charge-state distribution off the platform baseline

Trigger: The fully-tryptic 2+ fraction drifts from the rolling baseline. Mechanism: a tryptic peptide carries two basic sites (C-terminal K/R + N-terminus) so 2+ dominates; excess 3+/4+ comes from internal basic residues left by missed cleavages, excess 1+ from poor ionization, short peptides, or contaminants. Symptom: raised high-charge fraction (a digestion/chemistry signal) or raised 1+ fraction (an ionization/spray signal). Fix: read the charge distribution together with the missed-cleavage rate (high charge co-moving with missed cleavages = under-digestion) to separate a chemistry problem from a spray problem a raw protein-count drop cannot resolve alone.

TMT ratio compression from co-isolation

Trigger: Isobaric quant with a wide isolation window. Mechanism: near-isobaric co-eluting precursors are co-isolated and add their own reporters across all channels, a uniform pedestal. Symptom: every fold-change compressed toward 1 (a real 10:1 reads ~5:1). Fix: filter isolation interference < 50%; prefer SPS-MS3 (McAlister 2014) / FAIMS / narrow windows; run a TKO or empty-channel control to measure the floor.

Quantitative Thresholds

Common Errors

References

• Bielow C, Mastrobuoni G, Kempa S. Proteomics Quality Control: Quality Control Software for MaxQuant Results. J Proteome Res 2016;15(3):777-787.
• Kovalchik KA, Colborne S, Spencer SEP, et al. RawTools: Rapid and Dynamic Interrogation of Orbitrap Data Files for Mass Spectrometer System Management. J Proteome Res 2019;18(2):700-708.
• Morgenstern D, Barzilay R, Levin Y. RawBeans: A Simple, Vendor-Independent, Raw-Data Quality-Control Tool. J Proteome Res 2021;20(4):2098-2104.
• Kockmann T, Panse C. The rawrr R Package: Direct Access to Orbitrap Data and Beyond. J Proteome Res 2021;20(4):2028-2034.
• Trachsel C, Panse C, Kockmann T, et al. rawDiag: An R Package Supporting Rational LC-MS Method Optimization for Bottom-up Proteomics. J Proteome Res 2018;17(8):2908-2914.
• Ma ZQ, Polzin KO, Dasari S, et al. QuaMeter: Multivendor Performance Metrics for LC-MS/MS Proteomics Instrumentation. Anal Chem 2012;84(14):5845-5850.
• Bereman MS, Beri J, Sharma V, et al. An Automated Pipeline to Monitor System Performance in Liquid Chromatography-Tandem Mass Spectrometry Proteomic Experiments. J Proteome Res 2016;15(12):4763-4769.
• Kocher T, Swart R, Mechtler K. Ultra-High-Pressure RPLC Hyphenated to an LTQ-Orbitrap Velos Reveals a Linear Relation between Peak Capacity and Number of Identified Peptides. Anal Chem 2011;83(7):2699-2704.
• Tyanova S, Temu T, Sinitcyn P, et al. The Perseus computational platform for comprehensive analysis of (prote)omics data. Nat Methods 2016;13(9):731-740.
• Brenes AJ. Calculating and Reporting Coefficients of Variation for DIA-Based Proteomics. J Proteome Res 2024;23(12):5274-5278.
• McAlister GC, Nusinow DP, Jedrychowski MP, et al. MultiNotch MS3 Enables Accurate, Sensitive, and Multiplexed Detection of Differential Expression across Cancer Cell Line Proteomes. Anal Chem 2014;86(14):7150-7158.
• Huang T, Choi M, Tzouros M, et al. MSstatsTMT: Statistical Detection of Differentially Abundant Proteins in Experiments with Isobaric Labeling and Multiple Mixtures. Mol Cell Proteomics 2020;19(10):1706-1723.
• Neely BA, Palmblad M, et al. Quality Control in the Mass Spectrometry Proteomics Core: A Practical Primer. J Biomol Tech 2024;35(3).

Related Skills

• data-import - Load search-engine output and intensity matrices before QC
• quantification - Normalization and imputation mechanics that QC mandates running AFTER inspection
• differential-abundance - The moderated statistical test QC gates
• dia-analysis - DIA q-value/FDR internals behind the protein-count QC
• data-visualization/dimensionality-reduction-plots - PCA/MDS projection plotting
• workflows/proteomics-pipeline - End-to-end pipeline placing QC before differential testing