maxquant-proteomics

MaxQuant + Perseus — 蛋白质组学分析流程

概述

MaxQuant 是用于无标记定量（LFQ）和 SILAC 蛋白质组学的社区标准软件。它会从原始 LC-MS/MS 文件中执行数据库检索、蛋白质分组和基于强度的定量，并生成 proteinGroups.txt 作为主要输出。下游统计分析——筛选、归一化、插补、差异丰度检验和可视化——使用 pandas、scipy 和 matplotlib/seaborn 在 Python 中完成，在可复现的脚本环境中复现 Perseus 工作流。

适用场景

• 对多个生物学条件下的蛋白质执行无标记定量（LFQ）——MaxQuant 的 MaxLFQ 算法是社区基准
• 运行采用轻/重或三标记设计的 SILAC（稳定同位素标记）实验
• 通过 MaxQuant 的报告离子定量处理 iTRAQ 或 TMT 等压标记实验
• 在需要广泛引用的 MaxQuant 输出格式（proteinGroups.txt）以便与已发表数据集比较时，鉴定和定量蛋白质
• 无需安装 Perseus（仅提供 Windows 图形界面）即可对 MaxQuant 输出执行差异丰度统计分析
• 在可复现的 Python 工作流中，基于蛋白质组学数据生成适用于发表的火山图和 GO 富集分析
• 处理需要仪器原生处理或 Sequest HT 的 Thermo 原始文件时，改用 Proteome Discoverer
• 需要 GPU 加速数据库检索（快 3–10×）或处理 DIA（数据非依赖采集）数据时，改用 FragPipe/MSFragger
• 差异丰度分析或 GSEA 后，使用 omics-plotting SKILL 生成适用于发表的图表

前提条件

• MaxQuant：Windows 软件；从 https://maxquant.org/ 下载（v2.4+）；需要 .NET 6 运行时
• Python 包：pandas、numpy、scipy、matplotlib、seaborn、statsmodels、gseapy
• 数据要求：Thermo .raw 文件或转换为 mzML 的文件；FASTA 蛋白质数据库（UniProt reviewed + 污染物数据库）
• 环境：MaxQuant 在 Windows 上运行（图形界面或命令行）；Python 分析可跨平台运行

pip install pandas numpy scipy matplotlib seaborn statsmodels gseapy

# Install pyMaxQuant for programmatic mqpar.xml configuration
pip install pymaxquant

预先访谈

编写任何分析代码前，先与用户确认以下事项。

decisions:
  - id: D1
    param: quantificationStrategy
    kind: required
    source: user
    ask: "How were samples quantified - label-free, SILAC, or isobaric tags?"
    default: null

  - id: D2
    param: proteinSequenceDatabase
    kind: required
    source: user
    ask: "Which organism's protein FASTA, and which release, should spectra be searched against?"
    default: null

  - id: D3
    param: digestionEnzyme
    kind: required
    source: user
    ask: "Which protease was used, and how many missed cleavages should be allowed?"
    default: "trypsin, up to 2 missed cleavages"

  - id: D4
    param: modifications
    kind: required
    source: user
    ask: "Which modifications are fixed, and which variable - and is there an enrichment such as phospho to search for?"
    default: "fixed carbamidomethyl on cysteine; variable oxidation and N-terminal acetylation"

  - id: D5
    param: falseDiscoveryRates
    kind: required
    source: user
    ask: "What false-discovery rate should peptide and protein identifications be held to?"
    default: "1% at both levels"

  - id: D6
    param: matchBetweenRuns
    kind: required
    source: user
    depends_on: [D1]
    ask: "Should identifications be transferred between runs by retention time, recovering more proteins at the cost of some transferred errors?"
    default: "off"

  - id: D7
    param: lfqMinRatioCount
    kind: optional
    source: user
    depends_on: [D1]
    ask: "How many shared peptides must two samples have before their abundances are normalized against each other?"
    default: 2
    skip_if: "labelled quantification"

  - id: D8
    param: missingValueImputation
    kind: required
    source: user
    ask: "How should proteins missing from a sample be treated - left missing, or imputed as below the detection limit?"
    default: "imputed from a downshifted distribution"

  - id: D9
    param: differentialThresholds
    kind: required
    source: user
    ask: "What significance and fold-change cuts define a changed protein?"
    default: "FDR 0.05, log2 fold change 1.0"

