Tooluniverse Proteomics Analysis

Proteomics Analysis

RULE ZERO — Check for pre-computed results FIRST

Before following any instruction below, scan the data folder for:

• *_executed.ipynb → read with tu run read_executed_notebook '{"data_folder":"<path>","search":"<keyword>"}' and cite its cell outputs as the authoritative answer
• Pre-computed result files (CSV/TSV with names like *results*, *deseq*, *enrich*, *stats*, *_simplified.csv) → read directly and report the requested value
• Canonical analysis scripts (analysis.R, run_*.py, find_*.R, *.Rmd) → execute as-is and read the output

Only follow this skill's re-analysis recipe below if none of the above exist. Re-running from raw data produces different numbers than the published answer and is much slower (often 5-10× turn count).

---

Comprehensive analysis of mass spectrometry-based proteomics data from protein identification through quantification, differential expression, post-translational modifications, and systems-level interpretation.

When to Use This Skill

Triggers: User has proteomics MS output files, asks about protein abundance/expression, differential protein expression, PTM analysis, protein-RNA correlation, multi-omics integration involving proteomics, protein complex/interaction analysis, or proteomics biomarker discovery.

COMPUTE, DON'T DESCRIBE

When analysis requires computation (statistics, data processing, scoring, enrichment), write and run Python code via Bash. Don't describe what you would do — execute it and report actual results. Use ToolUniverse tools to retrieve data, then Python (pandas, scipy, statsmodels, matplotlib) to analyze it.

Core Capabilities

• Data Import: MaxQuant, Spectronaut, DIA-NN, Proteome Discoverer, FragPipe outputs
• Quality Control: Missing value analysis, intensity distributions, sample clustering
• Normalization: Median, quantile, TMM, VSN — choice depends on experimental design (see Interpretation Framework)
• Imputation: MinProb (MNAR), KNN (MAR), QRILC for missing values
• Differential Expression: Limma, DEP, MSstats for statistical testing
• PTM Analysis: Phospho-site localization, PTM enrichment, kinase prediction
• Protein-RNA Integration: Correlation analysis, translation efficiency
• Pathway Enrichment: Over-representation and GSEA for protein sets
• PPI Analysis: Protein complex detection, interaction networks via STRING/IntAct

Workflow Overview

Input: MS Proteomics Data
    |
Phase 1: Data Import & QC
Phase 2: Preprocessing (filter, impute, normalize)
Phase 3: Differential Expression Analysis
Phase 4: PTM Analysis (if applicable)
Phase 5: Functional Enrichment (GO, KEGG, Reactome)
Phase 6: Protein-Protein Interactions (STRING networks)
Phase 7: Multi-Omics Integration (optional, protein-RNA correlation)
Phase 8: Generate Report

See PHASE_DETAILS.md for detailed procedures per phase.

Integration with ToolUniverse

Quantified Minimums

• At least 500 proteins quantified (human: 3,000+ is reasonable; 10,000+ is deep coverage)
• At least 3 biological replicates per condition (non-negotiable for reliable statistics)
• Filter to proteins with 2+ unique peptides (single-peptide IDs are not reported as DE)
• Statistical test: limma or t-test with Benjamini-Hochberg multiple testing correction
• Pathway enrichment: at least one method (GO, KEGG, or Reactome)
• Report must include: QC summary, DE results with volcano plot, pathway analysis, visualizations

Interpretation Framework

Starting Point: Experimental Design

Quantitative proteomics compares protein abundance. LOOK UP DON'T GUESS — always verify the experimental method, platform, and replicate count before choosing an analysis strategy.

Quantification strategy decision tree:

