41 lines
1.1 KiB
Python
41 lines
1.1 KiB
Python
import numpy as np
|
||
from andromeda import andromeda_score
|
||
|
||
def lp_score(masses, spec_masses, spec_intensities, tol=0.1, p=0.06):
|
||
"""
|
||
计算局部化概率（LP）分数。
|
||
|
||
参数:
|
||
- masses: 理论质谱峰列表
|
||
- spec_masses: 实验质谱峰列表
|
||
- spec_intensities: 实验峰对应的强度
|
||
- tol: 质量容差
|
||
- p: 匹配概率
|
||
|
||
返回:
|
||
- LP分数
|
||
"""
|
||
# 初始化匹配计数
|
||
k = 0
|
||
for mass in masses:
|
||
# 找到最近的实验峰并检查匹配
|
||
closest_idx = np.abs(np.array(spec_masses) - mass).argmin()
|
||
if np.abs(spec_masses[closest_idx] - mass) <= tol:
|
||
k += 1
|
||
|
||
# 计算Andromeda得分
|
||
n = len(masses) # 理论峰总数
|
||
base_score = andromeda_score(k, n, p)
|
||
|
||
# 计算LP归一化
|
||
lp = np.exp((base_score - 100) * np.log(10) / 10)
|
||
lp /= np.sum(lp) # 归一化操作
|
||
|
||
return lp
|
||
|
||
# 示例计算
|
||
masses = [100, 105, 110, 115, 120] # 示例理论峰
|
||
spec_masses = [101, 106, 111, 116, 121] # 示例实验峰
|
||
spec_intensities = [10, 20, 30, 40, 50] # 实验峰强度
|
||
lp_score(masses, spec_masses, spec_intensities)
|