代谢流分析（13C标记）¶

# === 同位素标记数据处理（MID校正与分析） ===
import numpy as np
import pandas as pd

# 1. 原始质谱MID数据（以丙酮酸为例，3个碳原子）
# M+0, M+1, M+2, M+3 的原始强度
raw_mid = np.array([0.15, 0.25, 0.35, 0.25])  # 未校正的MID

# 2. 天然同位素校正（简化版）
# 13C天然丰度为0.011
def correct_natural_isotope(mid, n_carbon):
    """简化的天然同位素校正（矩阵法）"""
    from scipy.special import comb
    n = len(mid)
    # 构建校正矩阵
    correction_matrix = np.zeros((n, n))
    p_natural = 0.011  # 13C天然丰度
    for i in range(n):
        for j in range(i + 1):
            # 二项分布计算天然同位素贡献
            correction_matrix[i, j] = comb(n_carbon - j, i - j) * \
                                       p_natural**(i-j) * (1-p_natural)**(n_carbon-i)
    # 求解校正后的MID
    corrected = np.linalg.solve(correction_matrix[:n,:n], mid[:n])
    corrected = corrected / corrected.sum()  # 归一化到1
    return corrected

corrected_mid = correct_natural_isotope(raw_mid, n_carbon=3)
print(f"校正前 MID: {raw_mid}")
print(f"校正后 MID: {np.round(corrected_mid, 4)}")

# 3. 计算13C标记富集度
# 平均标记度 = Σ(i × MID_i) / n_carbon
n_carbon = 3
avg_enrichment = sum(i * corrected_mid[i] for i in range(len(corrected_mid))) / n_carbon
print(f"平均13C富集度: {avg_enrichment:.3f}")  # 0=无标记, 1=完全标记

# 4. 批量处理多个代谢物的MID
metabolites = {
    "Pyruvate": {"mid": [0.15, 0.25, 0.35, 0.25], "n_c": 3},
    "Citrate": {"mid": [0.10, 0.15, 0.20, 0.25, 0.18, 0.08, 0.04], "n_c": 6},
    "Succinate": {"mid": [0.20, 0.30, 0.30, 0.15, 0.05], "n_c": 4},
}

results = []
for name, info in metabolites.items():
    mid = np.array(info["mid"])
    mid = mid / mid.sum()  # 先归一化
    corrected = correct_natural_isotope(mid, info["n_c"])
    enrichment = sum(i * corrected[i] for i in range(len(corrected))) / info["n_c"]
    results.append({"Metabolite": name, "Enrichment": enrichment})

df_results = pd.DataFrame(results)
print(df_results)

# === 使用IsoCor做天然同位素校正 ===
# pip install isocor

# IsoCor命令行使用
isocor \
  --input raw_mid_data.tsv \     # 输入原始MID数据
  --output corrected_mid.tsv \   # 输出校正后的MID
  --tracer "C" \                 # 示踪元素是碳
  --resolution 60000 \           # 质谱分辨率
  --mz_tol 0.005                 # m/z容忍度