D8 是蛋白质组学区别于 RNA-seq 的地方，也是结果最常因未经审视的默认设置而发生变化的地方。某个蛋白质在所有对照样本中都未检出、却在所有处理样本中都存在，可能是最有力的结果，也可能是检测限造成的伪象；插补方式决定了火山图会呈现哪一种情况。D6 与之相互影响：转移的鉴定结果恰好会填补原本需要插补的数值。

快速开始

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()}")

工作流

步骤 1：通过 mqpar.xml 配置 MaxQuant 参数

MaxQuant 通过 XML 参数文件（mqpar.xml）控制。运行检索前，以编程方式编辑该文件，设置文件路径、酶、修饰和定量类型。

import xml.etree.ElementTree as ET

def update_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"],
)

mqpar.xml 的关键参数（在模板中设置或直接编辑）：

<enzymes>
  <string>Trypsin/P</string>
</enzymes>
<maxMissedCleavages>2</maxMissedCleavages>
<variableModifications>
  <string>Oxidation (M)</string>
  <string>Acetyl (Protein N-term)</string>
</variableModifications>
<fixedModifications>
  <string>Carbamidomethyl (C)</string>
</fixedModifications>


<lfqMode>1</lfqMode>                     
<lfqMinRatioCount>2</lfqMinRatioCount>   
<matchBetweenRuns>True</matchBetweenRuns>


<peptideFdr>0.01</peptideFdr>
<proteinFdr>0.01</proteinFdr>

步骤 2：从命令行运行 MaxQuant（Windows）

MaxQuant 可通过捆绑的 MaxQuantCmd.exe 从 Windows 命令提示符以无界面方式运行。

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 log
tail -f combined/proc/#runningTimes.txt

步骤 3：加载并筛选 proteinGroups.txt

筛除反向诱饵、潜在污染物以及仅通过修饰位点鉴定的蛋白质。

import pandas as pd
import numpy as np

def load_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', ...]

步骤 4：对 LFQ 强度进行 Log2 转换和中位数归一化

将强度为零的值替换为 NaN（MaxQuant 将缺失值导出为 0），进行 log2 转换，然后按样本执行中位数中心化。

def prepare_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, ...

步骤 5：插补缺失值（MNAR 策略）

非随机缺失（MNAR）值来自低于检测限的蛋白质。从观测到的强度分布低端进行插补——这是标准的 Perseus 方法。

def impute_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

第 6 步：统计检验——带 FDR 校正的 t 检验

对不同条件下的每种蛋白质执行双样本 t 检验，然后进行 Benjamini-Hochberg FDR 校正。

from scipy import stats
from statsmodels.stats.multitest import multipletests

def differential_abundance(lfq: pd.DataFrame,
                            group_a: list[str],
                            group_b: list[str],
                            alpha: float = 0.05) -> pd.DataFrame:
    """
    Two-sample t-test + BH FDR correction for all proteins.

    Parameters
    ----------
    group_a, group_b : sample name lists for each condition
    alpha : FDR threshold
    """
    results = []
    for protein_id, row in lfq.iterrows():
        a_vals = row[group_a].dropna().values
        b_vals = row[group_b].dropna().values

        if len(a_vals) >= 2 and len(b_vals) >= 2:
            t_stat, p_val = stats.ttest_ind(a_vals, b_vals, equal_var=False)
            log2fc = b_vals.mean() - a_vals.mean()
        else:
            t_stat, p_val, log2fc = np.nan, np.nan, np.nan

        results.append({
            "protein_id": protein_id,
            "log2FC": log2fc,
            "pvalue": p_val,
            "t_stat": t_stat,
        })

    res_df = pd.DataFrame(results).set_index("protein_id")

    # BH FDR correction on valid p-values
    valid = res_df["pvalue"].notna()
    _, padj, _, _ = multipletests(res_df.loc[valid, "pvalue"], method="fdr_bh")
    res_df.loc[valid, "padj"] = padj

    # Add significance flag
    res_df["significant"] = (res_df["padj"] < alpha) & (res_df["pvalue"].notna())
    sig_count = res_df["significant"].sum()
    print(f"Significant proteins (FDR < {alpha}): {sig_count}")

    return res_df.sort_values("padj")

# Define sample groups
group_ctrl  = ["ctrl_1",  "ctrl_2",  "ctrl_3"]
group_treat = ["treat_1", "treat_2", "treat_3"]

