"""
Mass list of aminoacids, protenimodifications, glycan monomers, etc.
Contains also all functions for mass calculation.
Source of mass values from:
http://web.expasy.org/findmod/findmod_masses.html
"""
#我需要知道那些支链也碎的mz，所以我需要知道怎么表明立体结构,以糖所处的位置为糖峰标号，新糖峰标记，比如以前的BY无法标支链碎裂峰，现在以糖来做下标
#pglyco结构转换为图结构
#可以做ficharge到percursor charge，丢水，丢NH3， ？三电荷的碎片有吗，我有看过，那还是保留吧 谱图预测可以快速匹配大量碎片
#对于糖来说，如果对称糖结构，那么不同的来源会有相同的mz，出现在一个峰上,可能只能强度平分
#所有的函数都汇聚在pepfragmass（）函数了，可以只用那个函数
#HCD 碎裂，by离子会把糖也全部丢掉，而且只碎糖的时候也会碎裂多根键 找一下HCD碎裂的诊断离子表
#所有可能的糖结构，需要能做成可以理解的模式。
# 。其实还是需要考虑一下，为什么要做谱图搜索，对于算法而言，糖产生一些细小的结构改变，只是去做一些分析。但是对谱图预测而言，需要重新去预测。那我现在最多只能拿已有的糖去做谱图预测。
#IgG每个位点有的糖就那么几十种，实在不行可以按照那个来。
#糖的峰要慢慢check，注意一下mz
#单糖可以通过多步反应，掉中间那个，可以不预测，只作为诊断离子，HCD模式如何计算碎裂多次
#用ETD模式，CID模式，和HCD模式，依次测试
#可以结合sequence searching 和spectral searching，先匹配预测的mz，再匹配预测的强度
#搜索做断两次+单糖+丢Fuc，预测做丢一次+丢Fuc
#有不同的碎裂模式，不同的碎裂模式要用相应的谱图验证
import re
import copy
import dgl
import torch
import ipdb
import itertools
# --------------------------- Other masses ----------------------------#
# http://www.sisweb.com/referenc/source/exactmas.htm
MASS = {}
MASS["H2O"] = 18.01057
MASS["NH3"] = 17.02655
MASS["H+"] = 1.00728
MASS["H"] = 1.007825
MASS["C"] = 12.0
MASS["N"] = 14.003074
MASS["O"] = 15.994915
MASS["S"] = 31.972072
MASS["P"] = 30.973763
MASS["Na"] = 22.989770
MASS["K"] = 38.963708
MASS["Mg"] = 23.985045
# --------------------------- Aminoacids ------------------------------#
AMINOACID = {}
AMINOACID["A"] = 71.03711
AMINOACID["R"] = 156.10111
AMINOACID["N"] = 114.04293
AMINOACID["J"] = 114.04293 #带糖的N
AMINOACID["D"] = 115.02694
AMINOACID["C"] = 103.00919
AMINOACID["E"] = 129.04259
AMINOACID["Q"] = 128.05858
AMINOACID["G"] = 57.02146
AMINOACID["H"] = 137.05891
AMINOACID["I"] = 113.08406
AMINOACID["L"] = 113.08406
AMINOACID["K"] = 128.09496
AMINOACID["M"] = 131.04049
AMINOACID["F"] = 147.06841
AMINOACID["P"] = 97.05276
AMINOACID["S"] = 87.03203
AMINOACID["T"] = 101.04768
AMINOACID["W"] = 186.07931
AMINOACID["Y"] = 163.06333
AMINOACID["V"] = 99.06841
# ------------------------------ Glycans ------------------------------#
#Glycan的列表，Monoisotopic mass,库中的所有修饰保持大写，输入的大小写都可以兼容，因为可以转换
GLYCAN = {}
GLYCAN["DHEX"] = 146.0579
GLYCAN["HEX"] = 162.0528
GLYCAN["HEXNAC"] = 203.0794
GLYCAN["NEUAC"] = 291.0954
GLYCAN["NEUGC"] = 307.0903
GLYCAN["FUC"] = 146.0579
# --------------------------- Compositions ----------------------------#
# http://www.webqc.org/aminoacids.php
COMPOSITION = {}
COMPOSITION['A'] = {'C': 3, 'H': 7, 'N': 1, 'O': 2}
COMPOSITION['C'] = {'C': 3, 'H': 7, 'N': 1, 'O': 2, 'S': 1}
COMPOSITION['D'] = {'C': 4, 'H': 7, 'N': 1, 'O': 4}
COMPOSITION['E'] = {'C': 5, 'H': 9, 'N': 1, 'O': 4}
COMPOSITION['F'] = {'C': 9, 'H': 11, 'N': 1, 'O': 2}
COMPOSITION['G'] = {'C': 2, 'H': 5, 'N': 1, 'O': 2}
COMPOSITION['H'] = {'C': 6, 'H': 9, 'N': 3, 'O': 2}
COMPOSITION['I'] = {'C': 6, 'H': 13, 'N': 1, 'O': 2}
COMPOSITION['K'] = {'C': 6, 'H': 14, 'N': 2, 'O': 2}
COMPOSITION['L'] = {'C': 6, 'H': 13, 'N': 1, 'O': 2}
COMPOSITION['M'] = {'C': 5, 'H': 11, 'N': 1, 'O': 2, 'S': 1}
COMPOSITION['N'] = {'C': 4, 'H': 8, 'N': 2, 'O': 3}
COMPOSITION['J'] = {'C': 4, 'H': 8, 'N': 2, 'O': 3}
COMPOSITION['P'] = {'C': 5, 'H': 9, 'N': 1, 'O': 2}
COMPOSITION['Q'] = {'C': 5, 'H': 10, 'N': 2, 'O': 3}
COMPOSITION['R'] = {'C': 6, 'H': 14, 'N': 4, 'O': 2}
COMPOSITION['S'] = {'C': 3, 'H': 7, 'N': 1, 'O': 3}
COMPOSITION['T'] = {'C': 4, 'H': 9, 'N': 1, 'O': 3}
COMPOSITION['V'] = {'C': 5, 'H': 11, 'N': 1, 'O': 2}
COMPOSITION['W'] = {'C': 11, 'H': 12, 'N': 2, 'O': 2}
COMPOSITION['Y'] = {'C': 9, 'H': 11, 'N': 1, 'O': 3}
COMPOSITION["HEX"] = {'C': 6, 'H': 12, 'N': 0, 'O': 6}
COMPOSITION["DHEX"] = {'C': 6, 'H': 12, 'N': 0, 'O': 5}
COMPOSITION["HEXNAC"] = {'C': 8, 'H': 15, 'N': 1, 'O': 6}
COMPOSITION["NEUAC"] = {'C': 11, 'H': 19, 'N': 1, 'O': 9}
COMPOSITION["NEUGC"] = {'C': 11, 'H': 19, 'N': 1, 'O': 10}
COMPOSITION["H2O"] = {'H': 2, 'O': 1}
# --------------------------- Proteinmodifications---------------------#
PROTEINMODIFICATION = {}
PROTEINMODIFICATION["ACE"] = {"mass": 42.0106,
"composition":{'C': 2, 'H': 2, 'O': 1}
} # Acetylation
PROTEINMODIFICATION["AMI"] = {"mass": -0.9840,
"targets": {"CTERM"},
"composition":{'H': 1, 'O': -1, 'N': 1}
} # Amidation
# Cys_CAM Idoacetamide treatment (carbamidation)
PROTEINMODIFICATION["CAR"] = {"mass": 57.021464,
"targets": {"C","NTERM"},
"composition":{'C':2, 'H':3, 'O':1, 'N':1}
}
# Cys_CM, Iodoacetic acid treatment (carboxylation)
PROTEINMODIFICATION["CM"] = {"mass": 58.005479,
"targets": {"C"},
"composition":{'C': 2, 'H': 2, 'O': 2}
}
PROTEINMODIFICATION["DEA"] = {"mass": 0.9840,
"targets": {"N", "Q"},
"composition":{'H': -1, 'O': 1, 'N': -1}
} # Deamidation
PROTEINMODIFICATION["DEH"] = {"mass": -18.0106,
"targets": {"S", "T"},
"composition":{'H': -2, 'O': -1}
} # Dehydration Serine
PROTEINMODIFICATION["GUA"] = {"mass": 42.0218,
"targets": {"K"},
"composition":{'C': 1, 'H': 2,'N': 2}
} # K Guandination
PROTEINMODIFICATION["HYD"] = {"mass": 15.9949,
"composition":{'O': 1}
} # Hydroxylation
PROTEINMODIFICATION["MET"] = {"mass": 14.0157,
"composition":{'C': 1, 'H': 2}
} # Methylation
PROTEINMODIFICATION["OXI"] = {"mass": 15.9949,
"targets": {"M"},
"composition":{'O': 1}
} # MSO
PROTEINMODIFICATION["GLY"] = {"mass":0.0 ,
"targets": {"N", "S", "T"},
"composition":{}
}
#按照后修饰和氨基酸对PROTEINMODIFICATION里的后修饰进行排序
def getModificationNames():
""" Return a sorted list of Modification names.
Only contains modifications with targets declaration"""
def sort_names(name):
if "_" in name:
amino, mod = name.split("_")
return mod, amino
return name, " "
names = set()
for key in PROTEINMODIFICATION:
if "targets" in PROTEINMODIFICATION[key]:
names.add(key)
return sorted(names, key=sort_names)
# print("Modifications in list",getModificationNames())
# ------------------------------- Glycan graph ---------------------------#
#表示糖的立体结构，通过DFS计算碎裂一条边以后，剩下的糖部分
node2dict={0:"peptide",1:"HexNAc",2:"Hex",3:"Fuc",4:"NeuAc",5:"NeuGc"}
node2char={"P":"peptide","H" : "Hex", "N" : "HexNAc", "A" : "NeuAc", "G" : "NeuGc" , "F" :"Fuc"}
node2idx={"P":0,"H" : 2, "N" : 1, "A" : 4, "G" : 5 , "F" :3}
def peptide_process(peptide):
"""Create dict for input sequende.
Input:LSVECAJK_5.Car._6_2_(N(F)(N(H(H)(N)(H(N)))))
Output:{'sequence': 'LSVECAJK', 'charge': '2', 'modifications': [{'name': '(Car)1', 'amino': 'C', 'position': 4, 'type': 'mod'},
{'name': '(Hex)3(HexNAc)4(Fuc)1', 'amino': 'J', 'position': 6, 'structure': '(N(F)(N(H(H)(N)(H(N)))))', 'type': 'glyco'}]}
"""
a=peptide.split("_")
sequence=a[0]
# print(sequence)
charge=a[3]
mod_dict={}
mod_dictg={}
mod_list=[]
if a[1] !="None":
mod=a[1].rstrip(".").split(".")
# print(mod)
for i in range(0,len(mod),2):
mod_dict=copy.deepcopy(mod_dict)
mod_dict["name"]="("+mod[i+1]+")1"
mod_dict["amino"]=sequence[int(mod[i])-1]
mod_dict["position"]=int(mod[i])-1
mod_dict["type"]="mod"
mod_list.append(mod_dict)
# print("mod_list",mod_list)
glycos = a[4].split(';') # ['(N(A)(H(A)))', '(N(A)(H(A)))']
positions = a[2].split(';') # ['1', '14']
for glyco, position in zip(glycos, positions):
composition = {
"H": glyco.count("H"),
"N": glyco.count("N"),
"A": glyco.count("A"),
"G": glyco.count("G"),
"F": glyco.count("F")
}
comp_str = ""
for k in composition:
if composition[k] > 0:
comp_str += f"({node2char[k]}){composition[k]}"
position = int(position)
# 每次循环创建新的mod_dictg
mod_dictg = {
"name": comp_str,
"amino": peptide[position],
"position": position,
"structure": glyco,
"type": "glyco"
}
mod_list.append(mod_dictg)
peptide_dict={"sequence":sequence,"charge":charge,"modifications":mod_list}
return peptide_dict
def glyco_process_0(glyco:str):
glyco_ind=re.sub("\w","",glyco)
glyco_cha=re.sub("\W","",glyco)
# ipdb.set_trace()
if len(glyco_ind) != 2*len(glyco_cha):
print("alert:len(glyco_ind) != 2*len(glyco_cha)")
raise AssertionError
node1_index=-1
node2_index=0
node1="P"
node_group=[]
# ipdb.set_trace()
for i in range(len(glyco_ind)):
if glyco_ind[i] =="(":
node2=glyco_cha[node2_index]
node1_index+=1
node2_index+=1
node=[node1_index,node2_index]
node_group.append(node)
# ipdb.set_trace()
# print(node_group)
# print("(",node1_index)
# print("(",node2_index)
node1=node2
node1_index=node2_index-1
if glyco_ind[i]==")":
node1_index-=1
node1=glyco_cha[node1_index]
# print(")",node1_index)
# ipdb.set_trace()
return node1_index,node2_index,node_group,glyco_cha
def glyco_process(glyco:str):
glyco_ind=re.sub("\w","",glyco)
glyco_cha=re.sub("\W","",glyco)
# ipdb.set_trace()
if len(glyco_ind) != 2*len(glyco_cha):
print("alert:len(glyco_ind) != 2*len(glyco_cha)")
raise AssertionError
##双指针,同时需要dict指明str的index和node index的关系
node1_ptr=0
node2_ptr=0
node1=0##多余的起点
indict={-1:0,}
node_count=0
edge_index=[]
stack=[-2]
# ipdb.set_trace()
while node2_ptr < len(glyco):
if glyco[node2_ptr]=="(":
# print("(")
node1_ptr=stack[-1]
stack.append(node2_ptr)
elif glyco[node2_ptr]==")":
# print(")")
node1_ptr=stack.pop()
else:##碰到字母
# print("add")
node_count+=1
indict[node2_ptr]=node_count
edge_index.append([indict[node1_ptr+1],indict[node2_ptr]])
# print(node2_ptr)
# ipdb.set_trace()
node2_ptr+=1
# print(edge_index)
return node1_ptr,node2_ptr,edge_index,glyco_cha
# print(glyco_process("(N(N(H(H)(H(H))))(F))"))
# glyco_process_0("(N(N(H(H)(H(H))))(F))")
# glyco_process_0("(N(F)(N(H(H(N(H)))(H(N(H(G)))))))")
# glyco_process_0("(N(N(H(H(H(H)))(H(H(H(H)))(H(H(H(H))))))))")
# glyco_process_0("(N(N(H(H(N(H(A)))(N(H)))(H(H))))(F))")
# glyco_process("(N(N(H(H(N(H(A)))(N(H)))(H(H))))(F))")
# ipdb.set_trace()
# global diff
# diff=0
# def glycoprocess_diff(glyco):
# glyco=glyco["H,N,A,G,F"]
# graphedge_0=glyco_process_0(glyco)[2]
# graphedge_1=glyco_process(glyco)[2]
# if graphedge_0!= graphedge_1:
# # print(glyco)
# # print(graphedge_1)
# # print(graphedge_0)
# global diff
# diff+=1
# # print(diff)
# # ipdb.set_trace()
# return diff
# # ipdb.set_trace()
# import pandas as pd
# mouse_gdb=pd.read_csv("/remote-home/yxwang/test/zzb/DeepGlyco/DeepSweet_v1/code/task_processing/AHGF/pGlyco-N-Mouse.gdb",sep="\t")
# mouse_gdb.apply(glycoprocess_diff,axis=1)
# ipdb.set_trace()
def GlycanFrag_struc(glyco_graph):
edge_index=glyco_graph[2]
# print("edge_index",edge_index)
#eg.[[0, 1], [1, 2], [1, 3], [3, 4], [4, 5], [4, 6], [6, 7], [7, 8], [5, 9], [9, 10], [10, 11]]
nodef="P"+glyco_graph[3]
nodef=[node2idx[i] for i in nodef]
# print("nodef",nodef)
#eg. nodef [0, 1, 3, 1, 2, 1, 2, 1, 2, 2, 1, 3]
g=dgl.graph(edge_index)
g.ndata["mononer"]=torch.Tensor(nodef).to(int)
return g
def GlycanFrag(glyco_graph):
edge_index=glyco_graph[2]
# print("edge_index",edge_index)
#eg.[[0, 1], [1, 2], [1, 3], [3, 4], [4, 5], [4, 6], [6, 7], [7, 8], [5, 9], [9, 10], [10, 11]]
nodef="P"+glyco_graph[3]
nodef=[node2idx[i] for i in nodef]
# print("nodef",nodef)
#eg. nodef [0, 1, 3, 1, 2, 1, 2, 1, 2, 2, 1, 3]
g=dgl.graph(edge_index)
g.ndata["mononer"]=torch.Tensor(nodef).to(int)
glycanfrag={}
for eid in range(len(g.edges()[0])):
graphe=dgl.remove_edges(g,eid)
removed_edge=g.edges()[0][eid].item(),g.edges()[1][eid].item()
#c从这里可以看出，removed_edge是按照eid的顺序从0到len(g.edges()[0])从图上裂掉的。所以removed_edge的顺序就是eid的顺序
eid_res=list(dgl.dfs_edges_generator(graphe,0))
noderes=[0]
graphe_edges=graphe.edges()
for eid2 in eid_res:
node1,node2=graphe_edges[0][eid2],graphe_edges[1][eid2]
noderes.append(node1.item())
noderes.append(node2.item())
noderes=list(set(noderes))
res_monomer=g.ndata["mononer"][noderes].numpy().tolist()
res_sugar=[node2dict[i] for i in res_monomer if i != 0]
# print("removed_edge:",removed_edge,"res_index:",res_monomer,"res_sugar:",res_sugar)
glycanfrag[removed_edge]=res_sugar
return glycanfrag
def GlycanFrag_HCD(glyco_graph):
edge_index=glyco_graph[2]
# print("edge_index",edge_index)
nodef="P"+glyco_graph[3]
nodef=[node2idx[i] for i in nodef]
# print("nodef",nodef)
g=dgl.graph(edge_index)
g.ndata["mononer"]=torch.Tensor(nodef).to(int)
glycanfrag={}
# print("g.edges",g.edges())
for eid in range(len(g.edges()[0])):
for eid1 in range(eid,len(g.edges()[0])):
graphe=dgl.remove_edges(g,[eid,eid1])
removed_edge1=g.edges()[0][eid].item(),g.edges()[1][eid].item()
removed_edge2=g.edges()[0][eid1].item(),g.edges()[1][eid1].item()
removed_edge=[removed_edge1,removed_edge2]
# ------------------- turn to bidirectional -------------------
u , v = graphe.edges(order = "eid")
graphe.add_edges(v , u) # bidirect
# print("removed_edge_list",removed_edge)
three_splits=[]
for dfs_initial_node in [removed_edge1[0],removed_edge1[1],removed_edge2[1]]:
# print("dfs_initial_node:",dfs_initial_node)
noderes=[dfs_initial_node]
eid_res=list(dgl.dfs_edges_generator(graphe,dfs_initial_node))
#很奇怪，这里可以接受0，但是不能接受removed_edge1[0] removed_edge1[1] removed_edge2[2]
# print("removed_edge1[0]",removed_edge1[0])
# print("eid_res",eid_res)
graphe_edges=graphe.edges()
# print("graphe_edges",graphe_edges)
for eid2 in eid_res:
# print("eid_res_eid2",eid2)
node1,node2=graphe_edges[0][eid2],graphe_edges[1][eid2]
noderes.append(node1.item())
noderes.append(node2.item())
noderes=list(set(noderes))
# ipdb.set_trace()
res_monomer=g.ndata["mononer"][noderes].numpy().tolist()
res_sugar=[node2dict[i] for i in res_monomer if i != 0]
three_splits.append( {"res_index":res_monomer,"res_sugar":res_sugar})
# print("removed_edge:",removed_edge,"three_splits:", three_splits)
removed_str=""
for n in removed_edge:
for n1 in n:
removed_str+=str(n1)
glycanfrag[removed_str]=three_splits
return glycanfrag
# import ipdb
# ipdb.set_trace()
# glyco_graph=glyco_process("(N(F)(N(H(H)(H))))")
# glycanfrag=GlycanFrag_HCD(glyco_graph)
# print(glycanfrag)
#对得到的糖碎片计算mz，算好了放到最后一个函数，做整合
# ------------------------------- Functions ---------------------------#
#计算肽段部分质量
def calcPeptideMass(peptide):
mass = 0
for s in peptide["sequence"]:
try:
mass += AMINOACID[s]
except KeyError:
print("Unknown Aminoacid '"+s+"'!")
raise
return mass
# print("Peptide",peptide)
# peptidemass=calcPeptideMass(peptide)
# print("Mass for peptide without modifications",peptidemass)
#计算后修饰的总质量
def calcModificationMass(peptide):
mass = 0
if peptide["modifications"] is None:
return mass
if len(peptide["modifications"]) ==0:
return mass
# TODO: Checks for modification consistency
# a) amount of amino acids must be in sequence
# b) a given position (except -1) can only be
for mod in peptide["modifications"]:
name=mod["name"]
for part in re.findall(r"\(.+?\)-?\d+", name.upper()):
monomer, amount = part.split(")")
monomer = monomer[1:]
amount = int(amount)
if monomer in GLYCAN:
mass += GLYCAN[monomer]*amount
elif monomer in MASS:
mass += MASS[monomer]*amount
elif monomer in PROTEINMODIFICATION:
mass += PROTEINMODIFICATION[monomer]["mass"]*amount
else:
raise Exception("cannot find monomer {} in {}".format(monomer, name))
return mass
# modificationmass=calcModificationMass(peptide)
# print("mass for modifications",modificationmass)
#算修饰的mz ver:算糖基质量，可以算加和物质
def calcModpepMass(peptide):
peptidemass = calcPeptideMass(peptide)
modificationmass = calcModificationMass(peptide)
ModpepMass=peptidemass+modificationmass
return ModpepMass
# ModpepMass=calcModpepMass(peptide)
# print("Mass for peptide with modifications",ModpepMass)
def pepfragmass(input,mode,maxcharge=3):
peptide=peptide_process(input)
# print("peptide",peptide)
for modification in peptide["modifications"]:
name=modification["name"].replace(" ", "")
if re.match(r"^(\(.+?\)-?\d+)+$", name) == None:
print("name",name)
raise Exception(r"""Input string '{}' doesn't follow the
regex '^(\(.+?\)-?\d+)+$'
Please surrond monomers by brackets followed by the amount of
the monomer, e.g. '(NeuAc)1(H2O)-1(H+)1'""".format(name))
sequence=peptide["sequence"]
length=len(sequence)
FragCharge=list(range(1,min((int(peptide["charge"])+1),maxcharge)))#FragCharge最大值调小了
b_ions=[]
y_ions=[]
B_ions=[]
Y_ions=[]
# 计算b/y离子
#检查一下糖的碎片模式，碎片电荷以及中性丢失的情况
#M+H-H2O，M+H， M+2H-H2O， M+2H ...
#计算B/Y ions
#多糖基化肽段另一侧应该完全碎裂，还是应该保留
if "HCD_1" in mode:
ModpepMass = calcModpepMass(peptide) + MASS["H2O"]
peptide_copy = copy.deepcopy(peptide) # 原始备份，最后还原
frag_index_offset = 0 # 新增：控制多个糖之间的 frag_index 编号递增
for i, mod in enumerate(peptide["modifications"]):
if mod["type"] != "glyco":
continue
structure = mod["structure"]
glyco_graph = glyco_process(structure)
glycanfrags = GlycanFrag(glyco_graph)
for local_index, (eid, glycan) in enumerate(glycanfrags.items()):
frag_index = frag_index_offset + local_index # ✅ 正确编号
composition = {
"H": glycan.count("Hex"),
"N": glycan.count("HexNAc"),
"A": glycan.count("NeuAc"),
"G": glycan.count("NeuGc"),
"F": glycan.count("Fuc")
}
comp_str = ''.join(
f"({node2char[k]}){v}" for k, v in composition.items() if v > 0
)
# 修改当前糖基的 name 字段，其它糖保持原样
peptide["modifications"][i]["name"] = comp_str
FragYMass = calcModpepMass(peptide) + MASS["H2O"]
FragBMass = ModpepMass - FragYMass
# print(peptide)
# ipdb.set_trace()
for ficharge in FragCharge:
base_label = f"{eid[0]}{eid[1]}_{frag_index}_{ficharge}"
#如果有两个糖，第二个糖的frag_index不从0开始，从前一个糖后面编，比如前一个糖frag_index最大到了4，后一个从5开始数
for loss_type in ["noloss", "loss_H2O", "loss_NH3", "loss_FUC"]:
loss_mass = {
"noloss": 0,
"loss_H2O": MASS["H2O"],
"loss_NH3": MASS["NH3"],
"loss_FUC": GLYCAN["FUC"] if "Fuc" in glycan else 0
}[loss_type]
# 计算 Y ion
FragTypeY = f"Y{base_label}"
if loss_type != "noloss" and loss_mass > 0:
FragTypeY += f"_{loss_type}"
if loss_mass==0 and loss_type=="loss_FUC":
pass
else:
Y_ions.append({FragTypeY: round((FragYMass - loss_mass) / ficharge + MASS["H+"], 5)})
# 计算 B ion
if loss_type!="loss_FUC":
FragTypeB = f"B{base_label}"
if loss_type != "noloss" and loss_mass > 0:
FragTypeB += f"_{loss_type}"
B_ions.append({FragTypeB: round((FragBMass - loss_mass) / ficharge + MASS["H+"], 5)})
frag_index_offset += len(glycanfrags)
peptide = copy.deepcopy(peptide_copy) # 还原
# print("mark",peptide)
# ipdb.set_trace()
# print(Y_ions)
if "HCD_by" in mode:
#HCD的b,y碎片需要脱糖，或者有一个HexNac
glyco_states = ["(HexNAc)0", "(HexNAc)1"]
glyco_mods = [
(i, mod) for i, mod in enumerate(peptide["modifications"]) if mod["type"] == "glyco"
]
gly_sites = [mod["position"] for _, mod in glyco_mods]
# 所有糖基的脱糖状态组合，如 [(0,0), (0,1), (1,0), (1,1)]
for state_combo in itertools.product(glyco_states, repeat=len(glyco_mods)):
# 修改所有糖基的 name 字段
for (i, _), gly_state in zip(glyco_mods, state_combo):
peptide["modifications"][i]["name"] = gly_state
ModpepMass = calcModpepMass(peptide) + MASS["H2O"]
for FrgNumC in range(1, len(sequence)):
Fragpepb = sequence[:FrgNumC]
Frgmodification = [
mod for mod in peptide["modifications"] if mod["position"] <= FrgNumC - 1
]
Fragb = {"sequence": Fragpepb, "modifications": Frgmodification}
FragbMass = calcModpepMass(Fragb)
FragyMass = ModpepMass - FragbMass
for ficharge in FragCharge:
for loss_type in ["noloss", "loss_H2O", "loss_NH3"]:
loss_mass = {
"noloss": 0,
"loss_H2O": MASS["H2O"],
"loss_NH3": MASS["NH3"]
}[loss_type]
frag_label_b = f"b{FrgNumC}_{ficharge}"
frag_label_y = f"y{len(sequence) - FrgNumC}_{ficharge}"
if loss_type != "noloss":
frag_label_b += f"_{loss_type}"
frag_label_y += f"_{loss_type}"
frag_label_b += "_" + "_".join(state_combo)
frag_label_y += "_" + "_".join(state_combo)
# 判断糖是否要保留
if not (frag_label_b.startswith("b") and any(FrgNumC <= s for s in gly_sites) and "(HexNAc)1" in state_combo):
b_ions.append({frag_label_b: round((FragbMass - loss_mass) / ficharge + MASS["H+"], 5)})
if not (frag_label_y.startswith("y") and any(FrgNumC > s for s in gly_sites) and "(HexNAc)1" in state_combo):
y_ions.append({frag_label_y: round((FragyMass - loss_mass) / ficharge + MASS["H+"], 5)})
#HCD_1 mode意味着肽段全保留，糖部分碎一刀,另外如果有Fuc可以丢Fuc
#B/Y ions糖可以碎裂两刀或者一刀,虽然可以碎三次以上，但是碎太多，搜索空间过大，而且不具有结构区分性
#另外加上单糖
#有很多m/z重复的，需要去重
return b_ions,y_ions,B_ions,Y_ions
# input="AEAVGETLTLPGLVSADJGTYTCEAANK_23.Car._63002_3_(N(F)(N(H(H(H)(H))(H(N(H)(F))))))"
# input="EDGMLPAJR_None_7_3_(N(F)(N(H(H(H))(H))))"
# input="GGJGTICDNQR_7.Car._2_2_(N(F)(N(H(H(N))(H(N)(N)))))"
# mz_calc=pepfragmass(input,["HCD_by"])
# print(mz_calc)
# import ipdb
# ipdb.set_trace()
# print("Fragment MZ for modifided peptides",pepfragmass(peptide))
# input="EDGMLPAJR_None_7_3_(N(F)(N(H(H(H))(H))))"
# input="GGJGTICDNQR_7.Car._2_2_(N(F)(N(H(H(N))(H(N)(N)))))"
# input="TNLPYSQDKAQPGTTNYQHHHHHH_None_0;13_5_(N(A));(N(A)(H(A)))"
# mz_calc=pepfragmass(input,["HCD_1"])
# print(mz_calc)
# import ipdb
# ipdb.set_trace()