<jupyter_start><jupyter_text>Author: <NAME>Date Created: 12 Jul 2018Distributed under the terms of the MIT License.  Visualize properties of a random subset of reaction systems in a directory for comparison- check:    - reaction components with Database website    - molecules with PUBCHEM/CHEBI    - molecule properties with PUBCHEM    - sequence with Database website    - sequence properties with UniPROT<jupyter_code>import glob
from ercollect import molecule as mol
from ercollect.molecule import molecule
from ercollect import rxn_syst
from ercollect.rxn_syst import reaction
import numpy as np
import random
from rdkit.Chem import Draw
from rdkit.Chem import AllChem as Chem
from rdkit.Chem.Draw import IPythonConsole
from rdkit.Chem.Draw import rdMolDraw2D
from rdkit.Chem.Draw.MolDrawing import MolDrawing, DrawingOptions
from IPython.display import clear_output
directory = '/home/atarzia/psp/screening_results/'
directory += 'kegg_111218/'
# directory += 'atlas_111218/'
# directory += 'new_reactions_kegg_atlas/'
# directory += 'new_reactions_sabio_wcharge/'
# directory += 'biomin_search_sabio_wcharge/'<jupyter_output><empty_output><jupyter_text>collect a random subset of X reactions<jupyter_code>no_ = 300
rs_to_test = []
no_rxns_in_total = len(glob.glob(directory+"*sRS*.gpkl"))
idx = np.random.randint(no_rxns_in_total, size=no_) # + 10000
generator = rxn_syst.yield_rxn_syst(output_dir=directory)
smiles_list = []
n_list = []
for i, rs in enumerate(generator):
    if i not in idx:
        continue
    if '-2_' not in rs.pkl:
        continue
    print('index:', i)
    print('pkl', rs.pkl)
    print('EC', rs.EC)
#     try:
#         print('reversible?', rs.reversible)
#     except AttributeError:
#         print('reversible?', 'unknown')
    if rs.skip_rxn is True:
        print('should be skipped?')
        try:
            print(rs.skip_reason)
        except AttributeError:
            pass
    else:
        print('----------------------------------------------')
        # check component properties
        smiles_list = []
        n_list = []
        for m in rs.components:
            print('-----------')
            print(m.name, '--', m.role)
            print('iupac:', m.iupac_name)
            if m.SMILES is not None:
                n_list.append(m.DB+' - '+m.name)
                smiles_list.append(m.SMILES)
#                 print('CHEBI ID:', m.chebiID)
                print('SMILES:', m.SMILES)
                print('PUBCHEM XlogP:', m.XlogP)
                print('PUBCHEM complexity:', m.complexity)
                print('RDKIT logP:', m.logP)
                print('RDKIT Synthetic accessibility:', m.Synth_score)   
                print('size:', m.mid_diam, 'angstrom')
        print('----------------------------------------------')
        print('change in complexity:', rs.delta_comp)
        print('change in synthetic accessibility:', rs.delta_sa)
        # check sequence properties
        print('----------------------------------------------')
        try:
            if rs.sequence is not None:
                print(rs.sequence)
                try:
                    print('uniprotID:', rs.UniprotID)
                except AttributeError:
                    pass
                print('add other sequence IDs for other DBs')
                print('sequence length:', len(rs.sequence))
                print('pI:', rs.pI)
                print('GRAVY:', rs.GRAVY)
                print('A index:', rs.A_index)
                print('I index:', rs.I_index)
                print('TM index:', rs.TM_index)
        except AttributeError:
            pass
    input('done?')
    clear_output()
m_list = [Chem.MolFromSmiles(i) for i in smiles_list]
Draw.MolsToGridImage(m_list, legends=n_list)<jupyter_output><empty_output><jupyter_text>Visualize a specific reaction system- includes search functions<jupyter_code>from ercollect.rxn_syst import reaction, get_RS, yield_rxn_syst
from rdkit.Chem import AllChem as Chem
from rdkit.Chem import Draw, Descriptors
directory = '/home/atarzia/psp/screening_results/'
directory += 'kegg_111218/'
# directory += 'atlas_111218/'
# directory += 'new_reactions_kegg_atlas/'
# directory += 'new_reactions_sabio_wcharge/'
# directory += 'biomin_search_sabio_wcharge/'
# get a list of RS with max_comp_size < XX and sequence != None
XX = 6.5
for rs in yield_rxn_syst(output_dir=directory):
    if rs.skip_rxn is True:
        continue
    if rs.max_comp_size is None:
        continue
    if rs.max_comp_size < XX:
        try:
            if rs.sequence is not None:
                print(rs.pkl, rs.TM_index, rs.A_index)
        except AttributeError:
            pass
pkl_name = 'sRS-2_10_1_1-KEGG-R09735.gpkl'
rs = get_RS(directory+pkl_name, output_dir=directory, verbose=True)
rs.__dict__
smiles_list = []
n_list = []
for m in rs.components:
    print(m.name, m.mid_diam, m.logP, m.role, m.Synth_score, m.complexity, m.pkl)
    print(m.name, m.DB)
    print('---', Descriptors.ExactMolWt(Chem.MolFromSmiles(m.SMILES)))
    n_list.append(m.name)
    smiles_list.append(m.SMILES)
    MOL = Chem.MolFromSmiles(m.SMILES)
