<jupyter_start><jupyter_text>Introduction NetworKIN Result Files- See Code/PreprocessingPredictionData/Running_NetworKIN.md for instructions on installing and running NetworKIN for preparing the initial Raw NetworKIN results file About the NetworKIN Raw Results File- Version: Last Run in KinPred Nov 2020- NetworKIN result files contain the following columns:    - Name    - Position    - Tree    - NetPhorest Group    - Kinase/Phosphatase/Phospho-binding domain    - NetworKIN score    - NetPhorest probability    - STRING score    - Target STRING ID    - Kinase/Phosphatase/Phospho-binding domain    - STRING ID    - Target description    - Kinase/Phosphatase/Phospho-binding domain description    - Target Name    - Kinase/Phosphatase/Phospho-binding domain Name    - Peptide sequence window    - Intermediate nodes         Preprocessing NetworKIN Raw1. Filtering Date: only keep the results of kinase prediction (TREE == 'KIN')2. Mapping accessions:    - **Mapping substrate accessions:**  Get the UniprotID from the sequence identifiers of the original fasta files that submitted in NetworKIN for prediction.     - **Mapping kinase accessions:** get uniprotID for the kinases using 'Kinase/Phosphatase/Phospho-binding domain description' column         - correcting duplicated names found above3. **Mapping sites:**    - formatting pep seq (e.g. ----MsGSKSV --> ____MSGSKSV)    - update 'site' if the orignal sequence submitted for perdiction is different from the reference sequence due to changes causing shift in position**Output Files Dataframe:**- **substrate_id:** unique IDs for the substrate phosphorylation site (substrate_acc + position)- **substrate_name:** substrate gene name used in NetworKIN- **substrate_acc:** substrate uniprotID- **site:** aa + position in protein sequence- **Position:** position in protein sequence- **pep:** +/- 5 AA- **kinase_name:** kinase name used in NetworKIN- **kinase_acc:** kinase uniprotID- **score:** NetworKIN score     Creating Resource Files1.  **'globalKinaseMap.csv':**     - creat a new or add unique kinases from NetworKIN to the globel substrate resource file.    - get and add the Kinase Name that would use across all perdictors for the kinases to the result files Standard Formatted NetworKIN**'NetworKIN_formatted.csv':** Standardize the preprocessed file with following columns:- **substrate_id** - unique IDs for the substrate phosphorylation site (substrate_acc + position)- **substrate_name** - gene name for the substrates- **substrate_acc** - mapped UniprotIDs for the substrates- **site** - phosphorylation  site- **pep** - +/- 5 AA peptide sequence around the site- **score**- **Kinase Name**<jupyter_code># IMPORTS
import pandas as pd
import os
import re
import glob

import humanProteomesReference, networKin_convert, getUniprotID, checkSite

#only need when testing the code
import time
# DEFINE FILE NAMES/DIRs
##################
# Version (Date) #
##################
version = '2019-12-11'

##################
# File Location  #
##################
# local (../../)
base = '../../'

##################################################
# For Prepare Fasta Files to Submit in NetworKIN #
##################################################

# Human Proteome fasta file
HP_fasta = base + 'Data/Raw/HumanProteome/humanProteome_' + version + '.fasta'
# Dir for splited Human Proteome fasta files
HP_dir = base + 'Data/Raw/HumanProteome/'

# human proteome referece file 
HP_csv = base + 'Data/Map/humanProteome_' + version + '.csv'

####################################################
# For Preprocessing NetworKIN Prediction Results   #
#--------------------------------------------------#
# . The files submitted for NetworKIN predictor is #
#   NOT the up-to-date human proteom sequences     #
# . There has been an update in human proteomes    #
#   from the time the perdiction results were got  #
#   to running the preporcessing steps.            #
####################################################

