<jupyter_start><jupyter_text>Install a few more packages required to run this notebook<jupyter_code>!pip install tqdm seaborn mols2grid sklearn pingouin
import useful_rdkit_utils as uru
import pandas as pd
from rdkit import Chem
from tqdm.auto import tqdm
import seaborn as sns
import mols2grid
from sklearn.model_selection import train_test_split
from sklearn.ensemble import ExtraTreesRegressor
from sklearn.metrics import r2_score
import numpy as np
import pingouin<jupyter_output><empty_output><jupyter_text>0. Preliminary Setup<jupyter_code>tqdm.pandas()
mol_buff = """Nc1ccc(S(=O)(=O)Nc2ccccn2)cc1 32842
CCC(C)C1(CC)C(=O)[N-]C(=O)NC1=O.[Na+] 674933
Cl.N=C(N)N 674679
CC1=CC(=O)c2ccccc2C1=O 19344
Cn1c(=O)c2[nH]cnc2n(C)c1=O.Cn1c(=O)c2[nH]cnc2n(C)c1=O.NCCN 794445
C[C@]1(O)CC[C@H]2[C@@H]3CCC4=CC(=O)CC[C@]4(C)[C@H]3CC[C@@]21C 328400
C=C1CC[C@H](O)C/C1=C/C=C1\CCC[C@@]2(C)[C@H]1CC[C@@H]2[C@H](C)/C=C/[C@H](C)C(C)C 430124
Cc1cc(OP(=O)([O-])[O-])c2ccccc2c1OP(=O)([O-])[O-].O.O.O.O.O.O.[Na+].[Na+].[Na+].[Na+] 2197599
CC/C(=C(/CC)c1ccc(O)cc1)c1ccc(O)cc1 860
Nc1ccc(S(=O)(=O)Nc2ncccn2)cc1 2131"""
mol_buff_list = [x.split() for x in mol_buff.split("\n")]
mol_df = pd.DataFrame(mol_buff_list,columns=["SMILES","Name"])
mol_df<jupyter_output><empty_output><jupyter_text>1. Molecular Geometry<jupyter_code>mol_df['mol'] = mol_df.SMILES.apply(Chem.MolFromSmiles)<jupyter_output><empty_output><jupyter_text>1.1 Generate a 3D structure<jupyter_code>mol_df.mol = mol_df.mol.apply(uru.gen_3d)<jupyter_output><empty_output><jupyter_text>1.2 View a 3D structure<jupyter_code>uru.MolTo3DView(mol_df.mol[0])<jupyter_output><empty_output><jupyter_text>1.3 Find the geometric center of a molecule<jupyter_code>mol_df[['x','y','z']] = mol_df.mol.apply(uru.get_center).to_list()
mol_df.head()<jupyter_output><empty_output><jupyter_text>1.4 Calculate and plot the first two shape moments for a set of moleculesRead a SMILES file with drugs from ChEMBL into a Pandas dataframe<jupyter_code>drug_df = pd.read_csv("https://raw.githubusercontent.com/PatWalters/datafiles/main/chembl_drugs.smi",
                     names=["SMILES","Name"],sep=" ")<jupyter_output><empty_output><jupyter_text>Add a molecule column to the dataframe<jupyter_code>drug_df['mol'] = drug_df.SMILES.progress_apply(Chem.MolFromSmiles)<jupyter_output><empty_output><jupyter_text>Strip salts, etc. by extracting the largest fragment from each molecule<jupyter_code>drug_df.mol = drug_df.mol.progress_apply(uru.get_largest_fragment)<jupyter_output><empty_output><jupyter_text>We want to limit this analysis to molecules with molecular weight < 500 so well calculate rule of 5 properties<jupyter_code>ro5_calc = uru.Ro5Calculator()
drug_df[ro5_calc.names] = drug_df.mol.progress_apply(ro5_calc.calc_mol).to_list()<jupyter_output><empty_output><jupyter_text>Remove molecules with MW > 500<jupyter_code>drug_df = drug_df.query("MolWt <= 500")<jupyter_output><empty_output><jupyter_text>Remove molecules with non-standard atom types<jupyter_code>weird_atom_smarts = "[!#1;!C;!O;!N;!S;!P;!F;!Cl;!Br;!I;!c;!o;!n;!s;!p;!Na;!K;!Mg;!Ca;!Li]"
weird_atom_pat = Chem.MolFromSmarts(weird_atom_smarts)
drug_df['odd'] = [x.HasSubstructMatch(weird_atom_pat) for x in drug_df.mol]
drug_df = drug_df.query("not odd")<jupyter_output><empty_output><jupyter_text>Add a 3D structure to the molecules<jupyter_code>drug_df.mol = drug_df.mol.progress_apply(uru.gen_3d)<jupyter_output><empty_output><jupyter_text>Remove molecules without a 3D structure<jupyter_code>drug_df.dropna(subset=['mol'],inplace=True)<jupyter_output><empty_output><jupyter_text>Define the pricipal moments<jupyter_code>drug_df[['NPR1','NPR2']] = drug_df.mol.apply(uru.get_shape_moments).to_list()<jupyter_output><empty_output><jupyter_text>Plot the principal moments<jupyter_code>sns.set_context('talk')
ax = sns.jointplot(x='NPR1', y='NPR2', data=drug_df, kind="hex",height=10)
ax.ax_marg_x.set_xlim(0, 1)
ax.ax_marg_y.set_ylim(0.45, 1.05);<jupyter_output><empty_output><jupyter_text>2. Clustering  2.1 Generating molecular fingerprints<jupyter_code>drug_df['fp'] = drug_df.mol.progress_apply(uru.mol2morgan_fp)<jupyter_output><empty_output><jupyter_text>2.2 Clustering using the Taylor-Butina algorithm<jupyter_code>drug_df['cluster'] = uru.taylor_butina_clustering(drug_df.fp.values)
drug_df.cluster.value_counts()
mols2grid.display(drug_df.query('cluster==0'),subset=["img"])<jupyter_output><empty_output><jupyter_text>3. Building a machine learning modelWe'll build two machine learning models, one using RDKit Morgan fingerprints as descriptors and another using Morgan fingerprints + descriptors.  After that we'll compare the models.   3.1 Building the modelsRead a dataset with SMILES, Name, and pIC50<jupyter_code>df = pd.read_csv("https://raw.githubusercontent.com/PatWalters/yamc/main/data/A2a.smi",
                 sep=" ",
                 names=["SMILES","Name","pIC50"])<jupyter_output><empty_output><jupyter_text>Add a molecule column to the dataframe<jupyter_code>df['mol'] = df.SMILES.apply(Chem.MolFromSmiles)<jupyter_output><empty_output><jupyter_text>Add a fingerprint column as a numpy array<jupyter_code>df['fp'] = df.mol.apply(uru.mol2numpy_fp)<jupyter_output><empty_output><jupyter_text>Calculate RDKit properties<jupyter_code>rdkit_props = uru.RDKitProperties()
