"""
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
# --------------------------- 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)
glyco=a[4]
# print("glyco",glyco)
position=int(a[2])
composition={"H":str(glyco.count("H")),"N":str(glyco.count("N")),
"A":str(glyco.count("A")),"G":str(glyco.count("G")),
"F":str(glyco.count("F"))}
comp_str=""
for k in list(composition.keys()):
if composition[k]!="0":
comp_str+="("+node2char[k]+")"
comp_str+=str(composition[k])
mod_dictg["name"]=comp_str
mod_dictg["amino"]=peptide[position]
# print(mod_dictg["amino"])
mod_dictg["position"]=position
mod_dictg["structure"]=glyco
mod_dictg["type"]="glyco"
mod_list.append(mod_dictg)
peptide_dict={"sequence":sequence,"charge":charge,"modifications":mod_list}
# print(peptide)
# print(peptide_dict)
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 ...
if "ETD" in mode: #ETD模式意味着，随着位置糖全部保留，然后碎裂肽段,ETD是cz碎裂注意一下
ModpepMass = calcModpepMass(peptide)+MASS["H2O"]
for FrgNumC in range(1,len(sequence)):
Frgmodification=[]
for modification in peptide["modifications"]:
if modification["position"]<=FrgNumC-1:
Frgmodification.append(modification)
Fragpepb=sequence[:FrgNumC]
Fragb={"sequence":Fragpepb,"modifications":Frgmodification}
#mass = MASS["H2O"] check一下 percursor, b/y ions, B/Y ions的质量
FragbMass = calcModpepMass(Fragb)+MASS["N"]+MASS["H"]
FragyMass=ModpepMass-FragbMass
for ficharge in FragCharge:
for loss_type in ["noloss","loss_H2O","loss_NH3"]:
if loss_type=="noloss":
FragTypeb= "c"+str(FrgNumC)+"_"+str(ficharge)
bion={FragTypeb:FragbMass/ficharge+MASS["H+"]}
b_ions.append(bion)
# print(b_ions)
FragTypey = "z" + str(len(sequence)-FrgNumC) + "_" + str(ficharge)
yion = {FragTypey: FragyMass / ficharge+MASS["H+"]}
y_ions.append(yion)
if loss_type=="loss_H2O":
FragTypeb= "c"+str(FrgNumC)+"_"+str(ficharge)+"_loss_H2O"
bion={FragTypeb:(FragbMass-MASS["H2O"])/ficharge+MASS["H+"]}
b_ions.append(bion)
# print(b_ions)
FragTypey = "z" + str(len(sequence)-FrgNumC) + "_" + str(ficharge)+"_loss_H2O"
yion = {FragTypey: (FragyMass-MASS["H2O"]) / ficharge+MASS["H+"]}
y_ions.append(yion)
if loss_type=="loss_NH3":
FragTypeb= "c"+str(FrgNumC)+"_"+str(ficharge)+"_loss_NH3"
bion={FragTypeb:(FragbMass-MASS["NH3"])/ficharge+MASS["H+"]}
b_ions.append(bion)
# print(b_ions)
FragTypey = "z" + str(len(sequence)-FrgNumC) + "_" + str(ficharge)+"_loss_NH3"
yion = {FragTypey: (FragyMass-MASS["NH3"]) / ficharge+MASS["H+"]}
y_ions.append(yion)
#计算B/Y ions
if "HCD_1" in mode: #HCD_1 mode意味着肽段全保留，糖部分碎一刀,另外如果有Fuc可以丢Fuc
ModpepMass = calcModpepMass(peptide)+MASS["H2O"]
for i in peptide["modifications"]:
if i["type"]=="glyco":
structure=i["structure"]
glyco_graph=glyco_process(structure)
# import ipdb
# ipdb.set_trace()
glycanfrag=GlycanFrag(glyco_graph)
keys=list(glycanfrag.keys())
values=list(glycanfrag.values())
# print(values)
# print(keys)
peptide_copy=copy.deepcopy(peptide)
# import ipdb
# ipdb.set_trace()
# print("FragY",peptide)
# print("FragY2",peptide_copy)
for k in range(len(keys)):
#将k认为是eid的列表，可以考虑是不是对
Glycan=values[k]
# print("Glycan",Glycan)
composition={"H":Glycan.count("Hex"),"N":Glycan.count("HexNAc"),
"A":Glycan.count("NeuAc"),"G":Glycan.count("NeuGc"),"F":Glycan.count("Fuc")}
comp_str=""
for m in list(composition.keys()):
if composition[m]!=0:
comp_str+="("+node2char[m]+")"
comp_str+=str(composition[m])
i["name"]=comp_str #这里改变了peptide,存了peptide的原始备份peptide_copy
# print("comp_str",comp_str)
FragYMass = calcModpepMass(peptide)+MASS["H2O"] #h2o需要核实
# print(FragYMass)
FragBMass=ModpepMass-FragYMass
for ficharge in FragCharge:
for loss_type in ["noloss","loss_H2O","loss_NH3","loss_FUC"]:
if loss_type=="noloss":
FragTypeY = "Y" + str(keys[k][0]) + str(keys[k][1])+"_" + str(k) + "_"+ str(ficharge)
Yion = {FragTypeY: round(FragYMass / ficharge+MASS["H+"],5)}
Y_ions.append(Yion)
FragTypeB= "B"+ str(keys[k][0]) + str(keys[k][1])+"_" + str(k) + "_" + str(ficharge)
Bion = {FragTypeB: round(FragBMass / ficharge+MASS["H+"],5)}
B_ions.append(Bion)
if loss_type=="loss_H2O":
FragTypeY= "Y" + str(keys[k][0]) + str(keys[k][1])+"_" + str(k) + "_" + str(ficharge)+"_loss_H2O"
Yion={FragTypeY: round((FragYMass-MASS["H2O"]) / ficharge+MASS["H+"],5)}
Y_ions.append(Yion)
FragTypeB = "B"+ str(keys[k][0]) + str(keys[k][1])+"_" + str(k) + "_" + str(ficharge)+"_loss_H2O"
Bion = {FragTypeB: round((FragBMass-MASS["H2O"]) / ficharge+MASS["H+"],5)}
B_ions.append(Bion)
if loss_type=="loss_NH3":
FragTypeY= "Y" + str(keys[k][0]) + str(keys[k][1])+"_" + str(k) + "_" + str(ficharge)+"_loss_NH3"
Yion={FragTypeY: round((FragYMass-MASS["NH3"]) / ficharge+MASS["H+"],5)}
Y_ions.append(Yion)
FragTypeB = "B"+ str(keys[k][0]) + str(keys[k][1])+"_" + str(k) + "_" + str(ficharge)+"_loss_NH3"
Bion = {FragTypeB: round((FragBMass-MASS["NH3"]) / ficharge+MASS["H+"],5)}
B_ions.append(Bion)
if len(Glycan)>1 and Glycan[1]== "Fuc" and loss_type=="loss_FUC":
FragTypeY= "Y" + str(keys[k][0]) + str(keys[k][1])+"_" + str(k) + "_" + str(ficharge)+"_loss_FUC"
Yion={FragTypeY: round((FragYMass-GLYCAN["FUC"]) / ficharge+MASS["H+"],5)}
Y_ions.append(Yion)
peptide=peptide_copy
# print(Y_ions)
if "HCD_by" in mode:
#HCD的b,y碎片需要脱糖，或者有一个HexNac
copy_num=0
for gly_HCD in ["(HexNAc)0","(HexNAc)1"]:
for modification in peptide["modifications"]:
if modification["type"]=="glyco":
if copy_num==0:
glyco_mod=copy.deepcopy(modification["name"])
glysite=modification["position"]
# print("glyco_mod",glyco_mod)
copy_num+=1
# assert modification["structure"][1] =="N"
modification["name"]=gly_HCD
glyco_aftermod=copy.deepcopy(modification)
ModpepMass = calcModpepMass(peptide)+MASS["H2O"]
for FrgNumC in range(1,len(sequence)):
Frgmodification=[]
for modification in peptide["modifications"]:
if modification["position"]<=FrgNumC-1:
Frgmodification.append(modification)
Fragpepb=sequence[:FrgNumC]
Fragb={"sequence":Fragpepb,"modifications":Frgmodification}
#mass = MASS["H2O"] check一下 percursor, b/y ions, B/Y ions的质量
FragbMass = calcModpepMass(Fragb)
FragyMass=ModpepMass-FragbMass
for ficharge in FragCharge:
for loss_type in ["noloss","loss_H2O","loss_NH3"]:
if loss_type=="noloss":
FragTypeb= "b"+str(FrgNumC)+"_"+str(ficharge)+"_"+gly_HCD
if FragTypeb.startswith("b") and int(FrgNumC)<=int(glysite) and gly_HCD=="(HexNAc)1":
pass
else:
bion={FragTypeb:round(FragbMass/ficharge+MASS["H+"],5)}
b_ions.append(bion)
# print(b_ions)
# import ipdb
# ipdb.set_trace()
FragTypey = "y" + str(length-FrgNumC) + "_" + str(ficharge)+"_"+gly_HCD
if FragTypey.startswith("y") and FrgNumC>int(glysite) and gly_HCD=="(HexNAc)1":
pass
else:
yion = {FragTypey: round(FragyMass / ficharge+MASS["H+"],5)}
y_ions.append(yion)
if loss_type=="loss_H2O":
FragTypeb= "b"+str(FrgNumC)+"_"+str(ficharge)+"_loss_H2O_"+gly_HCD
if FragTypeb.startswith("b") and int(FrgNumC)<=int(glysite) and gly_HCD=="(HexNAc)1":
pass
else:
bion={FragTypeb:round((FragbMass-MASS["H2O"])/ficharge+MASS["H+"],5)}
b_ions.append(bion)
# print(b_ions)
FragTypey = "y" + str(len(sequence)-FrgNumC) + "_" + str(ficharge)+"_loss_H2O_"+gly_HCD
if FragTypey.startswith("y") and FrgNumC>int(glysite) and gly_HCD=="(HexNAc)1":
pass
else:
yion = {FragTypey: round((FragyMass-MASS["H2O"]) / ficharge+MASS["H+"],5)}
y_ions.append(yion)
if loss_type=="loss_NH3":
FragTypeb= "b"+str(FrgNumC)+"_"+str(ficharge)+"_loss_NH3_"+gly_HCD
if FragTypeb.startswith("b") and int(FrgNumC)<=int(glysite) and gly_HCD=="(HexNAc)1":
pass
else:
bion={FragTypeb:round((FragbMass-MASS["NH3"])/ficharge+MASS["H+"],5)}
b_ions.append(bion)
# print(b_ions)
FragTypey = "y" + str(len(sequence)-FrgNumC) + "_" + str(ficharge)+"_loss_NH3_"+gly_HCD
if FragTypey.startswith("y") and FrgNumC>int(glysite) and gly_HCD=="(HexNAc)1":
pass
else:
yion = {FragTypey: round((FragyMass-MASS["NH3"]) / ficharge+MASS["H+"],5)}
y_ions.append(yion)
#HCD_1 mode意味着肽段全保留，糖部分碎一刀,另外如果有Fuc可以丢Fuc
#B/Y ions糖可以碎裂两刀或者一刀,虽然可以碎三次以上，但是碎太多，搜索空间过大，而且不具有结构区分性
#另外加上单糖
#有很多m/z重复的，需要去重
if "HCD_BY_2" in mode:
ModpepMass = calcModpepMass(peptide)+MASS["H2O"]
for i in peptide["modifications"]:
if i["type"]=="glyco":
structure=i["structure"]
# ipdb.set_trace()
p = re.compile(r'[()](.*?)[)]', re.S)
# ipdb.set_trace()
glycan_all=re.findall(p,i["name"])
# print("glycan_all",glycan_all)
# ipdb.set_trace()
# FragCharge=list(range(1,min((int(peptide["charge"])+1),4)))
for mono_sugar in glycan_all:
for ficharge in FragCharge:
for addition_type in ["","_add_H2O","_loss_H2O"]:
FragTypeB = "B_" + str(ficharge)+"_"+mono_sugar+addition_type
if addition_type =="":
FragBMass=GLYCAN[mono_sugar.upper()]
if addition_type =="_add_H2O":
FragBMass=GLYCAN[mono_sugar.upper()]+MASS["H2O"]
if addition_type =="_loss_H2O":
FragBMass=GLYCAN[mono_sugar.upper()]-MASS["H2O"]
Bion = {FragTypeB: FragBMass / ficharge+MASS["H+"]}
B_ions.append(Bion)
# ipdb.set_trace()
glyco_graph=glyco_process(structure)
glycanfrag=GlycanFrag_HCD(glyco_graph)
keys=list(glycanfrag.keys())
values=list(glycanfrag.values())
# print(values)
# print(keys)
# ipdb.set_trace()
peptide_copy=copy.deepcopy(peptide)
for k in range(len(keys)):
Glycan=values[k]
# print("Glycan",Glycan)
# print("bond",keys[k])
# ipdb.set_trace()
frg_num=0
for res in Glycan:
res_index=res["res_index"]
if 0 in res_index: #Y_ions
frg_num+=1
FragNumY=keys[k]+"_"+str(frg_num)
# print("FragNumY",FragNumY)
res_sugar=res["res_sugar"]
# print("res_sugar",res_sugar)
composition={"H":res_sugar.count("Hex"),"N":res_sugar.count("HexNAc"),
"A":res_sugar.count("NeuAc"),"G":res_sugar.count("NeuGc"),"F":res_sugar.count("Fuc")}
comp_str=""
for m in list(composition.keys()):
if composition[m]!=0:
comp_str+="("+node2char[m]+")"
comp_str+=str(composition[m])
i["name"]=comp_str #这里改变了peptide,存了peptide的原始备份peptide_copy
# ipdb.set_trace()
# print("comp_str",comp_str)
FragYMass = calcModpepMass(peptide)+MASS["H2O"] #h2o需要核实
# print(FragYMass)
for ficharge in FragCharge:
for loss_type in ["noloss","loss_H2O","loss_NH3","loss_FUC"]:
if loss_type=="noloss":
FragTypeY = "Y" + FragNumY+"_" + str(ficharge)
Yion = {FragTypeY: round(FragYMass / ficharge+MASS["H+"],5)}
Y_ions.append(Yion)
if loss_type=="loss_H2O":
FragTypeY= "Y" + FragNumY+"_" + str(ficharge)+"_loss_H2O"
Yion={FragTypeY: round((FragYMass-MASS["H2O"]) / ficharge+MASS["H+"],5)}
Y_ions.append(Yion)
if loss_type=="loss_NH3":
FragTypeY= "Y" + FragNumY+"_" + str(ficharge)+"_loss_NH3"
Yion={FragTypeY: round((FragYMass-MASS["NH3"]) / ficharge+MASS["H+"],5)}
Y_ions.append(Yion)
if len(res_sugar)>1 and res_sugar[1]== "Fuc" and loss_type=="loss_FUC":
FragTypeY= "Y" + FragNumY+"_" + str(ficharge) + "_" + str(ficharge)+"_loss_FUC"
Yion={FragTypeY: round((FragYMass-GLYCAN["FUC"]) / ficharge+MASS["H+"],5)}
Y_ions.append(Yion)
else:
frg_num+=1
FragNumB=keys[k]+"_"+str(frg_num)
# print(FragNumB)
res_sugar=res["res_sugar"]
composition={"H":res_sugar.count("Hex"),"N":res_sugar.count("HexNAc"),
"A":res_sugar.count("NeuAc"),"G":res_sugar.count("NeuGc"),"F":res_sugar.count("Fuc")}
FragBMass=0
for m in list(composition.keys()):
if composition[m]!=0:
# print(node2char[m].upper())
FragBMass+=GLYCAN[node2char[m].upper()]*composition[m]
# print("FragBMass",FragBMass)
for ficharge in FragCharge:
for loss_type in ["noloss","loss_H2O","loss_NH3","loss_FUC"]:
if loss_type=="noloss":
FragTypeB= "B"+ FragNumB+"_" + str(ficharge)
Bion = {FragTypeB: round(FragBMass / ficharge+MASS["H+"],5)}
B_ions.append(Bion)
if loss_type=="loss_H2O":
FragTypeB = "B"+ FragNumB+"_" + str(ficharge)+"_loss_H2O"
Bion = {FragTypeB: round((FragBMass-MASS["H2O"]) / ficharge+MASS["H+"],5)}
B_ions.append(Bion)
if loss_type=="loss_NH3":
if "HexNAc" in res_sugar or "NEUAC" in res_sugar or "NEUGC" in res_sugar:
# print("res_sugar",res_sugar)
FragTypeB = "B"+ FragNumB+"_" + str(ficharge)+"_loss_NH3"
Bion = {FragTypeB: round((FragBMass-MASS["NH3"]) / ficharge+MASS["H+"],5)}
B_ions.append(Bion)
if len(res_sugar)>1 and res_sugar[1]== "Fuc" and loss_type=="loss_FUC":
# print("res_sugar",res_sugar)
FragTypeB = "B"+ FragNumB+"_" + str(ficharge)+"_loss_FUC"
Bion = {FragTypeB: round((FragBMass-GLYCAN["FUC"]) / ficharge+MASS["H+"],5)}
B_ions.append(Bion)
peptide=peptide_copy
# print(Y_ions)
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))