DeepGPO / DeepGPO_code / 1_dataset_format_NO.py
1_dataset_format_NO.py
Raw
import pandas as pd
import os
import numpy as np
from pathlib import Path
import masses
import mgf_processing
from weights import *
# --------------------------- argparse ---------------------#
import argparse
def parsering():
    parser = argparse.ArgumentParser()
    # Training parameter
    parser.add_argument('--datafold', type=str, 
                        default="/remote-home1/yxwang/test/zzb/DeepGlyco/DeepSweet_v1/data/mouse/PXD005411/",
                        help='datafold ')
    parser.add_argument('--dfname', type=str, 
                        default="pGlycoDB-GP-FDR-Pro_PXD005411.txt",
                        help='pglyco3 crude result ')
    parser.add_argument('--mgfdatafold', type=str, 
                        default="MSConvert_mgf_PXD005411/" , 
                        help='mgf data fold')
    parser.add_argument('--output_name', type=str, 
                        default="PXD005411_MouseBrain_data_1st.csv", help='outputfile name')
    parser.add_argument('--only_duplicated', type=str,default="Drop_duplicated", help='Duplicated/Drop_duplicated/Retained_all')
    parser.add_argument('--mgfsourceorign', type=str,default="pGlyco3", help='Please ensure the tool for producing mgf (MsConvert or pGlyco3)')
    parser.add_argument('--fragmentation', type=str,default="HCD", help='HCD/EThCD/ETD')
    parser.add_argument('--enzyme', type=str,default="None", help='protease used')
    parser.add_argument('--filter_jsonname', type=str,default="SStruespectrum_filtered_O_", help='')
    parser.add_argument('--not_use_weights', action='store_true', help='')
    args = parser.parse_args()
    return args
args=parsering()
DFNAME=args.datafold+args.dfname
mgfdatafold=args.datafold+args.mgfdatafold
output_name=args.datafold+args.output_name
only_duplicated=args.only_duplicated
mgfsourceorign=args.mgfsourceorign
assert mgfsourceorign in ["pGlyco3","MsConvert"], "mgfsourceorign not in [pGlyco3,MsConvert]"
fragmentation=args.fragmentation
assert fragmentation in ["HCD","ETD","EThCD"], "fragmentation not in [HCD,ETD,EThCD]"
Enzyme=args.enzyme
if args.not_use_weights:
    # Code when not using weights
    print("Not using weights")
else:
    # Code when using weights
    print("Using weights")
# --------------------------- hyper paramaters ---------------------#
FRAG_AVA=["ETD","HCD_1","HCD_by","HCD_BY_2"]
if fragmentation=="HCD":
    FRAG_INDEX=[1,2] #"HCD_1" for BY prediction. "HCD_by" for by prediction
if fragmentation=="ETD":
    FRAG_INDEX=[0]
if fragmentation=="EThCD":
    FRAG_INDEX=[0,1,2]
FRAG_MODE=[x for x in FRAG_AVA if FRAG_AVA.index(x) in FRAG_INDEX]
print(f"FRAG_MODE: {FRAG_MODE}")
jsonfold= os.path.join(mgfdatafold, "json/")
jsonname="SStruespectrum.json"
filter_jsonname=args.filter_jsonname+args.only_duplicated+".json" #相比于SStruespectrum.json提取数据中有的scan,来减少搜索范围
TOLER=20
# --------------------------- pglyco3 result processing---------------------#
def pglyco3_result(DFNAME):
    df=pd.read_csv(DFNAME,sep="\t")
    df.reset_index(inplace=True,drop=True)
    df_column=list(df.columns)
    print(f"Columns of df {df_column}",end = "\n\n")
    print(f"df rank should all be 1.. Please check!!: {list(df['Rank'].drop_duplicates())}",end = "\n\n")
    assert list(df['Rank'].drop_duplicates())==[1]
    df["Peptide"]=df["Peptide"].str.replace("J","N")
    if "LocalizedSiteGroups" in df_column:
        if args.not_use_weights:
            use_weights=False
        else:
            use_weights=True
        print("weight is considered: ", use_weights)
    else:
        use_weights=False
    if use_weights:
        df=df[["RawName","Scan", 'Charge',"Peptide","Mod",
           "PlausibleStruct",'GlySite',"RT","PrecursorMZ",'TotalFDR',"LocalizedSiteGroups"]]
    else:
        df=df[["RawName","Scan", 'Charge',"Peptide","Mod",
           "PlausibleStruct",'GlySite',"RT","PrecursorMZ",'TotalFDR']]
    print(f"Row number of df {len(df)}",end = "\n\n")
    df.drop_duplicates(inplace=True)
    print(f"Row number of df after drop_duplicates {len(df)}",end = "\n\n")
    return df,use_weights

