<jupyter_start><jupyter_text>In this notebook we explore using random forest to determine activity values of molecules (using regression) and predicting an active/inactive molecule (using classification).  First we load the fingerprints from RDkit.<jupyter_code>import pandas as pd
import itertools
from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor
#from sklearn.model_selection import train_test_split
from sklearn.inspection import permutation_importance
from sklearn.metrics import classification_report, confusion_matrix
import matplotlib.pyplot as plt
import numpy as np

data = pd.read_pickle('../processed_data/combined_dataset.pkl')
N=len(data.loc[data['target'] == 'TMPRSS2']) #number of compounds tested against TMPRSS2 itself
fp_dark=pd.read_pickle('../processed_data/DarkChemicalMatter_processed.pkl.gz')
fp_dark=np.stack(fp_dark['morgan_fingerprint'])<jupyter_output><empty_output><jupyter_text>We have two sets of activity values since not all compound assays were performed against TMPRSS2 (the protease of interest in this study).  Thus, some compounds' activity values have to be scaled via a correlation of the protease they were originally purposed for, against TMPRSS2.  We save the binary values of whether a compound was deemed active/inactive as "activity", and the numeric activity value as "ac_merged".  Also we note if the values are taken from the "scaled" dataset.  Ac values are converted to negative log10.  Finally we convert everything to numpy array to input into random forest. Now we split the training and the test data.  To test how well the random forest works, we split the compounds with known TMPRSS2 activity (that is, their recorded ac value is from an assay against TMPRSS2) between the training and test sets.  Here I've used a 2/3 training 1/3 test split, which can be changed by changing variable "n_test".  The indices assigned by the split is randomly shuffled each time the code is run.<jupyter_code>fps_merged=np.stack(data[data['acvalue_target'].notnull()]['morgan_fingerprint'])
fps_merged_scaled=np.stack(data[data['acvalue_scaled_to_tmprss2'].notnull()]['morgan_fingerprint'])
ac_merged_scaled=np.asarray(data[data['acvalue_scaled_to_tmprss2'].notnull()]['acvalue_scaled_to_tmprss2'])
act_list=list(data[data['acvalue_target'].notnull()]['activity_target'])
act_list_2=list(data[data['acvalue_scaled_to_tmprss2'].notnull()]['activity_target'])
activity_scaled=[]
activity=[]
for i in range(len(act_list)):
    if act_list[i] == 'Active':
        activity.append(1)
    else:
        activity.append(0)
for i in range(len(act_list_2)):
    if act_list_2[i] == 'Active':
        activity_scaled.append(1)
    else:
        activity_scaled.append(0)
#alternatively
ac_merged=np.asarray(data[data['acvalue_target'].notnull()]['acvalue_target'])
activity=np.array(activity)
ac_merged=-np.log10(ac_merged)
activity_scaled=np.array(activity_scaled)
ac_merged_scaled=-np.log10(ac_merged_scaled)
#split train and test 
#split the tmprss2 active compounds into training and test
import random
s=np.arange(N)
random.seed(33)#set seed for reproducibility
random.shuffle(s)
n_test=round(N/3) #percentage to withold for testing
cut=n_test
test=s[0:cut]
train=s[cut::]<jupyter_output><empty_output><jupyter_text>Here we generate training and test sets for the classifier, which incorporates the "dark chemical matter" which is known to be inactive.  We split the dark chemical matter evenly between the test and training sets (same amount of data points as from the assays).  You can skip the following code block if not running the classifier.  The purpose of the classifier is to do a preliminary "filtering" of active/inactive molecules for the regressor.<jupyter_code>s=np.arange(len(fp_dark))
random.seed(33)
random.shuffle(s)
test_dark=s[0:cut]
train_dark=s[cut:len(activity_scaled)]
X_test=np.concatenate((fps_merged_scaled[test], fp_dark[test_dark]))
X_train=np.concatenate((fps_merged_scaled[train], fps_merged_scaled[N::], fp_dark[train_dark]))
y_train=np.concatenate((activity_scaled[train], activity_scaled[N::], np.zeros([len(train_dark)])))
y2_train=np.concatenate((ac_merged_scaled[train], ac_merged_scaled[N::], np.zeros([len(train_dark)])))
y_test=np.concatenate((activity_scaled[test], np.zeros([len(test_dark)])))
y2_test=np.concatenate((ac_merged_scaled[test],np.zeros([len(test_dark)])))