# NetworKIN results dir
NW_dir = base + 'Data/Raw/NetworKIN/'
NW_update_dir = base + 'Data/Raw/NetworKIN/updated/'
# NetworKIN temp dir
NW_temp_dir_acc = base + 'Data/Temp/NetworKIN/mappedAcc/'
NW_temp_dir_site = base + 'Data/Temp/NetworKIN/mappedSite/'
NW_temp_dir_acc_update = base + 'Data/Temp/NetworKIN/mappedAcc/updated/'
NW_temp_dir_site_update = base + 'Data/Temp/NetworKIN/mappedSite/updated/'

# Resource Files
HK_org = base + 'Data/Raw/HumanKinase/globalKinaseMap.txt'                  # orginal manually created kinase file
KinaseMap = base + 'Data/Map/globalKinaseMap.csv'                           # add all unique kinase in HPRD to the global file

# Standard formatted output file
NW_formatted = base + 'Data/Formatted/NetworKIN/NetworKIN_formatted_' + version + '.csv'       # preprocessed file with cloumns: substrate_id/substrate/substrate_acc/kinase/site/pep/score<jupyter_output><empty_output><jupyter_text>Preprocessing NetworKIN Raw In the NetworKIN raw file, most of the names in the 'Kinase/Phosphatase/Phospho-binding domain' column can not use to retrieve the uniprotID for the kinase, but all the names in the 'Kinase/Phosphatase/Phospho-binding domain description' column can. However, there are different 'Kinase/Phosphatase/Phospho-binding domain' with the same 'Kinase/Phosphatase/Phospho-binding domain description'.  We need to identify those and correct them if need.   Get unique kinases in NetworKIN<jupyter_code>all_results = glob.glob(NW_dir + '*.tsv')
#create empty df to store unique kinases
df_unique_kin = pd.DataFrame()
for filename in all_results:
    df = pd.read_csv(filename, usecols = ['Tree','Kinase/Phosphatase/Phospho-binding domain', 'Kinase/Phosphatase/Phospho-binding domain description'], sep = '\t')
    # the only type of perdiction we are intreseted in is 'KIN'
    df = df[df.Tree == 'KIN']
    df = df[['Kinase/Phosphatase/Phospho-binding domain', 'Kinase/Phosphatase/Phospho-binding domain description']].drop_duplicates()
    # append unique kinases found in each result files
    df_unique_kin = df_unique_kin.append(df)
    
