<jupyter_start><jupyter_text>Setting up and reading data files<jupyter_code>import numpy as np
import pandas as pd
from uniqify import uniqify
from unlistify import unlistify
import pickle
#-------------------------------------------------------------------------
# Loading all relevant data files. All relevant files are saved in the
# data_metabolism folder.
#
# NOTE: please unzip the file: data_metabolism/kegg_files.zip into the same 
# directory and make sure all the following files can be read by the program
# in order for the analysis to work.
#-------------------------------------------------------------------------
stoich_matrix = np.load('data_metabolism/akshit_stoich_matrix.npy')
names = [e[:-1] for e in open('data_metabolism/akshit_names.txt', 'r').readlines()]
rxns = ''.join(open('data_metabolism/akshit_rxns.txt', 'r').readlines()).split()
cpds = ''.join(open('data_metabolism/akshit_mets.txt', 'r').readlines()).split()

met_map = kegg_to_id = {cpds[e] : e for e in range(len(cpds))}
inv_met_map = id_to_kegg = {value: key for key, value in met_map.items()}
rxn_map = rxn_kegg_to_id = {rxns[e] : e for e in range(len(rxns))}
inv_rxn_map = rxn_id_to_kegg = {value: key for key, value in rxn_map.items()}
cpd_string_dict = {cpds[e] : names[e] for e in range(len(names))}

#-------------------------------------------------------------------------
# Parsing some KEGG reactions correctly.
#-------------------------------------------------------------------------
stoich_matrix[rxn_kegg_to_id['R11307']][kegg_to_id['C00760']] = -1
stoich_matrix[rxn_kegg_to_id['R00803']][kegg_to_id['C00103']] = 0
stoich_matrix[rxn_kegg_to_id['R00028']][kegg_to_id['C00031']] = 0
stoich_matrix[rxn_kegg_to_id['R00028']][kegg_to_id['C00267']] = 1
stoich_matrix[rxn_kegg_to_id['R00010']][kegg_to_id['C00031']] = 0
stoich_matrix[rxn_kegg_to_id['R00010']][kegg_to_id['C00267']] = 1
stoich_matrix[rxn_kegg_to_id['R00839']][kegg_to_id['C00031']] = 0

#-------------------------------------------------------------------------
# Initializing and setting the currency molecules.
#-------------------------------------------------------------------------
Currency = ''.join(open('data_metabolism/ncurr.txt', 'r').readlines()).split()

#-------------------------------------------------------------------------
# Defining reactant, product and reactant number vectors for the scope 
# expansion algorithm.
#-------------------------------------------------------------------------
rxnMat = ((stoich_matrix.clip(max = 0.0)) != 0) * 1
prodMat = (stoich_matrix.clip(min = 0.0) != 0) * 1
sumRxnVec = np.sum(rxnMat, axis = 1)
sumProdVec = np.sum(prodMat, axis = 1)<jupyter_output><empty_output><jupyter_text>Defining all 16 carbon sources used in the experiment.<jupyter_code>martina_cpd_to_kegg = dict( pd.read_csv( 'data_metabolism/compound_to_kegg.csv' ).values.tolist( ) )
martina_cpd_to_kegg[ 'glucose' ] = 'C00221'
martina_cpd_to_kegg[ 'cellobiose' ] = 'C04534'
martina_cpd_to_kegg<jupyter_output><empty_output><jupyter_text>Reading media<jupyter_code>martina_media_df = pd.read_excel( 'data_metabolism/media_from_martina_modified_names.xlsx', index_col=0 ).fillna( '' )
media_listoflist = list( martina_media_df.values.tolist( ) )
for idx, tm in enumerate( media_listoflist ):
    media_listoflist[ idx ] = list( filter( None, tm ) )
