{"cells":[{"metadata":{},"cell_type":"markdown","source":"# General information\n\nThe competition provide us several files with additional molecular features, such as dipole moments, magnetic shielding, mulliken charges and potential energy, which have information only for train and not for test datasets. \n\nSo, in order to use such attributes we would have to first create models from the training molecules and infer values for each test molecule.\n\n### Dealing With Molecules\n\nThe problem is that estimating properties relative to the position of the atoms in the molecule, it is useless to deal only with 1, 2 or 3 atoms, because they are often influenced by several (or all) neighboring atoms of a molecule.   \nThus, the ideal model must consider the entire molecular structure.   \nThe competition files provide the xyz positions of each atom, but using all the positions of all the atoms in a machine learning model can be a very difficult task.   \nTo overcome this problem, the concept of fingerprint of the molecule arises, where the molecular structure can be represented by an array, calculated from the spatial relationship (in three dimensions) between atoms in the molecule. Based on this, in theory you can create models for any molecular property based on the spatial structure.    \nThe most famous technique for defining molecular fingerprint is **Morgan**'s fingerprinting, for more information see the paper:    \n*Morgan, H. L. The Generation of a Unique Machine Description for Chemical Structures - A Technique Developed at Chemical Abstracts Service. J. Chem. Doc. 1965, 5: 107-112.*   \nThe process of creating the molecule's fingerprint is complex because it takes into account the entire 3D structure of the molecule, including the identity of the atoms, the type of chemical bonds, the degree of hybridization, conformation, relative angles, and so on.\nFor this task we will use the incredible DEEPCHEM library, which abstracts most of the complexity and requires only a few lines of code.\n\n### Main goal\nThe main purpose of this kernel is to teach how to use the Stanford DeepChem library to use the additional molecular features in this competition. No knowledge of chemistry is needed, just pure python !!\nIn this kernel we will use the potential energy data as an example, however the same concept can be applied to the other files provided.\n\n> **If you want to see only the DEEPCHEM library application, you can skip the introduction and preprocessing and go straight to the title: Deepchem to provide the potential energy, at the end of this kernel.**"},{"metadata":{},"cell_type":"markdown","source":"### Instaling packages, libraries, and so on..."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"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","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"!pip install deepchem","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"%%bash -e\nif ! [[ -f ./xyz2mol.py ]]; then\n  wget https://raw.githubusercontent.com/jensengroup/xyz2mol/master/xyz2mol.py\nfi","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":false},"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);","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"file_folder = '../input/champs-scalar-coupling' if 'champs-scalar-coupling' in os.listdir('../input/') else '../input'\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.read_csv(f'{file_folder}/train.csv')[0-60000:]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"import sys\n!conda install --yes --prefix {sys.prefix} -c rdkit rdkit","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Dealing with different file formats\nTo be able to use the deepchem library satisfactorily we need to have our molecules encoded in SMILES format. The files of the molecules supplied by Kaggle are in xyz formats.   \nSo we will modify them. This transformation will be in 3 steps:\n1. Convert XYZ to MOL format.\n2. Transform MOL to SMILES format.\n3. Create a csv file.\n\n## Read molecules in RDKit using XYZ2MOL package\n### XYZ file format\nThe XYZ file format is a chemical file format. Normally XYZ format specifies the molecule geometry by giving the number of atoms with Cartesian coordinates that will be read on the first line, a comment on the second, and the lines of atomic coordinates in the following lines.[1] The file format is used in computational chemistry programs for importing and exporting geometries. The units are generally in ångströms.(https://en.wikipedia.org/wiki/XYZ_file_format)    \nNormally, we find in a XYZ file:  \n-First line: total number of atoms (optional)  \n-Second line: molecule name or comment (optional)  \n-All other lines: element symbol or atomic number, x, y, and z coordinates, separated by spaces, tabs, or commas  \nExample:    \nglucose from 2gbp     \nC  35.884  30.895  49.120   \nC  36.177  29.853  50.124    \nC  37.296  30.296  51.074    \nC  38.553  30.400  50.259     \nC  38.357  31.290  49.044     \nC  39.559  31.209  48.082     \nO  34.968  30.340  48.234  \nO  34.923  29.775  50.910  \nO  37.441  29.265  52.113  \nO  39.572  30.954  51.086  \nO  37.155  30.858  48.364  \nO  39.261  32.018  46.920\n\n### Mol file format\nAn MDL Molfile is a file format for holding information about the atoms, bonds, connectivity and coordinates of a molecule.(https://en.wikipedia.org/wiki/Chemical_table_file)  \n\n### xyz2mol\nThis is a library that converts an xyz file to an RDKit mol object. The code is based on this paper Yeonjoon Kim and Woo Youn Kim \"Universal Structure Conversion Method for Organic Molecules: From Atomic Connectivity to Three-Dimensional Geometry\" Bull. Korean Chem. Soc. 2015, Vol. 36, 1769-1777 DOI: 10.1002/bkcs.10334. (https://github.com/jensengroup/xyz2mol)  \n\n### Smiles Format\nThe simplified molecular-input line-entry system (SMILES) is a specification in the form of a line notation for describing the structure of chemical species using short ASCII strings. SMILES strings can be imported by most molecule editors for conversion back into two-dimensional drawings or three-dimensional models of the molecules.(https://en.wikipedia.org/wiki/Simplified_molecular-input_line-entry_system)    \nIn terms of a graph-based computational procedure, SMILES is a string obtained by printing the symbol nodes encountered in a depth-first tree traversal of a chemical graph. The chemical graph is first trimmed to remove hydrogen atoms and cycles are broken to turn it into a spanning tree. Where cycles have been broken, numeric suffix labels are included to indicate the connected nodes. Parentheses are used to indicate points of branching on the tree.(https://en.wikipedia.org/wiki/Simplified_molecular-input_line-entry_system, http://opensmiles.org)    \nA wide variety of SMILES strings are acceptable as input. For example, all of the following represent ethanol:   \n\nCCO\t\nOCC\t\nC(O)C\t\n[CH3][CH2][OH]\t\n[H][C]([H])([H])C([H])([H])[O][H]\t\n\n### RDKit\nRDKit is a collection of cheminformatics and machine-learning software written in C++ and Python, created by Greg Landrum.(https://github.com/rdkit/rdkit)     \nIt is an amazing library for chemistry, computational biology and bioinformatics studies involving python programming.\nWe use it to convert mol format to Smiles format."},{"metadata":{"trusted":true},"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\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()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"## Testing formatting xyz files into smiles format"},{"metadata":{},"cell_type":"markdown","source":"Let´s test the conversion to smiles format."},{"metadata":{"trusted":true},"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","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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 is wrong!!')\n    pass ","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Visualize 3d molecules\n\nHow can we guarantee that the format conversion is correct?\nLet's look at some sample in 3D for verification."},{"metadata":{"trusted":true},"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!!'''","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Using DeepChem Library"},{"metadata":{},"cell_type":"markdown","source":"DeepChem aims to provide a high quality open-source toolchain that democratizes the use of deep-learning in drug discovery, materials science, quantum chemistry, and biology.\nThe library was originally created by Bharath Ramsundar with encouragement and guidance from Vijay Pande.\n\nIt started as a Pande group project at Stanford, and is now developed by many academic and industrial collaborators. \n\nI really recommend you take a look at these sites to get deeper into the tool.    \nhttps://github.com/deepchem/deepchem      \nhttps://deepchem.io    \nhttps://deepchem.io/docs/notebooks/index.html    "},{"metadata":{},"cell_type":"markdown","source":"## Prepare dataframe for using with DeepChem library"},{"metadata":{},"cell_type":"markdown","source":"In this example we will use the potential_energy file.\nThe idea is create a model using the potential energy from the train dataset in order to predict it into the test dataframe. "},{"metadata":{},"cell_type":"markdown","source":"### Potential Energy\n\nPotential energy is a concept that can be used to understand any type of change in a system. Any system (a collection of molecules, atoms, or electrons) has a potential energy associated with it. The origin of this energy and what it represents is not important, all that you need to know we can use the potential energy of a molecule to understand how and why it will undergo a change. These changes could include the formation or breaking of a bond, the gain or loss of an electron, or the change in physical orientation of the atoms in the molecule. Each physical change has an associated change in potential energy.(https://sop4cv.com/chapters/PotentialEnergyThermodynamicsAndKinetics.html)    \n\n\npotential_energy.csv - contains the potential energy of the molecules. The first column (molecule_name) contains the name of the molecule, the second column (potential_energy) contains the potential energy of the molecule.(https://www.kaggle.com/c/champs-scalar-coupling/data)   "},{"metadata":{"trusted":true},"cell_type":"code","source":"potential_energyDF = pd.read_csv(f'{file_folder}/potential_energy.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"potential_energyDF.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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()\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"#### Create smiles format from molecules"},{"metadata":{"trusted":true},"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]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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         ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df_temp.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"potential_energyDF_Smiles.shape","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Deepchem to provide the potential energy"},{"metadata":{"trusted":true},"cell_type":"code","source":"from rdkit import Chem\nimport random\nfrom deepchem.feat import CircularFingerprint\nimport deepchem as dc\nimport numpy as np","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Dataset file\ndataset_file = '../potential_energyDF_Smiles.csv'","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"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)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"After the aforementioned featurization, we will create our potential energy model. We will use random forest."},{"metadata":{"trusted":true},"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)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"R2 Score obtained was not bad for a first model and the small dataset used, of course you can test with other models and parameters to improve it."},{"metadata":{"trusted":true},"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()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now you can use your model to predict the potential energy of your test dataframe.    "},{"metadata":{},"cell_type":"markdown","source":"As we said earlier, the same concept can be applied to the other files provided, such as **dipole moments, mulliken charges**, and so on.\n\n**So, if you find it useful and liked this code what do you think of giving me an upvote? Thanks a lot.    \nIf you have any suggestions or questions please contact me.   \nKeep in touch because I want to update it with more tips.   \nThank you very much.**"}],"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.6.4","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat":4,"nbformat_minor":1}