results = differential_abundance(lfq_imputed, group_ctrl, group_treat)
print(results[results["significant"]].head(10))
# Significant proteins (FDR < 0.05): 312

第 7 步：火山图可视化

绘图：计算 DEG 表，然后阅读 skills/data-visualization/omics-plotting/SKILL.md 并遵循其中的「火山图」配方，处理导出的 CSV（→ figures/volcano_plot.png）。

# Volcano input: log2FC vs padj + protein labels; render with the omics-plotting SKILL (`skills/data-visualization/omics-plotting/SKILL.md`) "Volcano" recipe.
volcano_df = results[["log2FC", "padj"]].copy()
volcano_df["gene"] = pg["Gene names"].reindex(results.index)
volcano_df.dropna(subset=["log2FC", "padj"]).to_csv("volcano_input.csv")
print(f"Volcano input: {int(volcano_df['padj'].notna().sum())} proteins -> volcano_input.csv")
# -> figures/volcano_plot.png via omics-plotting (map log2FC -> log2FoldChange)

第 8 步：显著蛋白质的 GO/通路富集分析

使用 gseapy 的 Enrichr API，对显著上调和下调的蛋白质执行过度代表分析（ORA）。

import gseapy as gp

def run_enrichment(results: pd.DataFrame,
                   gene_names: pd.Series,
                   gene_sets: list[str] | None = None,
                   top_n: int = 20) -> dict[str, pd.DataFrame]:
    """
    ORA enrichment for up- and down-regulated proteins via Enrichr.

    Parameters
    ----------
    gene_sets : Enrichr gene set libraries (default: GO BP + KEGG)
    top_n     : top results to display per direction
    """
    if gene_sets is None:
        gene_sets = ["GO_Biological_Process_2023", "KEGG_2021_Human"]

    results_out = {}
    gene_map = gene_names.reindex(results.index).fillna("Unknown")

    for direction in ("up", "down"):
        sig = results[
            (results["significant"]) &
            (results["log2FC"] > 0 if direction == "up" else results["log2FC"] < 0)
        ]
        gene_list = gene_map.reindex(sig.index).tolist()
        print(f"{direction.capitalize()}-regulated: {len(gene_list)} proteins")

        if len(gene_list) < 5:
            print(f"  Too few proteins for enrichment (n={len(gene_list)}), skipping")
            continue

        enr = gp.enrichr(
            gene_list=gene_list,
            gene_sets=gene_sets,
            organism="Human",
            outdir=f"enrichment_{direction}",
            cutoff=0.05,
        )
        top_results = enr.results.sort_values("Adjusted P-value").head(top_n)
        results_out[direction] = top_results
        print(top_results[["Term", "Adjusted P-value", "Overlap"]].to_string(index=False))

    return results_out

enr_results = run_enrichment(results, pg["Gene names"])

关键参数

参数	默认值	范围 / 选项	作用
matchBetweenRuns	False	True / False	通过保留时间匹配，在不同运行之间转移鉴定结果；使定量蛋白质数量增加 10–30%
lfqMinRatioCount	2	1–5	LFQ 归一化所需的最少肽段对数量；较低的值可提高覆盖率，但会降低准确性
maxMissedCleavages	2	0–4	允许的胰蛋白酶漏切数；消化效果较差的样本可调高此值
peptideFdr / proteinFdr	0.01	0.001–0.05	肽段和蛋白质鉴定的 FDR 阈值
MNAR downshift	1.8	1.5–2.5	以列标准差为单位，将插补分布下移至检测限以下；值越大，插补越保守
MNAR width	0.3	0.1–0.5	插补分布相对于列标准差的宽度
t 检验 alpha	0.05	0.01–0.1	差异丰度的 FDR 显著性阈值
火山图 fc_threshold	1.0	0.5–2.0	火山图中标记为「显著」的 log2 倍数变化截断值

关键概念

MaxQuant 输出文件

文件	内容
proteinGroups.txt	主要输出：每个蛋白质组一行，包含 LFQ/SILAC 强度、肽段计数和序列覆盖率
peptides.txt	肽段级定量，包含电荷态和修饰
evidence.txt	单条 MS/MS 鉴定结果（每个肽段-谱图匹配一行）
msms.txt	完整的 MS/MS 扫描数据，包含碎片离子和评分
summary.txt	每个原始文件的统计信息：鉴定结果、MS/MS 计数和校准信息

