604 宏蛋白组分析流程（Metaproteomics）¶

DNA（宏基因组）  → "菜单上有什么菜"（功能潜力）
RNA（宏转录组）  → "厨师正在炒哪些菜"（正在表达的基因）
蛋白（宏蛋白组） → "桌上端出来了什么菜"（实际产出的蛋白） ← 最接近功能
代谢物（代谢组） → "客人吃了什么"（最终的代谢产物）

样本采集 → 蛋白提取 → 酶解（胰蛋白酶） → 肽段纯化
    → LC-MS/MS质谱检测 → 数据库搜索 → 蛋白鉴定/定量
    → 物种分类 + 功能注释 → 差异分析

========== 实验设计要点 ==========
1. 蛋白提取方法影响极大：推荐SDS裂解 + FASP超滤法
2. 酶解：胰蛋白酶（Trypsin），37°C过夜消化
3. 质谱平台：Orbitrap（Thermo）或timsTOF（Bruker）
4. 采集模式：DDA（传统）或DIA（2024趋势，覆盖更全）
5. 建议测序深度：每样本2h梯度，产生约2-5万张MS2谱图
6. 生物学重复：≥3个
7. 技术重复：建议2-3个（质谱重复性不如测序）

# ========== 步骤1：构建蛋白搜索数据库 ==========
# 宏蛋白组最大的挑战：数据库太大 → 搜索慢 + 假阳性高
# 解决策略：用匹配的宏基因组数据构建样本特异性数据库

# 策略1：从宏基因组组装结果预测蛋白（推荐）
# 假设已经有宏基因组的组装结果
prodigal \
    -i metagenome_assembly.fa \  # 输入：宏基因组组装的contigs
    -a predicted_proteins.faa \  # 输出：预测的蛋白序列（FASTA格式）
    -p meta \                    # 宏基因组模式
    -f gff                       # 输出GFF格式注释

# 去冗余（减小数据库大小）
cd-hit \
    -i predicted_proteins.faa \  # 输入：预测的蛋白
    -o proteins_nr.faa \         # 输出：去冗余后的蛋白
    -c 0.95 \                    # 95%序列一致性阈值
    -n 5 \                       # 词长
    -M 16000 \                   # 内存限制16G
    -T 16                        # 16个线程

echo "原始蛋白数:" && grep -c ">" predicted_proteins.faa  # 统计原始数量
echo "去冗余后:" && grep -c ">" proteins_nr.faa           # 统计去冗余后数量

# 策略2：使用公共数据库（无宏基因组数据时的备选）
# 下载IGC（肠道微生物基因组综合目录）
wget http://meta.genomics.cn/meta/dataTools  # IGC数据库下载页面
# 或使用UniRef50/90子集

# 策略3：MetaPep核心肽段数据库（2023发表，速度更快）
# MetaPep是预构建的肠道微生物核心肽段库
# 搜索速度比传统MAG搜索快数倍

# 添加污染物数据库（必须！）
cat proteins_nr.faa \                        # 微生物蛋白
    /db/contaminants/cRAP.fasta \            # 常见污染物蛋白（角蛋白、胰蛋白酶等）
    /db/human/human_proteome.fasta \         # 人类蛋白组（去宿主蛋白用）
    > search_database.fasta                  # 合并后的搜索数据库

echo "最终数据库蛋白数:" && grep -c ">" search_database.fasta

# ========== 步骤2：使用MetaLab进行自动化分析 ==========
# MetaLab是一站式宏蛋白组分析平台
# 内部调用MaxQuant做蛋白鉴定，整合物种分类和功能注释

# MetaLab工作流程：
# 1. 数据库构建（自动迭代搜索，生成样本特异性数据库）
# 2. 肽段鉴定和蛋白定量（调用MaxQuant）
# 3. 物种分类（LCA算法）
# 4. 功能注释（COG/KEGG/GO）

# MetaLab是GUI软件，主要通过图形界面操作：
# 1. 加载.raw质谱文件
# 2. 选择/上传蛋白数据库
# 3. 设置搜索参数
# 4. 点击"Run"开始分析

# MetaLab的关键参数设置参考：
# - 酶：Trypsin/P
# - 最大漏切：2
# - 固定修饰：Carbamidomethyl (C)  半胱氨酸烷基化
# - 可变修饰：Oxidation (M), Acetyl (Protein N-term)
# - 母离子质量容差：±20 ppm（first search），±4.5 ppm（main search）
# - 碎片离子质量容差：±0.5 Da（ion trap）或 ±20 ppm（Orbitrap）
# - FDR：1%（肽段和蛋白两个层面）

# ========== 如果不用MetaLab，可以手动用MaxQuant ==========

