<jupyter_start><jupyter_text>Querying the modelThis notebook provides an illustration of a workflow for querying the model with a compound of interest, to obtain likely promiscuous enzymes that can function on it. Once model development is at a satisfactory place in terms of performance, this workflow could be used as a framework for building production code for interacting with the model Inputs:- MASTER_DF.csv    - Should contain the full database of enzymes included in the model, along with paired compound information for each enzyme- Query compound    - Right now this needs to be a PubChem SID number     Outputs:- A ranked list of potential promiscuous enzymes and modeled probability values<jupyter_code># imports

from metamoles import *

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

import sklearn
from sklearn import linear_model
from sklearn.model_selection import train_test_split
from sklearn.metrics import roc_auc_score
from sklearn.metrics import roc_curve
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Step 1 Read in master dataframeThis dataframe has been curated, pre-structured, and stored for use with this workflow.<jupyter_code># read in the master_df

master_df = pd.read_csv('../datasets/MASTER_DF.csv') 
print(master_df.shape)
master_df.head()<jupyter_output>(6879, 24)<jupyter_text>Step 2  Get query SMILES string & pair query compound with each unique enzyme in the master DataFrameExample: PubChem SID 3480<jupyter_code># pair_query_compound_with_enzymes()

updated_df = pair_query_compound(master_df, 'entry', 'PubChem', 'SMILES', '3480')
print(updated_df.shape)
updated_df.head()<jupyter_output>(7395, 3)<jupyter_text>Step 3 Calculate molecular distances between products of the same enzymeThe model is based on the fact that we expect these distances to be closer for reactive enzyme-product pairs than for non-reactive enzyme-product pairs<jupyter_code># calculate_dist() with appended dataset

distance_df = calculate_dist(updated_df)
print(distance_df.shape)
distance_df.head()<jupyter_output>(7395, 6)<jupyter_text>Step 4 Get dummy variables to represent enzyme classWe expect that many enzyme properties could be predictive features for this model. Enzyme class should encapsulate many of these features at a high level.<jupyter_code># remove any rows that are not the query compound

reduced_df = distance_df[distance_df['PubChem'] == '3480']
# binarize_enzyme_class()

query_df = binarize_enzyme_class(reduced_df, 'entry')
query_df = query_df.reset_index(drop=True)
print(query_df.shape)
query_df.head()<jupyter_output>(516, 13)<jupyter_text>Step 5 Add in compound features with RDKitThis step uses the RDKit packages to generate descriptive features of the reaction product compounds<jupyter_code>query_df = create_cpd_info(query_df)
print(query_df.shape)
query_df.head()<jupyter_output>(516, 22)<jupyter_text>Step 6 Re-Instantiate model<jupyter_code># create new dataframe with desired features for training model
# current version of the model uses the following features
feature_df = master_df[['dist', 'enzyme_class_1', 'enzyme_class_2', 'enzyme_class_3',
       'enzyme_class_4', 'enzyme_class_5', 'enzyme_class_6', 'enzyme_class_7',
       'n_O', 'n_N', 'n_S', 'n_X', 'DoU']]

# shape array for regression
features = np.array(feature_df) 
reactions = list(master_df['reacts'])

feature_train, feature_test, reaction_train, reaction_test = train_test_split(features, reactions,
                                                  test_size=0.20, random_state=42)

model_1 = linear_model.LogisticRegression(solver='liblinear', penalty='l1', random_state=1, class_weight='balanced')
model_1.fit(feature_train, np.ravel(reaction_train))
# Test the model

logit_roc_auc = roc_auc_score(reaction_test, model_1.predict(feature_test))

fpr, tpr, thresholds = roc_curve(reaction_test, model_1.predict_proba(feature_test)[:,1])

plt.figure()
plt.plot(fpr, tpr, label='Logistic Regression (area = %0.2f)' % logit_roc_auc)
plt.plot([0, 1], [0, 1],'r--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver operating characteristic')
plt.legend(loc="lower right")
plt.show()