kegg_media = [ [ martina_cpd_to_kegg[ e ] for e in tm ] for tm in media_listoflist ]
CS = [id_to_kegg[kegg_to_id[e]] for e in uniqify( unlistify( kegg_media ) )]<jupyter_output><empty_output><jupyter_text>Checking that we have 75 defined media in total.<jupyter_code>len(kegg_media)<jupyter_output><empty_output><jupyter_text>The scope expansion algorithm<jupyter_code>def giveScope(rxnMat, prodMat, seedVec, sumVec):
    """
    Takes in the stoichiometric matrix in the form of a reaction 
    and product matrix, along with a sum vector mentioning the 
    number of reactants in each reaction in the matrix. Also takes
    in a vector of initially seeded metabolites (in seedVec).

    Expands, beginning from the seeded metabolites, a set of reactions
    and metabolites that can be reached iteratively from the former.
    

    RETURNS:
    Returns a set of 'scope-expanded' metabolites and reactions in the 
    stoichiometric matrix form.

    currScopeMets is the set of metabolites, with their usual IDs.
    rxnProc is the set of reactions with my personal IDs 
    """
    currScopeMets = np.copy(seedVec)
    scopeSize = sum(currScopeMets)
    prevScopeSize = 0.0

    while scopeSize > prevScopeSize:
        prevScopeSize = scopeSize

        # This generates the rho vector.
        rxnProc = ((np.dot(rxnMat, currScopeMets) - sumVec) == 0) * 1
        # This generates the phi vector.
        currScopeMets = (np.dot(np.transpose(prodMat), rxnProc) + currScopeMets > 0) * 1

        # Recalculating scope size.
        scopeSize = sum(currScopeMets)

    return currScopeMets, rxnProc<jupyter_output><empty_output><jupyter_text>Preparing all reaction maps from KEGG to use  Getting all amino acid reaction maps<jupyter_code>aa03 = ''.join(open('data_metabolism/metabolic_maps_kegg/aa_maps/map00250.txt', 'r').readlines()).split()
aa04 = ''.join(open('data_metabolism/metabolic_maps_kegg/aa_maps/map00260.txt', 'r').readlines()).split()
aa05 = ''.join(open('data_metabolism/metabolic_maps_kegg/aa_maps/map00270.txt', 'r').readlines()).split()
aa06 = ''.join(open('data_metabolism/metabolic_maps_kegg/aa_maps/map00280.txt', 'r').readlines()).split()
aa07 = ''.join(open('data_metabolism/metabolic_maps_kegg/aa_maps/map00290.txt', 'r').readlines()).split()
aa08 = ''.join(open('data_metabolism/metabolic_maps_kegg/aa_maps/map00300.txt', 'r').readlines()).split()
aa09 = ''.join(open('data_metabolism/metabolic_maps_kegg/aa_maps/map00310.txt', 'r').readlines()).split()
aa10 = ''.join(open('data_metabolism/metabolic_maps_kegg/aa_maps/map00330.txt', 'r').readlines()).split()
aa11 = ''.join(open('data_metabolism/metabolic_maps_kegg/aa_maps/map00340.txt', 'r').readlines()).split()
aa12 = ''.join(open('data_metabolism/metabolic_maps_kegg/aa_maps/map00350.txt', 'r').readlines()).split()
aa13 = ''.join(open('data_metabolism/metabolic_maps_kegg/aa_maps/map00360.txt', 'r').readlines()).split()
aa14 = ''.join(open('data_metabolism/metabolic_maps_kegg/aa_maps/map00380.txt', 'r').readlines()).split()
canaa = aa03 + aa04 + aa05 + aa06 + aa07 + aa08 + aa09 + aa10 + aa11 + aa12 + aa13 + aa14<jupyter_output><empty_output><jupyter_text>Getting all carbohydrate and carbon metabolism maps<jupyter_code>can1 = ''.join(open('data_metabolism/metabolic_maps_kegg/rn00330.txt', 'r').readlines()).split()
can2 = ''.join(open('data_metabolism/metabolic_maps_kegg/map01200.txt', 'r').readlines()).split()
can3 = ''.join(open('data_metabolism/metabolic_maps_kegg/map00051.txt', 'r').readlines()).split()
can4 = ''.join(open('data_metabolism/metabolic_maps_kegg/map00500.txt', 'r').readlines()).split()
can5 = ''.join(open('data_metabolism/metabolic_maps_kegg/map00040.txt', 'r').readlines()).split()
can6 = ''.join(open('data_metabolism/metabolic_maps_kegg/map00561.txt', 'r').readlines()).split()
can = can1 + can2 + can3 + can4 + can5 + can6 + canaa
can = list((set(can)) & set(rxns))<jupyter_output><empty_output><jupyter_text>Removing a manually curated list of reactions which are not specific to soil bacteria.<jupyter_code>fwrxns = [rxn_kegg_to_id[e] for e in can if e not in [ 'R03292', 'R03293', 'R11819', 'R04447', 'R04446', 'R03295', 'R03290', 'R03291', 'R02727', 'R00801', 'R05140', 'R11790', 'R04194', 'R00806', 'R02120', 'R01823', 'R00178', 'R11791', 'R11309', 'R00948', 'R00872', 'R09945', 'R04111', 'R01555',  'R01557', 'R01295', 'R11540', 'R00341', 'R00345', 'R00354', 'R00479', 'R00878', 'R09994', 'R02110']]
mp = int(len(stoich_matrix)/2)
rvrxns = [ mp + 295, mp + 286, mp + rxn_kegg_to_id['R01251'], mp + rxn_kegg_to_id['R00352'], mp + rxn_kegg_to_id['R01082'], mp + rxn_kegg_to_id['R01518'], mp + rxn_kegg_to_id['R09477'], mp + rxn_kegg_to_id['R07676'], mp + rxn_kegg_to_id['R01082'], mp + rxn_kegg_to_id['R00658'] ]
avrxns = fwrxns + rvrxns + [286, 189]

