蛋白质翻译后修饰（PTM）分析¶

# === 使用MSFragger进行开放搜索发现未知PTM ===
# MSFragger是FragPipe工具套件的核心引擎

# 1. 下载FragPipe: https://github.com/Nesvilab/FragPipe
# 包含MSFragger + Philosopher + PTM-Shepherd

# 2. 命令行运行开放搜索
java -jar MSFragger.jar \
  --database_name human_uniprot.fasta \  # 蛋白序列库
  --precursor_mass_lower -150 \          # 母离子质量偏移下限(-150 Da)
  --precursor_mass_upper 500 \           # 母离子质量偏移上限(+500 Da)
  --precursor_mass_units 1 \             # 单位：Da
  --fragment_mass_tolerance 20 \         # 碎片离子容忍度20ppm
  --output_format tsv \                  # 输出TSV格式
  sample1.mzML sample2.mzML             # 输入文件

# === PTM分析下游统计 ===
import pandas as pd
import numpy as np

# 1. 读取PTM搜索结果（以MaxQuant乙酰化位点为例）
acetyl = pd.read_csv("Acetyl (K)Sites.txt", sep="\t")

# 2. 基本过滤
acetyl = acetyl[acetyl["Reverse"] != "+"]  # 去反库
acetyl = acetyl[acetyl["Localization prob"] > 0.75]  # 高置信位点

# 3. 计算修饰化学计量比（Stoichiometry）
# 化学计量比 = 修饰肽段强度 / (修饰肽段强度 + 未修饰肽段强度)
# 反映实际修饰程度
if "Occupancy" in acetyl.columns:
    acetyl_occupied = acetyl[acetyl["Occupancy"] > 0]
    print(f"平均乙酰化占有率: {acetyl_occupied['Occupancy'].mean():.2%}")

# 4. PTM motif分析
from collections import Counter

# 提取磷酸化位点周围的序列motif（±7个氨基酸）
def extract_motif(sequence_window):
    """统计修饰位点周围的氨基酸频率"""
    if pd.isna(sequence_window) or len(sequence_window) < 15:
        return None
    return sequence_window  # 15个字符的窗口序列

motifs = acetyl["Sequence window"].apply(extract_motif).dropna()
# 统计中心位点（修饰位点）两侧的氨基酸偏好
center = 7  # 窗口中心位置
position_freq = {}
for pos in range(-3, 4):  # 分析修饰位点±3位的偏好
    aa_at_pos = [m[center + pos] for m in motifs if len(m) > center + pos]
    position_freq[pos] = Counter(aa_at_pos).most_common(3)

print("修饰位点周围氨基酸偏好:")
for pos, freq in position_freq.items():
    print(f"  位置{pos:+d}: {freq}")