# drop any duplicated kinases
df_unique_kin = df_unique_kin.drop_duplicates()
# get the  kinase(s) with the same domain description
duplicateRowsDF = df_unique_kin[df_unique_kin.duplicated(['Kinase/Phosphatase/Phospho-binding domain description'],keep=False)]
duplicateRowsDF.sort_values(['Kinase/Phosphatase/Phospho-binding domain description','Kinase/Phosphatase/Phospho-binding domain'])<jupyter_output><empty_output><jupyter_text>After manually check for the above kinases, MST4 should associate with STK26 not STK3. Create a dictionary to correct that in the result files.<jupyter_code>correct_kinase = {'MST4' : 'STK26'}<jupyter_output><empty_output><jupyter_text>Mapping Accessions (UniprotID) and Site1. Filtering Date: only keep the results of kinase prediction (TREE == 'KIN')2. Mapping accessions:    - **Mapping substrate accessions:**  Get the UniprotID from the sequence identifiers of the original fasta files that submitted in NetworKIN for prediction.     - **Mapping kinase accessions:** get uniprotID for the kinases using 'Kinase/Phosphatase/Phospho-binding domain description' column         - correcting duplicated names found above3. **Mapping sites:**    - formatting pep seq (e.g. ----MsGSKSV --> ____MSGSKSV)    - update 'site' if the orignal sequence submitted for perdiction is different from the reference sequence due to changes causing shift in position**Output Files Dataframe:**- **substrate_id:** unique IDs for the substrate phosphorylation site (substrate_acc + position)- **substrate_name:** substrate gene name used in NetworKIN- **substrate_acc:** substrate uniprotID- **site:** aa + position in protein sequence- **Position:** position in protein sequence- **pep:** +/- 5 AA- **kinase_name:** kinase name used in NetworKIN- **kinase_acc:** kinase uniprotID- **score:** NetworKIN score **Mapping Accessions**<jupyter_code># convert substrate_acc and kinase_acc
convert_type = 'acc'
networKin_convert.kin_convert_directory(NW_dir, 'na', NW_temp_dir_acc, convert_type)<jupyter_output>reading  19.tsv
getting unique sub
getting sub_acc
merge
done 2.626702308654785
getting unique kin
getting kin_acc
merge
done 240.26136994361877
saving
Done 30.140185832977295
reading  18.tsv
getting unique sub
getting sub_acc
merge
done 3.2055108547210693
getting unique kin
getting kin_acc
merge
done 251.04621005058289
saving
Done 38.56568384170532
reading  20.tsv
getting unique sub
getting sub_acc
merge
done 3.681636095046997
getting unique kin
getting kin_acc
merge
done 326.4341878890991
saving
Done 40.38838195800781
reading  21.tsv
getting unique sub
getting sub_acc
merge
done 1.4215641021728516
getting unique kin
getting kin_acc
merge
done 316.3874821662903
saving
Done 17.19068193435669
reading  8.tsv
getting unique sub
getting sub_acc
merge
done 2.9008588790893555
getting unique kin
getting kin_acc
merge
done 276.8457570075989
saving
Done 27.926394939422607
reading  9.tsv
getting unique sub
getting sub_acc
merge
done 2.908118724822998
getting unique kin
getting kin_acc
merge
done[...]<jupyter_text>**Mapping Site**<jupyter_code># map the site to the updated (new) human proteome reference
convert_type = 'site'
networKin_convert.kin_convert_directory(NW_temp_dir_acc, HP_csv, NW_temp_dir_site, convert_type)<jupyter_output>Set input file dir...
done
read the Human Proteome df...
done
processing  9 .tsv
Get unique substrate sites in  9 .tsv
P30988 :   corrected:  (seqNotFound)S4  site:  S4  seq len:  474 pep:  MQFSGEKIS
P30988 :   corrected:  (seqNotFound)S9  site:  S9  seq len:  474 pep:  SGEKISGQRDL
P30988 :   corrected:  (seqNotFound)S17  site:  S17  seq len:  474 pep:  RDLQKSKMRFT
P30988 :   corrected:  (seqNotFound)T22  site:  T22  seq len:  474 pep:  SKMRFTFTSRC
P30988 :   corrected:  (seqNotFound)T195  site:  T195  seq len:  474 pep:  FFRKLTTIFPL
P30988 :   corrected:  (seqNotFound)T196  site:  T196  seq len:  474 pep:  FRKLTTIFPLN
P30988 :   corrected:  (seqNotFound)Y204  site:  Y204  seq len:  474 pep:  PLNWKYRKALS
P30988 :   corrected:  (seqNotFound)S209  site:  S209  seq len:  474 pep:  YRKALSLGCQR
Q9BSG1 :   corrected:  (seqNotFound)S161  site:  S161  seq len:  425 pep:  LRRRRSALSRE
Saving  9 .tsv
9 .tsv Done  256.96863174438477
processing  8 .tsv
Get unique substrate sites in  8 .tsv
Q86VQ6 :[...]<jupyter_text>**Remove unmapped/outdated results**- Check if there is any unmapped substrate/site due to outdated uniprot sequence records    - get a list of outdated uniprotID<jupyter_code>all_results = glob.glob(NW_temp_dir_site + '*.csv')
    
updatedSub_li = []
for filename in all_results:
    start = time.time() 
    df_unmapped = pd.read_csv(filename, usecols = ['substrate_id','substrate_acc', 'site'])
    df_unmapped = df_unmapped[~(df_unmapped['site'].str.contains('S|T|Y', na=False)) | (df_unmapped['substrate_id'] == 'outdated')]
    df_unmapped = df_unmapped.substrate_acc.drop_duplicates()
    unmapped_li = df_unmapped.values.tolist()
    updatedSub_li.extend(unmapped_li)
    end = time.time()
    print (f"Time\t{(end-start):.3f}")
