MaxQuant + Perseus proteomics pipeline: run MaxQuant for LFQ and SILAC; parse proteinGroups.txt in Python; filter contaminants/decoys; log2 + median-normalize; impute MNAR; t-test with FDR; volcano plot; GO/pathway enrichment. Use Proteome Discoverer for Thermo-native processing; FragPipe/MSFragger for GPU-accelerated DB search.
MaxQuant + Perseus proteomics pipeline: run MaxQuant for LFQ and SILAC; parse proteinGroups.txt in Python; filter contaminants/decoys; log2 + median-normalize; impute MNAR; t-test with FDR; volcano plot; GO/pathway enrichment. Use Proteome Discoverer for Thermo-native processing; FragPipe/MSFragger for GPU-accelerated DB search.
license
Apache-2.0
MaxQuant + Perseus — Proteomics Analysis Pipeline
Overview
MaxQuant is the community-standard software for label-free quantification (LFQ) and SILAC proteomics. It performs database search, protein grouping, and intensity-based quantification from raw LC-MS/MS files, producing proteinGroups.txt as the primary output. Downstream statistical analysis — filtering, normalization, imputation, differential abundance testing, and visualization — is performed in Python using pandas, scipy, and matplotlib/seaborn, mirroring the Perseus workflow in a reproducible scripting environment.
When to Use
Performing label-free quantification (LFQ) of proteins across multiple biological conditions — MaxQuant's MaxLFQ algorithm is the community benchmark
Running SILAC (stable isotope labeling) experiments with light/heavy or triple-label designs
Processing iTRAQ or TMT isobaric labeling experiments via MaxQuant's reporter ion quantification
Identifying and quantifying proteins when you need the widely-cited MaxQuant output format (proteinGroups.txt) for comparison with published datasets
Performing statistical differential abundance analysis on MaxQuant outputs without installing Perseus (GUI-only, Windows)
Generating publication-quality volcano plots and GO enrichment from proteomics data in a reproducible Python workflow
Use instead when working with Thermo raw files requiring instrument-native processing or Sequest HT
Proteome Discoverer
Use FragPipe/MSFragger instead for GPU-accelerated database search (3–10× faster) or when processing DIA (data-independent acquisition) data
Prerequisites
MaxQuant: Windows software; download from https://maxquant.org/ (v2.4+); requires .NET 6 runtime
# Install pyMaxQuant for programmatic mqpar.xml configuration
pip install pymaxquant
Quick Start
import pandas as pd
import numpy as np
# Load MaxQuant output
df = pd.read_csv("combined/txt/proteinGroups.txt", sep="\t", low_memory=False)
print(f"Raw protein groups: {len(df)}")
# Filter contaminants, reverse decoys, only-by-site
mask = (
(df["Potential contaminant"] != "+") &
(df["Reverse"] != "+") &
(df["Only identified by site"] != "+")
)
df = df[mask].copy()
print(f"After filtering: {len(df)} protein groups")
# Extract LFQ intensity columns
lfq_cols = [c for c in df.columns if c.startswith("LFQ intensity ")]
print(f"LFQ columns: {lfq_cols}")
# Log2-transform (0 → NaN)
lfq = df[lfq_cols].replace(0, np.nan)
lfq = np.log2(lfq)
print(f"Valid values per sample:\n{lfq.notna().sum()}")
Workflow
Step 1: Configure MaxQuant Parameters via mqpar.xml
MaxQuant is controlled by an XML parameter file (mqpar.xml). Edit it programmatically to set file paths, enzyme, modifications, and quantification type before running the search.
import xml.etree.ElementTree as ET
defupdate_mqpar(template_path: str, output_path: str,
raw_files: list[str], fasta_path: str,
experiment_names: list[str]) -> None:
"""Update mqpar.xml with sample-specific file paths."""
tree = ET.parse(template_path)
root = tree.getroot()
# Set raw file paths
file_paths_node = root.find(".//filePaths")
file_paths_node.clear()
for rf in raw_files:
elem = ET.SubElement(file_paths_node, "string")
elem.text = rf
# Set experiment names (maps files to conditions)
experiments_node = root.find(".//experiments")
experiments_node.clear()
for name in experiment_names:
elem = ET.SubElement(experiments_node, "string")
elem.text = name
# Set FASTA database
fasta_node = root.find(".//fastaFiles/FastaFileInfo/fastaFilePath")
fasta_node.text = fasta_path
tree.write(output_path, xml_declaration=True, encoding="utf-8")
print(f"Written: {output_path}")
# Example usage
raw_files = [
r"C:\Data\ctrl_rep1.raw",
r"C:\Data\ctrl_rep2.raw",
r"C:\Data\treat_rep1.raw",
r"C:\Data\treat_rep2.raw",
]
update_mqpar(
template_path="mqpar_template.xml",
output_path="mqpar.xml",
raw_files=raw_files,
fasta_path=r"C:\Databases\human_uniprot_contaminants.fasta",
experiment_names=["ctrl", "ctrl", "treat", "treat"],
)
Key mqpar.xml parameters (set in template or edit directly):
MaxQuant can be run headlessly from the Windows command prompt using the bundled MaxQuantCmd.exe.
REM Windows Command Prompt — run MaxQuant with configured mqpar.xml
REM Adjust path to match your MaxQuant installation directory
set MQ_PATH=C:\Program Files\MaxQuant\bin\MaxQuantCmd.exe
set MQPAR=C:\Projects\proteomics\mqpar.xml
"%MQ_PATH%" "%MQPAR%"
REM For specific workflow steps only (useful for reruns):
REM Step IDs: 0=write tables, 1=feature detection, 7=peptide identification
"%MQ_PATH%" "%MQPAR%" --steps 1,7,11
# Cross-platform: run MaxQuant under Wine on Linux/macOS (CI/server use)
wine MaxQuantCmd.exe mqpar.xml
# Monitor progress logtail -f combined/proc/#runningTimes.txt
Step 3: Load and Filter proteinGroups.txt
Filter out reverse decoys, potential contaminants, and proteins only identified by modification site.
import pandas as pd
import numpy as np
defload_protein_groups(path: str) -> pd.DataFrame:
"""Load MaxQuant proteinGroups.txt with quality filters applied."""
df = pd.read_csv(path, sep="\t", low_memory=False)
print(f"Total protein groups: {len(df)}")
# Remove reverse decoys, contaminants, and only-by-site hits
n_before = len(df)
df = df[
(df.get("Reverse", pd.Series("")) != "+") &
(df.get("Potential contaminant", pd.Series("")) != "+") &
(df.get("Only identified by site", pd.Series("")) != "+")
].copy()
print(f"After quality filter: {len(df)} ({n_before - len(df)} removed)")
# Parse gene names (take first entry for multi-gene groups)
df["Gene names"] = df["Gene names"].fillna("Unknown").str.split(";").str[0]
# Set unique index on majority protein ID
df = df.set_index("Majority protein IDs")
return df
# Load output
pg = load_protein_groups("combined/txt/proteinGroups.txt")
# Identify LFQ intensity columns
lfq_cols = [c for c in pg.columns if c.startswith("LFQ intensity ")]
print(f"LFQ samples ({len(lfq_cols)}): {lfq_cols}")
# Output: LFQ samples (6): ['LFQ intensity ctrl_1', 'LFQ intensity ctrl_2', ...]
Step 4: Log2 Transform and Median Normalize LFQ Intensities
Replace zero intensities with NaN (missing values in MaxQuant are exported as 0), log2-transform, then apply per-sample median centering.
defprepare_lfq_matrix(df: pd.DataFrame, lfq_cols: list[str]) -> pd.DataFrame:
"""Extract, transform, and normalize LFQ intensity matrix."""# Extract and rename columns (strip 'LFQ intensity ' prefix)
lfq = df[lfq_cols].copy()
lfq.columns = [c.replace("LFQ intensity ", "") for c in lfq_cols]
# Replace 0 with NaN (MaxQuant encodes missing as 0)
lfq = lfq.replace(0, np.nan)
# Log2 transform
lfq = np.log2(lfq)
# Median centering per sample (subtract per-column median of valid values)
col_medians = lfq.median(axis=0)
global_median = col_medians.median()
lfq = lfq.subtract(col_medians, axis=1).add(global_median)
print(f"Matrix shape: {lfq.shape}")
print(f"Missing values per sample:\n{lfq.isna().sum()}")
print(f"Valid values per sample:\n{lfq.notna().sum()}")
return lfq
lfq_matrix = prepare_lfq_matrix(pg, lfq_cols)
# Matrix shape: (3241, 6)# Missing values per sample: ctrl_1: 421, ctrl_2: 389, ...
Step 5: Impute Missing Values (MNAR Strategy)
Missing-not-at-random (MNAR) values arise from proteins below the detection limit. Impute from the low end of the observed intensity distribution — the standard Perseus approach.
defimpute_mnar(lfq: pd.DataFrame,
width: float = 0.3,
downshift: float = 1.8,
random_state: int = 42) -> pd.DataFrame:
"""
Impute MNAR missing values from a downshifted Gaussian.
Parameters
----------
width : std of imputation distribution (fraction of sample std)
downshift : downshift in units of sample std below mean
random_state : for reproducibility
"""
rng = np.random.default_rng(random_state)
lfq_imp = lfq.copy()
for col in lfq_imp.columns:
col_data = lfq_imp[col].dropna()
col_mean = col_data.mean()
col_std = col_data.std()
n_missing = lfq_imp[col].isna().sum()
if n_missing > 0:
imputed = rng.normal(
loc=col_mean - downshift * col_std,
scale=width * col_std,
size=n_missing,
)
lfq_imp.loc[lfq_imp[col].isna(), col] = imputed
print(f"Imputed {lfq.isna().sum().sum()} missing values")
return lfq_imp
lfq_imputed = impute_mnar(lfq_matrix)
# Imputed 2847 missing values
Step 6: Statistical Testing — t-test with FDR Correction
Perform two-sample t-tests for each protein between conditions, then apply Benjamini-Hochberg FDR correction.
LFQ (Label-Free Quantification): MaxLFQ algorithm normalizes intensities across samples based on razor+unique peptide ratios. Use for cross-sample comparisons (fold changes). Stored in LFQ intensity <sample> columns.
iBAQ (intensity-Based Absolute Quantification): Divides summed peptide intensities by the number of theoretically observable peptides. Use for estimating copy numbers and comparing absolute abundance between proteins within a sample. Stored in iBAQ column.
SILAC ratio: Direct H/L ratio from isotope-labeled pairs. More accurate than LFQ for small fold changes.
Perseus Equivalent Operations in Python
Perseus step
Python equivalent
Filter rows by categorical column
df[df["Reverse"] != "+"]
Replace 0 with NaN
df.replace(0, np.nan)
Log2 transform
np.log2(df)
Median normalization
df.subtract(df.median()).add(global_median)
MNAR imputation (normal distribution)
impute_mnar() function above
Two-sample t-test
scipy.stats.ttest_ind() + multipletests()
Volcano plot
matplotlib.pyplot scatter + threshold lines
Hierarchical clustering
seaborn.clustermap()
Common Recipes
Recipe: SILAC Ratio Analysis
When to use: SILAC experiments with H/L or H/M/L labeling instead of LFQ.
import pandas as pd
import numpy as np
# Load proteinGroups.txt for SILAC experiment
df = pd.read_csv("combined/txt/proteinGroups.txt", sep="\t", low_memory=False)
# Filter contaminants and decoys
df = df[(df["Reverse"] != "+") & (df["Potential contaminant"] != "+")].copy()
# Extract H/L ratio columns (log2-transformed)
ratio_cols = [c for c in df.columns if c.startswith("Ratio H/L ") and"normalized"in c.lower()]
ifnot ratio_cols:
# Fall back to non-normalized
ratio_cols = [c for c in df.columns if c.startswith("Ratio H/L")]
print(f"SILAC ratio columns: {ratio_cols}")
ratios = df[ratio_cols].copy().replace(0, np.nan)
# Log2 transform ratios
log2_ratios = np.log2(ratios)
log2_ratios.columns = [c.replace("Ratio H/L normalized ", "") for c in ratio_cols]
# Summary statistics per sampleprint(log2_ratios.describe().round(3))
Recipe: Hierarchical Clustering Heatmap
When to use: visualizing patterns across all significant proteins simultaneously.
Differential abundance table: log2FC, pvalue, padj, significant per protein
volcano_plot.pdf
Volcano plot with up/down-regulated proteins colored and top proteins labeled
heatmap.pdf
Hierarchical clustering heatmap of top significant proteins (Z-score normalized)
enrichment_up/
gseapy output directory: GO/KEGG enrichment for up-regulated proteins
enrichment_down/
gseapy output directory: GO/KEGG enrichment for down-regulated proteins
Troubleshooting
Problem
Cause
Solution
MaxQuant produces 0 protein identifications
Wrong FASTA database or enzyme settings; raw file path not found
Verify .raw file paths in mqpar.xml are absolute Windows paths; confirm enzyme matches experiment (Trypsin/P vs Trypsin); check summary.txt for identification rate
All LFQ intensities are 0 after filtering
matchBetweenRuns off + sparse data, or wrong column selection
Check combined/txt/proteinGroups.txt directly; use pg.filter(like="LFQ intensity") to confirm column names; lower lfqMinRatioCount to 1
Too many missing values after log2 transform
Insufficient replicates, inconsistent sample loading, or undetected peptides
Enable matchBetweenRuns; verify equal protein loading (Bradford/BCA); consider stricter valid-value filter (require 3/3 per group) before imputation
Memory error loading proteinGroups.txt
File is large (>500 MB for DDA with many samples)
Use pd.read_csv(..., low_memory=False, usecols=[...]) to select only needed columns; or use pd.read_csv(..., chunksize=...)
gseapy Enrichr returns empty results
Gene symbols unrecognized or network timeout
Ensure gene list uses HGNC symbols (not UniProt IDs); check internet connectivity; use gp.enrichr(..., timeout=60)
Volcano plot: all proteins in "ns"
FDR threshold too stringent or padj not calculated
Verify multipletests returned valid FDR values; try relaxing alpha to 0.1; check sample group assignments are correct
MaxQuant run hangs at "Feature detection"
Low memory (MaxQuant needs 4–8 GB RAM per 3–4 raw files)
Process files in smaller batches; increase system RAM; close other applications
Imputation inflates false positives
Imputing too aggressively (low downshift)
Increase downshift to 2.0–2.5; alternatively, filter to proteins with ≥ 2 valid values per group before testing
