蛋白组Label-free定量¶

# === MaxQuant命令行运行 ===
# 下载MaxQuant: https://www.maxquant.org/

# 1. 准备参数文件mqpar.xml（通常用GUI配置后导出）
# 关键参数设置：
# - Label-free quantification: LFQ
# - Match between runs: 勾选
# - Min. ratio count: 2（至少2个肽段用于定量）
# - FASTA文件: uniprot_human.fasta

# 2. 命令行运行MaxQuant
dotnet MaxQuantCmd.exe mqpar.xml  # 运行分析（可能需要几小时到几天）

# === MaxQuant结果下游分析（Perseus风格） ===
import pandas as pd
import numpy as np
from scipy.stats import ttest_ind
from statsmodels.stats.multitest import multipletests

# 1. 读取MaxQuant的proteinGroups.txt输出
pg = pd.read_csv("combined/txt/proteinGroups.txt", sep="\t")

# 2. 基本过滤
pg = pg[pg["Reverse"] != "+"]       # 去除反库匹配（decoy）
pg = pg[pg["Potential contaminant"] != "+"]  # 去除常见污染物
pg = pg[pg["Only identified by site"] != "+"]  # 去除仅通过修饰位点鉴定的
print(f"过滤后蛋白数: {len(pg)}")

# 3. 提取LFQ强度列
lfq_cols = [c for c in pg.columns if c.startswith("LFQ intensity")]  # 找到所有LFQ列
lfq_data = pg[["Protein IDs", "Gene names"] + lfq_cols].copy()
lfq_data.set_index("Protein IDs", inplace=True)

# 4. 替换0为NaN，做log2变换
lfq_numeric = lfq_data[lfq_cols].replace(0, np.nan)  # MaxQuant用0表示缺失
lfq_log2 = np.log2(lfq_numeric)  # log2变换

# 5. 过滤：至少在一个组中70%样本有值
group1_cols = lfq_cols[:3]  # 前3个样本为组1
group2_cols = lfq_cols[3:6]  # 后3个样本为组2
valid_g1 = lfq_log2[group1_cols].notna().sum(axis=1) >= 2  # 组1至少2/3有值
valid_g2 = lfq_log2[group2_cols].notna().sum(axis=1) >= 2  # 组2至少2/3有值
lfq_filtered = lfq_log2[valid_g1 & valid_g2]

# 6. 差异蛋白分析
results = []
for idx in lfq_filtered.index:
    g1 = lfq_filtered.loc[idx, group1_cols].dropna().values  # 组1数据
    g2 = lfq_filtered.loc[idx, group2_cols].dropna().values  # 组2数据
    fc = g1.mean() - g2.mean()  # log2 fold change
    _, pval = ttest_ind(g1, g2)
    results.append({"Protein": idx, "log2FC": fc, "pvalue": pval})

df_res = pd.DataFrame(results)
df_res["padj"] = multipletests(df_res["pvalue"], method="fdr_bh")[1]  # FDR校正
sig_proteins = df_res[(df_res["padj"] < 0.05) & (abs(df_res["log2FC"]) > 1)]
print(f"显著差异蛋白: {len(sig_proteins)}")