df['desc'] = df.mol.apply(rdkit_props.calc_mol)
df['fp_plus_desc'] = [np.append(a,b) for a,b in df[['fp','desc']].values]<jupyter_output><empty_output><jupyter_text>Perform 10 cycles of cross validation and save the R^2 values for each cycle.<jupyter_code>r2_list = []
for i in tqdm(range(0,10)):
    train, test = train_test_split(df)
    train_X_fp = np.stack(train.fp)
    train_X_fp_plus_desc = np.stack(train.fp_plus_desc)
    train_y = train.pIC50.values
    test_X_fp = np.stack(test.fp)
    test_X_fp_plus_desc = np.stack(test.fp_plus_desc)
    test_y = test.pIC50.values
    model_fp = ExtraTreesRegressor()
    model_fp_plus_desc = ExtraTreesRegressor()
    model_fp.fit(train_X_fp, train_y)
    model_fp_plus_desc.fit(train_X_fp_plus_desc, train_y)
    pred_fp_r2 = r2_score(test_y,model_fp.predict(test_X_fp))
    pred_fp_plus_desc_r2 = r2_score(test_y,model_fp_plus_desc.predict(test_X_fp_plus_desc))
    r2_list.append([i,pred_fp_r2, pred_fp_plus_desc_r2])
r2_df = pd.DataFrame(r2_list,columns=["cycle","fp","fp_plus_desc"])<jupyter_output><empty_output><jupyter_text>3.2 Comparing the modelsTransform the data for plotting<jupyter_code>melt_df = r2_df.melt(id_vars="cycle")
melt_df.columns = ["cycle","method","R2"]
melt_df<jupyter_output><empty_output><jupyter_text>I to use the paired plots from the pigouin Python package as an inutitive means of comparing methods.<jupyter_code>ax = pingouin.plot_paired(data=melt_df,dv="R2",within="method",subject="cycle");
ax.set_xlabel("Method")
ax.set_ylabel("$R^2$");<jupyter_output><empty_output><jupyter_text>In order to have a more robust comparison we'll calculatee the effect size.  Since we have a small sample size we'll use the Hedges method in place of Cohen's. As reminder, here's a table summarizing the meaning of the effect size. * Small Effect Size: d=0.20* Medium Effect Size: d=0.50* Large Effect Size: d=0.80<jupyter_code>eff_size = pingouin.compute_effsize(r2_df.fp,r2_df.fp_plus_desc,eftype="cohen")
eff_size.round(2)<jupyter_output><empty_output><jupyter_text>We can also use the equation below to translate Cohen's d into a probability between 0 and 1. ![image.png](attachment:image.png)As we can see the effect size is small and probability of one method being better than the other is low.  Of course, in practice, we would do this comparison over many more datasets.<jupyter_code>(0.25 * abs(eff_size) + 0.5).round(2)<jupyter_output><empty_output><jupyter_text>4. Ring features and structure display<jupyter_code>spiro_mol = Chem.MolFromSmiles("C1CCC2(C1)Cc1ccccc1C2")
spiro_mol<jupyter_output><empty_output><jupyter_text>4.1 Labeling atomsLabel atoms by index<jupyter_code>uru.label_atoms(spiro_mol,[x.GetIdx() for x in spiro_mol.GetAtoms()])<jupyter_output><empty_output><jupyter_text>4.2 Identify spiro atoms<jupyter_code>spiro_atoms = uru.get_spiro_atoms(spiro_mol)
spiro_atoms<jupyter_output><empty_output><jupyter_text>4.3 Tagging spiro atoms<jupyter_code>uru.tag_atoms(spiro_mol,spiro_atoms)<jupyter_output><empty_output><jupyter_text>4.4 Identifying the largest ring<jupyter_code>uru.max_ring_size(spiro_mol)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1344737467331908, 'fraction_numerical': 0.026847232121644096, 'mean_word_length': 8.417225950782997, 'pattern_counts': {'":': 0, '<': 131, '<?xml version=': 0, '>': 130, 'https://': 2, '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': '32', 'dir': 'jupyter-structured-clean-dedup', 'id': '305675', 'license': 'MIT', 'n_tokens_mistral': 3617, 'n_tokens_neox': 3386, 'n_words': 523, 'path': 'notebooks/demo_useful_rdkit_utils.ipynb', 'repo_name': 'goraj/useful_rdkit_utils', 'size': '8418'}