updatedSub_li<jupyter_output>Time	8.223
Time	9.399
Time	9.811
Time	10.685
Time	9.016
Time	8.493
Time	9.990
Time	5.003
Time	12.716
Time	9.222
Time	8.385
Time	9.441
Time	8.693
Time	8.373
Time	8.131
Time	8.660
Time	11.187
Time	10.031
Time	8.585
Time	8.566
Time	9.941<jupyter_text>- Remove any record with unmapped substrate_acc/site in '*_mappedSite.csv'<jupyter_code>all_results = glob.glob(NW_temp_dir_site + '*.csv')

for filename in all_results:
    start = time.time()
    df_mapSite = pd.read_csv(filename)
    # remove the outdated records from df_subMap
    df_update = df_mapSite[~df_mapSite['substrate_acc'].isin(updatedSub_li)]
    df_update.to_csv(filename, chunksize=100000, index=False)
    end = time.time()
    print (f"chunk time\t{(end-start):.3f}")<jupyter_output>chunk time	32.552
chunk time	34.708
chunk time	37.731
chunk time	42.370
chunk time	35.697
chunk time	32.222
chunk time	37.086
chunk time	19.179
chunk time	45.215
chunk time	36.513
chunk time	33.216
chunk time	37.718
chunk time	33.076
chunk time	33.881
chunk time	32.461
chunk time	33.712
chunk time	39.761
chunk time	37.249
chunk time	32.760
chunk time	32.690
chunk time	37.350<jupyter_text>Update NetworKIN results- The next 6 cells is only for NetworKIN results from outdated Human Proteomes sequences (input sequences for NetworKIN prediction an earlier version than the referece human proteome sequence) 1. download the sequence fasta file of the above UniprotID (updatedSub_li) from Uniprot.org, save it as '../Data/Raw/HumanProteome/NetworKIN_updateSub.fasta'. Submit the 'NetworKIN_updateSub.fasta' in NetworKIN again2. run mapAcc and mapSite function for the result file from 'NetworKIN_updateSub.fasta'<jupyter_code># convert substrate_acc and kinase_acc
convert_type = 'acc'
networKin_convert.kin_convert_directory(NW_update_dir, 'na', NW_temp_dir_acc_update, convert_type)
# map the site to the updated (new) human proteome reference
convert_type = 'site'
networKin_convert.kin_convert_directory(NW_temp_dir_acc_update, HP_csv, NW_temp_dir_site_update, convert_type)<jupyter_output>Set input file dir...
done
read the Human Proteome df...
done
processing  2020 .tsv
Get unique substrate sites in  2020 .tsv
P55073 :   corrected:  (seqNotFound)S167  site:  S167  seq len:  304 pep:  VLNFGSCTPPF
P55073 :   corrected:  (seqNotFound)T169  site:  T169  seq len:  304 pep:  NFGSCTPPFMA
P59797 :   corrected:  (seqNotFound)Y269  site:  Y269  seq len:  346 pep:  LIRVTYCGLSY
P59797 :   corrected:  (seqNotFound)S275  site:  S275  seq len:  346 pep:  CGLSYSLRYIL
P59797 :   corrected:  (seqNotFound)S273  site:  S273  seq len:  346 pep:  TYCGLSYSLRY
P59797 :   corrected:  (seqNotFound)T268  site:  T268  seq len:  346 pep:  VLIRVTYCGLS
P59797 :   corrected:  (seqNotFound)Y274  site:  Y274  seq len:  346 pep:  YCGLSYSLRYI
Q9NZV6 :   corrected:  (seqNotFound)S92  site:  S92  seq len:  116 pep:  PKPGQSRFIFS
Q9NZV6 :   corrected:  (seqNotFound)S97  site:  S97  seq len:  116 pep:  SRFIFSSSLKF
Q9NZV6 :   corrected:  (seqNotFound)S98  site:  S98  seq len:  116 pep:  RFIFSSSLKFV
Q9NZV6 :   [...]<jupyter_text>3. save '2020_mappedSite.csv' under the same dir as other *_mappedSite.csv files    - NetworKIN ignorgs the 'U'(Selenocysteine) in some substrate sequences, this causes frame shift of the downstream sequences    - those unmapped sites will not be included<jupyter_code>start = time.time()
df_mapUpdateSite = pd.read_csv(NW_temp_dir_site_update + '2020_mappedSite.csv')
df_mapUpdateSite = df_mapUpdateSite[(df_mapUpdateSite['site'].str.contains('S|T|Y', na=False)) & (df_mapUpdateSite['substrate_id'] != 'outdated')]
df_mapUpdateSite.to_csv(NW_temp_dir_site + '2020_mappedSite.csv')
end = time.time()
print (f"Time\t{(end-start):.3f}")
df_mapUpdateSite<jupyter_output>Time	3.122<jupyter_text>Get the Gene Name of the Substrates from the Reference Human Proteome- get and add the Gene Name that would use across all perdictors for the substrates to the result files<jupyter_code># get the protein gene names and accessions from the Reference Human Proteome
df_unique_sub =  pd.read_csv(HP_csv, usecols = ['UniprotID','Gene Name'], sep = '\t')
df_unique_sub
# add the Gene Name that would use across all perdictors for the substrates to the result files
start = time.time()
all_results = glob.glob(NW_temp_dir_site + '*.csv')