X_test_u=np.concatenate((fps_merged[test], fp_dark[test_dark]))
X_train_u=np.concatenate((fps_merged[train], fps_merged[N::], fp_dark[train_dark]))
y_train_u=np.concatenate((activity[train], activity[N::], np.zeros([len(train_dark)])))
y2_train_u=np.concatenate((ac_merged[train], ac_merged[N::], np.zeros([len(train_dark)])))
y_test_u=np.concatenate((activity[test], np.zeros([len(test_dark)])))
y2_test_u=np.concatenate((ac_merged[test],np.zeros([len(test_dark)])))
#generate big training set
#X_train_all=np.concatenate((fps_merged_scaled, fp_dark[train_dark]))
X_train_all=np.concatenate((fps_merged, fp_dark[train_dark]))
y2_train_all=np.concatenate((ac_merged_scaled, np.zeros([len(train_dark)])))
#y_train_all=np.concatenate((activity_scaled, np.zeros([len(train_dark)])))
y_train_all=np.concatenate((activity, np.zeros([len(train_dark)])))<jupyter_output>425
403
394<jupyter_text>Here we run the classifier for the test case (where our test set has some compounds with known TMPRSS2 activity).<jupyter_code>###YOU CAN SKIP THIS STEP IF NOT RUNNING THE CLASSIFIER###
rf=RandomForestClassifier(verbose=2, n_estimators=100, random_state=111)
rf.fit(np.asarray(X_train_u), np.asarray(y_train_u))
y_pred=rf.predict(X_test_u)
print(rf.score(X_test_u, y_test_u))<jupyter_output>1.0<jupyter_text>and we see how well it does in terms of precision and recall.<jupyter_code>#get the index of the features (of RDkit fingerprint) that were important
importances = rf.feature_importances_
featurenums = np.array([str(x).zfill(2) for x in range(len(importances))])
indices = np.argsort(importances)[::-1][0:25]#get the 25 most important features
plt.title('Feature Importances (train set)')
plt.bar(range(len(indices)), importances[indices], align='center')
plt.ylabel('Relative Importance')
plt.xticks(range(len(indices)), featurenums[indices], rotation=90)
plt.show()
#output metrics
print(classification_report(y_test, y_pred))<jupyter_output><empty_output><jupyter_text>Now we train the classifier on the full dataset.  Note that I'm using the unscaled dataset here to keep more possible actives.<jupyter_code>rf=RandomForestClassifier(verbose=2, n_estimators=100, random_state=111)
rf.fit(np.asarray(X_train_all), np.asarray(y_train_all))<jupyter_output>[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:    0.0s remaining:    0.0s<jupyter_text>REGRESSIONHere we again take the same random splits of training and test set (set by the variable "n_test") as a proportion of compounds in the ground truth TMPSSR2 assay data, but we don't add the dark data, which can distort the results of regression. We train both the "scaled" version with scaled data and the version with unscaled data.<jupyter_code>import random
from sklearn.metrics import mean_squared_error
n_bootstrap = 50
mse_unscaled=np.zeros(n_bootstrap)
mse_scaled=np.zeros(n_bootstrap)
random.seed(33)
for i in range(n_bootstrap):
        s=np.arange(N)
        #set seed for reproducibility
        random.shuffle(s)
        n_test=round(N/3) #percentage to withold for testing
        cut=n_test
        test=s[0:cut]
        train=s[cut::]
        X_test_reg=fps_merged[test]
        X_train_reg=np.concatenate((fps_merged[train], fps_merged[N::]))
        y2_train_reg=np.concatenate((ac_merged[train], ac_merged[N::]))
        y2_test_reg=ac_merged[test]
        nan_array=np.isnan(y2_train_reg)
        X_train_reg=X_train_reg[~nan_array]
        y2_train_reg=y2_train_reg[~nan_array]
        #and generate for scaled arrays
        #don't use the dark data for regressor, hence 'pos only'
        X_train_reg_scaled=np.concatenate((fps_merged_scaled[train], fps_merged_scaled[92::]))
        y2_train_reg_scaled=np.concatenate((ac_merged_scaled[train], ac_merged_scaled[92::]))
        rg_pos_only=RandomForestRegressor(verbose=0, n_estimators=100, random_state=111)
        rg_pos_only_scaled=RandomForestRegressor(verbose=0, n_estimators=100, random_state=111)
        rg_pos_only.fit(np.asarray(X_train_reg), np.asarray(y2_train_reg))
        rg_pos_only_scaled.fit(np.asarray(X_train_reg_scaled), np.asarray(y2_train_reg_scaled))
        y_pred_regr_pos=rg_pos_only.predict(X_test_reg)
        mse_unscaled[i]=mean_squared_error(y2_test_reg, y_pred_regr_pos)
        y_pred_regr_scaled_pos=rg_pos_only_scaled.predict(X_test_reg)
        mse_scaled[i]=mean_squared_error(y2_test_reg, y_pred_regr_scaled_pos)<jupyter_output><empty_output><jupyter_text>Examine mean squared error (MSE) for both scaled and unscaled dataset, with 50 bootstrap random resamplings of the training set.  The scaled version appears to do slightly better (lower mean MSE) so we will use that.  Can someone please check if "%store -r mse_unscaled" will retrieve the variable on their computer?<jupyter_code>plt.hist(mse_unscaled, bins=33, label='unscaled')
plt.hist(mse_scaled, bins=33, label='scaled')
plt.legend()
print(np.mean(mse_unscaled))
print(np.mean(mse_scaled))<jupyter_output>0.6522259685038808
0.34818155556849
Stored 'mse_unscaled' (ndarray)
Stored 'mse_scaled' (ndarray)<jupyter_text>Plotting the predicted vs actual regression values (-log ac) from the validation set.  The red line denotes 1-to-1 correlation.  For smaller activity values (large negative log) the regressor appears to overpredict the AC values somewhat.<jupyter_code>plt.plot(y2_test_reg, y_pred_regr_pos, '*')
plt.plot(y2_test_reg, y2_test_reg,'r--')
plt.set_xlabel('')<jupyter_output><empty_output><jupyter_text>Here we select the 25 most important features from the morgan fingerprints.<jupyter_code>importances = rg_pos_only_scaled.feature_importances_
featurenums = np.array([str(x).zfill(2) for x in range(len(importances))])
indices = np.argsort(importances)[::-1][0:25]#get the 25 most important features
plt.title('Feature Importances (train set)')
plt.bar(range(len(indices)), importances[indices], align='center')
plt.ylabel('Relative Importance')
plt.xticks(range(len(indices)), featurenums[indices], rotation=90)
plt.show()<jupyter_output><empty_output><jupyter_text>When we're satisfied with tuning the random forest regressor, we can build the regressor with all of the data combined to run on an unknown screening set.  First we load the screening dataset<jupyter_code>screening_fp=pd.read_pickle('../processed_data/screening_data_processed.pkl')
screening_fp.head()
molnames_all=list(screening_fp['name'])
molnames_all=np.asarray(molnames_all)
screening=np.stack(screening_fp['morgan_fingerprint'])
sources=list(screening_fp['source'])
bitinfo=screening_fp['bitinfo']<jupyter_output><empty_output><jupyter_text>Then we run the classifier and truncate the screening dataset to classified actives.<jupyter_code>y_pred=rf.predict(screening)
print(len(y_pred[y_pred>0])) #number of active compounds
screening_fp['active'] = y_pred
screening_fp_trunc=screening_fp.loc[screening_fp['active']>0.0]<jupyter_output>[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:    0.0s remaining:    0.0s<jupyter_text>Then build the regressor based on the scaled dataset<jupyter_code>rg_all=RandomForestRegressor(verbose=2, n_estimators=100, random_state=111)
rg_all.fit(np.asarray(fps_merged_scaled), np.asarray(ac_merged_scaled))<jupyter_output>[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:    0.0s remaining:    0.0s<jupyter_text>Run regression on screening dataset that was pre-filtered by the lassifier.<jupyter_code>predicted_activity=rg_all.predict(screening[y_pred>0])
screening_fp_trunc['predicted_activity(log10)'] = predicted_activity
importances = rg_all.feature_importances_
featurenums = np.array([str(x).zfill(2) for x in range(len(importances))])
indices = np.argsort(importances)[::-1][0:25]#get the 25 most important features
plt.title('Feature Importances (train set)')
plt.bar(range(len(indices)), importances[indices], align='center')
plt.ylabel('Relative Importance')
plt.xticks(range(len(indices)), featurenums[indices], rotation=90)
plt.show()<jupyter_output><empty_output><jupyter_text>Predict the activities from screening library and print the top 20 hits.  We also plot the predicted activity values color coded by screening library.<jupyter_code>from seaborn import violinplot 
violinplot(x='source', y='predicted_activity(log10)', data=screening_fp_trunc)<jupyter_output><empty_output><jupyter_text>Drop the duplicate tautomers in the top hits, keeping the "first" or lowest predicted activity value.<jupyter_code>from rdkit.Chem import MolToSmiles
from indigo import *
indigo = Indigo()
canonical_smiles=[]
for smile in screening_fp_trunc['smiles']:
    mol1 = indigo.loadMolecule(smile)
    mol1.aromatize()
    canonical_smiles.append(mol1.canonicalSmiles())
test_smiles=screening_fp_trunc['smiles']
screening_fp_trunc['smiles']=canonical_smiles
screening_fp_trunc=screening_fp_trunc.sort_values(by='predicted_activity(log10)', ascending=False)
screening_fp_trunc['name'] =screening_fp_trunc['name'].str.upper() 
screening_fp_trunc=screening_fp_trunc.drop_duplicates(subset=['name'])
screening_fp_trunc=screening_fp_trunc.drop_duplicates(subset=['smiles'])
screening_fp_trunc['RMSE']=np.sqrt(np.mean(mse_scaled))

screening_fp_trunc.loc[:, ["name", "smiles", "predicted_activity(log10)", "RMSE"]].to_csv('../results/random_forest_results.csv')
screening_fp_trunc.loc[:, ["name", "smiles", "predicted_activity(log10)"]].head(25)<jupyter_output><empty_output><jupyter_text>Information in the bits appears to be set to 1 when the structure is present in one of the compound's atoms, 0 if not.  Many of the features deemed "important" by the random forest are actually lack of structures (0 in a bit location), so here we try to determine what positive and/or negative structures were important by retrieving these from the bit info of our top hits.<jupyter_code>bitinfo=screening_fp_trunc['bitinfo']
keys1=[bitinfo[key] for key in bitinfo.keys()]
rdkit_mols_draw=[]
bitinfos_draw=[]
compound_names=[]
ilocs=[]
for featurenum in featurenums[indices]:
    flag=0        
    location=0
    while flag == 0:
        ibit=screening_fp_trunc.iloc[location,:]['bitinfo']
        if int(featurenum) in list(ibit.keys()):
            compound_names.append(screening_fp_trunc.iloc[location,:]['name'])
            rdkit_mols_draw.append(screening_fp_trunc.iloc[location,:]['rdkit_mol'])
            bitinfos_draw.append(screening_fp_trunc.iloc[location,:]['bitinfo'])
            ilocs.append(location)
            flag=1
        else:
            location=location+1
          #  print('still searching for '+featurenum)<jupyter_output><empty_output><jupyter_text>Now we draw the snippets of the relevant "top hits".<jupyter_code>from rdkit.Chem import Draw
from IPython.display import SVG, display
for i,featurenum in enumerate(featurenums[indices]):
    print('found the bit in compound ' + str(ilocs[i])+ ' aka '+compound_names[i])
    if int(ilocs[i])<= 10:
        print('it is likely it is a positive feature.')
    else:
        print('it is likely to be a negative feature.')
    svg = Draw.DrawMorganBit(rdkit_mols_draw[i], int(featurenum), bitinfos_draw[i], useSVG=True)
    display(svg)<jupyter_output>found the bit in compound 67 aka N-ALPHA-(2-NAPHTHYLSULFONYL)-N(3-AMIDINO-L-PHENYLALANINYL)ISOPIPECOLINIC ACID METHYL ESTER
it is likely to be a negative feature.
{'fraction_non_alphanumeric': 0.11175616835994194, 'fraction_numerical': 0.015334132643402536, 'mean_word_length': 6.892430278884462, 'pattern_counts': {'":': 0, '<': 67, '<?xml version=': 0, '>': 69, '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': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '317773', 'license': 'MIT', 'n_tokens_mistral': 5673, 'n_tokens_neox': 5352, 'n_words': 1229, 'path': 'notebooks/random_forest_clean.ipynb', 'repo_name': 'srensi/tmprss2', 'size': '15847'}