• Cell culture, high accuracy needed → SILAC (ratios within same MS run, most accurate, but culture-only)
• Multiple conditions, multiplexing needed → TMT/iTRAQ (up to 18-plex in one run; TMM/VSN normalization; beware ratio compression artifact that reduces observed fold-changes)
• Discovery study, flexible design → Label-free (LFQ) (intensity-based; median/quantile normalization; more missing values; wider dynamic range)
• Replicates: n < 3 = unreliable fold changes (variance cannot be estimated). Minimum n = 3 biological replicates; n >= 4 preferred for clinical samples. Never report significance from duplicates.
• FDR cutoff: Benjamini-Hochberg correction mandatory. FDR < 0.05 standard; FDR < 0.01 stringent. Never report unadjusted p-values alone.

Protein Identification Reasoning

Protein identification from MS data follows a logical chain. LOOK UP DON'T GUESS — search UniProt and STRING for protein annotation rather than inferring function from name alone.

1. Peptide mass fingerprinting (PMF): Intact protein digested → measured peptide masses compared against theoretical digest of all database proteins. A match requires >=4 peptides covering >=15% of the protein sequence. Single-peptide hits are unreliable (could match multiple proteins). 2. Tandem MS (MS/MS): Fragment ion spectra matched to peptide sequences via search engines (Andromeda, SEQUEST, X!Tandem). Each peptide-spectrum match (PSM) scored; only PSMs above FDR threshold count. Unique peptides (mapping to one protein) are essential — shared peptides cannot distinguish between protein isoforms. 3. Protein inference: Multiple peptides → protein groups. When peptides are shared between homologs, report the protein group (not individual proteins). Use proteins_api_search or UniProt_search to resolve ambiguous protein groups. 4. Coverage matters: 2+ unique peptides is the minimum for a confident protein ID. Proteins identified by a single unique peptide should be flagged as tentative.

Post-Translational Modification (PTM) Analysis Reasoning

PTMs (phosphorylation, ubiquitination, acetylation, glycosylation) add biological complexity beyond protein abundance.

1. Site localization: A phospho-site is confidently localized only if the localization probability > 0.75 (MaxQuant) or ptmRS score > 75 (Proteome Discoverer). Ambiguous sites should not be reported as specific residue modifications. 2. Enrichment is required: Without phospho-enrichment (TiO2, IMAC), only the most abundant phosphopeptides are detected (~1% of phosphoproteome). An absence of a phospho-site in non-enriched data does not mean it is absent biologically. 3. Kinase prediction: If phospho-sites are identified, predict upstream kinases using motif analysis. Cross-reference with OpenTargets_get_target_safety_profile_by_ensemblID for kinase-disease associations. LOOK UP kinase-substrate relationships in PhosphoSitePlus rather than guessing from sequence motif alone. 4. Stoichiometry: A protein can be 5% or 95% phosphorylated at a given site — this matters enormously for function but is hard to measure. Report whether data supports stoichiometry estimation or only site identification.

Differential Expression Thresholds

• Strong: padj < 0.01, FC > 2.0, ≥5 unique peptides, <20% missing
• Moderate: padj 0.01-0.05, FC 1.5-2.0, 2-5 peptides, 20-50% missing
• Weak/unreliable: padj 0.05-0.1, FC 1.2-1.5, 1-2 peptides (single-peptide proteins are unreliable), >50% missing (imputation needed)

Evidence Grading

• T1: Validated by orthogonal method (Western blot, PRM) + functional study
• T2: Significant DE (padj < 0.05, FC > 1.5) in 2+ biological replicates
• T3: Significant DE in 1 experiment, or significant but low FC
• T4: Identified but not quantified, or single peptide identification

Synthesis Questions

1. How many proteins are differentially expressed? (>500 DE proteins suggests global perturbation; <50 suggests targeted effect) 2. Are key pathway proteins concordantly regulated? (all subunits of a complex changing = high confidence) 3. Do proteomics results correlate with transcriptomics? (low correlation is common — post-translational regulation) 4. Are PTM changes driving the phenotype? (check phosphoproteomics if available) 5. What is the coverage relative to the expected proteome? (human: ~10K quantified is good; <3K is limited)

---

Limitations