for filename in all_results:
    df = pd.read_csv(filename, usecols = ['substrate_id', 'substrate_acc', 'substrate_name', 'site', 'Position', 'pep', 'score', 'kinase_name', 'kinase_acc'])
    # merge df_subsMap with df_unique_sub to add the common substrate gene name to the df
    df = df.merge(df_unique_sub, left_on=['substrate_acc'], right_on=['UniprotID'], how = 'left')

    df = df.drop(columns = ['UniprotID'])
    df.to_csv(filename, index=False)
df<jupyter_output><empty_output><jupyter_text>Creating Resource Files**globalKinaseMap**- creat a new or add unique kinases from NetworKIN to the globel kinase resource file.- get and add the Kinase Name that would use across all perdictors for the kinases to the result files<jupyter_code># get unique kinases in the NetworKIN result files

all_results = glob.glob(NW_temp_dir_site + '*.csv')

df_unique_kin = pd.DataFrame()

for filename in all_results:
    df = pd.read_csv(filename, usecols = ['kinase_acc', 'kinase_name'])
    df = df.drop_duplicates()
    df_unique_kin = df_unique_kin.append(df, ignore_index=True)
    
df_unique_kin = df_unique_kin.drop_duplicates()
df_unique_kin
# start = time.time()
unmapped_list = pd.DataFrame()

if os.path.isfile(KinaseMap): 
    df_humanKinase = pd.read_csv(KinaseMap)
# if globalkinaseMap.csv file does not exist, create an new df using orginal human kinase map
else:
    df_humanKinase = pd.read_csv(HK_org, usecols = ['Kinase Name', 'Preferred Name', 'UniprotID', 'Type', 'description'], sep = '\t')
    df_humanKinase['description'].replace(regex=True,inplace=True,to_replace=r'\[Source.+\]',value=r'')

# add unique kinases from NetworKIN to the globel kinases resource file
for index, row in df_unique_kin.iterrows():
    kinase = df_unique_kin.at[index, 'kinase_acc']
    # if the kinase/other enzyme already in the globalKinaseMap.csv file
    if any(df_humanKinase.UniprotID == kinase):
        # get the index of the substrate in the globalKinaseMap.csv file 
        idx = df_humanKinase.index[df_humanKinase.UniprotID == kinase].values[0] 

        df_humanKinase.at[idx, 'NetworKIN_kinase_name'] = df_unique_kin.at[index, 'kinase_name']
        
    # if the kinase is not in the globalKinaseMap.csv file, we need a list to check annotations manullay
    else:
        unmapped_list = unmapped_list.append(row,sort=False).reset_index(drop=True)
        
print (unmapped_list)<jupyter_output>kinase_acc kinase_name
0     Q16654        PDK4
1     Q15120        PDK3
2     Q15119        PDK2
3     Q86VQ0        LCA5<jupyter_text>- Manually check the above unmapped kinase. Create a dictionay for the one(s) that are kinases to add in the globalKinaseMap.csv.  Create a list of the one(s) that are not kinases to drop records from the NetworKIN result files.<jupyter_code>new_kinase = {'EPHB6':'Ephrin type-B receptor 6'}