# MaxQuant也是GUI软件，命令行可以这样调用：
# （需要先在GUI中配置好参数并保存为mqpar.xml）
dotnet MaxQuantCmd.exe mqpar.xml  # 运行MaxQuant命令行模式

# MaxQuant核心输出文件：
# txt/proteinGroups.txt  — 蛋白组定量结果（最重要）
# txt/peptides.txt       — 肽段鉴定结果
# txt/evidence.txt       — 详细鉴定证据
# txt/msms.txt           — MS2谱图匹配

# 也可以用命令行工具替代MaxQuant：

# 用MSFragger（速度更快）
java -jar MSFragger.jar \
    --database search_database.fasta \     # 蛋白数据库
    --num_threads 16 \                     # 线程数
    *.mzML                                 # 输入：质谱数据（mzML格式）

# RAW文件转mzML格式
# 用ThermoRawFileParser（Linux可用）
ThermoRawFileParser \
    -i sample.raw \                        # 输入：Thermo RAW文件
    -o converted/ \                        # 输出目录
    -f 2                                   # 格式2=mzML

# ========== 步骤3：基于鉴定的肽段做物种分类 ==========
import pandas as pd  # 导入pandas

# 读取MaxQuant的肽段鉴定结果
peptides = pd.read_csv("txt/peptides.txt", sep="\t")  # 读取肽段表

# MetaLab会自动做LCA（最低共同祖先）分类
# 手动实现LCA的思路：
# 1. 每个肽段可能匹配多个蛋白（多物种共有肽段）
# 2. 找到这些蛋白的最低共同祖先
# 3. 物种特异性肽段 → 直接分到该物种
# 4. 共有肽段 → 分到上一级分类（属/科/目...）

# 用Unipept做在线物种分类（基于胰蛋白酶肽段）
# https://unipept.ugent.be/
# 上传肽段列表，自动返回物种分类和功能注释

# 从proteinGroups.txt提取物种信息
proteins = pd.read_csv("txt/proteinGroups.txt", sep="\t")  # 读取蛋白组结果

# 过滤：去除污染物、反库命中、只有1个肽段的蛋白
proteins_clean = proteins[
    (proteins["Potential contaminant"] != "+") &   # 不是污染物
    (proteins["Reverse"] != "+") &                  # 不是反库（decoy）
    (proteins["Peptides"] >= 2)                     # 至少2个肽段
].copy()

print(f"过滤后蛋白组数: {len(proteins_clean)}")  # 打印过滤后的蛋白组数

# 提取LFQ强度（Label-Free Quantification）
lfq_cols = [c for c in proteins_clean.columns if c.startswith("LFQ intensity")]  # 找LFQ列
lfq_data = proteins_clean[lfq_cols]  # 提取LFQ数据
lfq_data.index = proteins_clean["Protein IDs"]  # 用蛋白ID作为索引

# ========== 步骤4：功能注释 ==========

# 方法1：eggNOG-mapper注释（推荐，最全面）
emapper.py \
    -i identified_proteins.faa \     # 输入：鉴定到的蛋白序列
    -o eggnog_annotation \           # 输出前缀
    --cpu 16 \                       # CPU数
    -m diamond \                     # 比对方法
    --data_dir /db/eggnog/ \         # eggNOG数据库
    --go_evidence non-electronic \   # GO证据类型
    --tax_scope auto                 # 自动选择分类范围

# eggNOG-mapper输出中重要的列：
# COG_category — COG功能分类
# KEGG_ko      — KEGG直系同源基因编号
# KEGG_Pathway — KEGG通路
# GO_terms     — GO功能注释

# 方法2：InterProScan注释（更慢但更准确）
interproscan.sh \
    -i identified_proteins.faa \     # 输入蛋白序列
    -o interpro_results.tsv \        # 输出
    -f tsv \                         # TSV格式
    -cpu 16 \                        # CPU数
    -goterms \                       # 输出GO注释
    -pa                              # 输出通路注释

# ========== 步骤5：差异蛋白分析 ==========
library(limma)    # 线性模型差异分析
library(ggplot2)  # 可视化

# 读取LFQ定量数据
protein_data <- read.delim("proteinGroups.txt",  # 读取MaxQuant蛋白组结果
                            sep = "\t")

# 提取LFQ强度列
lfq <- protein_data[, grep("^LFQ", colnames(protein_data))]  # 提取LFQ列
rownames(lfq) <- protein_data$Protein.IDs  # 蛋白ID作为行名

# 过滤：去除污染物和反库
keep <- protein_data$Potential.contaminant != "+" &  # 非污染物
        protein_data$Reverse != "+" &                 # 非反库
        protein_data$Peptides >= 2                    # 至少2个肽段
lfq <- lfq[keep, ]  # 应用过滤

# log2转换
lfq_log <- log2(lfq)  # log2转换