• Platform-specific: Optimized for MS-based proteomics (not Western blot quantification)
• Missing values: High missing rate (>50% per protein) limits statistical power
• PTM analysis: Requires enrichment protocols for comprehensive PTM profiling
• Absolute quantification: Relative abundance only (unless TMT/SILAC used)
• Protein isoforms: Typically collapsed to gene level
• Dynamic range: MS has limited dynamic range vs mRNA sequencing

References

Methods: MaxQuant (doi:10.1038/nbt.1511), Limma for proteomics (doi:10.1093/nar/gkv007), DEP workflow (doi:10.1038/nprot.2018.107)

Databases: STRING, PhosphoSitePlus, CORUM

Reference Files

• PHASE_DETAILS.md - Detailed procedures for each analysis phase, including report template

Code Reference: Proteomics Analysis

Data Loading

MaxQuant proteinGroups.txt

def load_maxquant_proteins(protein_groups_file):
    """Load MaxQuant proteinGroups.txt file."""
    import pandas as pd

    df = pd.read_csv(protein_groups_file, sep='\t')

    # Extract intensity columns (LFQ or raw)
    intensity_cols = [col for col in df.columns if 'LFQ intensity' in col or 'Intensity ' in col]

    intensity_matrix = df[intensity_cols].copy()
    intensity_matrix.columns = [col.replace('LFQ intensity ', '').replace('Intensity ', '')
                                 for col in intensity_cols]

    metadata = df[['Protein IDs', 'Gene names', 'Fasta headers',
                   'Peptides', 'Sequence coverage [%]']].copy()

    return intensity_matrix, metadata

Quality Control

Missing Value Assessment

def assess_missing_values(intensity_matrix):
    """Calculate percentage of missing values per protein and sample."""
    missing_per_protein = (intensity_matrix == 0).sum(axis=1) / intensity_matrix.shape[1]
    missing_per_sample = (intensity_matrix == 0).sum(axis=0) / intensity_matrix.shape[0]
    return missing_per_protein, missing_per_sample

Intensity Distribution

def plot_intensity_distributions(intensity_matrix):
    """Plot log10 intensity distributions per sample."""
    import matplotlib.pyplot as plt
    import numpy as np

    log_intensities = np.log10(intensity_matrix.replace(0, np.nan))
    log_intensities.plot(kind='box')
    plt.ylabel('log10 Intensity')
    plt.title('Intensity Distribution per Sample')

Sample Correlation

def plot_sample_correlation(intensity_matrix):
    """Calculate and visualize sample-sample correlation."""
    import numpy as np
    import seaborn as sns

    log_data = np.log2(intensity_matrix.replace(0, np.nan))
    corr_matrix = log_data.corr(method='pearson')
    sns.heatmap(corr_matrix, annot=True, cmap='RdYlBu_r', vmin=0.8, vmax=1.0)

PCA

def perform_pca(intensity_matrix, sample_groups):
    """Principal component analysis for sample clustering."""
    import numpy as np
    import matplotlib.pyplot as plt
    from sklearn.decomposition import PCA

    log_data = np.log2(intensity_matrix.replace(0, np.nan))
    imputed = log_data.fillna(log_data.min().min())

    pca = PCA(n_components=2)
    pca_result = pca.fit_transform(imputed.T)

    plt.scatter(pca_result[:, 0], pca_result[:, 1], c=sample_groups)
    plt.xlabel(f'PC1 ({pca.explained_variance_ratio_[0]:.1%})')
    plt.ylabel(f'PC2 ({pca.explained_variance_ratio_[1]:.1%})')

Preprocessing

Filtering

def filter_proteins(intensity_matrix, metadata, min_valid=3):
    """Filter out low-confidence proteins."""
    valid_proteins = metadata['Peptides'] >= 2
    n_detected = (intensity_matrix > 0).sum(axis=1)
    valid_detection = n_detected >= min_valid
    is_contaminant = metadata['Protein IDs'].str.contains('CON__', na=False)
    is_reverse = metadata['Protein IDs'].str.contains('REV__', na=False)
    keep = valid_proteins & valid_detection & ~is_contaminant & ~is_reverse
    return intensity_matrix[keep], metadata[keep]