avoid_rxns = []
for trxn in avrxns:
    crs = set([kegg_to_id[e] for e in Currency])
    if (set(np.nonzero(rxnMat[trxn])[0]) & crs) == set(np.nonzero(rxnMat[trxn])[0]):
        avoid_rxns.append(trxn)

avoid_rxns = [rxn_id_to_kegg[e] for e in avoid_rxns]

for tind, trow in enumerate(stoich_matrix):
    if trow[kegg_to_id['C00718']] == 1 and tind < mp:
        avoid_rxns.append(rxn_id_to_kegg[tind])

fwrxns = [rxn_kegg_to_id[e] for e in can if e not in avoid_rxns + [ 'R03292', 'R03293', 'R11819', 'R04447', 'R04446', 'R03295', 'R03290', 'R03291', 'R02727', 'R00801', 'R05140', 'R11790', 'R04194', 'R00806', 'R02120', 'R01823', 'R00178', 'R11791', 'R11309', 'R00948', 'R00872', 'R09945', 'R04111', 'R01555',  'R01557', 'R01295', 'R11540', 'R00341', 'R00345', 'R00354', 'R00479', 'R00878', 'R09994', 'R02110']]
mp = int(len(stoich_matrix)/2)
rvrxns = [ mp + 295, mp + 286, mp + rxn_kegg_to_id['R01251'], mp + rxn_kegg_to_id['R00352'], mp + rxn_kegg_to_id['R01082'], mp + rxn_kegg_to_id['R01518'], mp + rxn_kegg_to_id['R09477'], mp + rxn_kegg_to_id['R07676'], mp + rxn_kegg_to_id['R01082'], mp + rxn_kegg_to_id['R00658'] ]
avrxns = fwrxns + rvrxns + [286, 189]
len( avrxns )<jupyter_output><empty_output><jupyter_text>Performing scope expansion<jupyter_code>from tqdm import tqdm 
calc_scopes = {}
rxn_scopes = {}
for media_id in tqdm( range( 1, len( kegg_media ) + 1 ) ):
    seedVec = np.zeros(len(rxnMat.T))
    media = kegg_media[ media_id - 1 ]
    seedVec[[kegg_to_id[e] for e in Currency + media]] = 1
    scopeMets, rxnProc = giveScope( rxnMat[ avrxns ], prodMat[ avrxns ], seedVec, sumRxnVec[ avrxns ])
    transrxn = np.zeros(len(stoich_matrix))
    transrxn[avrxns] = np.copy(rxnProc)
    rxnProc = np.copy(transrxn)
    calc_scopes[ media_id ] = scopeMets
    rxn_scopes[ media_id ] = rxnProc[:]<jupyter_output>100%|██████████████████████████████████████████████████████████████████████████████████| 75/75 [00:32<00:00,  2.28it/s]<jupyter_text>Calculating the scope of currency molecules.<jupyter_code>seedVec = np.zeros(len(rxnMat.T))
