maxquant-proteomics

MaxQuant + Perseus — プロテオミクス解析パイプライン

概要

MaxQuantは、ラベルフリー定量（LFQ）とSILACプロテオミクスのコミュニティ標準ソフトウェアです。raw LC-MS/MSファイルからデータベース検索、タンパク質グループ化、強度ベースの定量を実行し、主要な出力としてproteinGroups.txtを生成します。その後の統計解析（フィルタリング、正規化、代入、差次的存在量検定、可視化）は、再現可能なスクリプト環境でPerseusのワークフローを再現しながら、pandas、scipy、matplotlib/seabornを使用してPythonで実行します。

使用する場面

• 複数の生物学的条件間でタンパク質のラベルフリー定量（LFQ）を行う場合。MaxQuantのMaxLFQアルゴリズムはコミュニティのベンチマークです
• light/heavyまたはトリプルラベル設計でSILAC（安定同位体標識）実験を行う場合
• MaxQuantのレポーターイオン定量を使用して、iTRAQまたはTMT同重体標識実験を処理する場合
• 公開データセットとの比較のために、広く引用されているMaxQuantの出力形式（proteinGroups.txt）が必要な場合に、タンパク質を同定・定量する場合
• Perseus（GUIのみ、Windows）をインストールせずに、MaxQuantの出力に対する統計的な差次的存在量解析を行う場合
• 再現可能なPythonワークフローで、プロテオミクスデータから出版品質のボルケーノプロットとGOエンリッチメントを生成する場合
• Thermoのrawファイルに対して、機器ネイティブの処理またはSequest HTが必要な場合は、代わりに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 + contaminant database）
• 環境：MaxQuantはWindowsで実行します（GUIまたはCLI）。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）

同梱のMaxQuantCmd.exeを使用して、WindowsコマンドプロンプトからMaxQuantをヘッドレスで実行できます。

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 検定

条件間で各タンパク質について 2 標本 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 テーブルを計算してから、エクスポートした CSV に対して skills/data-visualization/omics-plotting/SKILL.md を読み、「Volcano」の手順に従ってください（→ 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	主な出力：タンパク質グループごとに 1 行。LFQ／SILAC 強度、ペプチド数、配列カバレッジを含みます
peptides.txt	電荷状態と修飾を含む、ペプチドレベルの定量データ
evidence.txt	個々の MS/MS 同定結果（ペプチド・スペクトル一致ごとに 1 行）
msms.txt	フラグメントイオンとスコアを含む、MS/MS スキャンの全データ
summary.txt	raw ファイルごとの統計情報：同定数、MS/MS 数、キャリブレーション

LFQ 強度と iBAQ

• LFQ（ラベルフリー定量）：MaxLFQ アルゴリズムは、razor ペプチドとユニークペプチドの比に基づき、サンプル間の強度を正規化します。サンプル間の比較（フォールドチェンジ）に使用します。LFQ intensity <sample> 列に格納されます。
• iBAQ（強度ベース絶対定量）：合計ペプチド強度を、理論上観測可能なペプチド数で割ります。コピー数の推定や、サンプル内でのタンパク質間の絶対存在量の比較に使用します。iBAQ 列に格納されます。
• SILAC 比：同位体標識ペアから直接求める H/L 比です。小さなフォールドチェンジでは LFQ よりも正確です。

Python での Perseus 相当の操作

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() 関数
2 標本 t 検定	scipy.stats.ttest_ind() + multipletests()
ボルケーノプロット	matplotlib.pyplot の散布図 + 閾値線
階層的クラスタリング	seaborn.clustermap()

よく使うレシピ

レシピ：SILAC 比解析

使用する場合：LFQ の代わりに H/L または H/M/L 標識を使用する 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 を読み、「Clustered expression heatmap」の手順に従ってください（→ 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	rawファイルごとのQC統計：同定率、MS/MS数
results_differential.csv	差次的存在量の表：タンパク質ごとのlog2FC、pvalue、padj、significant
volcano_plot.pdf	上方／下方制御されたタンパク質を色分けし、上位のタンパク質にラベルを付けたボルケーノプロット
heatmap.pdf	上位の有意なタンパク質を対象とする階層的クラスタリングのヒートマップ（Zスコア正規化済み）
enrichment_up/	gseapyの出力ディレクトリ：上方制御されたタンパク質のGO/KEGGエンリッチメント
enrichment_down/	gseapyの出力ディレクトリ：下方制御されたタンパク質のGO/KEGGエンリッチメント

トラブルシューティング

問題	原因	解決策
MaxQuantでタンパク質が1つも同定されない	FASTAデータベースまたは酵素設定が誤っている、rawファイルのパスが見つからない	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」で停止する	メモリ不足（MaxQuantはrawファイル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 による濃縮解析 — 共通レシピで使用するタンパク質ネットワークの濃縮解析エンドポイント