Missing Value Imputation

def impute_missing_values(intensity_matrix, method='MinProb'):
    """Impute missing protein intensities."""
    import numpy as np
    import pandas as pd

    if method == 'MinProb':
        min_val = intensity_matrix[intensity_matrix > 0].min().min()
        width, shift = 0.3, 1.8
        imputed = intensity_matrix.copy()
        missing_mask = imputed == 0
        n_missing = missing_mask.sum().sum()
        random_vals = np.random.normal(loc=min_val - shift, scale=width, size=n_missing)
        imputed.values[missing_mask.values] = random_vals
        return imputed
    elif method == 'KNN':
        from sklearn.impute import KNNImputer
        imputer = KNNImputer(n_neighbors=5)
        imputed = pd.DataFrame(
            imputer.fit_transform(intensity_matrix.replace(0, np.nan)),
            index=intensity_matrix.index, columns=intensity_matrix.columns
        )
        return imputed

Normalization

def normalize_intensities(intensity_matrix, method='median'):
    """Normalize protein intensities across samples."""
    if method == 'median':
        medians = intensity_matrix.median(axis=0)
        global_median = medians.median()
        norm_factors = global_median / medians
        return intensity_matrix * norm_factors
    elif method == 'quantile':
        from sklearn.preprocessing import quantile_transform
        import pandas as pd
        return pd.DataFrame(
            quantile_transform(intensity_matrix, axis=1),
            index=intensity_matrix.index, columns=intensity_matrix.columns
        )

Differential Expression

def differential_expression_limma(log2_intensities, group1_samples, group2_samples):
    """Perform differential expression using limma-like approach."""
    import numpy as np
    import pandas as pd
    from scipy import stats
    from statsmodels.stats.multitest import multipletests

    results = []
    for protein in log2_intensities.index:
        group1 = log2_intensities.loc[protein, group1_samples]
        group2 = log2_intensities.loc[protein, group2_samples]
        log2fc = group2.mean() - group1.mean()
        t_stat, p_value = stats.ttest_ind(group1, group2, equal_var=False)
        results.append({
            'protein': protein, 'log2FC': log2fc,
            'mean_group1': group1.mean(), 'mean_group2': group2.mean(),
            'p_value': p_value, 't_statistic': t_stat
        })

    results_df = pd.DataFrame(results)
    results_df['adj_p_value'] = multipletests(results_df['p_value'], method='fdr_bh')[1]
    results_df['significant'] = (
        (results_df['adj_p_value'] < 0.05) & (np.abs(results_df['log2FC']) > 1.0)
    )
    return results_df

Volcano Plot

def plot_volcano(de_results, title='Volcano Plot'):
    """Visualize differential expression results."""
    import matplotlib.pyplot as plt
    import numpy as np

    plt.figure(figsize=(8, 6))
    non_sig = de_results[~de_results['significant']]
    plt.scatter(non_sig['log2FC'], -np.log10(non_sig['p_value']), c='gray', alpha=0.5, s=10)
    sig = de_results[de_results['significant']]
    plt.scatter(sig['log2FC'], -np.log10(sig['p_value']), c='red', alpha=0.7, s=20)
    plt.axhline(-np.log10(0.05), color='blue', linestyle='--', label='p=0.05')
    plt.axvline(-1, color='blue', linestyle='--')
    plt.axvline(1, color='blue', linestyle='--', label='|log2FC|=1')
    plt.xlabel('log2 Fold Change')
    plt.ylabel('-log10(p-value)')
    plt.title(title)
    plt.legend()

PTM Analysis

def analyze_phosphosites(phospho_sites_file, intensity_matrix):
    """Analyze phosphorylation site changes from MaxQuant Phospho (STY)Sites.txt."""
    import pandas as pd

    phospho = pd.read_csv(phospho_sites_file, sep='\t')
    phospho_confident = phospho[phospho['Localization prob'] > 0.75]
    phospho_confident['site'] = (
        phospho_confident['Gene names'] + '_' +
        phospho_confident['Amino acid'] +
        phospho_confident['Position'].astype(str)
    )
    return phospho_confident