def combine_iden_pep(instance):
    a=instance["Peptide"]
    b=instance["Mod"]
    e=""
    if not pd.isna(b):
        b=b.rstrip(";")
        for i in b.split(";"):
            for k in i.split(","):
                k=k[:3]+"."
                e+=k
        b=e
    else:
        b=None
    c=instance["GlySite"]-1  #GlySite 是从1开始的,会比index J 大一
    d=instance["Charge"]
    e=instance["PlausibleStruct"]
    return str(a)+"_"+str(b)+"_"+str(c)+"_"+str(d)+"_"+str(e)

def pglyco3_processing(df,
                    only_duplicated="Drop_duplicated"):
    """Create required columns.
    Args:
    duplicated: True or False: whether or not only peak duplicated columns.
    True: only duplicated row are retained for repeatability test.
    False: only rows with lowest totalFDR for duplicated columns or unique columns are retained for training.
    
    """
    # ipdb.set_trace()
    df["iden_pep"]=df.apply(combine_iden_pep,axis=1) #eg. JASQNQDNVYQGGGVCLDCQHHTTGINCER_16.Car.19.Car.28.Car._0_4_(N(N(H(H(H))(H(H)))))
    if only_duplicated=="Duplicated":
        df1=df[["iden_pep"]].loc[df["iden_pep"].duplicated()].drop_duplicates()
        df=df.loc[df["iden_pep"].isin(df1["iden_pep"])]
    print("Waiting to process multiply glycopeptides")
    if only_duplicated == "Drop_duplicated":
        df.sort_values(by='TotalFDR',ascending=True,inplace=True)
        # ipdb.set_trace()
        df.drop_duplicates(subset=['iden_pep'],inplace=True)
        df.reset_index(drop=True,inplace=True)
    if only_duplicated == "Retained_all":
        pass
    return df
# --------------------------- spectrum filtration---------------------#
#从json中找到相应的谱图,缩小搜索空间
def json_extraction(jsonfold=jsonfold,
                    jsonname=jsonname,
                    filename=filter_jsonname,
                    mgfsourceorign=mgfsourceorign):
    datalis=pd.read_json(os.path.join(jsonfold, jsonname))
    datalis["title"]=datalis["SourceFile"].map(str) + "-" + datalis["Spectrum"].map(str)
    datalis=datalis.loc[datalis["title"].isin(df["GlySpec"])]
    print("Please ensure the Spectrum numbers of MsConvert json files match those of the pGlyco3 result!")
    datalis.reset_index(inplace=True, drop=True)
    datalis.to_json(os.path.join(jsonfold, filename))
    return datalis
# ----------------------- ions picking ------------------------------#
def fragment_training(instance):
    spectrum=instance["GlySpec"]
    datalis_1=datalis.loc[datalis["title"]==spectrum]
    datalis_1=datalis_1.reset_index(drop=True)
    iden_pep=instance["iden_pep"]
    mz_calc=masses.pepfragmass(iden_pep,FRAG_MODE,3) #iden_pep已经改成了glysite,避免多J的可能
    ppm=TOLER
    FragmentMz=[]
    for mz in mz_calc:
        for ion in mz:
            FragmentMz.append(list(ion.values())[0])
    FragmentMz=list(set(FragmentMz))
    mass={"FragmentMz":FragmentMz}
    #FragmentMz:所有算出来的理论质荷比
    mzdict=mgf_processing.putTintensity(ppm, mass, datalis_1)
    for k in list(mzdict.keys()):
        if mzdict[k]==0:
            del mzdict[k]
    mzdict_1={}
    #补上mzdict的碎裂类型
    for i in mz_calc:
        for a in i:
            mz_calc_1=list(a.values())[0]
            if mz_calc_1 in list(mzdict.keys()):
                # print("a",a)
                # print("mzdict[mz_calc_1]",mzdict[mz_calc_1])
                type=list(a.keys())[0]
                intensity=mzdict[mz_calc_1]
                if not mz_calc_1 in mzdict_1.keys():
                    type_list=[]
                    type_list.append(type)
                    ions=(type_list,intensity)
                    mzdict_1[mz_calc_1]=ions
                else:
                    type_list=mzdict_1[mz_calc_1][0]
                    type_list.append(type)
                    ions=(type_list,intensity)
                    mzdict_1[mz_calc_1]=ions
    return mzdict_1