LFQ 强度与 iBAQ

• LFQ（非标记定量）：MaxLFQ 算法基于 razor 肽段和唯一肽段的比值，对不同样本的强度进行归一化。用于跨样本比较（倍数变化）。存储在 LFQ intensity <sample> 列中。
• iBAQ（基于强度的绝对定量）：将肽段强度总和除以理论上可观测的肽段数量。用于估算拷贝数，以及比较同一样本中不同蛋白质的绝对丰度。存储在 iBAQ 列中。
• SILAC 比值：来自同位素标记配对的直接 H/L 比值。对于较小的倍数变化，其准确性高于 LFQ。

Perseus 在 Python 中的等效操作

Perseus 步骤	Python 等效操作
按分类列筛选行	df[df["Reverse"] != "+"]
将 0 替换为 NaN	df.replace(0, np.nan)
Log2 转换	np.log2(df)
中位数归一化	df.subtract(df.median()).add(global_median)
MNAR 插补（正态分布）	上面的 impute_mnar() 函数
双样本 t 检验	scipy.stats.ttest_ind() + multipletests()
火山图	使用 matplotlib.pyplot 散点图和阈值线
层次聚类	seaborn.clustermap()

常用配方

配方：SILAC 比值分析

适用场景：使用 H/L 或 H/M/L 标记而非 LFQ 的 SILAC 实验。

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()]
if not 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 sample
print(log2_ratios.describe().round(3))

配方：层次聚类热图

适用场景：同时可视化所有显著蛋白质的模式。计算下面的 z 分数矩阵，然后阅读 skills/data-visualization/omics-plotting/SKILL.md 并遵循其中的「聚类表达热图」配方（→ figures/heatmap.pdf）。

import pandas as pd

def prep_heatmap_matrix(lfq_imputed: pd.DataFrame,
                        results: pd.DataFrame,
                        gene_names: pd.Series,
                        top_n: int = 50,
                        out_path: str = "heatmap_matrix.csv") -> pd.DataFrame:
    """Z-scored expression matrix of the top significant proteins, for the omics-plotting heatmap."""
    sig_proteins = results[results["significant"]].nsmallest(top_n, "padj").index
    heatmap_data = lfq_imputed.loc[sig_proteins].copy()
    heatmap_data.index = gene_names.reindex(sig_proteins).fillna(sig_proteins)
    # Z-score per row
    heatmap_z = heatmap_data.subtract(heatmap_data.mean(axis=1), axis=0).divide(
        heatmap_data.std(axis=1).replace(0, 1), axis=0
    )
    heatmap_z.to_csv(out_path)
    return heatmap_z

mat = prep_heatmap_matrix(lfq_imputed, results, pg["Gene names"])
print(f"Heatmap matrix: {mat.shape} -> heatmap_matrix.csv")
# Render with the omics-plotting SKILL (`skills/data-visualization/omics-plotting/SKILL.md`) "Clustered expression heatmap" recipe -> figures/heatmap.pdf

配方：对显著蛋白进行 STRING-db 网络富集分析

适用场景：无需在本地下载基因集，即可进行蛋白质-蛋白质相互作用网络分析和富集分析。

import requests
import pandas as pd

def string_enrichment(gene_list: list[str],
                      species: int = 9606,
                      fdr_threshold: float = 0.05) -> pd.DataFrame:
    """Query STRING /enrichment endpoint for GO/KEGG enrichment."""
    url = "https://string-db.org/api/json/enrichment"
    params = {
        "identifiers": "\r".join(gene_list),
        "species": species,
        "caller_identity": "maxquant_proteomics_skill",
    }
    response = requests.post(url, data=params)
    response.raise_for_status()

    enr_df = pd.DataFrame(response.json())
    if enr_df.empty:
        print("No enrichment results returned")
        return enr_df

    enr_df = enr_df[enr_df["fdr"].astype(float) < fdr_threshold]
    enr_df = enr_df.sort_values("fdr")
    print(f"Enriched terms (FDR < {fdr_threshold}): {len(enr_df)}")
    print(enr_df[["category", "term", "description", "fdr", "number_of_genes"]].head(15).to_string(index=False))
    return enr_df

# Significant up-regulated gene names
up_genes = pg.loc[
    results[(results["significant"]) & (results["log2FC"] > 1)].index, "Gene names"
].tolist()

string_enr = string_enrichment(up_genes)

