{"cells":[{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","collapsed":true,"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":false},"cell_type":"markdown","source":"# Data Exploration\nThe goal of the competition is to predict NMR coupling constants from molecule structures and other information.\n\nBriefly, coupling constants represent the strength of the interaction between different atomic nuclei in a molecule. The interactions occur through bonds, rather than directly through space. Therefore the number of bonds between the two atoms of interest will strongly affect their corresponding coupling constant, while the spatial separation between the two atoms may have only a weak correlation with coupling constant. Coupling constants through more than three bonds can not normally be detected, except in certain special cases. The coupling constant is also strongly dependent on the chemical functional group to which the atoms belong. For example, the three-bond coupling between two alkene protons (hydrogen atoms) may differ significantly from that of two alkyl protons. The geometry of the bonds has a great impact on the coupling; the strength of proton-proton three bond couplings are strongly dependent on the dihedral angle between the two proton-carbon/heteroatom bonds as described by the [Karplus equation](https://en.wikipedia.org/wiki/Karplus_equation).\n\nAdditional literature:\n- [Vicinal Coupling Constants and Conformation of Biomolecules, Altona](https://doi.org/10.1002/9780470034590.emrstm0587)\n- [Scalar Coupling Constants—Their Analysis and Their Application for the Elucidation of Structures, Eberstadt _et al._](https://doi.org/10.1002/anie.199516711)\n- <a href=\"https://doi.org/10.1016/0040-4020(80)80155-4\">The relationship between proton-proton NMR coupling constants and substituent electronegativities, Haasnoot _et al._</a>\n\nIt may be difficult to identify the functional group involved from the raw structural data. However, the hybridisation (essentially the number of neighbouring atoms) and the identities of the other atoms in the bond my serve as proxies. The interatomic bond length can also be strongly influenced by the functional group that it belongs to:\n\n- [Typical interatomic distances: organic compounds, Allen _et al._](https://doi.org/10.1107/97809553602060000621)\n\nBefore we begin to explore our data, therefore, some features which we may be interested in engineering are:\n- Bond lengths\n- Functional groups (perhaps bond lengths may serve as a proxy?)\n- Number of neighbouring atoms (also may serve as a proxy for functional group)\n- Angles between bonds (in the case of two-bond couplings)\n- Dihedral angles (for three-bond couplings)\n- Mean coupling constants (possibly for coupling types which don't vary much)\n- Identities of the other atoms involved in two- or three-bond couplings\n- Electronegativities of the atoms involved\n\nLet's first import the necessary libraries, our training, and test data, and look at the information provided."},{"metadata":{"trusted":true},"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nimport seaborn as sns\nimport matplotlib.pyplot as plt\n\ntrain = pd.read_csv('../input/train.csv')\ntest = pd.read_csv('../input/train.csv')\ntrain.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We have a molecule name, allowing us to match the training data entry to a structure, and the indices of the two atoms that are coupled. Note that multiple coupling constants correspond to the same structure. Let's write a function to import the structural data and see what form it takes:"},{"metadata":{"trusted":true},"cell_type":"code","source":"def read_xyz(path, filename):\n    return pd.read_csv(path+filename, skiprows = 2, header = None, sep = ' ', usecols=[0, 1,2,3], names=['atom', 'x', 'y', 'z'])\n\npath = '../input/structures/'\nfilename = 'dsgdb9nsd_000001.xyz'\n\nread_xyz(path, filename)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"We have what appears to be a molecule of methane, with the positions of each of the atoms represented with cartesian coordinates. This molecule contains only carbon and hydrogen, but how many different atom types are present in our training data?"},{"metadata":{"trusted":true},"cell_type":"code","source":"# This is the code, but it is quite time consuming to run so I'll just provide the answer below\n\"\"\"\natom_list = []\nfor filename in os.listdir(\"../input/structures\"):\n    atom_list = atom_list + list(read_xyz(path, filename)['atom'])\natom_list = set(atom_list)\nprint(atom_list)\n\"\"\"\nprint(\"{'O', 'H', 'C', 'F', 'N'}\")\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"All the molecules are made up only of oxygen, hydrogen, carbon, fluorine, and nitrogen. "},{"metadata":{"trusted":true},"cell_type":"code","source":"x_list = []\ny_list = []\nz_list = []\nfor filename in os.listdir(\"../input/structures\"):\n    x_list = x_list + list(read_xyz(path, filename)['x'])\n    y_list = y_list + list(read_xyz(path, filename)['y'])\n    z_list = z_list + list(read_xyz(path, filename)['z'])\ndimfig, dimaxes = plt.subplots(3, 1, figsize = (6, 6))\nsns.distplot(x_list, ax=dimaxes[0])\nsns.distplot(y_list, ax=dimaxes[1])\nsns.distplot(z_list, ax=dimaxes[2])\nprint(\"x max: \" + str(np.max(x_list)) + \" x min : \" + str(np.min(x_list)))\nprint(\"y max: \" + str(np.max(y_list)) + \" y min : \" + str(np.min(y_list)))\nprint(\"z max: \" + str(np.max(z_list)) + \" z min : \" + str(np.min(z_list)))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The coordinates almost all appear to lie within 10 and -10 Å, and the great majority lie within 5 Å of the origin. Let's next look at how many distinct coupling types there are in the data."},{"metadata":{"trusted":true},"cell_type":"code","source":"coupling_types = set(train['type'])\nprint(coupling_types)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"So we have eight types of coupling, all of which are between a proton and another proton or heteroatom. We also have one, two and three bond couplings. Let's look at these coupling types and the distribution of their coupling constants in turn."},{"metadata":{"trusted":true},"cell_type":"code","source":"coupling_types = list(coupling_types)\ntotals = [np.sum(train['type'] == x) for x in coupling_types]\n\nsubsets = dict()\nfor x in coupling_types:\n    subsets[x] = train.loc[train['type'] == x]\n\nbar_fig, bar_axis = plt.subplots()\n\nsns.barplot(coupling_types, totals, ax = bar_axis)\n\ndist_fig, dist_axes = plt.subplots(len(subsets), 1, figsize = (6, 12))\n\nfor (x, y) in zip(dist_axes, coupling_types):\n    sns.distplot(subsets[y]['scalar_coupling_constant'], ax=x)\n    x.set_title(y)\n\ndist_fig.tight_layout()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The distributions of the coupling constants look quite complex \n\n# Feature Engineering\n\nWe will write some functions to allow us to engineer some features and see what relationship they have with the coupling constants. The algorithm for determining the dihedral angle was based on [this Stack Exchange thread.](https://math.stackexchange.com/questions/47059/how-do-i-calculate-a-dihedral-angle-given-cartesian-coordinates)"},{"metadata":{"trusted":true},"cell_type":"code","source":"def length(data, index1, index2):\n    \"\"\"Takes an xyz file imported by read_xyz and calculates the distance between two points\"\"\"\n    return np.sqrt(np.sum(np.square(data[['x', 'y', 'z']].loc[index1]-data[['x', 'y', 'z']].loc[index2])))\n\ndef neighbours(data, index):\n    \"\"\"Takes an xyz file imported by read_xyz and calculates the number of neighbours within sqrt(3) Å of the indexed atom\"\"\"\n    l2 = np.array([np.sum(np.square(data[['x', 'y', 'z']].loc[index]-data[['x', 'y', 'z']].loc[x])) for x in range(len(data))])\n    return np.sum(l2 < 3) - 1\n\ndef nearest(data, index):\n    \"\"\"Takes an xyz file imported by read_xyz and finds the index of the nearest atom\"\"\"\n    #data['index'] = data.index\n    point = data.loc[index][['x', 'y', 'z']]\n    data = data[data['atom'] != 'H'][['x', 'y', 'z']]\n    data[['x', 'y', 'z']] = data[['x', 'y', 'z']] - point\n    data[['x', 'y', 'z']] = np.square(data[['x', 'y', 'z']])\n    data = np.sum(data, axis = 1)\n    if index in data.index: data[index] = 999\n    return np.argmin(data)\n\ndef magnitude(vector):\n    \"\"\"Calculates the magnitude of a vector\"\"\"\n    return np.sqrt(np.sum(np.square(vector)))\n    \ndef dihedral(point1, point2, point3, point4):\n    \"\"\"Calculates the dihederal angle between two bonds\"\"\"\n    b1 = point1-point2\n    b2 = point2-point3\n    b3 = point3-point4\n    n1 = np.cross(b1, b2)\n    n1 = n1/magnitude(n1)\n    n2 = np.cross(b2, b3)\n    n2 = n2/magnitude(n2)\n    m1 = np.cross(n1, b2/magnitude(b2))\n    x = np.dot(n1, n2)\n    y = np.dot(m1, n2)\n    return np.arctan2(x, y)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The next step is to generate the features appropriate to each coupling type. We will do this for one thousand training examples for each coupling type for the purposes of visualisation:"},{"metadata":{"trusted":true},"cell_type":"code","source":"def single_bond(coupling_type):    \n    feature_list = []\n    \n    for x in range(1000):#len(subsets[coupling_type])):\n        current = subsets[coupling_type].iloc[x]\n        index0 = current['atom_index_0']\n        index1 = current['atom_index_1']\n        filename = current['molecule_name'] + '.xyz'\n        data = read_xyz(path, filename)\n        feature_list.append((length(data, index0, index1), neighbours(data, index1), current['scalar_coupling_constant']))\n    \n    return pd.DataFrame(feature_list, columns = ['length', 'hybrid', 'coupling'])\n\ndef two_bond(coupling_type):\n    feature_list = []\n    for x in range(1000):\n        current = subsets[coupling_type].iloc[x]\n        data = read_xyz(path, current['molecule_name'] + '.xyz')\n        index_0 = current['atom_index_0']\n        index_1 = current['atom_index_1']\n        shared = nearest(data, index_0)\n        length1 = length(data, index_0, shared)\n        length2 = length(data, index_1, shared)\n        vector1 = data[['x', 'y', 'z']].loc[index_0]-data[['x', 'y', 'z']].loc[shared]\n        vector2 = data[['x', 'y', 'z']].loc[index_1]-data[['x', 'y', 'z']].loc[shared]\n        cosine = np.dot(vector1, vector2)/(length1 * length2)\n        shared_hybrid = neighbours(data, shared)\n        carbon_hybrid = neighbours(data, index_1)\n        feature_list.append((length1, length2, cosine, data['atom'].iloc[shared], shared_hybrid, carbon_hybrid, current['scalar_coupling_constant']))\n    return pd.DataFrame(feature_list, columns = ['length1', 'length2', 'cosine', 'atom', 'hybrid1', 'hybrid2', 'coupling'])\n\ndef three_bond(coupling_type):\n    feature_list = []\n    for x in range(1000):\n        current = subsets[coupling_type].iloc[x]\n        data = read_xyz(path, current['molecule_name'] + '.xyz')\n        index_0 = current['atom_index_0']\n        index_1 = current['atom_index_1']\n        shared1 = nearest(data, index_0)\n        shared2 = nearest(data, index_1)\n        length1 = length(data, index_0, shared1)\n        length2 = length(data, index_1, shared2)\n        length_shared = length(data, index_0, index_1)\n        cosine = dihedral(data[['x', 'y', 'z']].loc[index_0], data[['x', 'y', 'z']].loc[shared1], data[['x', 'y', 'z']].loc[shared2], data[['x', 'y', 'z']].loc[index_1])\n        shared1_hybrid = neighbours(data, shared1)\n        shared2_hybrid = neighbours(data, shared2)\n        terminal_hybrid = neighbours(data, index_1)\n        feature_list.append((length1, length2, length_shared, cosine, data['atom'].iloc[shared1], data['atom'].iloc[shared2], shared1_hybrid, shared2_hybrid, terminal_hybrid, current['scalar_coupling_constant']))\n    return pd.DataFrame(feature_list, columns = ['length1', 'length2', 'length_shared', 'angle', 'atom1', 'atom2', 'hybrid1', 'hybrid2', 'terminal_hybrid', 'coupling'])\n\nfunction_dict = {'1': single_bond, '2': two_bond, '3': three_bond}\nengineered = {x:function_dict[x[0]](x) for x in coupling_types}","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Now let's look at the relationship between some of the engineered features and the coupling constant. Firstly, the dihedral angle anf the 3JHH couplings: "},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.scatterplot(engineered['3JHH']['angle'], engineered['3JHH']['coupling'], hue=engineered['3JHH']['length_shared'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There dihedral angle and the coupling constant appear to be related by a periodic function, as one might expect (Karplus equation). Interestingly, the dihedral angle also correlates with the length of the central bond, possibly due to the relationship between the relationship between hybridisation and bond length. Let's look at the relationships of all the engineered features with a pairplot, paying special attention to the last row - the plots of the engineered features against the couping constant."},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.pairplot(engineered['3JHH'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Does the same hold true for the other three-bond couplings? Next we will analyse the three-bond proton-carbon coupling"},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.pairplot(engineered['3JHC'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The relationship between dihedral angle and coupling constant is not so clear in this case. The relationship between the other engineered features and the coupling constant is also weak. The last of the three-bond couplings is proton to nitrogen:"},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.pairplot(engineered['3JHN'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here we see essentially similar patterns as in the proton-carbon couplings. Next come the two-bond couplings. Let's start with the proton-proton couplings."},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.pairplot(engineered['2JHH'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here we can see a clear correlation between the cosine of the angle between the two bonds and the coupling constant. Unsurprisingly (as both protons are bonded to the same atom) the lengths of both bonds are strongly correlated. Next, the two bond couplings between proton and carbon."},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.pairplot(engineered['2JHC'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Here we may be able to use the bond lengths, the hybridisation, and the cosine of the angle between bonds to predict coupling constants. The last of the two bond coupling constants in proton-nitrogen:"},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.pairplot(engineered['2JHN'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The correlations here appear rather weak. Letäs now look at the single bond couplings. These should be significantly simpler to analyse as we have calculated fewer engineered features. The only (non-exotic) molecule in which single-bond proton-proton interactions is possible is molecular hydrogen, and so we only have two-bond proton-nitrogen and proton-nitrogen couplings. Let's start with proton-carbon."},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.pairplot(engineered['1JHC'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Both the hybridisation of the carbon atom and the length of the proton-carbon bond appear to correlate with the coupling constant - these should be good features to use to make predictions. Finally, the single-bond proton-nitrogen couplings:"},{"metadata":{"trusted":true},"cell_type":"code","source":"sns.pairplot(engineered['1JHN'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"A very nice relationship between the length of the proton-nitrogen bond can be observed here. The hybridisation of the nitrogen also offers some predictive power - sp$_2$ hybridised nitrogen appears to couple to protons in a very narrow range around 35 Hz.\n\n# Conclusions\n\nThe number of features analysed here has been quite limited; information such as the identity of the functional groups which the coupled atoms belong can be expected to yield more information. However, several of the features engineered here have a significant correlation with the coupling constant and yield significant improvements in predictive ability if incorporated into machine learning models."}],"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}