Bio Proteomics Ptm Analysis

Version Compatibility

Reference examples tested with: MSstatsPTM 2.4+, pandas 2.2+, numpy 1.26+, scipy 1.12+

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.

PTM and Phosphoproteomics Analysis -- Three Inference Layers Stacked on a Biased Extraction

"Find the regulated phosphosites in my enriched samples" -> Localize each modification, then test whether its abundance change survives subtracting the protein-level change -- because a PTM result is three separate inferences (enrichment, localization, quantification) and each fails silently if the layer below is treated as solved.

• R: MSstatsPTM::groupComparisonPTM() for protein-adjusted site testing (the load-bearing tool)
• Python: pandas to expand MaxQuant Phospho (STY)Sites multiplicity and filter localization probability
• R: KSEAapp / PTM-SEA (ssGSEA2.0) for kinase-activity inference from the site fold-changes

Scope: this skill OWNS enrichment-chemistry framing, site localization and FLR, multiplicity-resolved site quant, protein-level adjustment, motif analysis, and kinase-activity inference. Peptide identification and open/variable-mod search route to peptide-identification; the underlying protein-level (unenriched) quant routes to quantification and differential-abundance; DIA acquisition mechanics route to dia-analysis. OUT OF SCOPE: intact-glycopeptide glycan-composition search (pGlyco3/MSFragger-Glyco) and absolute occupancy from three-ratio SILAC are noted but not implemented here.

The Single Most Important Modern Insight -- A PTM Result Is Three Inferences, Not One

1. Enrichment IS the experiment, and the chemistry is a filter confounded with biology. The data only contain what the chemistry captured. TiO2 and Fe-IMAC give partially-overlapping phosphoproteomes; anti-K-GG enriches ubiquitin, NEDD8, and ISG15 indistinguishably; a lectin reports only its cognate glycoforms. A between-method or between-lab "biological difference" must FIRST be excluded as a chemistry artifact before it is called biology.
2. Identifying a peptide is NOT localizing the modification. A phosphopeptide with two S/T and one phosphate has isobaric positional isomers of identical precursor mass and identical peptide-level score; the localization is a SECOND inference decided only by site-determining fragment ions, with its own error rate (false localization rate, FLR). Target-decoy peptide FDR cannot estimate FLR: a wrong localization is the correct sequence with the mod one residue over, not a decoy sequence (Fermin 2013). A 1% peptide FDR does NOT yield a 1% site FDR -- report them separately.
3. A change in phosphopeptide abundance is NOT a change in phosphorylation (the biggest quant trap). Observed PTM signal ~ (site occupancy) x (protein abundance) x (enrichment/ionization factor), so log2FC(PTM_observed) = log2FC(occupancy) + log2FC(protein). Without a paired global (unenriched) proteome run on the SAME samples to subtract log2FC(protein), every protein-abundance change masquerades as a regulated site. Because co-regulated proteins move together, the false positives are pathway-coherent and look biologically convincing -- the worst kind of artifact. This is the entire reason MSstatsPTM exists (Kohler 2023).
4. Most identified sites have no known function. Fewer than ~5% of phosphosites are functionally annotated; a fold-change alone says nothing about regulatory relevance. Functional triage (conservation, stoichiometry, Ochoa functional score, confident kinase assignment) is a separate fourth layer on top of the quant (Ochoa 2020).

Bottom line: report THREE numbers, not one -- peptide/PSM FDR, per-site localization probability with its threshold, and an empirically estimated global FLR -- and never call a site "regulated" from a phospho-only run without protein-level adjustment.

Tool Taxonomy

Enrichment chemistry (phospho)

Naming trap: Ti4+/Zr4+-IMAC (chelated ions) is DIFFERENT chemistry from TiO2/ZrO2 (bulk oxide). Glycolic acid is the modern additive standard (load 80% ACN / 5% TFA / 0.1 M glycolic acid).

Other-PTM enrichment and identity traps

Localization scoring

Site-FDR / FLR (the layer most pipelines skip)

Kinase-activity inference