配方：使用 pyMaxQuant 解析 MaxQuant 输出

适用场景：使用更高级别的 API 读取和筛选 MaxQuant 文本文件。

# pyMaxQuant provides typed accessors for MaxQuant output files
# Install: pip install pymaxquant
from maxquant.io import read_protein_groups

# Load with built-in contaminant filtering
pg_clean = read_protein_groups(
    "combined/txt/proteinGroups.txt",
    filter_invalid=True,      # removes reverse, contaminant, only-by-site
)
print(f"Loaded {len(pg_clean)} filtered protein groups")

# Access LFQ columns via helper
lfq_df = pg_clean.filter(like="LFQ intensity")
print(f"LFQ matrix: {lfq_df.shape}")

预期输出

文件	描述
combined/txt/proteinGroups.txt	MaxQuant 主要输出：包含 LFQ 强度、肽段计数、唯一肽段和 iBAQ 的蛋白组
combined/txt/peptides.txt	肽段水平的定量结果，包含修饰和电荷状态
combined/txt/summary.txt	每个原始文件的 QC 统计：鉴定率、MS/MS 计数
results_differential.csv	差异丰度表：每种蛋白对应的 log2FC、pvalue、padj、significant
volcano_plot.pdf	火山图：用颜色区分上调/下调蛋白，并标注排名靠前的蛋白
heatmap.pdf	显著性最高的蛋白的层次聚类热图（经 Z-score 归一化）
enrichment_up/	gseapy 输出目录：上调蛋白的 GO/KEGG 富集结果
enrichment_down/	gseapy 输出目录：下调蛋白的 GO/KEGG 富集结果

故障排除

问题	原因	解决方案
MaxQuant 未鉴定出任何蛋白	FASTA 数据库或酶设置错误；未找到原始文件路径	检查 mqpar.xml 中的 .raw 文件路径是否为绝对 Windows 路径；确认酶设置与实验相符（Trypsin/P 或 Trypsin）；查看 summary.txt 中的鉴定率
筛选后所有 LFQ 强度均为 0	matchBetweenRuns 未启用且数据稀疏，或列选择错误	直接检查 combined/txt/proteinGroups.txt；使用 pg.filter(like="LFQ intensity") 确认列名；将 lfqMinRatioCount 降至 1
log2 转换后缺失值过多	重复样本不足、样本上样量不一致，或肽段未被检测到	启用 matchBetweenRuns；确认蛋白上样量相等（Bradford/BCA）；在填补前考虑采用更严格的有效值筛选条件（要求每组 3/3 个有效值）
加载 proteinGroups.txt 时发生内存错误	文件较大（DDA 且样本较多时超过 500 MB）	使用 pd.read_csv(..., low_memory=False, usecols=[...]) 仅选取所需列；或使用 pd.read_csv(..., chunksize=...)
gseapy Enrichr 返回空结果	无法识别基因符号或网络超时	确保基因列表使用 HGNC 符号（而非 UniProt ID）；检查网络连接；使用 gp.enrichr(..., timeout=60)
火山图中所有蛋白都标为“ns”	FDR 阈值过严或未计算 padj	确认 multipletests 返回了有效的 FDR 值；尝试将 alpha 放宽至 0.1；检查样本组别分配是否正确
MaxQuant 运行卡在“Feature detection”阶段	内存不足（每 3–4 个原始文件需要 4–8 GB RAM）	分批处理文件；增加系统 RAM；关闭其他应用程序
填补导致假阳性增多	填补过于激进（downshift 值过低）	将 downshift 增至 2.0–2.5；也可以在检验前筛选出每组有 ≥ 2 个有效值的蛋白

参考资料

• MaxQuant 软件和文档 — MaxQuant 官方下载、手册和参数指南
• Perseus 计算平台 — MaxQuant 配套的统计分析工具；上述 Python 工作流的基础
• Cox et al. (2008) Nature Biotechnology — MaxQuant 论文 — 介绍 MaxLFQ 算法的 MaxQuant 原始论文
• Tyanova et al. (2016) Nature Methods — Perseus 论文 — Perseus 统计框架参考文献；上述所有工作流步骤均与 Perseus 操作对应
• pyMaxQuant GitHub — 用于解析 MaxQuant 输出文件的 Python 库
• gseapy 文档 — 第 8 步中使用的基因集富集分析库
• STRING API 富集分析 — 通用方案中使用的蛋白质网络富集分析接口