score = model_1.score(feature_test, reaction_test)
print('The model score is ' + str(round(score*100, 2)) + '%.')<jupyter_output><empty_output><jupyter_text>Step 7  Use model to predict reactivity of pairs<jupyter_code># select query features
# must match features for model training

query_feat_df = query_df[['dist', 'enzyme_class_1', 'enzyme_class_2', 'enzyme_class_3',
       'enzyme_class_4', 'enzyme_class_5', 'enzyme_class_6', 'enzyme_class_7',
       'n_O', 'n_N', 'n_S', 'n_X', 'DoU']]
# query the model

predictions = model_1.predict(query_feat_df) 

pred = model_1.predict_proba(query_feat_df)

prediction_values = pd.DataFrame(pred)
model_descriptive_df = pd.DataFrame()
model_descriptive_df['0']=prediction_values[0]
model_descriptive_df['1']=prediction_values[1]
# create a dataframe of model predictions

prediction_df = pd.merge(model_descriptive_df, query_df, left_index=True, right_index=True) 
print(prediction_df.shape)
prediction_df.head()
# sort the data so "Yes, it will react" predictions are at the top

prediction_df = prediction_df.sort_values(by=['1'], ascending=False)
prediction_df.head()<jupyter_output><empty_output><jupyter_text>Full Model ValidationThe bit of code below is used to test the model on all of the validation compounds. It runs each of the compounds set aside prior to model training through the model query workflow outlined above, and outputs the results to a csv for further analysis.<jupyter_code># read in validation data

validation_df = pd.read_csv('../datasets/validation_data/validation_data.csv')
# read in smiles string data
smiles_df = pd.read_csv('../datasets/PubChem_compounds_SMILES.csv', dtype='str')
# join on validation_df to get sid numbers
valid_df = pd.merge(validation_df, smiles_df, left_on='product', right_on='KEGG', how='left')
valid_df = valid_df.drop(columns=['Unnamed: 0_x', 'product', 'Unnamed: 0_y', 'entry_y', 'CID', 'SMILES'])
valid_df = valid_df.rename(columns={'entry_x': 'true_enzyme'})
valid_df = valid_df.drop_duplicates()
valid_df = valid_df.dropna(subset=['PubChem'])
valid_df = valid_df.reset_index(drop=True)
print(valid_df.shape)
valid_df.head()
# run each of these validation compounds through the query_model() function
# WARNING: this takes a very long time to run as is

headers = ['true_enzyme', 'KEGG', 'PubChem', 'prediction_rank', 'prediction_prob']
valid_pred_df = valid_df.reindex(columns=headers)

cpd_validation_list = valid_pred_df['PubChem'].unique().tolist()
# cpd_validation_list = ['135626334', '254741367']

for cpd in cpd_validation_list:
    try:
        prediction_df = query_model(master_df, cpd)
    except:
        pass
    handle = '../datasets/validation_data/{}.csv'.format(cpd)
    prediction_df.to_csv(handle, index=False)
    for index, row in valid_pred_df.loc[valid_pred_df.loc[:, 'PubChem'] == cpd].iterrows():
        enzyme = row['true_enzyme']
        for _, pred_row in prediction_df.loc[prediction_df.loc[:,'entry'] == enzyme].iterrows():
            valid_pred_df.loc[index, 'prediction_rank'] = pred_row['rank']
            valid_pred_df.loc[index, 'prediction_prob'] = pred_row['p_reacts']
    valid_pred_df.to_csv('validation_summary.csv', index=False)

valid_pred_df.head()
valid_pred_df.head()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10418332001065814, 'fraction_numerical': 0.015720756727950973, 'mean_word_length': 6.176864244741874, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 32, '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': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '290334', 'license': 'MIT', 'n_tokens_mistral': 2579, 'n_tokens_neox': 2410, 'n_words': 663, 'path': 'notebooks/Backup/model_query_workflow.ipynb', 'repo_name': 'metamoles/metamoles', 'size': '7506'}