not_kinase = ['PDK2','PDK3','PDK4','LCA5']<jupyter_output><empty_output><jupyter_text>- add the 1 protein kinases to the globalKinaseMap.csv using above dictionary<jupyter_code># get the length (last index) of the current df_humanKinase
len = df_humanKinase.UniprotID.count()

# add the 1 new kinase in the globalKinaseMap.csv
for key in new_kinase:
    # add the kinase name and description for the new kinase 
    df_humanKinase = df_humanKinase.append({'Kinase Name': key}, ignore_index=True)
    df_humanKinase.at[len,'Preferred Name'] = key
    df_humanKinase.at[len,'description'] = new_kinase[key]
    # get the index of where the new kinase is in the df_unique_kinase
    i = df_unique_kin.index[df_unique_kin['kinase_name'] == key].values[0] 
    # add the uniprotID for the new kinase 
    df_humanKinase.at[len,'UniprotID'] = df_unique_kin.at[i, 'kinase_acc']
    # add the accs used in NetworKIN for the new kinase
    df_humanKinase.at[len,'NetworKIN_kinase_name'] = df_unique_kin.at[i, 'kinase_name']

    len += 1

df_humanKinase.to_csv(KinaseMap,index = False) 
df_humanKinase
# get the new list kinase with common kinase name that would use across all referece and the uniprotID for these kinase
df_unique_kin = df_humanKinase[['Kinase Name','UniprotID']]
df_unique_kin<jupyter_output><empty_output><jupyter_text>- remove the 4 enzyme records that are not protein kinase from the NetworKIN result files- add the Kinase Name that would use across all perdictors for the kinase to the result files<jupyter_code>all_results = glob.glob(NW_temp_dir_site + '*.csv')

for filename in all_results:
    df = pd.read_csv(filename)
    # remove the 4 enzyme records that are not protein kinase from the df_kinaseMap
    df = df[~df['kinase_name'].isin(not_kinase)]

    # merge with df_unique_kinase to add the common kinase name to the df
    df = df.merge(df_unique_kin, left_on='kinase_acc', right_on='UniprotID', how = 'left')
    # drop the duplicated uniprot id for kinases
    df = df.drop(columns = 'UniprotID')

    df.to_csv(filename,index=False)  

df<jupyter_output><empty_output><jupyter_text>Standard Formatted NetWorKIN 'NetworKIN_formatted.csv'Standardize the preprocessed file with following columns:- substrate_id - unique IDs for the substrate phosphorylation site (substrate_acc + position)- substrate_name - gene name for the substrates- substrate_acc - mapped UniprotIDs for the substrates- site - phosphorylation site- pep - +/- 5 AA peptide sequence around the site- score - perdiction score- Kinase Name - Kinase name<jupyter_code>all_results = glob.glob(NW_temp_dir_site + '*.csv')
NW = []
for filename in all_results:
    df = pd.read_csv(filename, usecols = ['substrate_id','substrate_acc','Gene Name', 'site','pep', 'score', 'Kinase Name'])
    df = df.rename(columns={'Gene Name': 'substrate_name'})
    NW.append(df)
        
df_final = pd.concat(NW)
df_final = df_final.drop_duplicates()
df_final.to_csv(NW_formatted, chunksize = 1000000, index = False)
df_final<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10504300418292052, 'fraction_numerical': 0.03962024721530291, 'mean_word_length': 4.710681696188942, 'pattern_counts': {'":': 0, '<': 54, '<?xml version=': 0, '>': 56, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '63228', 'license': 'MIT', 'n_tokens_mistral': 7782, 'n_tokens_neox': 6944, 'n_words': 2123, 'path': 'Code/PreprocessingPredictionData/FormattingNetworKIN.ipynb', 'repo_name': 'NaegleLab/KinPred', 'size': '21277'}