Benchmark result (Mueller-Dott 2025): across ~19 methods, simple z-score (KSEA/RoKAI) matched or beat sophisticated methods. Performance is PRIOR-limited, not algorithm-limited; all methods inherit PhosphoSitePlus curation bias toward CK2/CDK1/PKA/MAPK, and the dark kinome is structurally invisible. Spend effort on the substrate prior, not the estimator.

Decision Tree by Scenario

Default when uncertain: localize with the search engine's probability (class I >=0.75), expand MaxQuant multiplicity, run MSstatsPTM with a paired global proteome, and call only ADJUSTED.Model hits regulated.

Expand the MaxQuant Site Table Before Any Quant

Goal: Produce a long, multiplicity-resolved, class-I-filtered phosphosite intensity matrix from Phospho (STY)Sites.txt.

Approach: Each site row spreads its quant across Intensity___1/___2/___3 (singly/doubly/triply-phospho forms, THREE underscores); the collapsed base Intensity mixes phospho-states and can fake dephosphorylation. Drop Reverse/contaminant, filter Localization prob, then melt the per-multiplicity columns into rows.

import pandas as pd
import numpy as np

# Filename has a SPACE in the modification name; accept either form.
phospho = pd.read_csv('Phospho (STY)Sites.txt', sep='\t', low_memory=False)

# Newer MaxQuant uses 'Potential contaminant'; older uses 'Contaminant'.
contaminant_col = 'Potential contaminant' if 'Potential contaminant' in phospho.columns else 'Contaminant'
phospho = phospho[(phospho['Reverse'] != '+') & (phospho[contaminant_col] != '+')]

CLASS_I_PROB = 0.75  # Olsen 2006 class-I convention; comparability standard, not a calibrated FLR
phospho = phospho[phospho['Localization prob'] >= CLASS_I_PROB].copy()

gene = phospho['Gene names'].where(phospho['Gene names'].notna(), phospho['Protein'])
phospho['site_id'] = gene.str.split(';').str[0] + '_' + phospho['Amino acid'] + phospho['Position'].astype(int).astype(str)

# Multiplicity columns carry THREE underscores: collapsing them mixes phospho-states.
mult_cols = [c for c in phospho.columns if '___' in c and c.split('___')[-1] in {'1', '2', '3'} and c.startswith('Intensity')]
long = phospho.melt(id_vars=['site_id', 'Amino acid', 'Position', 'Localization prob'], value_vars=mult_cols, var_name='run_multiplicity', value_name='intensity')
long['multiplicity'] = long['run_multiplicity'].str.split('___').str[-1]
long['run'] = long['run_multiplicity'].str.replace(r'___[123]$', '', regex=True).str.replace('Intensity ', '', regex=False)
long = long[long['intensity'] > 0]
long['log2_intensity'] = np.log2(long['intensity'])

Protein-Level Adjustment with MSstatsPTM

Goal: Decide whether each site change is real after subtracting the matched protein-abundance change.

Approach: MSstatsPTM carries TWO datasets -- a PTM dataset (enriched) and a PROTEIN dataset (global/unenriched). groupComparisonPTM fits independent linear models to each and returns a list of THREE: PTM.Model (unadjusted), PROTEIN.Model, and ADJUSTED.Model. The adjustment is dFC_adj = dFC_PTM - dFC_protein with SE_adj = sqrt(SE_PTM^2 + SE_protein^2), so adjustment ADDS uncertainty -- a site can be significant unadjusted yet lose significance after adjustment. A confident regulation call requires significance in ADJUSTED.Model.

library(MSstatsPTM)

# Converters are <Tool>toMSstatsPTMFormat and return a list with $PTM and $PROTEIN.
# MaxQtoMSstatsPTMFormat reads the MaxQuant 'evidence.txt' (NOT the Phospho (STY)Sites
# table -- the pandas multiplicity-expansion above is a SEPARATE workflow); the FASTA maps
# peptides back to site coordinates. Supply BOTH the enriched evidence and the global
# proteinGroups; without the protein dataset there is nothing to adjust against.
# Arg-name note: the FASTA argument is `fasta_path` in current MSstatsPTM; older builds may
# differ -- run `?MaxQtoMSstatsPTMFormat` to confirm before relying on it.
input <- MaxQtoMSstatsPTMFormat(
  evidence = read.table('evidence.txt', sep = '\t', header = TRUE, quote = ''),
  annotation = read.csv('annotation_ptm.csv'),
  fasta_path = 'uniprot_human.fasta',
  fasta_protein_name = 'uniprot_ac',
  proteinGroups = read.table('proteinGroups.txt', sep = '\t', header = TRUE, quote = ''),
  annotation_protein = read.csv('annotation_protein.csv'),
  mod_id = '\\(Phospho \\(STY\\)\\)',
  which_proteinid_ptm = 'Proteins',
  use_unmod_peptides = FALSE
)

