{"cells":[{"metadata":{},"cell_type":"markdown","source":"# How to find cis–trans isomerism\n\nCis indicates that the functional groups are on the same side of the carbon chain while trans conveys that functional groups are on opposing sides of the carbon chain."},{"metadata":{},"cell_type":"markdown","source":""},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true},"cell_type":"code","source":"!conda install -c openbabel openbabel -y","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"import openbabel\nimport numpy as np\nimport pandas as pd\nimport os\nfrom tqdm import tqdm\ntqdm.pandas()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Load data"},{"metadata":{"trusted":true},"cell_type":"code","source":"print(os.listdir('../input'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"base_path = '../input/champs-scalar-coupling'\n\nstructures = pd.read_csv(f'{base_path}/structures.csv')\n\n# data reduction\ntrain = pd.read_csv(f'{base_path}/train.csv')[::10]\ntest = pd.read_csv(f'{base_path}/test.csv')[::10]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"obConversion = openbabel.OBConversion()\nobConversion.SetInFormat(\"xyz\")\nxyz_path = f'{base_path}/structures/'","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Find cis-trans"},{"metadata":{"trusted":true},"cell_type":"code","source":"def cis_trans_bond_indices(molecule_name):\n mol = openbabel.OBMol()\n obConversion.ReadFile(mol, f'{xyz_path}/{molecule_name}.xyz')\n obs = openbabel.OBStereoFacade(mol)\n has_ct = [obs.HasCisTransStereo(n) for n in range(mol.NumBonds())]\n return [i for i, x in enumerate(has_ct) if x == True] if has_ct else []","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df = pd.DataFrame(structures.molecule_name.unique(), columns=['molecule_name'])\ndf.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df['bond_indices'] = df.molecule_name.progress_apply(lambda x: cis_trans_bond_indices(x))\ndf['len_bond_indices'] = df.bond_indices.progress_apply(lambda x:len(x))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df.len_bond_indices.unique()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"df[df['len_bond_indices']!=0].head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.merge(train, df, how='left', on='molecule_name')\ntest = pd.merge(test, df, how='left', on='molecule_name')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def is_cis_trans(molecule_name, bond_indices, atom_index_0, atom_index_1):\n if len(bond_indices) == 0:\n return pd.Series([0,0])\n\n mol = openbabel.OBMol()\n obConversion.ReadFile(mol, f'{xyz_path}/{molecule_name}.xyz')\n obs = openbabel.OBStereoFacade(mol)\n \n is_cis = [obs.GetCisTransStereo(i).IsCis(atom_index_0, atom_index_1) for i in bond_indices]\n is_trans = [obs.GetCisTransStereo(i).IsTrans(atom_index_0, atom_index_1) for i in bond_indices]\n return pd.Series([int(True in is_cis), int(True in is_trans)])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train[['is_cis','is_trans']] = train.progress_apply(lambda x: is_cis_trans(x.molecule_name,\n x.bond_indices,\n x.atom_index_0,\n x.atom_index_1), axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#test[['is_cis','is_trans']] = test.progress_apply(lambda x: is_Cis_Trans(x.molecule_name,\n# x.bond_indices,\n# x.atom_index_0,\n# x.atom_index_1), axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#train.to_csv('train_cis_trans.csv')\n#test.to_csv('test_cis_trans.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Visualize molecule\n\nReference
\nhttps://www.kaggle.com/mykolazotko/3d-visualization-of-molecules-with-plotly"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import plotly.graph_objs as go\nfrom plotly.offline import init_notebook_mode, iplot\nfrom sympy.geometry import Point3D\n\n\ninit_notebook_mode(connected=True)\n\ndef plot_molecule(molecule_name, structures_df):\n \"\"\"Creates a 3D plot of the molecule\"\"\"\n \n atomic_radii = dict(C=0.77, F=0.71, H=0.38, N=0.75, O=0.73) \n cpk_colors = dict(C='black', F='green', H='white', N='blue', O='red')\n\n molecule = structures_df[structures_df.molecule_name == molecule_name]\n coordinates = molecule[['x', 'y', 'z']].values\n x_coordinates = coordinates[:, 0]\n y_coordinates = coordinates[:, 1]\n z_coordinates = coordinates[:, 2]\n elements = molecule.atom.tolist()\n radii = [atomic_radii[element] for element in elements]\n \n def get_bonds():\n \"\"\"Generates a set of bonds from atomic cartesian coordinates\"\"\"\n ids = np.arange(coordinates.shape[0])\n bonds = dict()\n coordinates_compare, radii_compare, ids_compare = coordinates, radii, ids\n \n for _ in range(len(ids)):\n coordinates_compare = np.roll(coordinates_compare, -1, axis=0)\n radii_compare = np.roll(radii_compare, -1, axis=0)\n ids_compare = np.roll(ids_compare, -1, axis=0)\n distances = np.linalg.norm(coordinates - coordinates_compare, axis=1)\n bond_distances = (radii + radii_compare) * 1.3\n mask = np.logical_and(distances > 0.1, distances < bond_distances)\n distances = distances.round(2)\n new_bonds = {frozenset([i, j]): dist for i, j, dist in zip(ids[mask], ids_compare[mask], distances[mask])}\n bonds.update(new_bonds)\n return bonds \n \n def atom_trace():\n \"\"\"Creates an atom trace for the plot\"\"\"\n colors = [cpk_colors[element] for element in elements]\n markers = dict(color=colors, line=dict(color='lightgray', width=2), size=7, symbol='circle', opacity=0.8)\n trace = go.Scatter3d(x=x_coordinates, y=y_coordinates, z=z_coordinates, mode='markers', marker=markers,\n text=elements, name='')\n return trace\n\n def bond_trace():\n \"\"\"\"Creates a bond trace for the plot\"\"\"\n trace = go.Scatter3d(x=[], y=[], z=[], hoverinfo='none', mode='lines',\n marker=dict(color='grey', size=7, opacity=1))\n for i, j in bonds.keys():\n trace['x'] += (x_coordinates[i], x_coordinates[j], None)\n trace['y'] += (y_coordinates[i], y_coordinates[j], None)\n trace['z'] += (z_coordinates[i], z_coordinates[j], None)\n return trace\n \n bonds = get_bonds()\n \n zipped = zip(range(len(elements)), x_coordinates, y_coordinates, z_coordinates)\n annotations_id = [dict(text=num, x=x, y=y, z=z, showarrow=False, yshift=15)\n for num, x, y, z in zipped]\n \n annotations_length = []\n for (i, j), dist in bonds.items():\n p_i, p_j = Point3D(coordinates[i]), Point3D(coordinates[j])\n p = p_i.midpoint(p_j)\n annotation = dict(text=dist, x=float(p.x), y=float(p.y), z=float(p.z), showarrow=False, yshift=15)\n annotations_length.append(annotation) \n \n updatemenus = list([\n dict(buttons=list([\n dict(label = 'Atom indices',\n method = 'relayout',\n args = [{'scene.annotations': annotations_id}]),\n dict(label = 'Bond lengths',\n method = 'relayout',\n args = [{'scene.annotations': annotations_length}]),\n dict(label = 'Atom indices & Bond lengths',\n method = 'relayout',\n args = [{'scene.annotations': annotations_id + annotations_length}]),\n dict(label = 'Hide all',\n method = 'relayout',\n args = [{'scene.annotations': []}])\n ]),\n direction='down',\n xanchor = 'left',\n yanchor = 'top'\n ), \n ])\n \n data = [atom_trace(), bond_trace()]\n axis_params = dict(showgrid=False, showticklabels=False, zeroline=False, titlefont=dict(color='white'),\n showbackground=False)\n layout = dict(scene=dict(xaxis=axis_params, yaxis=axis_params, zaxis=axis_params, annotations=annotations_id), \n margin=dict(r=0, l=0, b=0, t=0), showlegend=False, updatemenus=updatemenus)\n\n fig = go.Figure(data=data, layout=layout)\n iplot(fig)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## cis "},{"metadata":{"trusted":true},"cell_type":"code","source":"cis_mol = train[train.is_cis == 1].iloc[0]\nprint(cis_mol[['atom_index_0','atom_index_1']])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_molecule(cis_mol.molecule_name, structures)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## trans"},{"metadata":{"trusted":true},"cell_type":"code","source":"trans_mol = train[train.is_trans == 1].iloc[0]\nprint(trans_mol[['atom_index_0','atom_index_1']])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"plot_molecule(trans_mol.molecule_name, structures)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# The effects of cis–trans isomerism on scalar coupling constant"},{"metadata":{},"cell_type":"markdown","source":"## Dihedrals Data\n\nkernel: https://www.kaggle.com/soerendip/calculate-angles-and-dihedrals-with-networkx/output
\ndataset: https://www.kaggle.com/soerendip/angle-and-dihedral-for-the-champs-structures"},{"metadata":{"trusted":true},"cell_type":"code","source":"angles = pd.read_csv(f'../input/calculate-angles-and-dihedrals-with-networkx/angles.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = pd.merge(train, \n angles[['molecule_name','atom_index_0','atom_index_1','dihedral']],\n how='left',\n on=['molecule_name','atom_index_0','atom_index_1'])","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Calculate distance\nReference
\nhttps://www.kaggle.com/seriousran/just-speed-up-calculate-distance-from-benchmark"},{"metadata":{"trusted":true},"cell_type":"code","source":"def map_atom_info(df, atom_idx):\n df = pd.merge(df, structures, how = 'left',\n left_on = ['molecule_name', f'atom_index_{atom_idx}'],\n right_on = ['molecule_name', 'atom_index'])\n \n df = df.drop('atom_index', axis=1)\n df = df.rename(columns={'atom': f'atom_{atom_idx}',\n 'x': f'x_{atom_idx}',\n 'y': f'y_{atom_idx}',\n 'z': f'z_{atom_idx}'})\n return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = map_atom_info(train, 0)\ntrain = map_atom_info(train, 1)\n\ntrain_p_0 = train[['x_0', 'y_0', 'z_0']].values\ntrain_p_1 = train[['x_1', 'y_1', 'z_1']].values\n\ntrain['dist'] = np.linalg.norm(train_p_0 - train_p_1, axis=1)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Plot 3JHC mol feature"},{"metadata":{"trusted":true},"cell_type":"code","source":"from mpl_toolkits.mplot3d import Axes3D\nfrom scipy import genfromtxt\nimport matplotlib.pyplot as plt\nfrom pylab import rcParams\nimport pylab","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_3JHC = train[train['type'] == '3JHC']\ntrain_normal = train_3JHC[(train_3JHC['is_cis'] == 0) & (train_3JHC['is_trans'] == 0)]\ntrain_cis = train_3JHC[train_3JHC['is_cis'] == 1]\ntrain_trans = train_3JHC[train_3JHC['is_trans'] == 1]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"rcParams['figure.figsize'] = 10,4\ntrain_normal['scalar_coupling_constant'].hist(bins=50)\npylab.suptitle(\"Not cis-trans\")\nplt.show()\ntrain_cis['scalar_coupling_constant'].hist(bins=30)\npylab.suptitle(\"Cis\")\nplt.show()\ntrain_trans['scalar_coupling_constant'].hist(bins=30)\npylab.suptitle(\"Trans\")\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"rcParams['figure.figsize'] = 10,7\n\nfig = plt.figure()\nax = Axes3D(fig)\n\nax.set_xlabel(\"dihedral\")\nax.set_ylabel(\"scalar_coupling_constant\")\nax.set_zlabel(\"dist\")\n\nax.plot(train_normal.dihedral, \n train_normal.dist, \n train_normal.scalar_coupling_constant,\n \"o\", color='gray', ms=1, mew=0.2)\nax.plot(train_cis.dihedral,\n train_cis.dist,\n train_cis.scalar_coupling_constant,\n \"o\", color='blue', ms=2, mew=2)\nax.plot(train_trans.dihedral,\n train_trans.dist,\n train_trans.scalar_coupling_constant,\n \"o\", color='red', ms=2, mew=2)\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Thank you for reading.
\nPlease let me know if I made a mistake."},{"metadata":{},"cell_type":"markdown","source":"# Add analysis"},{"metadata":{"trusted":true},"cell_type":"code","source":"train_3 = train[(train['type']=='3JHC') | (train['type']=='3JHH') | (train['type']=='3JHN')]\ntrain_normal = train_3[(train_3['is_cis'] == 0) & (train_3['is_trans'] == 0)]\ntrain_cis = train_3[train_3['is_cis'] == 1]\ntrain_trans = train_3[train_3['is_trans'] == 1]","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"There are very few cis-trans isomerism in dataset.
\nSo the score won't improve significantly."},{"metadata":{"trusted":true},"cell_type":"code","source":"rcParams['figure.figsize'] = 5,7\nplt.style.use('ggplot')\nleft = np.array(['not Cis-Trans', 'Cis', 'Trans'])\nheight = np.array([len(train_normal), len(train_cis), len(train_trans)])\nplt.bar(left, height)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train_normal.plot(kind='scatter', x='dist', y='scalar_coupling_constant',\n figsize=(10,4), title='not Cis and not Trans')\ntrain_cis.plot(kind='scatter', x='dist', y='scalar_coupling_constant',\n figsize=(10,4), title='is Cis')\ntrain_trans.plot(kind='scatter', x='dist', y='scalar_coupling_constant',\n figsize=(10,4), title='is Trans')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"When I tried, is-Cis feature slightly improved the score.
"}],"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}