PPI Network

def build_protein_network(protein_list, confidence=0.7):
    """Build PPI network using STRING database."""
    import networkx as nx
    from tooluniverse import ToolUniverse
    tu = ToolUniverse()

    interactions = tu.run_one_function({
        "name": "string_get_interactions",
        "arguments": {
            "proteins": ",".join(protein_list),
            "species": 9606, "score_threshold": int(confidence * 1000)
        }
    })

    G = nx.Graph()
    for interaction in interactions['data']:
        G.add_edge(interaction['protein1'], interaction['protein2'], score=interaction['score'])
    return G

def detect_protein_modules(network_graph):
    """Identify tightly connected protein modules."""
    from networkx.algorithms import community
    communities = community.greedy_modularity_communities(network_graph)
    return [{'module_id': i, 'proteins': list(comm), 'size': len(comm)}
            for i, comm in enumerate(communities)]

Protein-RNA Correlation

def correlate_protein_rna(protein_data, rna_data, common_samples):
    """Correlate protein and mRNA levels. Expected r ~ 0.4-0.6."""
    from scipy.stats import spearmanr

    common_genes = set(protein_data.index) & set(rna_data.index)
    correlations = {}
    for gene in common_genes:
        protein = protein_data.loc[gene, common_samples]
        rna = rna_data.loc[gene, common_samples]
        r, p = spearmanr(protein, rna)
        correlations[gene] = {'r': r, 'p': p}
    return correlations

Proteomics Analysis - Phase Details

Phase 1: Data Import & QC

Supported input formats:

• MaxQuant: proteinGroups.txt, evidence.txt, Phospho (STY)Sites.txt
• Spectronaut: *_Report.tsv
• DIA-NN: report.tsv, report.pr_matrix.tsv
• Proteome Discoverer: *_Proteins.txt, *_PSMs.txt

QC checks: 1. Missing value assessment per protein and per sample 2. Log10 intensity distribution boxplots (expect similar median/spread) 3. Sample-sample Pearson correlation (high within replicates) 4. PCA for sample clustering (expect separation by condition)

Phase 2: Preprocessing & Normalization

1. Filter: Keep proteins with 2+ unique peptides and detected in min 3 samples. Remove contaminants (CON__) and reverse (REV__) sequences. 2. Impute: MinProb for MNAR assumption (random low values), KNN for MAR assumption. 3. Normalize: Median normalization (divide by sample median, multiply by global median), or quantile normalization.

Phase 3: Differential Expression

1. For each protein, calculate log2 fold change between conditions 2. Welch's t-test (unequal variance) per protein 3. BH multiple testing correction 4. Classify significant: adj_p < 0.05 AND |log2FC| > 1.0 5. Generate volcano plot

Phase 4: PTM Analysis

1. Load modification-specific peptides (e.g., Phospho (STY)Sites.txt) 2. Filter by localization probability > 0.75 3. Construct site identifiers: GENE_AminoAcidPosition (e.g., AKT1_S473) 4. Differential analysis same as Phase 3 5. Predict upstream kinases using kinase-substrate databases

Phase 5: Functional Enrichment

Run enrichment for significant proteins using ToolUniverse:

• Enrichr_enrich with libraries: KEGG_2021_Human, Reactome_2022, GO_Biological_Process_2021
• Protein complex enrichment against CORUM database
• Tissue-specific enrichment if relevant

Phase 6: Protein-Protein Interactions

1. Query STRING for interaction networks (confidence > 0.7) 2. Build network graph (nodes = proteins, edges = interactions) 3. Detect modules/communities using greedy modularity 4. Identify hub proteins (high degree/betweenness centrality) 5. Annotate modules with enriched functions

Phase 7: Multi-Omics Integration

