#对预测得到的结果转化为谱图的形式 #目前改了by和BY，exp暂时还没改 import pandas as pd import numpy as np import ipdb from utils import * import torch import masses from masses import * # ----------------------- dict ------------------------------# MAXFICHARGE = int(num_col/12) if MAXFICHARGE == 2: ionname_by = "b1,b1n,b1o,b1h,b1nh,b1oh,y1,y1n,y1o,y1h,y1nh,y1oh,"\ "b2,b2n,b2o,b2h,b2nh,b2oh,y2,y2n,y2o,y2h,y2nh,y2oh".split( ',') ionname_BY = "B1,B1n,B1o,B1f,Y1,Y1n,Y1o,Y1f,"\ "B2,B2n,B2o,B2f,Y2,Y2n,Y2o,Y2f".split( ',') if MAXFICHARGE == 3: ionname_by = "b1,b1n,b1o,b1h,b1nh,b1oh,y1,y1n,y1o,y1h,y1nh,y1oh,"\ "b2,b2n,b2o,b2h,b2nh,b2oh,y2,y2n,y2o,y2h,y2nh,y2oh,"\ "b3,b3n,b3o,b3h,b3nh,b3oh,y3,y3n,y3o,y3h,y3nh,y3oh".split( ',') ionname_BY = "B1,B1n,B1o,B1f,Y1,Y1n,Y1o,Y1f,"\ "B2,B2n,B2o,B2f,Y2,Y2n,Y2o,Y2f,"\ "B3,B3n,B3o,B3f,Y3,Y3n,Y3o,Y3f".split( ',') lossdict_by = {"noloss": '', "H2O": 'o', "NH3": 'n', '(HexNAc)1': 'h', '1(NH3)1(+HexNAc)': 'nh', '1(H2O)1(+HexNAc)': 'oh'} lossdict_BY = {"noloss": '', "H2O": 'o', "NH3": 'n', 'FUC': 'f'} fields_byBY = "SourceFile,id,PEP.StrippedSequence,iden_pep,PlausibleStruct,PP.Charge,metric,FragmentMz,FI.FrgType,FI.LossType,FI.Charge,FI.Intensity,FI.FrgNum,Precurmass_cal".split(',') fields_exp = "SourceFile,id,PEP.StrippedSequence,iden_pep,PlausibleStruct,PP.Charge,FragmentMz,FI.FrgType,FI.LossType,FI.Charge,FI.Intensity,FI.FrgNum,Precurmass_cal".split(',') # ----------------------- func ------------------------------# def byBYextract(row): outputdataframe1 = pd.DataFrame(columns=fields_byBY) # StrippedSequence = "".join([vocab.to_word(i) for i in eval(row['repsequence'])]) StrippedSequence=row["ipeptide"] PPCharge=int(row["charge"]) metric=float(row["metric"]) id=row["id"] global id_glyspec global id_iden_pep glyspec=id_glyspec[id] ms2by=row["ms2by"] ms2BY=row["ms2BY"] iden_pep=id_iden_pep[id] # import ipdb # ipdb.set_trace() PlausibleStruct=iden_pep.split("_")[-1] glysite=iden_pep.split("_")[2] # Precurmass_exp=glyspec_pepmass_mgf[glyspec] Precurmass_cal= (calcModpepMass(peptide_process(iden_pep))+MASS["H2O"])/PPCharge+MASS["H+"]#根据iden_pep计算母离子质荷比 # deltamz=abs(Precurmass_exp-Precurmass_cal)/Precurmass_exp FRAG_MODE=["HCD_by","HCD_1"] mz_calc=masses.pepfragmass(iden_pep,FRAG_MODE) #得到iontype和m/z对 # ipdb.set_trace() #by jmax=ms2by.size()[1] imax=ms2by.size()[0] # ipdb.set_trace() assert jmax==num_col for j in range(0,jmax): for i in range(0,imax): intensity_pred=ms2by[i][j] #从ms2中提取对应位置的强度 if intensity_pred <=0: continue intensity_pred=float(intensity_pred) ionname=ionname_by[j] #从ionname列表中获得对应结果的离子类型 FIFrgType=ionname[0] FICharge=int(ionname[1]) FILossType=ionname[2:] FILossType=list(lossdict_by.keys())[list(lossdict_by.values()).index(FILossType)] assert FIFrgType in ["b","y"] # ipdb.set_trace() if FIFrgType=="b": FragmentNumber=i+1 masskey="b" + str(FragmentNumber) + "_" + str(FICharge) if FIFrgType=="y": FragmentNumber=imax - i masskey="y" + str(FragmentNumber) + "_" + str(FICharge) # import ipdb if FIFrgType=="b" and int(FragmentNumber)<=int(glysite) and ionname.endswith("h"): continue if FIFrgType=="y" and int(i+1)>int(glysite) and ionname.endswith("h"): continue if FILossType!="noloss": loss=ionname[2] if loss !="h": loss=list(lossdict_by.keys())[list(lossdict_by.values()).index(loss)] masskey+="_loss_"+FILossType if ionname.endswith("h"): masskey+="_"+"(HexNAc)1" # import ipdb # ipdb.set_trace() else: masskey+="_"+"(HexNAc)0" counter=0 for k in mz_calc: for dict in k: key=[str(key) for key in dict.keys()] assert len(key)==1 if masskey == key[0]: value=[float(value) for value in dict.values()] assert len(value)==1 FragmentMz=value[0] counter+=1 if counter>1: raise ValueError("loop number more than one") # ipdb.set_trace() if counter==1: outputdataframe1.loc[len(outputdataframe1)] = [glyspec,id, StrippedSequence,iden_pep,PlausibleStruct,\ PPCharge,metric,FragmentMz,\ FIFrgType, FILossType,FICharge,intensity_pred,\ FragmentNumber,Precurmass_cal] #BY jmax=ms2BY.size()[1] imax=ms2BY.size()[0] assert jmax==num_col*2/3 for j in range(0,jmax): for i in range(0,imax): intensity_pred=ms2BY[i][j] #从ms2中提取对应位置的强度 if intensity_pred <=0: continue intensity_pred=float(intensity_pred) ionname=ionname_BY[j] #从ionname列表中获得对应结果的离子类型 FIFrgType=ionname[0] FICharge=int(ionname[1]) FILossType=ionname[2:] if FILossType=="f" and "F" not in PlausibleStruct: continue if FILossType=="f" and FIFrgType=="B": continue FILossType=list(lossdict_BY.keys())[list(lossdict_BY.values()).index(FILossType)] FragmentNumber=i #"_" + str(k) + "_" + str(ficharge)+"_loss_H2O" masskey="_" + str(FragmentNumber) + "_" + str(FICharge) if FILossType!="noloss": masskey+="_loss_"+FILossType counter=0 # if FILossType=="FUC": # ipdb.set_trace() for k in mz_calc: for dict in k: key=[str(key) for key in dict.keys()] assert len(key)==1 if FIFrgType == key[0][0]: if key[0].endswith(masskey): value=[float(value) for value in dict.values()] assert len(value)==1 FragmentMz=value[0] counter+=1 if counter>1: raise ValueError("loop number more than one") # import ipdb # ipdb.set_trace() if counter==1: outputdataframe1.loc[len(outputdataframe1)] = [glyspec,id, StrippedSequence,iden_pep,PlausibleStruct,\ PPCharge,metric,FragmentMz,\ FIFrgType, FILossType,FICharge,intensity_pred,\ FragmentNumber,Precurmass_cal] # import ipdb # ipdb.set_trace() return outputdataframe1 def expbyBYextract(row): outputdataframe1 = pd.DataFrame(columns=fields_exp) StrippedSequence = row["peptide"] PPCharge=int(row["charge"]) id=row["_id"] glyspec=row["GlySpec"] # import ipdb # ipdb.set_trace() ms2by=row["ions_by"] ms2BY=row["ions_BY"] # iden_pep=id_iden_pep[id] # import ipdb # ipdb.set_trace() iden_pep=row["iden_pep"] PlausibleStruct=iden_pep.split("_")[-1] glysite=iden_pep.split("_")[2] # Precurmass_exp=glyspec_pepmass_mgf[glyspec] Precurmass_cal= (calcModpepMass(peptide_process(iden_pep))+MASS["H2O"])/PPCharge+MASS["H+"]#根据iden_pep计算母离子质荷比 # deltamz=abs(Precurmass_exp-Precurmass_cal)/Precurmass_exp FRAG_MODE=["HCD_by","HCD_1"] mz_calc=masses.pepfragmass(iden_pep,FRAG_MODE) #得到iontype和m/z对 # ipdb.set_trace() #by jmax=ms2by.size()[1] imax=ms2by.size()[0] # ipdb.set_trace() assert jmax==num_col for j in range(0,jmax): for i in range(0,imax): intensity_pred=ms2by[i][j] #从ms2中提取对应位置的强度 if intensity_pred <=0: continue intensity_pred=float(intensity_pred) ionname=ionname_by[j] #从ionname列表中获得对应结果的离子类型 FIFrgType=ionname[0] FICharge=int(ionname[1]) FILossType=ionname[2:] FILossType=list(lossdict_by.keys())[list(lossdict_by.values()).index(FILossType)] assert FIFrgType in ["b","y"] # ipdb.set_trace() if FIFrgType=="b": FragmentNumber=i+1 masskey="b" + str(FragmentNumber) + "_" + str(FICharge) if FIFrgType=="y": FragmentNumber=imax - i masskey="y" + str(FragmentNumber) + "_" + str(FICharge) # import ipdb if FIFrgType=="b" and int(FragmentNumber)<=int(glysite) and ionname.endswith("h"): continue if FIFrgType=="y" and int(i+1)>int(glysite) and ionname.endswith("h"): continue if FILossType!="noloss": loss=ionname[2] if loss !="h": loss=list(lossdict_by.keys())[list(lossdict_by.values()).index(loss)] masskey+="_loss_"+FILossType if ionname.endswith("h"): masskey+="_"+"(HexNAc)1" # import ipdb # ipdb.set_trace() else: masskey+="_"+"(HexNAc)0" counter=0 for k in mz_calc: for dict in k: key=[str(key) for key in dict.keys()] assert len(key)==1 if masskey == key[0]: value=[float(value) for value in dict.values()] assert len(value)==1 FragmentMz=value[0] counter+=1 if counter>1: raise ValueError("loop number more than one") # ipdb.set_trace() if counter==1: outputdataframe1.loc[len(outputdataframe1)] = [glyspec,id, StrippedSequence,iden_pep,PlausibleStruct,\ PPCharge,FragmentMz,\ FIFrgType, FILossType,FICharge,intensity_pred,\ FragmentNumber,Precurmass_cal] #BY jmax=ms2BY.size()[1] imax=ms2BY.size()[0] assert jmax==num_col*2/3 for j in range(0,jmax): for i in range(0,imax): intensity_pred=ms2BY[i][j] #从ms2中提取对应位置的强度 if intensity_pred <=0: continue intensity_pred=float(intensity_pred) ionname=ionname_BY[j] #从ionname列表中获得对应结果的离子类型 FIFrgType=ionname[0] FICharge=int(ionname[1]) FILossType=ionname[2:] if FILossType=="f" and "F" not in PlausibleStruct: continue if FILossType=="f" and FIFrgType=="B": continue FILossType=list(lossdict_BY.keys())[list(lossdict_BY.values()).index(FILossType)] FragmentNumber=i #"_" + str(k) + "_" + str(ficharge)+"_loss_H2O" masskey="_" + str(FragmentNumber) + "_" + str(FICharge) if FILossType!="noloss": masskey+="_loss_"+FILossType counter=0 # if FILossType=="FUC": # ipdb.set_trace() for k in mz_calc: for dict in k: key=[str(key) for key in dict.keys()] assert len(key)==1 if FIFrgType == key[0][0]: if key[0].endswith(masskey): value=[float(value) for value in dict.values()] assert len(value)==1 FragmentMz=value[0] counter+=1 if counter>1: raise ValueError("loop number more than one") # import ipdb # ipdb.set_trace() if counter==1: outputdataframe1.loc[len(outputdataframe1)] = [glyspec,id, StrippedSequence,iden_pep,PlausibleStruct,\ PPCharge,FragmentMz,\ FIFrgType, FILossType,FICharge,intensity_pred,\ FragmentNumber,Precurmass_cal] # import ipdb # ipdb.set_trace() return outputdataframe1 def byBYextract_reidentification(row): outputdataframe1 = pd.DataFrame(columns=["GlySpec","modelpep","FragmentMz",\ "FIFrgType", "FILossType","FICharge","intensity_pred",\ "FragmentNumber","FragmentGlycan"]) import ipdb GlySpec=row["GlySpec"] FRAG_MODE=["HCD_1"] # ipdb.set_trace() #BY for num in ["1","2"]: modelpep=row["modelpep"+num] ms2BY=row["maskmodel"+num] mz_calc=masses.pepfragmass(modelpep,FRAG_MODE) #得到iontype和m/z对 # ipdb.set_trace() jmax=ms2BY.size()[1] imax=ms2BY.size()[0] assert jmax==num_col*2/3 # ipdb.set_trace() for j in range(0,jmax): for i in range(0,imax): intensity_pred=ms2BY[i][j] #从ms2中提取对应位置的强度 if intensity_pred <=0: continue intensity_pred=float(intensity_pred) ionname=ionname_BY[j] #从ionname列表中获得对应结果的离子类型 FIFrgType=ionname[0] FICharge=int(ionname[1]) FILossType=ionname[2:] FILossType=list(lossdict_BY.keys())[list(lossdict_BY.values()).index(FILossType)] FragmentNumber=i masskey="_" + str(FragmentNumber) + "_" + str(FICharge) if FILossType!="noloss": masskey+="_loss_"+FILossType counter=0 # ipdb.set_trace() for k in mz_calc: for dict in k: key=[str(key) for key in dict.keys()] assert len(key)==1 if FIFrgType == key[0][0]: if key[0].endswith(masskey): value=[float(value) for value in dict.values()] assert len(value)==1 FragmentMz=value[0] fragnode=str(key[0].split("_")[0][1:]) # ipdb.set_trace() counter+=1 if counter>1: raise ValueError("loop number more than one") # import ipdb # ipdb.set_trace() if counter==1: outputdataframe1.loc[len(outputdataframe1)] = [GlySpec,modelpep,FragmentMz,\ FIFrgType, FILossType,FICharge,intensity_pred,\ FragmentNumber,fragnode] # import ipdb # ipdb.set_trace() return outputdataframe1 # ----------------------- execute ------------------------------# def postprocessing(savename): filename=savename+"_byBY_result.csv" df=pd.read_csv(filename) # import ipdb # ipdb.set_trace() # df=df.loc[df["metric"]==max(df["metric"])] #取一个进行测试 # df.to_csv(datafold+"del.csv",index=False) print(df.columns) df["ms2by"]=df["ms2by"].apply(lambda x:torch.Tensor(eval(x))) df["ms2BY"]=df["ms2BY"].apply(lambda x:torch.Tensor(eval(x))) #训练的时候ms2使用了log2（x+1)，这里变回来#检查一下 #'repsequence', 'charge', 'ipeptide', 'iPlausibleStruct', 'ms2', 'cos','id' print(f"The number of the lines of the input files {len(df)}") filename_expjson=savename+"_byBY_testdata.json" dfexp=pd.read_json(filename_expjson) # dfexp=dfexp.loc[dfexp["_id"]==2665] # dfexp.to_csv(datafold+"exp_6923_del.csv",index=False) print(dfexp.columns) # import ipdb # ipdb.set_trace() global id_glyspec global id_iden_pep id_glyspec = dict(zip(dfexp['_id'],dfexp['GlySpec'])) #产生_id,GlySpec字典 id_iden_pep=dict(zip(dfexp['_id'],dfexp['iden_pep'])) # import ipdb # ipdb.set_trace() # dfinput=pd.read_csv(dfinput) # ipdb.set_trace() outputdataframe_BY=pd.concat(df.apply(byBYextract,axis=1).to_list()) outputdataframe_BY.reset_index(inplace=True,drop=True) print(f"The number of the lines of the output files {len(outputdataframe_BY)}") outputdataframe_BY.to_csv(savename+"_predictmonomz.csv",index=False) # import ipdb # ipdb.set_trace() dfexp["ions_by"]=dfexp["ions_by"].apply(lambda x:torch.Tensor(x)) dfexp["ions_BY"]=dfexp["ions_BY"].apply(lambda x:torch.Tensor(x)) print(dfexp.columns) # dfexp.to_csv(datafold+"exp_6923_del.csv",index=False print(f"The number of the lines of the input files {len(dfexp)}") # outputdataframe_expBY = pd.DataFrame(columns=fields_exp) # import ipdb # ipdb.set_trace() outputdataframe_expBY=pd.concat(dfexp.apply(expbyBYextract,axis=1).to_list()) outputdataframe_expBY.reset_index(inplace=True,drop=True) print(f"The number of the lines of the output files {len(outputdataframe_expBY)}") outputdataframe_expBY.to_csv(savename+"_expmonomz.csv",index=False) def Precurmass_cal(innstance): iden_pep=innstance["iden_pep"] PPCharge=int(innstance["charge"]) Precurmass_cal=(calcModpepMass(peptide_process(iden_pep))+MASS["H2O"])/PPCharge+MASS["H+"] return Precurmass_cal def postprocessing_nointensity(savename,dfinputname): filename=savename+"_byBY_result.csv" df=pd.read_csv(filename) print(df.columns) print(f"The number of the lines of the input files {len(df)}") filename_expjson=savename+"_byBY_testdata.json" dfexp=pd.read_json(filename_expjson) print(dfexp.columns) print(f"The number of the lines of the input exp files {len(dfexp)}") # import ipdb # ipdb.set_trace() dfexp=dfexp[['_id','GlySpec',"iden_pep"]] df=df[[ 'charge', 'ipeptide', 'iPlausibleStruct','metric', "metricBY_cos",'id']] df.columns=['charge', 'PEP.StrippedSequence', 'PlausibleStruct','metric',"metricBY", '_id'] df=pd.merge(df,dfexp,on="_id",how="left") # import ipdb # ipdb.set_trace() dfinput=pd.read_csv(dfinputname) # ipdb.set_trace() dfinput=dfinput[['GlySpec','pepmass_mgf',"matching_ration","iden_pep"]].drop_duplicates() print(len(dfinput)) df=pd.merge(df,dfinput,on=["GlySpec","iden_pep"],how="left") # ipdb.set_trace() df["Precurmass_cal"]= df.apply(Precurmass_cal,axis=1) df["deltamzpercent"]=abs(df["pepmass_mgf"]-df["Precurmass_cal"])/df["Precurmass_cal"] df.to_csv(savename+"_nointensity_predictmonomz.csv",index=False) # if __name__ == '__main__': # savename="/remote-home/yxwang/test/zzb/DeepGlyco/DeepSweet_v1/data/human/IgG/PXD015360_IgG/test_replace_predict/test_byBY_PXD015360_IgG_data_1st_target_decoy_model" # dfinputname="/remote-home/yxwang/test/zzb/DeepGlyco/DeepSweet_v1/data/human/IgG/PXD015360_IgG/PXD015360_IgG_data_1st_redo1_nocoelute_target_decoy_Retained_all.csv" # postprocessing_nointensity(savename,dfinputname) # savename="/remote-home/yxwang/test/zzb/DeepGlyco/DeepSweet_v1/data/human/PXD009654/test_replace_predict/test_byBY_PXD009654_human_data_1st_target_decoy_model" # dfinputname="/remote-home/yxwang/test/zzb/DeepGlyco/DeepSweet_v1/data/human/PXD009654/PXD009654_data_1st_redo1_nocoelute_target_decoy_Retained_all.csv" # postprocessing_nointensity(savename,dfinputname) # savename="/remote-home/yxwang/test/zzb/DeepGlyco/DeepSweet_v1/data/human/IgG/PXD009716/test_replace_predict/PXD009716_IgG_data_1st_target_decoy_NAGF_model_0.930816" # dfinputname="/remote-home/yxwang/test/zzb/DeepGlyco/DeepSweet_v1/data/human/IgG/PXD009716/PXD009716_IgG_data_1st_target_decoy_Retained_all.csv" # postprocessing_nointensity(savename,dfinputname) # savename="/remote-home/yxwang/test/zzb/DeepGlyco/DeepSweet_v1/data/syn/PXD023980/test_replace_predict/test_byBY_PXD023980_syn_test_9_glycans_redo1_model_0.927153" # dfinputname="/remote-home/yxwang/test/zzb/DeepGlyco/DeepSweet_v1/data/syn/PXD023980/PXD023980_Syn_test_9_glycans_redo1_data_1st.csv" # postprocessing_nointensity(savename,dfinputname) # savename="/remote-home/yxwang/test/zzb/DeepGlyco/DeepSweet_v1/data/syn/PXD023980/test_replace_predict/test_byBY_PXD023980_syn_identification_task3" # dfinputname="/remote-home/yxwang/test/zzb/DeepGlyco/DeepSweet_v1/data/syn/PXD023980/Syn_PXD023980_task3_data_1st.csv" # postprocessing_nointensity(savename,dfinputname) # savename="/remote-home/yxwang/test/zzb/DeepGlyco/DeepSweet_v1/data/human/PXD016428/test_replace_predict/test_byBY_PXD016428_IgG_data_1st_target_decoy_model_0.933599" # dfinputname="/remote-home/yxwang/test/zzb/DeepGlyco/DeepSweet_v1/data/human/PXD016428/PXD016428_serum_data_1st_target_decoy_Retained_all.csv" # postprocessing_nointensity(savename,dfinputname) # savename="/remote-home/yxwang/test/zzb/DeepGlyco/DeepSweet_v1/data/human/PXD016428/test_replace_predict/test_byBY_PXD016428_IgG_data_1st_target_decoy_model_0.937599" # dfinputname="/remote-home/yxwang/test/zzb/DeepGlyco/DeepSweet_v1/data/human/PXD016428/PXD016428_serum_data_1st_target_decoy_Retained_all.csv" # postprocessing_nointensity(savename,dfinputname) # savename="/remote-home/yxwang/test/zzb/DeepGlyco/DeepSweet_v1/data/mouse/Five_tissues/test_replace_predict/test_byBY_PXD005411_model_0.927740" # postprocessing(savename)