{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## Predicting Molecular Properties based on the Stanford DeepChem Library","metadata":{}},{"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\npd.set_option('max_rows', 5)\nimport os\nimport matplotlib.pylab as plt\nimport seaborn as sns\nfrom sklearn import metrics\nimport warnings\nwarnings.simplefilter(action='ignore', category=FutureWarning)\nplt.style.use('ggplot')\ncolor_pal = [x['color'] for x in plt.rcParams['axes.prop_cycle']]\nimport os\nfrom IPython.display import Image, HTML, display","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install deepchem","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%bash -e\nif ! [[ -f ./xyz2mol.py ]]; then\n  wget https://raw.githubusercontent.com/jensengroup/xyz2mol/master/xyz2mol.py\nfi","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pip install py3Dmol # This is built on the object-oriented, webGL based JavaScript library for online molecular visualization 3Dmol.js (Rego & Koes, 2015);","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"file_folder = '../input/champs-scalar-coupling' if 'champs-scalar-coupling' in os.listdir('../input/') else '../input'","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train = pd.read_csv(f'{file_folder}/train.csv')[0-60000:]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train.head()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import sys\n!conda install --yes --prefix {sys.prefix} -c rdkit rdkit","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Few Snippets from https://www.kaggle.com/sunhwan/using-rdkit-for-atomic-feature-and-visualization\n# rdkit & xyz2mol\nfrom rdkit import Chem\nfrom rdkit.Chem import AllChem\nfrom rdkit.Chem.Draw import IPythonConsole #Needed to show molecules\nfrom rdkit.Chem import Draw\nfrom rdkit.Chem.Draw.MolDrawing import MolDrawing, DrawingOptions #Only needed if modifying defaults\nDrawingOptions.bondLineWidth=1.8\nfrom rdkit.Chem.rdmolops import SanitizeFlags\n\n# https://github.com/jensengroup/xyz2mol\nfrom xyz2mol import xyz2mol, xyz2AC, AC2mol, read_xyz_file\nfrom pathlib import Path\nimport pickle\n\nCACHEDIR = Path('./')\n\ndef chiral_stereo_check(mol):\n    # avoid sanitization error e.g., dsgdb9nsd_037900.xyz\n    Chem.SanitizeMol(mol, SanitizeFlags.SANITIZE_ALL - SanitizeFlags.SANITIZE_PROPERTIES)\n    Chem.DetectBondStereochemistry(mol,-1)\n    # ignore stereochemistry for now\n    #Chem.AssignStereochemistry(mol, flagPossibleStereoCenters=True, force=True)\n    #Chem.AssignAtomChiralTagsFromStructure(mol,-1)\n    return mol\n\ndef xyz2mol(atomicNumList,charge,xyz_coordinates,charged_fragments,quick):\n    AC,mol = xyz2AC(atomicNumList,xyz_coordinates)\n    new_mol = AC2mol(mol,AC,atomicNumList,charge,charged_fragments,quick)\n    new_mol = chiral_stereo_check(new_mol)\n    return new_mol\n\ndef MolFromXYZ(filename):\n    mol=''\n    charged_fragments = True\n    quick = True\n    cache_filename = CACHEDIR/f'{filename.stem}.pkl'\n    if cache_filename.exists():\n        return pickle.load(open(cache_filename, 'rb'))\n    else:\n        try:\n            atomicNumList, charge, xyz_coordinates = read_xyz_file(filename)\n            mol = xyz2mol(atomicNumList, charge, xyz_coordinates, charged_fragments, quick)\n            # commenting this out for kernel to work.\n            # for some reason kernel runs okay interactively, but fails when it is committed.\n            pickle.dump(mol, open(cache_filename, 'wb'))\n        except:\n            print(filename)\n    return mol\n#mol = MolFromXYZ(xyzfiles[1])\n#m = Chem.MolFromSmiles(Chem.MolToSmiles(mol, allHsExplicit=True)); m\n\nfrom multiprocessing import Pool\nfrom tqdm import *\nfrom glob import glob\n\ndef MolFromXYZ_(filename):\n    return filename.stem, MolFromXYZ(filename)\n\nmols = {}\nn_cpu = 4\nwith Pool(n_cpu) as p:\n    molecule_names = np.concatenate([train.molecule_name.unique()])\n    molecule_names = molecule_names[:400]\n    xyzfiles = [Path(f'../input/structures/{f}.xyz') for f in molecule_names]\n    n = len(xyzfiles)\n    with tqdm(total=n) as pbar:\n        for res in p.imap_unordered(MolFromXYZ_, xyzfiles):\n            mols[res[0]] = res[1]\n            pbar.update()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"molecule_names = np.concatenate([train.molecule_name.unique()])\nxyzfiles = [(f'{f}.xyz') for f in molecule_names]\nfew_molecule_names = molecule_names[:5]\nfew_molecule_names","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create mols format from xyz files\n# Create smiles format\n\ntry:\n    for molecule_name in few_molecule_names:\n        #print('Molecule: {}'.format(molecule_name))\n        m = MolFromXYZ(Path(f'../input/champs-scalar-coupling/structures/{molecule_name}.xyz'))\n        smile_fromMolecule = Chem.MolToSmiles(m)\n        print('Smile format from molecule: {}'.format(smile_fromMolecule))\n\nexcept:\n    print ('...something be goofy!!')\n    pass ","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Inspired by https://github.com/greglandrum/rdkit_blog/tree/master/notebooks\n\nimport py3Dmol # Amazing library for 3D visualization\nfrom rdkit import Chem\nfrom rdkit.Chem import AllChem\nfrom ipywidgets import interact, interactive, fixed\ndef drawit(m,p,confId=-1):\n    mb = Chem.MolToMolBlock(m,confId=confId)\n    p.removeAllModels()\n    p.addModel(mb,'sdf')\n    p.setStyle({'stick':{}})\n    p.setBackgroundColor('0xeeeeee')\n    p.zoomTo()\n    return p.show()\np = py3Dmol.view(width=400,height=400)\n\n# now construct the 3d view:\nm = MolFromXYZ(Path(f'../input/champs-scalar-coupling/structures/dsgdb9nsd_129764.xyz'))\nsmile_fromMolecule = Chem.MolToSmiles(m)\nm = Chem.MolFromSmiles(smile_fromMolecule)\n#print(smile_fromMolecule)\nm = Chem.AddHs(m)\nAllChem.EmbedMultipleConfs(m,randomSeed=0xf00d,useExpTorsionAnglePrefs=True, useBasicKnowledge=True)\ninteract(drawit, m=fixed(m),p=fixed(p));\n\n''' Trick: If you want click and drag on the image to spin the molecule, and scroll for zoom!!'''","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"potential_energyDF = pd.read_csv(f'{file_folder}/potential_energy.csv')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"potential_energyDF.head()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plot the distribution of potential_energy\npotential_energyDF['potential_energy'].plot(kind='hist', figsize=(25, 5), bins=500, title='Distribution of Molecular Potential Energy', color='g')\nplt.show()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import tensorflow as tf\nimport deepchem as dc\nimport numpy as np\nimport pandas as pd\n\n# Getting only few samples\nmolecule_names = molecule_names[0:400]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create smiles format\nxyzfiles = [(f'{f}.xyz') for f in molecule_names]\nn = len(xyzfiles)\ndf_temp = pd.DataFrame({'molecule_name':[],'smile_fromMolecule': []})\ntry:\n    for molecule_name in molecule_names:\n        #print('Molecule: {}'.format(molecule_name))\n        m = MolFromXYZ(Path(f'../input/champs-scalar-coupling/structures/{molecule_name}.xyz'))\n        smile_fromMolecule = Chem.MolToSmiles(m)\n        #print('Smile format from molecule: {}'.format(smile_fromMolecule))\n        df_temp = df_temp.append({'molecule_name': molecule_name,'smile_fromMolecule': smile_fromMolecule}, ignore_index=True)\nexcept:\n    print ('...something is wrong!!')\n    pass ","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_temp.head()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Merge two dataframes\npotential_energyDF_Smiles = pd.merge(df_temp, potential_energyDF, how = 'left', on='molecule_name')\npotential_energyDF_Smiles.dropna(inplace=True)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"potential_energyDF_Smiles.shape","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Our dataset file must contain a column with the SMILES sentence and another with our target (potential energy) in order to use with DeepChem\npotential_energyDF_Smiles[['potential_energy', 'smile_fromMolecule']].to_csv('../potential_energyDF_Smiles.csv', index=False)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from rdkit import Chem\nimport random\nfrom deepchem.feat import CircularFingerprint\nimport deepchem as dc\nimport numpy as np","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Dataset file\ndataset_file = '../potential_energyDF_Smiles.csv'","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Featurizer will create our fingerprint, and turn it into an array with 1024 bits.\nfeaturizer = dc.feat.CircularFingerprint(size=1024)\n\n# Prepare our dataset file\nloader = dc.data.CSVLoader(\n      tasks=[\"potential_energy\"], smiles_field=\"smile_fromMolecule\",\n      featurizer=featurizer)\n\ndataset = loader.featurize(dataset_file)\n\n# Split train, validation and test\nsplitter = dc.splits.ScaffoldSplitter(dataset_file)\ntrain_dataset, valid_dataset, test_dataset = splitter.train_valid_test_split(dataset)\ntrain_mols = [Chem.MolFromSmiles(compound)\n              for compound in train_dataset.ids]\nvalid_mols = [Chem.MolFromSmiles(compound)\n              for compound in valid_dataset.ids]\n# Normalize them\ntransformers = [dc.trans.NormalizationTransformer(transform_y=True, dataset=train_dataset)]\n\nfor dataset in [train_dataset, valid_dataset, test_dataset]:\n    for transformer in transformers:\n        dataset = transformer.transform(dataset)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.ensemble import RandomForestRegressor\n\n# Our regressor\nsklearn_model = RandomForestRegressor(n_estimators=100, criterion='mse', max_depth=20, \n                                      min_samples_split=2, min_samples_leaf=1, \n                                      min_weight_fraction_leaf=0.0, max_features=None,\n                                      max_leaf_nodes=None, min_impurity_decrease=0.0, \n                                      min_impurity_split=None, bootstrap=True, oob_score=False, \n                                      n_jobs=None, random_state=1, warm_start=False)\nmodel = dc.models.SklearnModel(sklearn_model)\nmodel.fit(train_dataset)\n\n# Evaluation\nfrom deepchem.utils.evaluate import Evaluator\nmetric = dc.metrics.Metric(dc.metrics.r2_score)\nevaluator = Evaluator(model, valid_dataset, transformers)\nr2score = evaluator.compute_model_performance([metric])\nprint(r2score)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# let’s plot the predicted R2 scores versus the true R2 scores for the constructed model.\nllim = -600\nulim = -300\ntask = \"measured potential_energy\"\npredicted_test = model.predict(test_dataset)\ntrue_test = test_dataset.y\nplt.scatter(predicted_test, true_test)\nplt.xlim((llim, ulim))\nplt.ylim((llim, ulim))\nplt.plot([llim, ulim], [llim, ulim])\nplt.xlabel('Predicted potential_energy')\nplt.ylabel('True potential_energy')\nplt.title(r'RF- predicted vs. true potential_energy')\nplt.show()","metadata":{},"execution_count":null,"outputs":[]}]}