1. Match protein and RNA data by gene name across common samples 2. Spearman correlation per gene (expected r ~ 0.4-0.6) 3. Classify regulation: transcriptional (r > 0.6), translational (high protein + low RNA, r < 0.2), degradation (low protein + high RNA) 4. Enrichment analysis on post-transcriptionally regulated genes

Phase 8: Report Generation

Report Structure

1. Dataset Summary: Platform, samples, conditions, proteins quantified 2. QC Results: Missing value heatmap, intensity distributions, PCA plot, replicate correlations 3. Differential Expression: Volcano plot, significant protein table (gene, log2FC, adj_p), MA plot 4. PTM Analysis (if applicable): Modified sites table, kinase activity predictions 5. Pathway Enrichment: Top enriched GO terms, KEGG/Reactome pathways with bar/dot plots 6. PPI Networks: STRING network visualization, hub proteins, module annotations 7. Multi-Omics (if applicable): Protein-RNA correlation scatter, translation-regulated genes 8. Conclusions: Key findings, limitations, suggested follow-up experiments

Report Template: Proteomics Analysis

# Proteomics Analysis Report

## Dataset Summary
- **Samples**: X (Y disease, Z control)
- **Proteins Identified**: N
- **Proteins Quantified**: N (at least 3 samples)
- **Platform**: [instrument], [software version]

## Quality Control
- **Missing Values**: X% average per protein
- **Sample Correlation**: X-X within groups
- **PCA**: [separation description] (PC1: X% variance)

## Differential Expression
- **Significant Proteins**: N (adj. p < 0.05, |log2FC| > 1)
  - Upregulated: N proteins
  - Downregulated: N proteins
- **Top upregulated**: [gene] (log2FC=X), [gene] (log2FC=X)
- **Top downregulated**: [gene] (log2FC=X), [gene] (log2FC=X)

## Phosphoproteomics
- **Phosphosites Quantified**: N
- **Differentially Phosphorylated**: N sites (p < 0.05)
- **Top Predicted Kinases**: [kinase1], [kinase2], [kinase3]

## Pathway Enrichment
### Top Pathways (Upregulated)
1. **[Pathway]** (p=X) - N proteins
2. **[Pathway]** (p=X) - N proteins

### Top Pathways (Downregulated)
1. **[Pathway]** (p=X) - N proteins
2. **[Pathway]** (p=X) - N proteins

## Protein Network Analysis
- **Network**: N nodes, N edges (STRING confidence > 0.7)
- **Modules Detected**: N functional modules
  - Module 1: [function] (N proteins)
  - Module 2: [function] (N proteins)

## Protein-RNA Correlation
- **Overall Correlation**: r = X (moderate, expected)
- **High Correlation**: N genes (r > 0.6) - transcriptional regulation
- **Low Correlation**: N genes (r < 0.2) - post-transcriptional regulation
- **Translation-Regulated**: N proteins (high protein, low RNA)

## Biological Interpretation
[Summary of key biological findings]

## Potential Biomarkers
Top proteins for classification (AUC=X):
1. [protein] (type)
2. [protein] (type)

Example Use Cases

Use Case 1: Cancer Proteomics

Question: "Analyze proteomics data from breast cancer vs normal tissue" Workflow: Load MaxQuant -> QC/filter -> Impute/normalize -> DE (432 sig) -> Pathway enrichment -> STRING network -> Integrate with RNA-seq -> Report with biomarkers

Use Case 2: Phosphoproteomics Signaling

Question: "What kinase signaling is activated in response to drug treatment?" Workflow: Load Phospho Sites -> Filter by localization -> Differential phosphorylation -> Kinase prediction -> MAPK/PI3K pathway enrichment -> Report

Use Case 3: Protein-RNA Integration

Question: "Which proteins are regulated post-transcriptionally?" Workflow: Load proteomics + RNA-seq -> Match samples -> Correlate per gene -> Classify low-correlation genes -> Enrichment for post-transcriptional regulators -> Report