seedVec[[kegg_to_id[e] for e in Currency]] = 1
currScope, currRxn = giveScope( rxnMat[ avrxns ], prodMat[ avrxns ], seedVec, sumRxnVec[ avrxns ])

transrxn = np.zeros(len(stoich_matrix))
transrxn[avrxns] = np.copy(currRxn)
currRxn = np.copy(transrxn)
rxn_curr = set(np.nonzero( currRxn )[0])

currMets = set(np.nonzero( currScope )[0])

# 3 reactions can be performed using only the currency metabolites.
print( len( rxn_curr ) )<jupyter_output>3<jupyter_text>Removing the scope of all currency molecules to get final scopes for each medium.<jupyter_code>filt_scope = {}
filt_scope_size = {}
for ind in tqdm( range( 1, len( kegg_media ) + 1 ) ):
    rxnProc = np.zeros_like( currRxn )
    rxn_tm = set(np.nonzero( rxn_scopes[ind] )[0])
    rxnProc[ list( rxn_tm - rxn_curr ) ] = 1
    rxnProc[ 189 ] = 1
    seedVec = np.zeros(len(rxnMat.T))
    media = kegg_media[ ind - 1 ]
    seedVec[[kegg_to_id[e] for e in Currency + media]] = 1
    currScopeMets = np.copy( seedVec )
    currScopeMets = (np.dot(np.transpose(prodMat[ avrxns ]), rxnProc[avrxns]) + currScopeMets > 0) * 1
    currScopeMets[ [ kegg_to_id[e] for e in Currency ] ] = 0
    filt_scope[ ind ] = currScopeMets
    filt_scope_size[ ind ] = np.count_nonzero( currScopeMets )<jupyter_output>100%|██████████████████████████████████████████████████████████████████████████████████| 75/75 [00:03<00:00, 24.61it/s]<jupyter_text>Results: number of predicted metabolites for each single resource media<jupyter_code>{ cpd_string_dict[kegg_media[i][0]] :filt_scope_size[i+1] for i in range(len(kegg_media)) if len(kegg_media[i]) == 1}<jupyter_output><empty_output><jupyter_text>Checking the scope for a particular medium (example for medium no. 46)<jupyter_code>{id_to_kegg[e]: cpd_string_dict[id_to_kegg[e]] for e in np.nonzero( filt_scope[ 46 ] )[0] if id_to_kegg[e] not in Currency }<jupyter_output><empty_output><jupyter_text>Plotting the predicted number of metabolites with the number of carbon sources<jupyter_code>num_cs_media_dict = dict( zip( list( range( 1, len( kegg_media ) + 1 ) ), list( map( len, kegg_media ) ) ) )
import matplotlib.pyplot as plt
%matplotlib inline
fig, ax = plt.subplots(figsize=(10,7))
plt.scatter( list( map( len, kegg_media ) ), list( filt_scope_size.values() ) )
plt.xlabel( 'Number of carbon sources' )
plt.ylabel( 'Number of predicted metabolites' )<jupyter_output><empty_output><jupyter_text>Saving metabolite matrix  Calculating which metabolites were produced across all media (we will only save those). There are only ~100 out of 1000s of metabolites, and we are only interested in saving those that were at least produced once.<jupyter_code>cpds_scope = np.sort(uniqify(unlistify([ [e for e in np.nonzero( filt_scope[ sn ] )[0]] for sn in range(1, len(kegg_media) + 1) ])))