#     Draw.MolToFile(MOL, fileName=m.name+'.svg')
m_list = [Chem.MolFromSmiles(i) for i in smiles_list]
Draw.MolsToGridImage(m_list, legends=n_list)<jupyter_output>C19848 0 -1.3538999999999997 reactant 5.567428995960881 1440 /home/atarzia/psp/molecule_DBs/atarzia/ATRS_3648.gpkl
C19848 KEGG
--- 724.0648315320001
C06232 5.593 -1.3541 reactant 6.068191459211439 81.3 /home/atarzia/psp/molecule_DBs/atarzia/ATRS_3782.gpkl
C06232 KEGG
--- 163.900716744
C18237 0 -0.24120000000000014 product 5.576142290786926 817 /home/atarzia/psp/molecule_DBs/atarzia/ATRS_3649.gpkl
C18237 KEGG
--- 522.8918991740001
C00020 8.859 -1.8630000000000009 product 3.8048894493072236 481 /home/atarzia/psp/molecule_DBs/atarzia/ATRS_73.gpkl
C00020 KEGG
--- 347.0630844179999
C00001 2.843 -0.8247 product 5.868476923076923 0 /home/atarzia/psp/molecule_DBs/atarzia/ATRS_7.gpkl
C00001 KEGG
--- 18.010564684<jupyter_text>Analyse all RS for a certain skip_reaction reason<jupyter_code>import os
directory = '/home/atarzia/psp/screening_results/'
directory += 'kegg_111218/'
# directory += 'atlas_111218/'
# directory += 'new_reactions_kegg_atlas/'
# directory += 'new_reactions_sabio_wcharge/'
# directory += 'biomin_search_sabio_wcharge/'
count = 0
count_total = 0
max_logP = 0
min_logS = 0
for rs in rxn_syst.yield_rxn_syst(output_dir=directory):
    count_total += 1
    try:
        if rs.reversible is True:
            count += 1
        max_logP = max([max_logP, rs.max_logP])
        min_logS = min([min_logS, rs.min_logS])
    except AttributeError:
        pass
#     if rs.skip_rxn is True:
#         if rs.skip_reason == 'one component has no SMILES':
#             count += 1
#             print(rs.pkl)
#             os.system('rm '+rs_dir+rs.pkl)
max_logP, min_logS
print(count, count_total, count/count_total * 100)<jupyter_output>0 12182 0.0<jupyter_text>Analyse the number of times each skip_reason is used<jupyter_code>from ercollect import rxn_syst
from rdkit.Chem import Descriptors
import matplotlib.pyplot as plt
reasons = {'one component failed resolution': 0, 'SABIO E-ID is for mutant': 0, 
           'SABIO R-ID not found': 0,  'DNA present - SABIO has a bug': 0,
           'CHEBI ID of a component not available': 0, 'a component is in skip_names': 0,
           'one component is ?': 0, 'No result for KEGG URL search - likely outdated': 0,
           'CHEBI ID not available for one component': 0,
           'one component has invalid SMILES': 0, 'one component has no SMILES': 0, 
           'one component has wildcard SMILES': 0,
           'one component has no molecule - rxn is incomplete or generic': 0,
           'one component could not have diameter calculated': 0,
           'KEGG ID could not be converted to MOL': 0,
           'KEGG ID gave generic structure': 0,
           'KEGG rxn includes polymeric species': 0}
directory = '/home/atarzia/psp/screening_results/'
directory += 'kegg_111218/'
# directory += 'new_reactions_atlas_111218/'
# directory += 'new_reactions_kegg_atlas/'
# directory += 'new_reactions_sabio_wcharge/'
# directory += 'biomin_search_sabio_wcharge/'
count_total = 0
count_skipped = 0
for rs in rxn_syst.yield_rxn_syst(output_dir=directory):
#     if 'KEGG' not in rs.pkl:
#         continue
    count_total += 1
    if rs.skip_rxn is True:
        count_skipped += 1
        reasons[rs.skip_reason] += 1
#         if rs.skip_reason == 'CHEBI ID not available for one component':
#             print(rs.__dict__)
#             input()
print(count_skipped, count_total, count_skipped/count_total * 100)
reasons
MWs =[]
for rs in rxn_syst.yield_rxn_syst(output_dir=directory):
    rs_MW = []
    if rs.skip_rxn is True:
        if rs.skip_reason == 'one component could not have diameter calculated':
            if rs.components is not None:
                for m in rs.components:
                    MW = Descriptors.MolWt(m.mol)
                    rs_MW.append(MW)
                    MWs.append(MW)
                if max(rs_MW) < 500:
                    print(rs.pkl)<jupyter_output>sRS-1_13_11_XX-KEGG-R09163.gpkl
