{"cells":[{"metadata":{},"cell_type":"markdown","source":"This kernel shows the importance of ring structure. This feature is proposed by  \nhttps://www.kaggle.com/sunhwan/using-rdkit-for-atomic-feature-and-visualization. I've just plotted the histograms. If this kernel is helpful, please not upvote on me, but upvote on him."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true,"_kg_hide-output":true,"collapsed":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list the files in the input directory\n\nimport os\nprint(os.listdir(\"../input\"))\n\n# Any results you write to the current directory are saved as output.","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","trusted":true,"_kg_hide-output":true,"_kg_hide-input":true,"collapsed":true},"cell_type":"code","source":"import sys\n!conda install --yes --prefix {sys.prefix} -c rdkit rdkit","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":true,"_kg_hide-input":true,"collapsed":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-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# only reading 10% of data for debug\ntrain = pd.read_csv('../input/train.csv')[::10]\ntest = pd.read_csv('../input/test.csv')[::10]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true,"collapsed":true},"cell_type":"code","source":"# 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    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(), test.molecule_name.unique()])\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,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# highlight the bond given in train/test\n# http://rdkit.blogspot.com/2015/02/new-drawing-code.html\nm = mols['dsgdb9nsd_002129']\natom_index_0 = int(train[train.molecule_name == 'dsgdb9nsd_002129'].iloc[0].atom_index_0)\natom_index_1 = int(train[train.molecule_name == 'dsgdb9nsd_002129'].iloc[0].atom_index_1)\nhighlight = [atom_index_0, atom_index_1]\n\nfrom IPython.display import SVG\nChem.rdDepictor.Compute2DCoords(m)\ndrawer = Draw.rdMolDraw2D.MolDraw2DSVG(400, 200)\ndrawer.DrawMolecule(m,highlightAtoms=highlight)\ndrawer.FinishDrawing()\nsvg = drawer.GetDrawingText().replace('svg:','')\n#SVG(svg)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# 3JHH\nrow = train[train.type == '3JHH'].iloc[1]\nm = mols[row.molecule_name]\natom_index_0 = int(row.atom_index_0)\natom_index_1 = int(row.atom_index_1)\nhighlight = [atom_index_0, atom_index_1]\n\nfrom IPython.display import SVG\nChem.rdDepictor.Compute2DCoords(m)\ndrawer = Draw.rdMolDraw2D.MolDraw2DSVG(400, 200)\ndrawer.DrawMolecule(m,highlightAtoms=highlight)\ndrawer.FinishDrawing()\nsvg = drawer.GetDrawingText().replace('svg:','')\n#SVG(svg)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"from pathlib import Path\nfrom sklearn import *\nPATH = Path('../input')\n\n# again, only using 10% of data for the committed kernel to work...\ntrain = pd.read_csv(PATH/'train.csv')[::10]\ntest = pd.read_csv(PATH/'test.csv')[::10]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# https://www.kaggle.com/jazivxt/all-this-over-a-dog\ntrain['atom1'] = train['type'].map(lambda x: str(x)[2])\ntrain['atom2'] = train['type'].map(lambda x: str(x)[3])\ntest['atom1'] = test['type'].map(lambda x: str(x)[2])\ntest['atom2'] = test['type'].map(lambda x: str(x)[3])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"lbl = preprocessing.LabelEncoder()\nfor i in range(4):\n    train['type'+str(i)] = lbl.fit_transform(train['type'].map(lambda x: str(x)[i]))\n    test['type'+str(i)] = lbl.transform(test['type'].map(lambda x: str(x)[i]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"structures = pd.read_csv(PATH/'structures.csv').rename(columns={'atom_index':'atom_index_0', 'x':'x0', 'y':'y0', 'z':'z0', 'atom':'atom1'})\ntrain = pd.merge(train, structures, how='left', on=['molecule_name', 'atom_index_0', 'atom1'])\ntest = pd.merge(test, structures, how='left', on=['molecule_name', 'atom_index_0', 'atom1'])\ndel structures","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"structures = pd.read_csv(PATH/'structures.csv').rename(columns={'atom_index':'atom_index_1', 'x':'x1', 'y':'y1', 'z':'z1', 'atom':'atom2'})\ntrain = pd.merge(train, structures, how='left', on=['molecule_name', 'atom_index_1', 'atom2'])\ntest = pd.merge(test, structures, how='left', on=['molecule_name', 'atom_index_1', 'atom2'])\ndel structures","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"def feature_atom(atom):\n    prop = {}\n    nb = [a.GetSymbol() for a in atom.GetNeighbors()] # neighbor atom type symbols\n    nb_h = sum([_ == 'H' for _ in nb]) # number of hydrogen as neighbor\n    nb_o = sum([_ == 'O' for _ in nb]) # number of oxygen as neighbor\n    nb_c = sum([_ == 'C' for _ in nb]) # number of carbon as neighbor\n    nb_n = sum([_ == 'N' for _ in nb]) # number of nitrogen as neighbor\n    nb_na = len(nb) - nb_h - nb_o - nb_n - nb_c\n    prop['degree'] = atom.GetDegree()\n    prop['hybridization'] = int(atom.GetHybridization())\n    prop['inring'] = int(atom.IsInRing()) # is the atom in a ring?\n    prop['inring3'] = int(atom.IsInRingSize(3)) # is the atom in a ring size of 3?\n    prop['inring4'] = int(atom.IsInRingSize(4)) # is the atom in a ring size of 4?\n    prop['inring5'] = int(atom.IsInRingSize(5)) # ...\n    prop['inring6'] = int(atom.IsInRingSize(6))\n    prop['inring7'] = int(atom.IsInRingSize(7))\n    prop['inring8'] = int(atom.IsInRingSize(8))\n    prop['nb_h'] = nb_h\n    prop['nb_o'] = nb_o\n    prop['nb_c'] = nb_c\n    prop['nb_n'] = nb_n\n    prop['nb_na'] = nb_na\n    return prop","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-output":true,"_kg_hide-input":true},"cell_type":"code","source":"# atom feature of dsgdb9nsd_002129 atom_index_0\nmolecule_name = 'dsgdb9nsd_002129'\nrow = train[train.molecule_name == molecule_name].iloc[0]\natom_index_0 = int(row.atom_index_0)\natom_index_1 = int(row.atom_index_1)\nm=mols[molecule_name]\na0 = m.GetAtomWithIdx(atom_index_0)\n#feature_atom(a0)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true,"collapsed":true},"cell_type":"code","source":"# extract some simple atomic feature for atom_index_0 and atom_index_1\n\n# use cached rdkit mol object to save memory\nif 'mols' in locals(): del mols\nimport gc\ngc.collect()\n\n# fix atom bonds\n# dsgdb9nsd_059827: hydrogen has is far apart\nnblist = {\n    'dsgdb9nsd_059827': {\n        13: 3\n    }\n}\n\ndef _features(args):\n    idx, row = args\n    molecule_name = row.molecule_name\n    atom_index_0 = int(row.atom_index_0)\n    atom_index_1 = int(row.atom_index_1)\n    \n    prop = {'molecule_name': molecule_name,\n            'atom_index_0': atom_index_0,\n            'atom_index_1': atom_index_1}\n\n    # atom_0 is always hydrogen\n    m = MolFromXYZ(PATH/f'structures/{molecule_name}.xyz') # less memory intensive in multiprocessing.Pool\n    a0 = m.GetAtomWithIdx(atom_index_0)\n\n    a1 = m.GetAtomWithIdx(atom_index_1)\n    a1_prop = feature_atom(a1)\n    prop.update({'a1_'+k: a1_prop[k] for k in a1_prop.keys()})\n\n    # skipping below for time constraint\n    # neighbor of atom_0\n    try:\n        a0_nb_idx = [a.GetIdx() for a in a0.GetNeighbors() if a.GetIdx() != a0].pop()\n    except:\n        if molecule_name in nblist and atom_index_0 in nblist[molecule_name]:\n            a0_nb_idx = nblist[molecule_name][atom_index_0]\n        else:\n            print(molecule_name)\n            print(row)\n\n    a0_nb = m.GetAtomWithIdx(a0_nb_idx)\n    a0_nb_prop = feature_atom(a0_nb)\n    for k in a0_nb_prop.keys():\n        prop['a0_nb_'+k] = a0_nb_prop[k]\n        \n    c = m.GetConformer()\n    #prop['dist_a0_a0_nb'] = np.linalg.norm(c.GetAtomPosition(atom_index_0) - c.GetAtomPosition(a0_nb_idx))\n    prop['x_a0_nb'] = c.GetAtomPosition(a0_nb_idx)[0]\n    prop['y_a0_nb'] = c.GetAtomPosition(a0_nb_idx)[1]\n    prop['z_a0_nb'] = c.GetAtomPosition(a0_nb_idx)[2]\n\n    # neighbor of atom_1\n    try:\n        a1_nb_idx = [a.GetIdx() for a in a1.GetNeighbors() if a.GetIdx() != a1].pop()\n    except:\n        if molecule_name in nblist and atom_index_1 in nblist[molecule_name]:\n            a1_nb_idx = nblist[molecule_name][atom_index_1]\n        else:\n            print(molecule_name)\n            print(row)\n    a1_nb = m.GetAtomWithIdx(a1_nb_idx)\n    a1_nb_prop = feature_atom(a1_nb)\n    for k in a1_nb_prop.keys():\n        prop['a1_nb_'+k] = a1_nb_prop[k]\n    prop['x_a1_nb'] = c.GetAtomPosition(a1_nb_idx)[0]\n    prop['y_a1_nb'] = c.GetAtomPosition(a1_nb_idx)[1]\n    prop['z_a1_nb'] = c.GetAtomPosition(a1_nb_idx)[2]\n    #prop['dist_a1_a1_nb'] = np.linalg.norm(c.GetAtomPosition(a1.GetIdx()) - c.GetAtomPosition(a1_nb.GetIdx()))\n    #prop['dist_a0_a1_nb'] = np.linalg.norm(c.GetAtomPosition(a0.GetIdx()) - c.GetAtomPosition(a1_nb.GetIdx()))\n    #prop['dist_a1_a0_nb'] = np.linalg.norm(c.GetAtomPosition(a1.GetIdx()) - c.GetAtomPosition(a0_nb.GetIdx()))\n    return prop\n\ndef features(df):\n    prop = []\n    n_cpu = 4\n    with Pool(n_cpu) as p:\n        n = len(df)\n        res = _features((0, df.iloc[0]))\n        keys = res.keys()\n        _df = df[['molecule_name', 'atom_index_0', 'atom_index_1']]\n        with tqdm(total=n) as pbar:\n            for res in p.imap_unordered(_features, _df.iterrows()):\n                # this is faster than using dict\n                prop.append([res[_] for _ in keys])\n                pbar.update()\n        del _df\n    \n    prop = pd.DataFrame.from_records(prop, columns=keys)\n    df = pd.merge(df, prop, how='left', on=['molecule_name', 'atom_index_0', 'atom_index_1'])\n    return df\n\ntrain = features(train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"#https://www.kaggle.com/artgor/molecular-properties-eda-and-models\ntrain_p0 = train[['x0', 'y0', 'z0']].values\ntrain_p1 = train[['x1', 'y1', 'z1']].values\ntrain['dist'] = np.linalg.norm(train_p0 - train_p1, axis=1)\ntrain['dist_to_type_mean'] = train['dist'] / train.groupby('type')['dist'].transform('mean')\ndel train_p0, train_p1","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Visualization\nLet's focus on the features of the ring size 1st atom belongs to. There are several kinds of cyclic compounds.\n'none(0)-cyclic', 'tri(3)-cyclic', 'tetra(4)-cyclic', 'penta(5)-cyclic', 'hexa(6)-cyclic', 'hepta(7)-cyclic' and 'opta(8)-cyclic'. I'd like to see the frequency distribution of 'scalar_coupling_constant' of each ring structure features."},{"metadata":{"trusted":true},"cell_type":"code","source":"import matplotlib.pyplot as plt\nmol_types=train[\"type\"].unique()\n\ndef plot_cyclic_vs_scalar_coupling(mol_type):\n    train_=train[train[\"type\"]==mol_type].copy()\n    colorlist = [\"r\", \"g\", \"b\", \"c\", \"m\", \"y\", \"k\", \"w\"]\n    \n    fig, ax = plt.subplots(figsize = (12, 6))\n    plt.scatter(train_.loc[train_['a0_nb_inring']== 0 , 'scalar_coupling_constant'],train_.loc[train_['a0_nb_inring'] == 0, 'dist'],  label=str(0)+\"cyclic\",c=\"black\",alpha=0.2);\n    for i in [3,4,5,6,7,8]:\n        plt.scatter(train_.loc[train_['a0_nb_inring'+str(i)]== 1 , 'scalar_coupling_constant'],train_.loc[train_['a0_nb_inring'+str(i)] == 1, 'dist'],  label=str(i)+\"cyclic\",c=colorlist[i-3],alpha=0.2);\n    lgnd = ax.legend(loc=\"upper right\", numpoints=1, fontsize=15)\n    for i in range(7):\n        lgnd.legendHandles[i]._sizes = [300]\n    plt.ylabel(\"distance between atoms\", fontsize=15)\n    plt.xlabel(\"scalar_coupling_const\", fontsize=15)\n    plt.title(f'{mol_type} a0_neighbor ring_size vs scalar_coupling_const vs distance', fontsize=15)\n\n    fig, ax = plt.subplots(figsize = (12, 6))\n    plt.hist(train_.loc[train_['a0_nb_inring']== 0 , 'scalar_coupling_constant'],bins=30, label=\"0 cyclic\",color=\"black\",alpha=0.2)\n    for i in [3,4,5,6,7,8]:\n        plt.hist(train_.loc[train_['a0_nb_inring'+str(i)]== 1 , 'scalar_coupling_constant'],bins=30, label=str(i)+\" cyclic\",color=colorlist[i-3],alpha=0.5)\n    lgnd = ax.legend(loc=\"upper right\", numpoints=1, fontsize=15)\n    for i in range(7):\n        lgnd.legendHandles[i]._sizes = [300]\n    plt.title(f'{mol_type} a0_neighbor ring_size vs scalar_coupling_const', fontsize=15)\n    plt.xlabel(\"scalar_coupling_const\", fontsize=15)\n\nfor mol_type in mol_types:\n    plot_cyclic_vs_scalar_coupling(mol_type)\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"For some molecule types, the distribution of 'scalar_coupling_constant' looks depending on the ring size the atom belongs to. So this would be effective features to predict scalar coupling constant."}],"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}