def mz_matching(instance):
    spectrum=instance["GlySpec"]
    datalis_1=datalis.loc[datalis["title"]==spectrum]
    datalis_1=datalis_1.reset_index(drop=True)
    iden_pep=instance["iden_pep"]
    mz_calc=masses.pepfragmass(iden_pep,["HCD_BY_2"],4) #iden_pep已经改成了glysite,避免多J的可能
    ppm=TOLER
    FragmentMz=[]
    for mz in mz_calc:
        for ion in mz:
            FragmentMz.append(list(ion.values())[0])
    FragmentMz=list(set(FragmentMz))
    FragmentMz.sort()
    mzexp=datalis_1["mz"][0]
    mzexp.sort()
    matchmz=[]
    for k in mzexp:
        i = (np.abs(np.array(FragmentMz) - k)).argmin()
        # ipdb.set_trace()
        if abs(FragmentMz[i] - k) < k * TOLER * 1 / 1000000:  #args.ppm=tolerance here,可以改回args版本
            matchmz.append(k)
    return {"matchmz":len(matchmz),"calc":len(FragmentMz),"mzexp":len(mzexp)}
# --------------------------- execution ---------------------#
if __name__=="__main__":
    DFNAME_path = Path(DFNAME)
    print(DFNAME_path)
    assert DFNAME_path.exists()
    # pglyco3 formatted result
    df_fp,use_weights=pglyco3_result(DFNAME)
    df=pglyco3_processing(df_fp,
                        only_duplicated=only_duplicated)
    # if mgfsourceorign=="MsConvert":
    df["GlySpec"]=df["RawName"].map(str) + "-" + df["Scan"].map(str)
    #json file
    json_path=Path(jsonfold,jsonname)
    if json_path.exists():
        print(f"{jsonname} exists.")
    else:
        print(f"{jsonname} does not exist. Begin mgf_process to produce required file...")
        datalis=mgf_processing.mgf_process(mgfdatafold=mgfdatafold,sourceorign=mgfsourceorign)
    #filtered json file
    file3_name_path = Path(jsonfold,filter_jsonname)
    # if file3_name_path.exists():
    #     print(f"{filter_jsonname} exists.")
    #     datalis=pd.read_json(os.path.join(jsonfold, filter_jsonname))
    # else:
    print(f"{file3_name_path} does not exist. Begin json_extraction to produce required file...")
    datalis=json_extraction(jsonfold=jsonfold,
                jsonname=jsonname,
                filename=filter_jsonname,
                mgfsourceorign=mgfsourceorign)
    datalis.drop_duplicates(subset="title",inplace=True)
    df=df[df["GlySpec"].isin(datalis["title"])]
    assert len(df["GlySpec"].drop_duplicates())==len(datalis["title"].drop_duplicates())
    if use_weights:
        df["weights"] = df.apply(lambda x: weights(x, Enzyme), axis=1)
        df=df[[ "GlySpec",'Charge',"RT","PrecursorMZ",'Peptide', 'Mod', 'PlausibleStruct', 'GlySite', 'iden_pep',"TotalFDR","weights"]]
    else:
        df=df[[ "GlySpec",'Charge',"RT","PrecursorMZ",'Peptide', 'Mod', 'PlausibleStruct', 'GlySite', 'iden_pep',"TotalFDR"]]
    df.drop_duplicates(subset=["GlySpec",'Charge',"RT","PrecursorMZ",'Peptide', 'Mod', 'PlausibleStruct', 'GlySite', 'iden_pep',"TotalFDR"],inplace=True)
    df.reset_index(drop=True,inplace=True)
    print("len(df_iden_pep.drop_duplicates())",len(df["iden_pep"].drop_duplicates()))
    print("len(df)",len(df))
    df["ions"]=df.apply(fragment_training,axis=1)
    
    df.to_csv(output_name,index=False)