sRS-1_13_11_XX-KEGG-R09206.gpkl
sRS-1_14_11_12-KEGG-R07184.gpkl
sRS-1_6_XX_XX-KEGG-R09218.gpkl
sRS-4_1_1_XX-KEGG-R09164.gpkl
sRS-4_1_1_XX-KEGG-R09207.gpkl
sRS-4_1_1_XX-KEGG-R09208.gpkl
sRS-4_2_3_127-KEGG-R10006.gpkl
sRS-4_2_3_42-KEGG-R06313.gpkl
sRS-4_2_3_80-KEGG-R10726.gpkl
sRS-XX_XX_XX_XX-KEGG-R06314.gpkl
sRS-XX_XX_XX_XX-KEGG-R06316.gpkl
sRS-XX_XX_XX_XX-KEGG-R06317.gpkl
sRS-XX_XX_XX_XX-KEGG-R08420.gpkl<jupyter_text>Analyse a list of RS- example given here of RS with max_comp_size < some threshold- using threshold 4.2 angstrom 10/12/18<jupyter_code>from ercollect.rxn_syst import reaction, get_RS, yield_rxn_syst
from rdkit.Chem import AllChem as Chem
from rdkit.Chem import Draw
RS_des = ['sRS-1_11_1_21-KEGG-R00009.gpkl', 'sRS-1_11_1_21-KEGG-R00602.gpkl',
          'sRS-1_11_1_8-KEGG-R02810.gpkl', 'sRS-1_13_11_49-KEGG-R05721.gpkl',
          'sRS-1_15_1_1-KEGG-R00275.gpkl', 'sRS-1_16_3_1-KEGG-R00078.gpkl',
          'sRS-1_1_3_13-KEGG-R00608.gpkl', 'sRS-1_2_2_4-KEGG-R00276.gpkl',
          'sRS-1_2_98_1-KEGG-R00614.gpkl', 'sRS-1_4_3_21-KEGG-R06154.gpkl',
          'sRS-1_7_1_14-KEGG-R09809.gpkl', 'sRS-1_7_3_6-KEGG-R10230.gpkl',
          'sRS-1_8_3_4-KEGG-R01851.gpkl', 'sRS-3_13_1_5-KEGG-R10534.gpkl',
          'sRS-3_13_1_5-KEGG-R10535.gpkl', 'sRS-3_13_1_5-KEGG-R10538.gpkl',
          'sRS-3_5_1_49-KEGG-R00524.gpkl', 'sRS-3_5_5_8-KEGG-R05780.gpkl',
          'sRS-3_5_5_XX-KEGG-R00152.gpkl', 'sRS-4_2_1_112-KEGG-R05380.gpkl',
          'sRS-4_2_1_66-KEGG-R01408.gpkl', 'sRS-4_99_1_2-KEGG-R09339.gpkl',
          'sRS-XX_XX_XX_XX-KEGG-R00793.gpkl', 'sRS-XX_XX_XX_XX-KEGG-R09094.gpkl',
          'sRS-XX_XX_XX_XX-KEGG-R09996.gpkl']
directory = '/home/atarzia/psp/screening_results/'
directory += 'kegg_111218/'
# directory += 'atlas_111218/'
# directory += 'new_reactions_kegg_atlas/'
# directory += 'new_reactions_sabio_wcharge/'
# directory += 'biomin_search_sabio_wcharge/'
count = 0
for pkl in RS_des:
    print(pkl)
    rs = get_RS(directory+pkl, output_dir=directory, verbose=False)
    smiles_list = []
    n_list = []
    for m in rs.components:
        if m.mid_diam == 0:
            print(m.name, m.mid_diam, m.SMILES)
            print(m.pkl)
    print('---------------------------------------------')
    count += 1
print(count)
m.mol<jupyter_output><empty_output><jupyter_text>Modify some attribute of all RS<jupyter_code>directory = '/home/atarzia/psp/screening_results/'
directory += 'kegg_111218/'
# directory += 'atlas_111218/'
# directory += 'new_reactions_kegg_atlas/'
# directory += 'new_reactions_sabio_wcharge/'
# directory += 'biomin_search_sabio_wcharge/'
rs.__dict__
for rs in rxn_syst.yield_rxn_syst(output_dir=rs_dir):
    rs.mol_collected = False
    rs.max_logP = None
    rs.max_XlogP = None
    rs.max_logS = None
    rs.min_logP = None
    rs.min_XlogP = None
    rs.min_logS = None
    rs.save_object(rs_dir+rs.pkl)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12765781418052452, 'fraction_numerical': 0.06717961446181848, 'mean_word_length': 2.8669225847728725, 'pattern_counts': {'":': 0, '<': 31, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 3, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '36110', 'license': 'MIT', 'n_tokens_mistral': 5079, 'n_tokens_neox': 4475, 'n_words': 891, 'path': 'archived_code/notebooks/analyse_RS.ipynb', 'repo_name': 'andrewtarzia/enzyme_screen', 'size': '12094'}