# 缺失值处理：用最小值填充（模拟检测限以下的低丰度蛋白）
min_val <- min(lfq_log[is.finite(as.matrix(lfq_log))], na.rm = TRUE)  # 找最小值
lfq_log[!is.finite(as.matrix(lfq_log))] <- min_val - 1  # NA和-Inf用最小值-1填充

# 设计矩阵
group <- factor(c(rep("Control", 5), rep("T2D", 5)))  # 分组
design <- model.matrix(~ 0 + group)  # 设计矩阵（无截距）
colnames(design) <- c("Control", "T2D")  # 列名

# limma差异分析
fit <- lmFit(lfq_log, design)  # 拟合线性模型

# 构建对比矩阵
contrast <- makeContrasts(T2D - Control,  # T2D vs Control
                          levels = design)
fit2 <- contrasts.fit(fit, contrast)  # 应用对比
fit2 <- eBayes(fit2)                  # 经验贝叶斯调节

# 提取结果
results <- topTable(fit2,        # 从拟合结果提取
                    number = Inf,  # 所有蛋白
                    sort.by = "p") # 按P值排序

# 筛选显著差异蛋白
sig_proteins <- results[results$adj.P.Val < 0.05 &  # 校正P值<0.05
                        abs(results$logFC) > 1, ]     # |log2FC|>1
cat("显著差异蛋白数:", nrow(sig_proteins), "\n")       # 打印数量

# 火山图
ggplot(results, aes(x = logFC, y = -log10(adj.P.Val))) +  # X=log2FC, Y=-log10P
  geom_point(aes(color = adj.P.Val < 0.05 & abs(logFC) > 1),  # 按显著性着色
             size = 1.5) +
  scale_color_manual(values = c("grey", "red")) +  # 灰色=不显著，红色=显著
  geom_hline(yintercept = -log10(0.05), linetype = "dashed") +  # P值阈值线
  geom_vline(xintercept = c(-1, 1), linetype = "dashed") +      # FC阈值线
  labs(title = "Metaproteomics: T2D vs Control",  # 标题
       x = "log2(Fold Change)",                    # X轴
       y = "-log10(adjusted P-value)") +           # Y轴
  theme_minimal()                                   # 简约主题

# ========== metaExpertPro：支持DIA的宏蛋白组分析 ==========
# 2024年发表在Molecular & Cellular Proteomics
# 比MetaLab多鉴定约22%的肽段

# metaExpertPro工作流程：
# 1. 构建样本特异性谱图库（Spectral Library）
# 2. DIA数据解析（用DIA-NN）
# 3. 蛋白定量和物种分类
# 4. 功能注释

# DIA模式的优势：
# - 覆盖更全面（不像DDA会漏掉低丰度蛋白）
# - 定量更准确（减少缺失值）
# - 重复性更好

# 安装DIA-NN（metaExpertPro的核心引擎）
# 从GitHub下载：https://github.com/vdemichev/DiaNN

# DIA-NN运行示例
diann \
    --f sample.mzML \                          # 输入：DIA质谱数据
    --lib spectral_library.tsv \               # 谱图库
    --fasta search_database.fasta \            # 蛋白数据库
    --threads 16 \                             # 线程数
    --out results.tsv \                        # 输出结果
    --qvalue 0.01 \                            # FDR阈值1%
    --matrices                                 # 输出定量矩阵

========== 宏蛋白组标准分析流程 ==========
RAW数据 → 格式转换(mzML) → 数据库构建(宏基因组预测蛋白)
    → 数据库搜索(MaxQuant/MSFragger) → 蛋白鉴定/定量
    → 物种分类(LCA) + 功能注释(eggNOG/KEGG)
    → 差异分析(limma) → 可视化

========== 核心工具 ==========
自动化平台：MetaLab（推荐）/ metaExpertPro（DIA模式）
搜索引擎  ：MaxQuant / MSFragger / Comet
DIA分析   ：DIA-NN / Spectronaut
格式转换  ：ThermoRawFileParser / msconvert
蛋白注释  ：eggNOG-mapper / InterProScan / Unipept
差异分析  ：limma / Perseus / MSstats
物种分类  ：Unipept / MetaLab内置LCA

========== 数据库选择 ==========
最推荐：宏基因组组装+Prodigal预测的蛋白（样本特异性）
备选1：IGC（肠道基因组综合目录，~10M基因）
备选2：UniRef50/90（通用蛋白数据库）
备选3：MetaPep（核心肽段库，速度快）
必加  ：cRAP污染物数据库 + 人类蛋白组

========== DDA vs DIA ==========
DDA：传统方式，工具成熟（MaxQuant），适合已建立流程的实验室
DIA：2024趋势，覆盖更全，定量更准（DIA-NN），推荐新项目使用