{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":67356,"databundleVersionId":8006601,"sourceType":"competition"},{"sourceId":174975929,"sourceType":"kernelVersion"},{"sourceId":174976137,"sourceType":"kernelVersion"},{"sourceId":175147484,"sourceType":"kernelVersion"},{"sourceId":175147555,"sourceType":"kernelVersion"}],"dockerImageVersionId":30698,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"!pip install rdkit","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-05-10T22:18:06.924569Z","iopub.execute_input":"2024-05-10T22:18:06.924980Z","iopub.status.idle":"2024-05-10T22:18:24.706600Z","shell.execute_reply.started":"2024-05-10T22:18:06.924947Z","shell.execute_reply":"2024-05-10T22:18:24.704765Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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\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 read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\nfrom rdkit import Chem\nimport os\nfor dirname, _, filenames in os.walk('../input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-05-10T22:18:24.709505Z","iopub.execute_input":"2024-05-10T22:18:24.710012Z","iopub.status.idle":"2024-05-10T22:18:25.310384Z","shell.execute_reply.started":"2024-05-10T22:18:24.709976Z","shell.execute_reply":"2024-05-10T22:18:25.309092Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## What's going on?\n\nThis notebook provides precomputed matches between building blocks and corresponding fragments found in the molecules. The matching procedure is described at [my main notebook](https://www.kaggle.com/code/latticetower/belka-is-nuts-train-preprocessing-eda#Fragment-analysis). It is a time-consuming computation, so all the data is splitted to several parts and processed separately. Also, I focus only on positive samples, since it it important to understand if there is a difference in building blocks connections and linker position or not.\n\n**This notebook will be updated, however, it is an addition to my main notebook with EDA.**\n\nIn all input files used here there are 3 additional columns: frag1, frag2, frag3. The naming is arbitrary, the important thing here is that `frag1` is the molecular fragment in `molecule_smiles` which corresponds to `buildingblock1_smiles`, `frag2` corresponds to `buildingblock2_smiles`, `frag3` corresponds to `buildingblock3_smiles`.\n\nBy looking at them we can understand how the building blocks are combine to form the molecules from our dataset and how the DNA linker was attached.\n\nThe main limitation for my code is that I've processed only the molecules which were connected with triazine core, so only the part of the test dataset was processed in this manner.","metadata":{}},{"cell_type":"code","source":"train_mixed_sEH_df1 = pd.read_parquet(\"../input/belka-is-nuts-precompute-fragments-for-train/train_all_mixed_with_fragments_sEH_p1.parquet.gzip\")\ntrain_mixed_sEH_df2 = pd.read_parquet(\"../input/belka-is-nuts-fragments-train-seh-p2/train_all_mixed_with_fragments_sEH_p2.parquet.gzip\")\ntrain_mixed_not_sEH_df = pd.read_parquet(\"../input/belka-is-nuts-precompute-fragments-for-train-2/train_all_mixed_with_fragments_not_sEH.parquet.gzip\")\ntrain_mixed_sEH_df1.shape, train_mixed_sEH_df2.shape, train_mixed_not_sEH_df.shape","metadata":{"execution":{"iopub.status.busy":"2024-05-10T22:18:25.317222Z","iopub.execute_input":"2024-05-10T22:18:25.317926Z","iopub.status.idle":"2024-05-10T22:18:28.419342Z","shell.execute_reply.started":"2024-05-10T22:18:25.317885Z","shell.execute_reply":"2024-05-10T22:18:28.418189Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fragment_columns = [f\"frag{i+1}\" for i in range(3)]","metadata":{"execution":{"iopub.status.busy":"2024-05-10T22:18:28.420713Z","iopub.execute_input":"2024-05-10T22:18:28.421025Z","iopub.status.idle":"2024-05-10T22:18:28.426277Z","shell.execute_reply.started":"2024-05-10T22:18:28.420999Z","shell.execute_reply":"2024-05-10T22:18:28.425127Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# train_ids1 = ~train_mixed_sEH_df1[fragment_columns].isnull().any(axis=1)\n# train_ids2 = ~train_mixed_sEH_df2[fragment_columns].isnull().any(axis=1)\n\ntrain_index1 = train_mixed_sEH_df1.index[:400000]\ntrain_index2 = train_mixed_sEH_df2.index[400000:]\n\ntrain_mixed_sEH_df = pd.concat([\n    train_mixed_sEH_df1.loc[train_index1],\n    train_mixed_sEH_df2.loc[train_index2]\n])","metadata":{"execution":{"iopub.status.busy":"2024-05-10T22:18:28.427738Z","iopub.execute_input":"2024-05-10T22:18:28.428280Z","iopub.status.idle":"2024-05-10T22:18:28.803769Z","shell.execute_reply.started":"2024-05-10T22:18:28.428225Z","shell.execute_reply":"2024-05-10T22:18:28.802472Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_mixed_df = pd.concat([train_mixed_sEH_df, train_mixed_not_sEH_df])\n\ntrain_mixed_df.to_parquet(\"train_mixed_with_fragments.parquet.gzip\", compression=\"gzip\")","metadata":{"execution":{"iopub.status.busy":"2024-05-10T22:18:28.805495Z","iopub.execute_input":"2024-05-10T22:18:28.805953Z","iopub.status.idle":"2024-05-10T22:18:37.162483Z","shell.execute_reply.started":"2024-05-10T22:18:28.805914Z","shell.execute_reply":"2024-05-10T22:18:37.161658Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# train_mixed_sEH_df = train_mixed_sEH_df[train_mixed_sEH_df.sEH == 1]\ntrain_all1_df = pd.read_parquet(\"../input/belka-is-nuts-precompute-fragments-for-train-2/train_all1_with_fragments.parquet\")\ntrain_all1_df.to_parquet(\"train_all1_with_fragments.parquet\")","metadata":{"execution":{"iopub.status.busy":"2024-05-10T22:18:37.164042Z","iopub.execute_input":"2024-05-10T22:18:37.164755Z","iopub.status.idle":"2024-05-10T22:18:37.184217Z","shell.execute_reply.started":"2024-05-10T22:18:37.164711Z","shell.execute_reply":"2024-05-10T22:18:37.183171Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# np.sum(train_mixed_sEH_df2[train_mixed_sEH_df2.sEH!=1][fragment_columns].isnull().sum(1) < 1)\ntrain_any1_df = pd.concat([train_mixed_df[train_all1_df.columns], train_all1_df], ignore_index=True)","metadata":{"execution":{"iopub.status.busy":"2024-05-10T22:18:37.185468Z","iopub.execute_input":"2024-05-10T22:18:37.186044Z","iopub.status.idle":"2024-05-10T22:18:37.689629Z","shell.execute_reply.started":"2024-05-10T22:18:37.186012Z","shell.execute_reply":"2024-05-10T22:18:37.688028Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"matched_fragments_ids = ~train_any1_df[fragment_columns].isnull().any(axis=1)\nprint(\"Molecules with no match between blocks and fragments: \", (~matched_fragments_ids).sum())\nfor i in range(3):\n    df = train_any1_df[matched_fragments_ids].groupby(f\"buildingblock{i+1}_smiles\", as_index=False).agg({\n        f\"frag{i+1}\": lambda x: np.unique(x)\n    })\n    print(df[f\"frag{i+1}\"].apply(len).value_counts())","metadata":{"execution":{"iopub.status.busy":"2024-05-10T22:18:37.693135Z","iopub.execute_input":"2024-05-10T22:18:37.693628Z","iopub.status.idle":"2024-05-10T22:18:42.662417Z","shell.execute_reply.started":"2024-05-10T22:18:37.693572Z","shell.execute_reply":"2024-05-10T22:18:42.660835Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"matched_fragments_ids = ~train_any1_df[fragment_columns].isnull().any(axis=1)\nprint(\"Molecules with no match between blocks and fragments: \", (~matched_fragments_ids).sum())\nfor i in range(3):\n    df = train_any1_df[matched_fragments_ids].groupby(f\"frag{i+1}\", as_index=False).agg({\n        f\"buildingblock{i+1}_smiles\": lambda x: np.unique(x)\n    })\n    print(df[f\"buildingblock{i+1}_smiles\"].apply(len).value_counts())","metadata":{"execution":{"iopub.status.busy":"2024-05-10T22:19:13.955990Z","iopub.execute_input":"2024-05-10T22:19:13.956351Z","iopub.status.idle":"2024-05-10T22:19:18.856963Z","shell.execute_reply.started":"2024-05-10T22:19:13.956322Z","shell.execute_reply":"2024-05-10T22:19:18.855698Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We know now that all the fragments appearing in binding small molecules can be matched with building blocks with my algo and each block corresponds to no more than 1 fragment. This means that every block is processed in one unique way and it corresponds to one particular set of atoms and bonds withing the binding molecule.\n\nWe won't run this computation on non-binding part of the dataset, we'll assume that in negative samples this rule holds and each block is used to extract the same specific fragment.","metadata":{}},{"cell_type":"code","source":"# from rdkit.Chem import rdFMCS\n\n# fragments_set = list(Chem.GetMolFrags(all_splits[1], asMols=True, sanitizeFrags=True))\n\n# mol_pairs = []\n# all_hit_atoms = []\n# all_hit_bonds = []\n# legends = []\n# for i in range(len(fragments_set)):\n#     fragment = Chem.AddHs(fragments_set[i])\n#     for j in range(3):\n#         bb_mol = Chem.AddHs(mol_bb[j])\n        \n#         match_result = rdFMCS.FindMCS(\n#             [fragment, bb_mol], \n#             matchValences=True,\n#             ringMatchesRingOnly=True,\n#             completeRingsOnly=True\n#         )\n#         pattern = Chem.MolFromSmarts(match_result.smartsString)\n        \n#         hit_atoms = list(fragment.GetSubstructMatch(pattern))\n#         hit_bonds = []\n#         for bond in pattern.GetBonds():\n#             aid1 = hit_atoms[bond.GetBeginAtomIdx()]\n#             aid2 = hit_atoms[bond.GetEndAtomIdx()]\n#             hit_bonds.append(fragment.GetBondBetweenAtoms(aid1,aid2).GetIdx())\n#         all_hit_atoms.append(hit_atoms)\n#         all_hit_bonds.append(hit_bonds)\n        \n#         hit_atoms = list(bb_mol.GetSubstructMatch(pattern))\n#         hit_bonds = []\n#         for bond in pattern.GetBonds():\n#             aid1 = hit_atoms[bond.GetBeginAtomIdx()]\n#             aid2 = hit_atoms[bond.GetEndAtomIdx()]\n#             hit_bonds.append(bb_mol.GetBondBetweenAtoms(aid1,aid2).GetIdx())\n#         all_hit_atoms.append(hit_atoms)\n#         all_hit_bonds.append(hit_bonds)\n        \n#         d = max(fragment.GetNumAtoms(), bb_mol.GetNumAtoms())\n#         has_linker1 = np.any([a.GetSymbol() == 'Dy' for a in fragment.GetAtoms()])\n#         if has_linker1:\n#             d = fragment.GetNumAtoms()\n#         has_linker2 = np.any([a.GetSymbol() == 'Dy' for a in bb_mol.GetAtoms()])\n#         if has_linker2:\n#             d = bb_mol.GetNumAtoms()\n#         score = match_result.numAtoms / d\n#         mol_pairs.extend([fragment, bb_mol])\n#         legends.extend([\"\", f\"fragment {i} and building block {j}: score {score:3.3f}\"])\n        \n\n# Chem.Draw.MolsToGridImage(\n#     mol_pairs, \n#     subImgSize=(300,300), \n#     legends=legends,\n#     molsPerRow=6,\n#     highlightAtomLists=all_hit_atoms\n# )","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-05-10T22:18:42.663971Z","iopub.execute_input":"2024-05-10T22:18:42.664298Z","iopub.status.idle":"2024-05-10T22:18:42.670078Z","shell.execute_reply.started":"2024-05-10T22:18:42.664271Z","shell.execute_reply":"2024-05-10T22:18:42.669170Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# from scipy.optimize import linear_sum_assignment\n# from rdkit import Chem\n# from rdkit.Chem import rdFMCS\n\n# def make_query(mol):\n#     p = Chem.AdjustQueryParameters.NoAdjustments()\n#     p.makeDummiesQueries = True\n#     mol = Chem.AdjustQueryProperties(mol, p)\n#     mol.UpdatePropertyCache()\n#     return mol\n\n\n# def remove_dummy_number(m):\n#     star = Chem.MolFromSmiles(\"*\")\n#     new_mol, = Chem.ReplaceSubstructs(m, star, star, replaceAll=True)\n#     return new_mol\n\n\n# def split_by_triazine_core(smiles, debug=False):\n#     molecule = Chem.MolFromSmiles(smiles)\n#     triazine_core = make_query(Chem.MolFromSmiles(\"c1ncncn1\"))\n#     for match in molecule.GetSubstructMatches(triazine_core):\n#         fragments = Chem.ReplaceCore(\n#             molecule, \n#             triazine_core, \n#             match, replaceDummies=True, \n#             # requireDummyMatch=True\n#             requireDummyMatch=False\n#         )\n#         fragments = [remove_dummy_number(x) for x in Chem.GetMolFrags(fragments, asMols=True, sanitizeFrags=True)]\n#         fragments = [Chem.MolToSmiles(x) for x in fragments]\n#         if len(fragments) != 3:\n#             if debug:\n#                 print(fragments)\n#             continue\n#         yield fragments\n\n\n# def score_pair(mol1, mol2, use_mcs=True, debug=False):\n#     if use_mcs:\n#         from rdkit.Chem import rdFMCS\n#         match_result = rdFMCS.FindMCS(\n#             [mol1, mol2],\n#             matchValences=True, \n#             ringMatchesRingOnly=True,\n#             completeRingsOnly=True,\n#             bondCompare=rdFMCS.BondCompare.CompareOrderExact\n#         )\n        \n#         d = max(mol1.GetNumAtoms(), mol2.GetNumAtoms())\n#         has_linker1 = np.any([a.GetSymbol() == 'Dy' for a in mol1.GetAtoms()])\n#         if has_linker1:\n#             d = mol1.GetNumAtoms()\n#         has_linker2 = np.any([a.GetSymbol() == 'Dy' for a in mol2.GetAtoms()])\n#         if has_linker2:\n#             d = mol2.GetNumAtoms()\n#         score = match_result.numAtoms/d\n#         if debug:\n#             print(repr(Chem.MolToSmiles(Chem.RemoveHs(mol1))), \n#                   \"   \", \n#                   repr(Chem.MolToSmiles(Chem.RemoveHs(mol2))), \n#                   \" - \", score)\n#         return score\n\n#     matches = mol2.GetSubstructMatches(mol1, useChirality=True)\n#     if len(matches) == 0:\n#         return 0.\n#     r = max([len(m) for m in matches])\n#     return r / d\n\n\n# def score_sets(frag_mols, mols_bb, debug=False):\n#     scores = [[\n#             score_pair(a, b, debug=debug) for b in mols_bb\n#         ]\n#         for a in frag_mols\n#     ]\n#     frag_ind, bb_ind = linear_sum_assignment(scores, maximize=True)\n#     if debug:\n#         print(\"\\nScores:\")\n#         print(scores)\n#         print(frag_ind, bb_ind)\n#     score = np.sum([scores[i][j] for i, j in zip(frag_ind, bb_ind)])\n#     if debug:\n#         print(\"\\nScore:\")\n#         print([scores[i][j] for i, j in zip(frag_ind, bb_ind)])\n#         print(score)\n#     # print(scores)\n#     # print(row_ind, col_ind)\n#     idx = np.argsort(bb_ind)\n#     bb_ind = bb_ind[idx]\n#     frag_ind = frag_ind[idx]\n#     # ordered molecular fragments corresponding to building blocks 0, 1, 2\n#     return (frag_ind, bb_ind), score \n\n\n# def match_fragments(smiles, smiles_bb, debug=False):\n#     mol_bb = [Chem.MolFromSmiles(bb, sanitize=True) for bb in smiles_bb]\n#     mol_bb = [Chem.AddHs(mol) for mol in mol_bb]\n#     max_score = 0\n#     indices = None\n#     assignment = None\n#     for fragments in split_by_triazine_core(smiles):\n#         if debug:\n#             print(\"split\", fragments)\n#         frag_mols = [Chem.MolFromSmiles(smiles) for smiles in fragments]\n#         frag_mols = [Chem.AddHs(mol) for mol in frag_mols]\n#         (frag_ind, bb_ind), score = score_sets(frag_mols, mol_bb, debug=debug)\n#         if debug:\n#             print(\"Score\", score)\n#         if score > max_score:\n#             max_score = score\n#             indices = (bb_ind, frag_ind)\n#             assignment = [fragments[i] for i in frag_ind]\n#     return assignment\n\n# match_fragments(smiles, smiles_bb, debug=True)","metadata":{"_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2024-05-10T22:18:42.671372Z","iopub.execute_input":"2024-05-10T22:18:42.671886Z","iopub.status.idle":"2024-05-10T22:18:42.955588Z","shell.execute_reply.started":"2024-05-10T22:18:42.671855Z","shell.execute_reply":"2024-05-10T22:18:42.954548Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# fragments_set = list(Chem.GetMolFrags(all_splits[1], asMols=True, sanitizeFrags=True))\n\n# mol_pairs = []\n# all_hit_atoms = []\n# all_hit_bonds = []\n# legends = []\n# for i in range(len(fragments_set)):\n#     fragment = fragments_set[i]\n#     for j in range(3):\n#         bb_mol = bb_molecules[j]\n        \n#         match_result = rdFMCS.FindMCS(\n#             [fragment, bb_mol], \n#             matchValences=True,\n#             ringMatchesRingOnly=True,\n#             completeRingsOnly=True\n#         )\n#         pattern = Chem.MolFromSmarts(match_result.smartsString)\n        \n#         hit_atoms = list(fragment.GetSubstructMatch(pattern))\n#         hit_bonds = []\n#         for bond in pattern.GetBonds():\n#             aid1 = hit_atoms[bond.GetBeginAtomIdx()]\n#             aid2 = hit_atoms[bond.GetEndAtomIdx()]\n#             hit_bonds.append(fragment.GetBondBetweenAtoms(aid1,aid2).GetIdx())\n#         all_hit_atoms.append(hit_atoms)\n#         all_hit_bonds.append(hit_bonds)\n        \n#         hit_atoms = list(bb_mol.GetSubstructMatch(pattern))\n#         hit_bonds = []\n#         for bond in pattern.GetBonds():\n#             aid1 = hit_atoms[bond.GetBeginAtomIdx()]\n#             aid2 = hit_atoms[bond.GetEndAtomIdx()]\n#             hit_bonds.append(bb_mol.GetBondBetweenAtoms(aid1,aid2).GetIdx())\n#         all_hit_atoms.append(hit_atoms)\n#         all_hit_bonds.append(hit_bonds)\n        \n#         d = max(fragment.GetNumAtoms(), bb_mol.GetNumAtoms())\n#         has_linker1 = np.any([a.GetSymbol() == 'Dy' for a in fragment.GetAtoms()])\n#         if has_linker1:\n#             d = fragment.GetNumAtoms()\n#         has_linker2 = np.any([a.GetSymbol() == 'Dy' for a in bb_mol.GetAtoms()])\n#         if has_linker2:\n#             d = bb_mol.GetNumAtoms()\n#         score = match_result.numAtoms / d\n#         mol_pairs.extend([fragment, bb_mol])\n#         legends.extend([\"\", f\"fragment {i} and building block {j}: score {score:3.3f}\"])\n        \n\n# Chem.Draw.MolsToGridImage(\n#     mol_pairs, \n#     subImgSize=(300,300), \n#     legends=legends,\n#     molsPerRow=6,\n#     highlightAtomLists=all_hit_atoms\n# )","metadata":{"_kg_hide-output":true,"_kg_hide-input":true,"execution":{"iopub.status.busy":"2024-05-10T22:18:42.957922Z","iopub.execute_input":"2024-05-10T22:18:42.958376Z","iopub.status.idle":"2024-05-10T22:18:42.977242Z","shell.execute_reply.started":"2024-05-10T22:18:42.958336Z","shell.execute_reply":"2024-05-10T22:18:42.976352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}