summarized <- dataSummarizationPTM(input, use_log_file = FALSE)
# LabelFree run: data.type = 'LF' (use 'TMT' for isobaric); contrast.matrix defaults to
# full pairwise. groupComparisonPTM has NO `model` argument -- it always fits independent
# PTM and PROTEIN models, then adjusts.
result <- groupComparisonPTM(summarized, data.type = 'LF')

# Three models; the adjusted one is the deliverable.
adjusted <- result$ADJUSTED.Model
regulated <- adjusted[!is.na(adjusted$adj.pvalue) & adjusted$adj.pvalue < 0.05 & abs(adjusted$log2FC) > 1, ]

# How much of each call was protein-driven: compare PTM.Model vs ADJUSTED.Model.

Motif Analysis with the Correct Background

Goal: Find kinase/writer motifs around the modified residue without rediscovering amino-acid composition bias.

Approach: Use the Sequence window (+/-15 residues, 31-mer) MaxQuant already provides, centered on the site. The background MUST be an experiment-matched S/T/Y set drawn from the identified proteins (or a central-residue-preserving shuffle), NOT the whole proteome or IUPAC-random -- those just report the composition of phospho-rich disordered regions. motif-x and MoMo p-values are only valid when the background is built this way.

from collections import Counter

# 'Sequence window' is a 31-mer (+/-15) centered on the modified residue.
WINDOW_HALF = 7  # +/-7 flanking is the standard kinase-motif window
foreground = [w[15 - WINDOW_HALF: 16 + WINDOW_HALF] for w in confident['Sequence window'].dropna() if len(w) >= 31]

# Background: same-residue windows from the matched dataset, NOT the whole proteome.
def position_frequencies(windows):
    counts = {i: Counter() for i in range(-WINDOW_HALF, WINDOW_HALF + 1)}
    for w in windows:
        for offset, aa in zip(range(-WINDOW_HALF, WINDOW_HALF + 1), w):
            if aa not in '_X':
                counts[offset][aa] += 1
    return counts

For a publication-grade enrichment logo, hand the foreground and a matched background to a dedicated tool (motif-x / MoMo) and render with data-visualization/sequence-logos.

A Note on Home-Grown Ascore

The function below is an ILLUSTRATIVE approximation, NOT real Ascore. Real Ascore (Beausoleil 2006) competes the best localization against the second-best, sweeps peak depth 1-10 per 100 Th, and restricts to site-determining ions -- none of which this captures. Use the search engine's own localization probability (MaxQuant Localization prob, ptmRS, PTMProphet) for real work, or pyOpenMS AScore (introspect the exact API before relying on it). The home-grown form is here only to show the binomial intuition.

import numpy as np
from scipy.stats import binom

def illustrative_localization_score(matched_site_ions, total_ions, depth_p=0.04):
    '''Binomial intuition only; NOT Ascore (no best-vs-second competition or depth sweep).'''
    if total_ions == 0 or matched_site_ions == 0:
        return 0.0
    p_random = 1 - binom.cdf(matched_site_ions - 1, total_ions, depth_p)
    return -10 * np.log10(p_random) if p_random > 0 else 100.0

Per-Method Failure Modes

Skipping protein-level adjustment

Trigger: Differential testing on a phospho-only run with no paired global proteome. Mechanism: log2FC(PTM_observed) = log2FC(occupancy) + log2FC(protein); the two terms are inseparable. Symptom: Pathway-coherent "regulated sites" that are pure protein-abundance changes (cyclins/histones in cell cycle, stabilized substrates under drug). Fix: Run a matched global proteome and adjust via MSstatsPTM; route the protein-level quant to quantification.