smap = {cpds_scope[e] : e+1 for e in range(len(cpds_scope))}
smap_inv = {value: key for key, value in smap.items()}<jupyter_output><empty_output><jupyter_text>There are 108 metabolites, in fact.<jupyter_code>len(smap)<jupyter_output><empty_output><jupyter_text>Merging custom KEGG metabolites Some KEGG metabolites, such as alpha and beta D glucose, are treated as separate metabolites by KEGG, but we want to merge them into one metabolite, so as not to double count them. The file merge_dict contains a manually curated mapping of which metabolites need to be merged into which others.<jupyter_code>mdf = pd.read_csv( 'data_metabolism/merge_dict.txt', sep='\t', header=None )
mdict = dict( mdf.values )
mdict<jupyter_output><empty_output><jupyter_text>First, we merge the metabolites and make sure if one of them is present in the scope, then it is counted. We do this by setting the first metabolite in the mapping to be present (equal to 1 in the scope).<jupyter_code>for mkey in mdict:
    for media_id in range( len( filt_scope ) ):
        if filt_scope[ media_id + 1 ][ smap_inv[ mkey ] ] == 1 and filt_scope[ media_id + 1 ][ smap_inv[ mdict[ mkey ] ] ] == 0:
            filt_scope[ media_id + 1 ][ smap_inv[ mdict[ mkey ] ] ] = 1<jupyter_output><empty_output><jupyter_text>Remove merged metabolites from the matrix and dicts Here, we are removing the second metabolite in the mapping from all scopes. This is to prevent double counting.<jupyter_code>pres_mets = { i+1 :np.nonzero(filt_scope[i+1])[0] for i in range(len(kegg_media))}

to_filt_out = [ smap_inv[ e ] for e in mdict ] 

filt_pres_mets = {}
for media_id in range( 1, len( filt_scope ) + 1 ):
    filt_pres_mets[ media_id ] = np.asarray( [ e for e in pres_mets[ media_id ] if e not in to_filt_out ] )<jupyter_output><empty_output><jupyter_text>Saving the final matrix after filtering The shape of the matrix is the number of media (75) by the number of metabolites (108). The first row and column are deliberately left blank.<jupyter_code>metm = np.zeros((len(kegg_media)+1, len(smap)+1))
for i in range(len(kegg_media)):
    to_pos = list(filt_pres_mets[i+1])
    for t in to_pos:
        metm[i+1][smap[t]] = 1
    if not (np.array([smap_inv[e] for e in np.where(metm[i+1])[0]]) == filt_pres_mets[i+1]).all():
        print('ERROR')
    
np.savetxt("data_metabolism/metabolic_matrix.csv", metm, delimiter=",")<jupyter_output><empty_output><jupyter_text>Creating a mapping between metabolite number and name for human readability<jupyter_code>scope_map_dict = {e+1:cpd_string_dict[id_to_kegg[smap_inv[e+1]]] for e in range(len(smap_inv))}
import csv

with open('data_metabolism/scope_map_table.csv', 'w') as f:
    w = csv.DictWriter(f, scope_map_dict.keys())
    w.writeheader()
    w.writerow(scope_map_dict)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1819625277865989, 'fraction_numerical': 0.05246109876151159, 'mean_word_length': 4.528792134831461, 'pattern_counts': {'":': 0, '<': 85, '<?xml version=': 0, '>': 85, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '613579', 'license': 'MIT', 'n_tokens_mistral': 6642, 'n_tokens_neox': 5908, 'n_words': 1309, 'path': 'metabolic_scope_expansion.ipynb', 'repo_name': 'mdalbello/Resource-diversity-relationship', 'size': '15745'}